convert_folium.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:dc8cf08d2051ebb2e6cae1082e8bf840c4ffeddc738fb94addb86b8a7d257ecf"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import HTML\n",
      "from lxml import etree\n",
      "import requests\n",
      "import pandas as pd\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import matplotlib.cm as cm\n",
      "from matplotlib.colors import Normalize, LinearSegmentedColormap\n",
      "from matplotlib.collections import PatchCollection\n",
      "import matplotlib.font_manager as fm\n",
      "from mpl_toolkits.basemap import Basemap\n",
      "from shapely.geometry import Point, Polygon, MultiPoint, MultiPolygon, shape, mapping\n",
      "from shapely.prepared import prep\n",
      "from pysal.esda.mapclassify import Natural_Breaks as nb\n",
      "from descartes import PolygonPatch\n",
      "import mpl_toolkits.basemap.pyproj as pyproj\n",
      "import fiona\n",
      "from itertools import chain"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# folium stuff\n",
      "import json\n",
      "import folium"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# ensure we have Helvetica available\n",
      "fm.findSystemFonts(fontpaths='/users/sth/fonts', fontext='ttf')\n",
      "prop = fm.FontProperties(fname='/users/sth/fonts/helvetica.ttf')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Parse an XML file into a dict"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# parse XML into dict\n",
      "tree = etree.parse(\"data/plaques_20140619.xml\")\n",
      "root = tree.getroot()\n",
      "\n",
      "output = dict()\n",
      "output['raw'] = []\n",
      "output['crs'] = []\n",
      "output['lon'] = []\n",
      "output['lat'] = []\n",
      "\n",
      "for each in root.xpath('/openplaques/plaque/geo'):\n",
      "    # check what we got back\n",
      "    output['crs'].append(each.get('reference_system', None))\n",
      "    output['lon'].append(each.get('longitude', None))\n",
      "    output['lat'].append(each.get('latitude', None))\n",
      "    # now go back up to plaque\n",
      "    r = each.getparent().xpath('inscription/raw')[0]\n",
      "    if isinstance(r.text, str):\n",
      "        output['raw'].append(r.text.lstrip().rstrip())\n",
      "    else:\n",
      "        output['raw'].append(None)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Set up a DataFrame, and perform type conversion(s)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df = pd.DataFrame(output)\n",
      "df = df.replace({'raw': 0}, None)\n",
      "df = df.dropna()\n",
      "df[['lon', 'lat']] = df[['lon', 'lat']].astype(float)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Get GeoJSON data, and create bounds, width and height"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "shp = fiona.open('data/wards.geojson')\n",
      "crs_data = shp.crs\n",
      "bds = shp.bounds\n",
      "shp.close()\n",
      "extra = 0.01\n",
      "ll = (bds[0], bds[1])\n",
      "ur = (bds[2], bds[3])\n",
      "coords = list(chain(ll, ur))\n",
      "w, h = coords[2] - coords[0], coords[3] - coords[1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Create a Basemap instance, and open the shapefile"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = Basemap(\n",
      "    projection='tmerc',\n",
      "    lon_0 = -2.,\n",
      "    lat_0 = 49.,\n",
      "    ellps = 'WGS84',\n",
      "    llcrnrlon=coords[0] - extra * w,\n",
      "    llcrnrlat=coords[1] - extra + 0.01 * h,\n",
      "    urcrnrlon=coords[2] + extra * w,\n",
      "    urcrnrlat=coords[3] + extra + 0.01 * h,\n",
      "    lat_ts=0,\n",
      "    resolution='i',\n",
      "    suppress_ticks=True)\n",
      "m.readshapefile(\n",
      "    'data/london_wards',\n",
      "    'london',\n",
      "    color='none',\n",
      "    zorder=2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "(649,\n",
        " 5,\n",
        " [-0.5103750689005356, 51.28676016315085, 0.0, 0.0],\n",
        " [0.3340155643740321, 51.691874116909894, 0.0, 0.0],\n",
        " <matplotlib.collections.LineCollection at 0x10dd58550>)"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARoAAADtCAYAAABgU7qvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAAzdJREFUeJzt2buNwzAQQEHKcEVqQXWrBbZEN+BfcO8MCDMhl8FGDyC4\nrbUGQOn26wWA6xMaICc0QE5ogJzQADmhAXL3d8PzPP19A187jmN7dv42NGOMse/7328DXM6c8+XM\n0wnICQ2QExogJzRATmiAnNAAOaEBckID5IQGyAkNkBMaICc0QE5ogJzQADmhAXJCA+SEBsgJDZAT\nGiAnNEBOaICc0AA5oQFyQgPkhAbICQ2QExogJzRATmiAnNAAOaEBckID5IQGyAkNkBMaICc0QE5o\ngJzQADmhAXJCA+SEBsgJDZATGiAnNEBOaICc0AA5oQFyQgPkhAbICQ2QExogJzRATmiAnNAAOaEB\nckID5IQGyAkNkBMaICc0QE5ogJzQADmhAXJCA+SEBsgJDZATGiAnNEBOaICc0AA5oQFyQgPkhAbI\nCQ2QExogJzRATmiAnNAAOaEBckID5IQGyAkNkBMaICc0QE5ogJzQADmhAXJCA+SEBsgJDZATGiAn\nNEBOaICc0AA5oQFyQgPkhAbICQ2QExogJzRATmiAnNAAOaEBckID5IQGyAkNkBMaICc0QE5ogJzQ\nADmhAXJCA+SEBsgJDZATGiAnNEBOaICc0AA5oQFyQgPkhAbICQ2QExogJzRATmiAnNAAOaEBckID\n5IQGyAkNkBMaICc0QE5ogJzQADmhAXJCA+SEBsgJDZATGiAnNEBOaICc0AA5oQFyQgPkhAbICQ2Q\nExogJzRATmiAnNAAOaEBckID5IQGyAkNkBMaICc0QE5ogJzQADmhAXJCA+SEBsgJDZATGiAnNEBO\naICc0AA5oQFyQgPkhAbICQ2QExogJzRATmiAnNAAOaEBckID5IQGyAkNkBMaICc0QE5ogJzQADmh\nAXJCA+SEBsgJDZATGiAnNEBOaICc0AA5oQFyQgPkhAbICQ2QExogJzRATmiAnNAAOaEBckID5IQG\nyAkNkBMaICc0QE5ogJzQADmhAXJCA+SEBsgJDZATGiAnNEBOaICc0AA5oQFyQgPkhAbICQ2QExog\nJzRATmiA3P3ThTnnf+wBXNi21vr1DsDFeToBOaEBckID5IQGyAkNkHsAuYYRQ7xy/dQAAAAASUVO\nRK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10e7e5810>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Create a DataFrame holding ward polygons and area"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# set up a map dataframe\n",
      "df_map = pd.DataFrame({\n",
      "    'poly': [Polygon(xy) for xy in m.london],\n",
      "    'ward_name': [w['NAME'] for w in m.london_info],\n",
      "})\n",
      "df_map['area_m'] = df_map['poly'].map(lambda x: x.area)\n",
      "df_map['area_km'] = df_map['area_m'] / 100000\n",
      "\n",
      "# Create Point objects in map coordinates from dataframe lon and lat values\n",
      "map_points = pd.Series(\n",
      "    [Point(m(mapped_x, mapped_y)) for mapped_x, mapped_y in zip(df['lon'], df['lat'])])\n",
      "plaque_points = MultiPoint(list(map_points.values))\n",
      "wards_polygon = prep(MultiPolygon(list(df_map['poly'].values)))\n",
      "# calculate points that fall within the London boundary\n",
      "ldn_points = filter(wards_polygon.contains, plaque_points)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Calculate plaque density data and bins"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# create a MultiPoint which we can check for set membership\n",
      "df_map['count'] = df_map['poly'].map(lambda x: int(len(filter(prep(x).contains, ldn_points))))\n",
      "df_map['density_m'] = df_map['count'] / df_map['area_m']\n",
      "df_map['density_km'] = df_map['count'] / df_map['area_km']\n",
      "# it's easier to work with NaN values when classifying\n",
      "df_map.replace(to_replace={'density_m': {0: np.nan}, 'density_km': {0: np.nan}}, inplace=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate Jenks natural breaks for density\n",
      "breaks = nb(\n",
      "    df_map[df_map['density_km'].notnull()].density_km.values,\n",
      "    initial=300,\n",
      "    k=6)\n",
      "# the notnull method lets us match indices when joining\n",
      "jb = pd.DataFrame({'jenks_bins': breaks.yb}, index=df_map[df_map['density_km'].notnull()].index)\n",
      "df_map = df_map.join(jb)\n",
      "df_map.jenks_bins.fillna(-1, inplace=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Use the [Folium](https://github.com/wrobstory/folium) library to impose a choropleth map onto an in-browser map (OSM, MapBox, Stamen etc.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# In order to use custom maps, pass any of the following as a tiles keyword:\n",
      "\"\"\"\n",
      "-\"OpenStreetMap\"\n",
      "-\"MapQuest Open\"\n",
      "-\"MapQuest Open Aerial\"\n",
      "-\"Mapbox Bright\" (Limited levels of zoom for free tiles)\n",
      "-\"Mapbox Control Room\" (Limited levels of zoom for free tiles)\n",
      "-\"Stamen Terrain\"\n",
      "-\"Stamen Toner\"\n",
      "-\"Cloudmade\" (Must pass API key)\n",
      "-\"Mapbox\" (Must pass API key)\n",
      "\"\"\"\n",
      "# we're binding the dataframe ward_name column to the GeoJSON feature.properties.NAME value\n",
      "# this will vary according to shapefile / GeoJSON\n",
      "# the data to be displayed as choropleth comes from the second 'columns' entry\n",
      "map_osm = folium.Map(\n",
      "    location=[51.5, 0.],\n",
      "    zoom_start=10,\n",
      "    )\n",
      "map_osm.geo_json(\n",
      "    geo_path='wards.geojson',\n",
      "    data=df_map,\n",
      "    columns=['ward_name', 'density_km'],\n",
      "    key_on='feature.properties.NAME',\n",
      "    threshold_scale = breaks.bins,\n",
      "    fill_color='BuPu', fill_opacity=0.75, line_opacity=0.2,\n",
      "    legend_name='Density (km2)')\n",
      "map_osm.create_map(path='ward_map.html')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "HTML('<iframe src=ward_map.html width=800 height=550></iframe>')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe src=ward_map.html width=800 height=550></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<IPython.core.display.HTML at 0x10e7131d0>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}