notebooks/04-gis-python-vectors.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import fiona\n",
    "import geopandas as gpd\n",
    "from osgeo import gdal\n",
    "from osgeo import ogr\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "## Fiona: reading vector data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Fiona is a minimalist python package for reading (and writing) vector\n",
    "data in python.  Fiona provides python objects (e.g. a dictionary for\n",
    "each record) to geospatial data in various formats.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading in vector files with Fiona is - *good practice* - done using the `with` environment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "with fiona.open('../data/deelbekkens/Deelbekken.shp') as deelbekkens:\n",
    "    feature = next(iter(deelbekkens)) # Just one checking the first\n",
    "    print(\"Bekken: \", feature['properties']['BEKNAAM'])\n",
    "    print(\"Vectortype: \", feature['geometry']['type'])\n",
    "    print(feature['geometry']['coordinates'][0][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is actually equivalent to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "deelbekkens = fiona.open('../data/deelbekkens/Deelbekken.shp') \n",
    "feature = next(iter(deelbekkens))\n",
    "print(\"Bekken: \", feature['properties']['BEKNAAM'])\n",
    "print(\"Stroomgebied: \", feature['properties']['STRMGEB'])\n",
    "print(\"Vectortype: \", feature['geometry']['type'])\n",
    "deelbekkens.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But people tend to forget to use the `close`(!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>REMEMBER</b>: \n",
    "    <ul>\n",
    "        <li> Instead of `open`-ing a file, do something and `close` again, use `with` (general Python advice)\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The projection information is interpreted by Fiona as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "with fiona.open('../data/deelbekkens/Deelbekken.shp') as deelbekkens:\n",
    "    #print(deelbekkens.crs_wkt)\n",
    "    print(deelbekkens.crs)\n",
    "    print(deelbekkens.crs_wkt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    <p>Actually, **GeoPandas relies on Fiona** and translates this information into a tabular information, compatible to Pandas:</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "df_deelbekkens = gpd.read_file('../data/deelbekkens/Deelbekken.shp')\n",
    "df_deelbekkens.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The capabilities to read file with `Fiona` is dependent on the available `drivers` on your computer. You can get an overview by asking fiona for the supported drivers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "fiona.supported_drivers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It provides also information about what you can do with the different types of data:\n",
    "* `r`: reading\n",
    "* `a`: appending (adding elements)\n",
    "* `w`: writing "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>: \n",
    "<ul>\n",
    "    <li>Open the `EUgrid10.geojson` file in the data folder using the Fiona package and check the `type` of `geometry` and the `cellcode` in the `properties` of the first feature in the dataset</li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# %load ../notebooks/_solutions/04-gis-python-vectors7.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>: \n",
    "<ul>\n",
    "    <li>Read again the EUgrid10.geojson file in the data folder, but use the GeoPandas library and show only the first feature</li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# %load ../notebooks/_solutions/04-gis-python-vectors8.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, what about GDAL?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "driver = ogr.GetDriverByName('ESRI Shapefile')\n",
    "deelbekken_gdal = driver.Open('../data/deelbekkens/Deelbekken.shp')\n",
    "deelbekken_gdal.GetLayerCount()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "deelbekken_layer = deelbekken_gdal.GetLayer()\n",
    "print(\"The number of features of the deelbekken shapefile is: \", deelbekken_layer.GetFeatureCount())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's why they call it `not pythonic`...\n",
    "\n",
    "However, `GDAL` is powerfull for other activities, we'll come back to that later. Moreover, it supports mostly more file formats as fiona (GeoPandas) provide:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr --formats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, the interaction with an Esri File GeoDataBase, is provided in [this](https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#get-all-layers-in-an-esri-file-geodatabase) example..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>REMEMBER</b>: \n",
    "    <ul>\n",
    "        <li> try `GeoPandas` first!\n",
    "        <li> If that does not work, check if your available GDAL drivers support the format.\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shapely: working with GEOMetries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Shapely is a Python library for geometric operations using the **GEOS** library.</p>\n",
    "\n",
    "<p>Shapely can perform:</p>\n",
    "<ul>\n",
    "<li>geometry validation</li>\n",
    "<li>geometry creation (e.g. collections)</li>\n",
    "<li>geometry operations</li>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import LineString, Polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "line = LineString([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])\n",
    "line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "dilated = line.buffer(0.5, cap_style=3)\n",
    "dilated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting back the coorindates as x, y arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "polygon = Polygon(LineString([(0.2, 0.2), (0.2, 1), (1, 1), (1, 0.2)]))\n",
    "polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "line.coords.xy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binary predicates to check if the line is contained by the polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "line.within(polygon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>: \n",
    "<ul>\n",
    "    <li>Check if the polygon intersects with the line object</li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# %load ../notebooks/_solutions/04-gis-python-vectors29.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    <p>Actually, **GeoPandas relies on shapely** and uses these geometric properties and operations</p>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "deelbekkens = gpd.read_file(\"../data/deelbekkens/Deelbekken.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "demer = deelbekkens[deelbekkens[\"BEKNAAM\"] == \"Demerbekken\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "demer.plot()\n",
    "demer.convex_hull.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## pyproj: handling projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pyproj handles the Spatial Reference Systems (SRS) transformations. Actually, the only items to remember are:\n",
    "* `pyproj.Proj`: define a SRS\n",
    "* `pyproj.transform`: execute a transformation between two defined SRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from pyproj import Proj, transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define SRS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The easiest way is mostly just using the [`EPSG` code](http://spatialreference.org/):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lambert 72:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "srs_lambert = Proj(init='epsg:31370')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "srs_lambert.srs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "WGS84:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "srs_wgs84 = Proj(init='epsg:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "srs_wgs84.is_latlong()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other methods to define the SRS are also available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "utm32 = Proj(proj=\"utm\", zone=\"32\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transform coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transform function takes the two defined SRS objects, together with an X/Y coordinate combination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "transform(srs_lambert, srs_wgs84, 98710.32800000161, 162573.7030000016)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    <p>Actually, **GeoPandas relies on pyproj** to perform SRS transformations using the EPSG code (but calls it `coordinate reference system (CRS)`...)</p>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "deelbekkens = gpd.read_file(\"../data/deelbekkens/Deelbekken.shp\")\n",
    "deelbekkens.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "deelbekkens.to_crs(epsg=\"4326\").head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, the function is useful in general to transform the X/Y combinations in any coordinate table. Consider the following example of species data with the coordinates in WGS84."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "inv_data = pd.read_excel(\"../data/invasive_extract.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "inv_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conversion to Lambert 72 is supported by the `pyproj` Package as follows. As I want to support future usage of this functionality as well, I'll write my own custom function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def transform_dfrow_projection(row, x_name, y_name, srs_from, srs_to):\n",
    "    \"\"\"\n",
    "    Converts the x and y coordinates of a given DataFrame row into a Series of the\n",
    "    longitude and latitude.\n",
    "    \"\"\"\n",
    "    from pyproj import transform\n",
    "    return pd.Series(transform(srs_from, srs_to, row[x_name], row[y_name]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a single row, this looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "transform_dfrow_projection(inv_data.loc[0], \"decimalLongitude\", \"decimalLatitude\", \n",
    "                           srs_wgs84, srs_lambert)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Apply`ing this to the entire table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "inv_data[[\"x\", \"y\"]] = inv_data.apply(transform_dfrow_projection, axis=1, \n",
    "                                      args=(\"decimalLongitude\", \"decimalLatitude\", \n",
    "                                            srs_wgs84, srs_lambert))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This feels too complicated?!?Understandable ,but let's just once take the effort to compile this into a reusable function..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def add_trf_coordinates(filename, x_name=\"decimalLongitude\", y_name=\"decimalLatitude\", \n",
    "                        srs_from='epsg:4326', srs_to='epsg:31370', \n",
    "                        new_x=\"X\", new_y=\"Y\"):\n",
    "    \"\"\"Add new X/Y columns with transformed SRS to a given Excel dataset with X/Y columns\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    filename : str\n",
    "        Path to excel file to read dat from\n",
    "    x_name : str\n",
    "        Column containing the x-coordinates (default decimalLongitude)\n",
    "    y_name : str\n",
    "        Column containing the y-coordinates (default decimalLatitude)  \n",
    "    srs_from: str\n",
    "        epsg:n code for the current SRS of the coordinates\n",
    "    srs_to: str\n",
    "        epsg:n code for the SRS of the added coordinates\n",
    "    new_x: str\n",
    "        name of the added x-coordinate column\n",
    "    new_y: str\n",
    "        name of the added y-coordinate column    \n",
    "    \n",
    "    returns\n",
    "    -------\n",
    "    df with 2 additional columns\n",
    "    \"\"\"\n",
    "    df = pd.read_excel(filename)\n",
    "    srs_from = Proj(init=srs_from)\n",
    "    srs_to = Proj(init=srs_to)\n",
    "    df[[new_x, new_y]] = df.apply(transform_dfrow_projection, axis=1, \n",
    "                                  args=(x_name, y_name, srs_from, srs_to))\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "df_enriched = add_trf_coordinates(\"../data/invasive_extract.xlsx\")\n",
    "df_enriched.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remark** that GDAL also has a utility [gdalsrsinfo](http://www.gdal.org/gdalsrsinfo.html) to get/check the SRS information of GIS file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "! gdalsrsinfo -o proj4 ../data/EUgrid10.geojson  #try also wiht -o all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or to execute the transformation of a dataset together with the transformation of the data type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "!ogr2ogr -t_srs \"EPSG:4326\" -f \"GeoJSON\" \"../scratch/deelbekken.geojson\" \"../data/deelbekkens/Deelbekken.shp\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(More on this in the next notebook)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geopandas: Pandas + geometries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, we got a whole sequence of data types (objects) to end up with GeoDataFrames:\n",
    "* lists/dictionaries -> standard Python library \n",
    "* Numpy arrays -> table of elements (usually numbers), all of the same type, indexed by positive integers\n",
    "* DataFrames -> Labeled Numpy arrays\n",
    "* GeoDataFrames -> DataFrames with geometry and SRS information attached to it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `df_enriched` is a Pandas DataFrame, so we can convert to a GeoDataFrame by converting the coordinates to a set of POINT geometries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import Point\n",
    "\n",
    "geometry = [Point(xy) for xy in zip(df_enriched[\"X\"], df_enriched[\"Y\"])]\n",
    "crs = {'init': 'epsg:31370'}\n",
    "geo_df = gpd.GeoDataFrame(df_enriched, crs=deelbekkens.crs, geometry=geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import mplleaflet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "geo_df.to_crs(epsg=\"4326\").plot(markersize=10, ax=ax) # on the fly conversion to WGS84\n",
    "mplleaflet.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>: \n",
    "<ul>\n",
    "    <li>Add the UIDN and BEKNAAM names from the `deelbekkens` data to the occurrence data `df_enriched` with a **spatial JOIN**\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# %load ../notebooks/_solutions/04-gis-python-vectors208.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>: \n",
    "<ul>\n",
    "    <li>Plot the `geo_df` together with a buffered version of 500m of each of the individual points, use mplleaflet to make the plot interactive\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# %load ../notebooks/_solutions/04-gis-python-vectors215.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For those interested, also check the following example from Joris Vandenbossche, explaining how to derive the distance to the nearest green area for all adresses in the centre of Antwerp:\n",
    "* https://gist.github.com/jorisvandenbossche/7b30ed43366a85af8626"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Nbtutor - export exercises",
  "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
  },
  "toc_position": {
   "height": "863px",
   "left": "0px",
   "right": "1663px",
   "top": "106px",
   "width": "212px"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}