notebooks/05-the-power-of-gdal.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import ogr\n",
    "import subprocess\n",
    "\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GDAL command line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: reprojection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "**GDAL** is a really powerful library for handling GIS data. It provides a number of functionalities to interact with spatial data. As a typical example, take the **reprojection** of a shapefile to another CRS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr ../scratch/deelbekkens_wgs84 -t_srs \"EPSG:4326\" ../data/deelbekkens/Deelbekken.shp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is this combination of commands?\n",
    "\n",
    "* `!` this is a jupyter notebook-thing, telling it we're running something on the command line instead of in Python\n",
    "* `../scratch/deelbekkens_wgs84` the output location of the created file\n",
    "* `-t_srs \"EPSG:4326\"` the CRS information for to which the data should be projected\n",
    "* `../data/deelbekkens/Deelbekken.shp` the source file location"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The documentation is a bit overwhelming:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but there are [great online resources](https://github.com/dwtkns/gdal-cheat-sheet) with good examples you can easily copy paste for your own applications..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Accessing online webservice data (Web Feature Service - WFS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A lot of expensive terminology... \n",
    "\n",
    "Let's illustrate this with an example: The information about municipalities is available as [open data on geopunt](http://www.geopunt.be/catalogus/webservicefolder/1/11c37274-f9db-526e-2067-6606-b324-23d1-f285dbe7) (coming from informatie Vlaanderen). The publication is provided as a [WFS service](http://docs.geoserver.org/stable/en/user/services/wfs/reference.html)... \n",
    "\n",
    "Take home message -> **GDAL can handle WFS web services** ;-)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Downloading the province boundaries from the WFS service provided by *informatie Vlaanderen/Geopunt* to a geojson file is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr -f 'Geojson' ../scratch/provinces.geojson WFS:\"https://geoservices.informatievlaanderen.be/overdrachtdiensten/VRBG/wfs\" Refprv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can start working with this date..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "provinces = gpd.read_file(\"../scratch/provinces.geojson\")\n",
    "provinces.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, GDAL can directly **query the WFS data**:\n",
    "\n",
    "Let's say I only need the province of `Antwerp`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr -f 'Geojson' ../scratch/antwerp_prov.geojson WFS:\"https://geoservices.informatievlaanderen.be/overdrachtdiensten/VRBG/wfs\" Refprv -where \"NAAM = 'Antwerpen'\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "antwerp = gpd.read_file(\"../scratch/antwerp_prov.geojson\")\n",
    "antwerp.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I do know that the [*Meetplaatsen Oppervlaktewaterkwaliteit*](http://www.geopunt.be/catalogus/webservicefolder/4/435c1aae-6619-0801-c317-2bf9-9f58-ce7f-e262a96c)  are also available as a WFS web service. However, I'm only interested in the locations for *fytoplankton*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr -f 'Geojson' ../scratch/metingen_fytoplankton.geojson WFS:\"https://geoservices.informatievlaanderen.be/overdrachtdiensten/MeetplOppervlwaterkwal/wfs\" Mtploppw -where \"FYTOPLANKT = '1'\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import mplleaflet\n",
    "fyto = gpd.read_file(\"../scratch/metingen_fytoplankton.geojson\")\n",
    "fyto.head()\n",
    "fyto.to_crs('+init=epsg:4326').plot(markersize=5)\n",
    "mplleaflet.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "fyto.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, the same type of subselection are also possible on shapefiles,... \n",
    "\n",
    "Extracting a specific *DEELBEKKEN* from the deelbekken shapefile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr ../scratch/subcat.shp ../data/deelbekkens/Deelbekken.shp -where \"DEELBEKKEN = '10-10'\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(If you're wondering how I know how to setup these commands and arguments, check the (draft) `introduction_webservices.ipynb` in the `scratch` folder. <br>\n",
    "I use the python interafce of GDAL/OGR and the package [`owslib`](https://geopython.github.io/OWSLib/) to find out how to setup the arguments.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GDAL command line, but inside Python..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No problem if this is still unclear... an example application!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clipping example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example we will use is to **clip raster data using a shapefile**. We use a data set from [natural earth](http://www.naturalearthdata.com/downloads/50m-natural-earth-1/50m-natural-earth-i-with-shaded-relief/), which we will unzip to start working on it (off course using Python itself):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import zipfile\n",
    "zip_ref = zipfile.ZipFile(\"../data/NE1_50m_SR.zip\", 'r')\n",
    "zip_ref.extractall(\"../scratch\")\n",
    "zip_ref.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **GDAL** function that support the clipping of a raster file, is called [`gdalwarp`](http://www.gdal.org/gdalwarp.html). Again, the documentation looks rather overwhelming... Let's start with an example execution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!gdalwarp ../scratch/NE1_50M_SR/NE1_50M_SR.tif ../scratch/cliptest.tif -cutline \"../scratch/subcat.shp\" -crop_to_cutline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* `!` this is a jupyter notebook-thing, telling it we're running something on the command line instead of in Python\n",
    "* `gdalwarp` is the GDAL command to use\n",
    "* `../scratch/NE1_50M_SR/NE1_50M_SR.tif` the source file location\n",
    "* `../scratch/cliptest.tif`  the output location of the created file\n",
    "* `-cutline \"../scratch/subcat.shp\"` the shape file to cut the raster with\n",
    "* `-crop_to_cutline` an additional argument to GDAL to make the clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "%matplotlib inline\n",
    "img=mpimg.imread('../scratch/cliptest.tif')\n",
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is off course a dummy example (to keep runtime low), but it illustrates the concept."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## the `subprocess` trick..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doing the same using pure Python code and `from osgeo import gdal` is actually not so beneficial, as the command above is reather straighforward... However, the dependency of the command line provides a switch of environment in any data analysis pipeline. I actually do want to have the **best of both worlds**:\n",
    "Using Python code, but running the command line version of GDAL...\n",
    "\n",
    "...therefore we need `subprocess`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import subprocess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doing the same as above, but actually using Python code to run the command with given variables as input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "inraster = '../scratch/NE1_50M_SR/NE1_50M_SR.tif'\n",
    "outraster = inraster.replace('.tif', '{}.tif'.format(\"_out\")) # same location, but adding _out to the output \n",
    "inshape = \"../scratch/subcat.shp\"\n",
    "subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, \n",
    "                     '-crop_to_cutline', '-overwrite'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remark** when GDAL provides a zero as return statement, this is a GOOD sign!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "%matplotlib inline\n",
    "img=mpimg.imread('../scratch/NE1_50M_SR/NE1_50M_SR_out.tif')\n",
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, the result is the same, but calling the command from Python. By **writing a Python function** for this routine, I do have a reusable functionality in my toolbox that I can load in any other Python script:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def clip_raster(inraster, outraster, invector):\n",
    "    \"\"\"clip a raster image with a vector file\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    inraster : GDAL compatible raster format\n",
    "    outraster : GDAL compatible raster format\n",
    "    invector : GDAL compatible vector format\n",
    "    \"\"\"\n",
    "    response = subprocess.call(['gdalwarp', inraster, outraster, '-cutline', \n",
    "                                invector, '-crop_to_cutline', '-overwrite'])\n",
    "    return(response)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "inraster = '../scratch/NE1_50M_SR/NE1_50M_SR.tif'\n",
    "outraster = inraster.replace('.tif', '{}.tif'.format(\"_out\")) # same location, but adding _out to the output \n",
    "inshape = \"../scratch/subcat.shp\"\n",
    "clip_raster(inraster, outraster, inshape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More advanced clipping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the data set of the provinces we called from the WFS server earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "provinces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can actually use a selection of the provinces data set to execute the clipping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "inraster = '../scratch/NE1_50M_SR/NE1_50M_SR.tif'\n",
    "outraster = inraster.replace('.tif', '{}.tif'.format(\"_OostVlaanderen\"))    \n",
    "invector = \"../scratch/provinces.geojson\"\n",
    "subprocess.call(['gdalwarp', inraster, outraster, '-cutline', invector, \n",
    "                 '-cwhere', \"NAAM='OOST-VLAANDEREN'\", \n",
    "                 '-crop_to_cutline', \n",
    "                 '-overwrite'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "%matplotlib inline\n",
    "img=mpimg.imread('../scratch/NE1_50M_SR/NE1_50M_SR_OostVlaanderen.tif')\n",
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By having it as a Python call, we can do the same action for each of the individual provinces in the dataset and **create for each of the provinces a clipped raster data set**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import ogr\n",
    "\n",
    "inraster = '../scratch/NE1_50M_SR/NE1_50M_SR.tif'\n",
    "invector = \"../scratch/provinces.geojson\"\n",
    "\n",
    "# GDAL magic...\n",
    "ds = ogr.Open(inshape)\n",
    "lyr = ds.GetLayer(0)\n",
    "\n",
    "lyr.ResetReading()\n",
    "ft = lyr.GetNextFeature()\n",
    "\n",
    "# clipping for each of the features (provincesin this case)\n",
    "while ft:\n",
    "\n",
    "    province_name = ft.GetFieldAsString('NAAM')\n",
    "    print(province_name)\n",
    "\n",
    "    outraster = inraster.replace('.tif', '_%s.tif' % province_name.replace('-', '_'))    \n",
    "    subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, \n",
    "                     '-crop_to_cutline', '-cwhere', \"NAAM='%s'\" %province_name])\n",
    "\n",
    "    ft = lyr.GetNextFeature()\n",
    "\n",
    "ds = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "%matplotlib inline\n",
    "img=mpimg.imread('../scratch/NE1_50M_SR/NE1_50M_SR_West_Vlaanderen.tif') # check also Antwerpen,...\n",
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.3"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}