{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import geopandas\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case study - Conflict mapping: mining sites in eastern DR Congo\n",
    "\n",
    "In this case study, we will explore a dataset on artisanal mining sites located in eastern DR Congo.\n",
    "\n",
    "**Note**: this tutorial is meant as a hands-on session, and most code examples are provided as exercises to be filled in. I highly recommend actually trying to do this yourself, but if you want to follow the solved tutorial, you can find this in the `_solved` directory.\n",
    "\n",
    "---\n",
    "\n",
    "#### Background\n",
    "\n",
    "[IPIS](http://ipisresearch.be/), the International Peace Information Service, manages a database on mining site visits in eastern DR Congo: http://ipisresearch.be/home/conflict-mapping/maps/open-data/\n",
    "\n",
    "Since  2009, IPIS has visited artisanal mining sites in the region during various data collection campaigns. As part of these campaigns, surveyor teams visit mining sites in the field, meet with miners and complete predefined questionnaires. These contain questions about the mining site, the minerals mined at the site and the armed groups possibly present at the site.\n",
    "\n",
    "Some additional links:\n",
    "\n",
    "* Tutorial on the same data using R from IPIS (but without geospatial aspect): http://ipisresearch.be/home/conflict-mapping/maps/open-data/open-data-tutorial/\n",
    "* Interactive web app using the same data: http://www.ipisresearch.be/mapping/webmapping/drcongo/v5/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Importing and exploring the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The mining site visit data\n",
    "\n",
    "IPIS provides a WFS server to access the data. We can send a query to this server to download the data, and load the result into a geopandas GeoDataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import json\n",
    "\n",
    "wfs_url = \"http://geo.ipisresearch.be/geoserver/public/ows\"\n",
    "params = dict(service='WFS', version='1.0.0', request='GetFeature',\n",
    "              typeName='public:cod_mines_curated_all_opendata_p_ipis', outputFormat='json')\n",
    "\n",
    "r = requests.get(wfs_url, params=params)\n",
    "data_features = json.loads(r.content.decode('UTF-8'))\n",
    "data_visits = geopandas.GeoDataFrame.from_features(data_features, crs={'init': 'epsg:4326'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the data is also provided in the tutorial materials as a GeoJSON file, so it is certainly available during the tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Read the GeoJSON file `data/cod_mines_curated_all_opendata_p_ipis.geojson` using geopandas, and call the result `data_visits`.</li>\n",
    "  <li>Inspect the first 5 rows, and check the number of observations</li>\n",
    " </ul> \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping3.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping4.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping5.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The provided dataset contains a lot of information, much more than we are going to use in this tutorial. Therefore, we will select a subset of the column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_visits = data_visits[['vid', 'project', 'visit_date', 'name', 'pcode', 'workers_numb', 'interference', 'armed_group1', 'mineral1', 'geometry']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_visits.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before starting the actual geospatial tutorial, we will use some more advanced pandas queries to construct a subset of the data that we will use further on: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take only the data of visits by IPIS\n",
    "data_ipis = data_visits[data_visits['project'].str.contains('IPIS') & (data_visits['workers_numb'] > 0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For those mining sites that were visited multiple times, take only the last visit\n",
    "data_ipis_lastvisit = data_ipis.sort_values('visit_date').groupby('pcode', as_index=False).last()\n",
    "data = geopandas.GeoDataFrame(data_ipis_lastvisit, crs=data_visits.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data on protected areas in the same region\n",
    "\n",
    "Next to the mining site data, we are also going to use a dataset on protected areas (national parks) in Congo. This dataset was downloaded from http://www.wri.org/our-work/project/congo-basin-forests/democratic-republic-congo#project-tabs and included in the tutorial repository: `data/cod_conservation.zip`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Extract the `data/cod_conservation.zip` archive, and read the shapefile contained in it. Assign the resulting GeoDataFrame to a variable named `protected_areas`.</li>\n",
    "  <li>Quickly plot the GeoDataFrame.</li>\n",
    " </ul> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping10.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping11.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conversion to a common Coordinate Reference System\n",
    "\n",
    "We will see that both datasets use a different Coordinate Reference System (CRS). For many operations, however, it is important that we use a consistent CRS, and therefore we will convert both to a commong CRS.\n",
    "\n",
    "But first, we explore problems we can encounter related to CRSs.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Goma](https://en.wikipedia.org/wiki/Goma) is the capital city of North Kivu province of Congo, close to the border with Rwanda. It's coordinates are 1.66°S 29.22°E.\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Create a single Point object representing the location of Goma. Call this `goma`.</li>\n",
    "  <li>Calculate the distances of all mines to Goma, and show the 5 smallest distances (mines closest to Goma).</li>\n",
    " </ul> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping12.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping13.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping14.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping15.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distances we see here in degrees, which is not helpful for interpreting those distances. That is a reason we will convert the data to another coordinate reference system (CRS) for the remainder of this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Make a visualization of the national parks and the mining sites on a single plot.</li>\n",
    " </ul> \n",
    " \n",
    "<p>Check the first section of the [04-more-on-visualization.ipynb](04-more-on-visualization.ipynb) notebook for tips and tricks to plot with GeoPandas.</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping16.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will notice that the protected areas and mining sites do not map to the same area on the plot. This is because the Coordinate Reference Systems (CRS) differ for both datasets. Another reason we will need to convert the CRS!\n",
    "\n",
    "Let's check the Coordinate Reference System (CRS) for both datasets.\n",
    "\n",
    "The mining sites data uses the [WGS 84 lat/lon (EPSG 4326)](http://spatialreference.org/ref/epsg/4326/) CRS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The protected areas dataset, on the other hand, uses a [WGS 84 / World Mercator (EPSG 3395)](http://spatialreference.org/ref/epsg/wgs-84-world-mercator/) projection (with meters as unit):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protected_areas.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will convert both datasets to a local UTM zone, so we can plot them together and that distance-based calculations give sensible results.\n",
    "\n",
    "To find the appropriate UTM zone, you can check http://www.dmap.co.uk/utmworld.htm or https://www.latlong.net/lat-long-utm.html, and in this case we will use UTM zone 35, which gives use EPSG 32735: https://epsg.io/32735\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Convert both datasets (`data` and `protected_areas`) to EPSG 32735. Name the results `data_utm` and `protected_areas_utm`.</li>\n",
    "  <li>Try again to visualize both datasets on a single map.</li>\n",
    " </ul> \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping19.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping20.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More advanced visualizations\n",
    "\n",
    "<p>For the following exercises, check the first section of the [04-more-on-visualization.ipynb](04-more-on-visualization.ipynb) notebook for tips and tricks to plot with GeoPandas.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Make a visualization of the national parks and the mining sites on a single plot.</li>\n",
    "  <li>Pay attention to the following details:\n",
    "     <ul>\n",
    "      <li>Make the figure a bit bigger.</li>\n",
    "      <li>The protected areas should be plotted in green</li>\n",
    "      <li>For plotting the mining sites, adjust the markersize and use an `alpha=0.5`.</li>\n",
    "      <li>Remove the figure border and x and y labels (coordinates)</li>\n",
    "     </ul> \n",
    "   </li>\n",
    " </ul> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping21.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping22.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " \n",
    " In addition to the previous figure:\n",
    " <ul>\n",
    "  <li>Give the mining sites a distinct color based on the `'interference'` column, indicating whether an armed group is present at the mining site or not.</li>\n",
    " </ul> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping23.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " \n",
    " In addition to the previous figure:\n",
    " <ul>\n",
    "  <li>Give the mining sites a distinct color based on the `'mineral1'` column, indicating which mineral is the primary mined mineral.</li>\n",
    " </ul> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping24.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Spatial operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " \n",
    " <ul>\n",
    "  <li>Access the geometry of the \"Kahuzi-Biega National park\".</li>\n",
    "  <li>Filter the mining sites to select those that are located in this national park.</li>\n",
    " </ul> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping25.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping26.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping27.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>: Determine for each mining site the \"closest\" protected area:\n",
    " \n",
    " <ul>\n",
    "  <li> PART 1 - do this for a single mining site:\n",
    "   <ul>\n",
    "    <li>Get a single mining site, e.g. the first of the dataset.</li>\n",
    "    <li>Calculate the distance (in km's) to all protected areas for this mining site</li>\n",
    "    <li>Get the index of the minimum distance (tip: `idxmin()`) and get the name of the protected are corresponding to this index.</li>\n",
    "   </ul> \n",
    "  </li>\n",
    "  <li> PART 2 - apply this procedure on each geometry:\n",
    "   <ul>\n",
    "    <li>Write the above procedure as a function that gets a single site and the protected areas dataframe as input and returns the name of the closest protected area as output.</li>\n",
    "    <li>Apply this function to all sites using the `.apply()` method on `data_utm.geometry`.</li>\n",
    "   </ul> \n",
    "  </li>\n",
    " </ul> \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping28.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping29.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping30.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping31.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping32.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Using spatial join to determine mining sites in the protected areas\n",
    "\n",
    "Based on the analysis and visualizations above, we can already see that there are mining sites inside the protected areas. Let's now do an actual spatial join to determine which sites are within the protected areas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mining sites in protected areas\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Add information about the protected areas to the mining sites dataset, using a spatial join:\n",
    "   <ul>\n",
    "    <li>Call the result `data_within_protected`</li>\n",
    "    <li>If the result is empty, this is an indication that the coordinate reference system is not matching. Make sure to re-project the data (see above).</li>\n",
    "      \n",
    "   </ul>\n",
    "  </li>\n",
    "  <li>How many mining sites are located within a national park?</li>\n",
    "  <li>Count the number of mining sites per national park (pandas tip: check `value_counts()`)</li>\n",
    "\n",
    " </ul> \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping33.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping34.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping35.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping36.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mining sites in the borders of protected areas\n",
    "\n",
    "And what about the borders of the protected areas? (just outside the park)\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    " <b>EXERCISE</b>:\n",
    " <ul>\n",
    "  <li>Create a new dataset, `protected_areas_borders`, that contains the border area (10 km wide) of each protected area:\n",
    "   <ul>\n",
    "    <li>Tip: one way of doing this is with the `buffer` and `difference` function.</li>\n",
    "    <li>Plot the resulting borders as a visual check of correctness.</li>\n",
    "   </ul>\n",
    "  </li>\n",
    "  <li>Count the number of mining sites per national park that are located within its borders</li>\n",
    "\n",
    " </ul> \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping37.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping38.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping39.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping40.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "clear_cell": true
   },
   "outputs": [],
   "source": [
    "# %load _solved/solutions/case-conflict-mapping41.py"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}