{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5k0c_GZZAlSc"
   },
   "source": [
    "This code chunk will process your raster and turn it into a shapefile, then re-import the shapefile and clip it to the extent of AOI (or catchment) and export that as a geojson as the final output for vectorization. The code beneath is another method to automatically process the raster as a geojson instead of as a shapefile then re-importing, clipping to extent, and exporting, but I have found the first method does take less time.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "50Ipejl4AicO"
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal, ogr, osr\n",
    "import geopandas as gpd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_ytn710YAfck"
   },
   "outputs": [],
   "source": [
    "raster = gdal.Open(r\"C:/Users/jtrum/world_bank/data/IHME_LMIC_WASH_2000_2017_S_OD_PERCENT_MEAN_2017_Y2020M06D02.tif\")\n",
    "band = raster.GetRasterBand(1) # can use whichever band, but 1 works fine in most cases\n",
    "band.ReadAsArray()\n",
    "\n",
    "proj = raster.GetProjection()\n",
    "shp_proj = osr.SpatialReference()\n",
    "shp_proj.ImportFromWkt(proj)\n",
    "\n",
    "output_file = 'C:/Users/jtrum/Desktop/wb_outputs/sanitation.shp' # change to your output path\n",
    "call_drive = ogr.GetDriverByName('ESRI Shapefile') # exports as shp\n",
    "create_shp = call_drive.CreateDataSource(output_file)\n",
    "shp_layer = create_shp.CreateLayer('pct', srs=shp_proj) # name of layer in shp\n",
    "new_field = ogr.FieldDefn(str('ID'), ogr.OFTInteger)# will become the column name in geodataframe that raster values are extrapolated from\n",
    "shp_layer.CreateField(new_field)\n",
    "\n",
    "gdal.Polygonize(band, None, shp_layer, 0, [], callback=None)\n",
    "create_shp.Destroy()\n",
    "raster = None\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "z0L1AXOnBLl1"
   },
   "outputs": [],
   "source": [
    "# re-import the shapefile as a geodataframe\n",
    "gdf = gpd.read_file(r\"C:/Users/jtrum/Desktop/wb_outputs/sanitation.shp\")\n",
    "# # load in aoi as geodataframe\n",
    "aoi = gpd.read_file('C:/Users/jtrum/Desktop/wb_outputs/aoiLuanda.geojson')\n",
    "# ensure matching crs\n",
    "gdf.crs = aoi.crs\n",
    "# clip the shapefile to the aoi\n",
    "gdf = gpd.clip(gdf, aoi)\n",
    "# export to geojson\n",
    "gdf.to_file(r\"C:/Users/jtrum/Desktop/wb_outputs/sanitation.geojson\", driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lKHT_j9pBEGY"
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal, ogr, osr\n",
    "\n",
    "raster = gdal.Open(r\"C:/Users/jtrum/world_bank/data/IHME_LMIC_WASH_2000_2017_S_OD_PERCENT_MEAN_2017_Y2020M06D02.tif\")\n",
    "band = raster.GetRasterBand(1)\n",
    "band_array = band.ReadAsArray()\n",
    "\n",
    "proj = raster.GetProjection()\n",
    "shp_proj = osr.SpatialReference()\n",
    "shp_proj.ImportFromWkt(proj)\n",
    "\n",
    "output_file = 'C:/Users/jtrum/Desktop/wb_outputs/sanitation.geojson'\n",
    "driver = ogr.GetDriverByName('GeoJSON')\n",
    "create_geojson = driver.CreateDataSource(output_file)\n",
    "geojson_layer = create_geojson.CreateLayer('pct', srs=shp_proj, geom_type=ogr.wkbPolygon)\n",
    "\n",
    "new_field = ogr.FieldDefn(str('ID'), ogr.OFTInteger)\n",
    "geojson_layer.CreateField(new_field)\n",
    "\n",
    "gdal.Polygonize(band, None, geojson_layer, 0, [], callback=None)\n",
    "\n",
    "create_geojson.Destroy()\n",
    "raster = None\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyOvOmHMaAqN6DGW32A2wtjW",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}