{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), the definitive guide to high-performance scientific computing and data science in Python.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 14.6. Manipulating geospatial data with Shapely and basemap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to run this recipe, you will need the following packages:\n",
    "\n",
    "* [GDAL/OGR](http://www.gdal.org/ogr/)\n",
    "* [fiona](http://toblerity.org/fiona/README.html)\n",
    "* [Shapely](http://toblerity.org/shapely/project.html)\n",
    "* [descartes](https://pypi.python.org/pypi/descartes)\n",
    "* [basemap](http://matplotlib.org/basemap/)\n",
    "\n",
    "On Windows, you can find binary installers for all those packages except descartes on Chris Gohlke's webpage. (http://www.lfd.uci.edu/~gohlke/pythonlibs/)\n",
    "\n",
    "Installing descartes is easy: `pip install descartes`.\n",
    "\n",
    "On other systems, you can find installation instructions on the projects' websites. GDAL/OGR is a C++ library that is required by fiona. The other packages are regular Python packages.\n",
    "\n",
    "Finally, you need to download the *Africa* dataset on the book's website. (http://ipython-books.github.io)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Let's import all the packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.collections as col\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import fiona\n",
    "import shapely.geometry as geom\n",
    "from descartes import PolygonPatch\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. We load the Shapefile dataset with fiona. This dataset notably contains the borders of all countries in the world."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# natural earth data\n",
    "countries = fiona.open(\"data/ne_10m_admin_0_countries.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. We select the countries in Africa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "africa = [c for c in countries \n",
    "          if c['properties']['CONTINENT'] == 'Africa']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Now, we create a basemap map showing the African continent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = Basemap(llcrnrlon=-23.03,\n",
    "            llcrnrlat=-37.72,\n",
    "            urcrnrlon=55.20,\n",
    "            urcrnrlat=40.58)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. We need to convert the geographical coordinates of the countries' borders in map coordinates, so that we can display then in basemap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def _convert(poly, m):\n",
    "    if isinstance(poly, list):\n",
    "        return [_convert(_, m) for _ in poly]\n",
    "    elif isinstance(poly, tuple):\n",
    "        return m(*poly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for _ in africa:\n",
    "    _['geometry']['coordinates'] = _convert(\n",
    "        _['geometry']['coordinates'], m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. The next step is to create matplotlib `PatchCollection` objects from the Shapefile dataset loaded with fiona. We use Shapely and descartes to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_patch(shape, **kwargs):\n",
    "    \"\"\"Return a matplotlib PatchCollection from a geometry \n",
    "    object loaded with fiona.\"\"\"\n",
    "    # Simple polygon.\n",
    "    if isinstance(shape, geom.Polygon):\n",
    "        return col.PatchCollection([PolygonPatch(shape, **kwargs)],\n",
    "                                   match_original=True)\n",
    "    # Collection of polygons.\n",
    "    elif isinstance(shape, geom.MultiPolygon):\n",
    "        return col.PatchCollection([PolygonPatch(c, **kwargs)\n",
    "                                    for c in shape],\n",
    "                                   match_original=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_patches(shapes, fc=None, ec=None, **kwargs):\n",
    "    \"\"\"Return a list of matplotlib PatchCollection objects\n",
    "    from a Shapefile dataset.\"\"\"\n",
    "    # fc and ec are the face and edge colors of the countries.\n",
    "    # We ensure these are lists of colors, with one element\n",
    "    # per country.\n",
    "    if not isinstance(fc, list):\n",
    "        fc = [fc] * len(shapes)\n",
    "    if not isinstance(ec, list):\n",
    "        ec = [ec] * len(shapes)\n",
    "    # We convert each polygon to a matplotlib PatchCollection\n",
    "    # object.\n",
    "    return [get_patch(geom.shape(s['geometry']), \n",
    "                      fc=fc_, ec=ec_, **kwargs) \n",
    "            for s, fc_, ec_ in zip(shapes, fc, ec)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. We also define a function to get countries colors depending on a specific field in the Shapefile dataset. Indeed, our dataset not only contains countries borders, but also a few administrative, economical and geographical properties for each country. Here, we will choose the color according to the population and GDP of the countries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_colors(field, cmap):\n",
    "    \"\"\"Return one color per country, depending on a specific\n",
    "    field in the dataset.\"\"\"\n",
    "    values = [country['properties'][field] for country in africa]\n",
    "    values_max = max(values)\n",
    "    return [cmap(v / values_max) for v in values]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "8. Finally, we display the maps. We display the coast lines with basemap, and the countries with our Shapefile dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6));\n",
    "# Display the countries color-coded with their population.\n",
    "ax = plt.subplot(121);\n",
    "m.drawcoastlines();\n",
    "patches = get_patches(africa, \n",
    "                      fc=get_colors('POP_EST', plt.cm.Reds), \n",
    "                      ec='k')\n",
    "for p in patches:\n",
    "    ax.add_collection(p)\n",
    "plt.title(\"Population\");\n",
    "# Display the countries color-coded with their population.\n",
    "ax = plt.subplot(122);\n",
    "m.drawcoastlines();\n",
    "patches = get_patches(africa, \n",
    "                      fc=get_colors('GDP_MD_EST', plt.cm.Blues), \n",
    "                      ec='k')\n",
    "for p in patches:\n",
    "    ax.add_collection(p)\n",
    "plt.title(\"GDP\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).\n",
    "\n",
    "> [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages)."
   ]
  }
 ],
 "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.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}