exercises/Working with Vectors in Python.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:680f715a61852e6bfe3fa34908b6fc2e66481f8915621ab4e63bcd71474f04e9"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Working with Vectors in Python\n",
      "\n",
      "## Points in Polygon Analysis\n",
      "\n",
      "For this first exercise, we're going to go through the steps of performing a simple 'points in polygon' analysis together, using NYC boroughs and public WiFi hotspots to demonstrate the simplicity of spatial and aspatial queries with GeoPandas."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import geopandas as gpd\n",
      "import sqlite3\n",
      "import os"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since I didn't want to have to rely on getting PostGIS, Spatialite, or some other spatial database up and running on everyone's laptop, we're going to fake it using SQLite.\n",
      "\n",
      "We'll use plain ol' Pandas to query the database (which contains `WKB` geometries), and then 'upgrade' the `DataFrame` to a `GeoDataFrame` for further processing\\*."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create connection to sqlite database\n",
      "conn = sqlite3.connect(os.path.join('..', 'data', 'wifi.sqlite'))\n",
      "sql = \"select * from wifi\"\n",
      "\n",
      "# GeoPandas has its own reference to Pandas for convinience\n",
      "wifi = gpd.pd.read_sql(sql, conn)\n",
      "wifi.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We could add `where` clauses and other queries to the above `SQL` (for example, to select only the 'Free' WiFi hotspots), but for the sake of this demonstration we'll do most of this stuff in Python instead.\n",
      "\n",
      "Now that we have a `DataFrame` with a `GEOMETRY` column, we need to convert the `WKB` to something that GeoPandas (and Shapely) understands. We'll create a helper function for converting from `WKB` to Shapely `geometry` objects:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from shapely import wkb\n",
      "from binascii import unhexlify\n",
      "from codecs import encode\n",
      "\n",
      "# Helper function\n",
      "def as_geom(hex):\n",
      "    return wkb.loads(unhexlify(encode(hex, \"hex\")))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that we have our helper function, let's create a new column in the data frame called `geom`, and then use our awesome `set_geometry` method to upgrade the `DataFrame` to a `GeoDataFrame`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "wifi[\"geom\"] = wifi[\"GEOMETRY\"].apply(as_geom)\n",
      "del wifi[\"GEOMETRY\"]  # Don't need this anymore\n",
      "\n",
      "# Coerse to GeoDataFrame... cool!\n",
      "wifi = wifi.set_geometry(\"geom\", crs={\"init\":\"epsg:2263\"})\n",
      "wifi.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'll also need a polygon layer to work with. For this, we'll use the NYC boroughs shapefile in your data directory:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "boros_file = os.path.join('..', 'data', 'nybb', 'nybb.shp')\n",
      "boros = gpd.read_file(boros_file)\n",
      "boros"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For our points in polygon analysis, we'll focus on the borough of Manhattan:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set the index to the borough names\n",
      "boros.set_index('BoroName', inplace=True)\n",
      "boros.sort()\n",
      "\n",
      "# We'll focus on Manhattan\n",
      "poly = boros.geometry[\"Manhattan\"]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In order to facilitate our point in polygon analysis, we're going to 'prepare' our selected polygon for efficient queries. Shapely geometries can be processed into a state that supports more efficient batches of operations by creating 'prepared geometries'. For more information on this, see the [Shapely documentation](http://toblerity.org/shapely/manual.html#prepared-geometry-operations).\n",
      "\n",
      "Note that this ignores any CRS differences, so if your two 'layers' are in different CRSs, you may get unexpected results. Best practice is to ensure that they are in the same CRS (in this case, they are)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from shapely.prepared import prep\n",
      "prep_poly = prep(poly)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "How many WiFi hotspots in manhattan?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mht_wifi = wifi.geometry.apply(poly.contains)\n",
      "print(sum(mht_wifi))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What percentage of WiFi hotspots in Manhattan are **Free**?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "free = wifi[wifi[\"type\"]==\"Free\"]\n",
      "import numpy as np\n",
      "\n",
      "mht_free = free.geometry.apply(poly.contains)\n",
      "print(float(sum(mht_free)) / float(sum(mht_wifi)) * 100.)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ok, let's plot *only* the **Free** WiFi hotspots in the borough of **Manhattan**:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "fig = plt.figure(figsize=(8, 8))\n",
      "ax = free[mht_free].plot()\n",
      "\n",
      "# Let's clean things up a bit...\n",
      "for i, l in enumerate(ax.lines):\n",
      "    l.set_markersize(16)\n",
      "    l.set_markeredgecolor(\"green\")\n",
      "    l.set_color(\"green\")\n",
      "ax.axis('off')  # Turn off axes\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, just because we can, let's make a 'web-map' of this 'analysis'. This time we'll use a non-default basemap. I like [Stamen](http://stamen.com/)'s maptiles, but if you have the tile url, you can use any maptiles you want!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import mplleaflet\n",
      "# Credit where credit is due...\n",
      "attr = 'Map tiles by <a href=\"http://stamen.com/\">Stamen</a> Data by <a href=\"http://openstreetmap.org\">OpenStreetMap</a>, under <a href=\"http://creativecommons.org/licenses/by-sa/3.0\">CC BY SA</a>.'\n",
      "tile = \"http://tile.stamen.com/%s/{z}/{x}/{y}.jpg\" % \"toner\" # Could also use \"terrain\" or \"watercolor\"\n",
      "mplleaflet.display(fig=ax.figure, crs=free.crs, tiles=(tile, attr))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Geology Summaries\n",
      "\n",
      "This second exercise is self-directed. Try to perform the following tasks using a shapefile of the geology of New Jersey that was downloaded from the USGS."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "geol_file = os.path.join(\"..\", \"data\", \"njgeol\", \"njgeol_poly_dd.shp\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 1a) Read in the `nygeol_poly_dd.shp` shapefile using GeoPandas"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 1b) What is the CRS for this 'layer'? Is it 'projected'?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 2a) Create a list of all the unique rock types in the data (in properties ROCKTYPE1 and ROCKTYPE2). How many are there?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 2b) Calculate the total area of each primary rocktype (ROCKTYPE1) by summing the AREA column.\n",
      "\n",
      "Note: The areas are in square degrees, not true area."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 2c) What would be a better way to compute these areas?\n",
      "\n",
      "Hint: The New Jersey State Plane EPSG code is 32111."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 2d) Recompute the areas for each ROCKTYPE1 using this method."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 3a) Plot the polygons (using any packages you want), coloring by primary rock type.\n",
      "\n",
      "Keep the plotting simple, and you can use random colors for the rock types."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Bonus) How would you share this 'map' with colleagues? It would be useful if they could also edit the map themselves.\n",
      "\n",
      "Warning: This is relatively large vector file (for the web) and may be extremely slow, maybe you should subset the data first!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "# Take a random subset\n",
      "rand = np.random.randint(0, geol.shape[0], 100)\n",
      "subset = geol.iloc[rand]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import geojsonio\n",
      "\n",
      "# Other commands go here"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If geojson.io isn't working (it might not work nicely with Safari?), go to the following url (yours will differ) to find the anonymous gist:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"https://gist.github.com/anonymous/%s\" % res[res.find(\"gist:/\")+6:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\\*In future releases of GeoPandas, this will work much more smoothly, but for now we're 'hacking it'."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}