{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#Drawing a Map from Pub Locations with the Matplotlib Basemap Toolkit\n",
    "\n",
    "**Author**: [Ramiro Gómez](http://ramiro.org/)\n",
    "\n",
    "In this notebook I show how you can draw a map of [Britain and Ireland](https://en.wikipedia.org/wiki/British_Isles) from location data using the [Matplotlib Basemap Toolkit](http://matplotlib.org/basemap/). The data points that will be drawn are pub locations extracted from [OpenStreetMap](http://www.openstreetmap.org/) and provided by [osm-x-tractor](http://osm-x-tractor.org/Data.aspx).\n",
    "\n",
    "When you download and extract the Points Of Interest (POI) dataset as a CSV file it has a file size of about 800 MB and more than 9 million entries for different types of locations. To filter out only the locations tagged as pubs you can use csvgrep, which is part of the [csvkit toolkit](http://csvkit.readthedocs.org/).\n",
    "\n",
    "    csvgrep -c amenity -m pub POIWorld.csv > POIWorld-pub.csv\n",
    "    \n",
    "The resulting CSV file is much smaller with a size of 7.7 MB and contains a little more than 120,000 entries for pub locations all over the world. Since the coverage varies across countries and regions, I chose to limit the map to Britain and Ireland where POI coverage seems quite comprehensive and there are a lot of pubs. Who could have thought?\n",
    "\n",
    "Next we load the required libraries and define a function that checks whether a given location tuple is within the given bounding box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "No module named 'mpl_toolkits.basemap'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-e3f914fb5dfc>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmpl_toolkits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbasemap\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mBasemap\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mImportError\u001b[0m: No module named 'mpl_toolkits.basemap'"
     ]
    }
   ],
   "source": [
    "# %load_ext signature\n",
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "\n",
    "def within_bbox(bbox, loc):\n",
    "    \"\"\"Determine whether given location is within given bounding box.\n",
    "    \n",
    "    Bounding box is a dict with ll_lon, ll_lat, ur_lon and ur_lat keys\n",
    "    that locate the lower left and upper right corners.\n",
    "    \n",
    "    The loction argument is a tuple of longitude and latitude values.\n",
    "    \"\"\"\n",
    "    \n",
    "    return bbox['ll_lon'] < loc[0] < bbox['ur_lon'] and bbox['ll_lat'] < loc[1] < bbox['ur_lat']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next statements load the pub dataset into a Pandas DataFrame, remove all columns that have missing values and print the length and the first few entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('csv/POIWorld-pub.csv')\n",
    "df.dropna(axis=1, how='any', inplace=True)\n",
    "print(len(df))\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restrict dataset to actual pubs\n",
    "\n",
    "Reddit user [kwwillett](https://www.reddit.com/user/kwwillett) pointed out that the dataset not only contains pubs [in a comment](https://www.reddit.com/r/MapPorn/comments/3erli3/britain_ireland_drawn_from_pubs_1144x1499_oc/ctibgig). The `csvgrep` call also matches strings containing it, see some examples below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.amenity.value_counts().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This changes number of pubs significantly so we need to perform additional filtering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = df[df.amenity.str.match(r'\\bpub\\b')]\n",
    "df.amenity.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The remaining location classifications are good to include, so now we have the actual pubs of the world as covered by OpenStreetMap in July 2015.\n",
    "\n",
    "## Limiting locations geographically\n",
    "\n",
    "To limit the dataset and the displayed map to Britain and Ireland, we now create a dict with appropriate coordinates from the [GeoPlanet Explorer](http://isithackday.com/geoplanet-explorer/index.php?woeid=24865730) and filter the longitude and latitude values from the DataFrame that are within the bounding box. We end up with roughly 28,000 records, meaning that almost a quarter of the pubs in the POI dataset are located on Britain and Ireland."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bbox = {\n",
    "    'lon': -5.23636,\n",
    "    'lat': 53.866772,\n",
    "    'll_lon': -10.65073,\n",
    "    'll_lat': 49.16209,\n",
    "    'ur_lon': 1.76334,\n",
    "    'ur_lat': 60.860699\n",
    "}\n",
    "\n",
    "locations = [loc for loc in zip(df.Longitude.values, df.Latitude.values) if within_bbox(bbox, loc)]\n",
    "len(locations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Drawing the map\n",
    "\n",
    "Now we get to the actual drawing part. We set the size of the image and some variables for [marker diplay](http://matplotlib.org/api/markers_api.html), create a `Basemap` object using the [Miller Cylindrical Projection](http://matplotlib.org/basemap/users/mill.html) and the coordinates looked up before to center the map and limit the area. Also we make sure that the background is white and no border is printed around the map image.\n",
    "\n",
    "The previously created list of location tuples is then unpacked into 2 lists of longitudes and latitudes, which are provided as arguments in the call to the Basemap instance. This converts the lon/lat values (in degrees) to x/y map projection coordinates (in meters). \n",
    "\n",
    "The location dots are drawn with the scatter method passing the x and y coordinate lists, and the marker display variables. The map is then annotated with a title that says what the map shows, a footer with information about the author and data source and finally displayed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(20, 30))\n",
    "markersize = 1\n",
    "markertype = ','  # pixel\n",
    "markercolor = '#325CA9'  # blue\n",
    "markeralpha = .8 #  a bit of transparency\n",
    "\n",
    "m = Basemap(\n",
    "    projection='mill', lon_0=bbox['lon'], lat_0=bbox['lat'],\n",
    "    llcrnrlon=bbox['ll_lon'], llcrnrlat=bbox['ll_lat'],\n",
    "    urcrnrlon=bbox['ur_lon'], urcrnrlat=bbox['ur_lat'])\n",
    "\n",
    "# Avoid border around map.\n",
    "m.drawmapboundary(fill_color='#ffffff', linewidth=.0)\n",
    "\n",
    "# Convert locations to x/y coordinates and plot them as dots.\n",
    "lons, lats = zip(*locations)\n",
    "x, y = m(lons, lats)\n",
    "m.scatter(x, y, markersize, marker=markertype, color=markercolor, alpha=markeralpha)\n",
    "\n",
    "# Set the map title.\n",
    "plt.annotate('Britain & Ireland\\ndrawn from pubs',\n",
    "             xy=(0, .87),\n",
    "             size=120, \n",
    "             xycoords='axes fraction',\n",
    "             color='#888888', \n",
    "             family='Gloria')\n",
    "\n",
    "# Set the map footer.\n",
    "plt.annotate('Author: Ramiro Gómez - ramiro.org • Data: OpenStreetMap - openstreetmap.org', \n",
    "             xy=(0, 0), \n",
    "             size=14, \n",
    "             xycoords='axes fraction',\n",
    "             color='#666666',\n",
    "             family='Droid Sans')\n",
    "\n",
    "plt.savefig('img/britain-ireland-drawn-from-pubs.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting map of Britain and Ireland shows the contour of the islands and unsurprisingly which areas are more or less densely populated. A large version of the map is available as prints from [Redbubble](http://www.redbubble.com/people/ramiro/works/15675223-britain-and-ireland-drawn-from-pubs-map-print) and [Zazzle](http://www.zazzle.com/britain_ireland_drawn_from_pubs_map_print-228420337846418460?rf=238355915198956003&tc=r)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%signature"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
