gis-tools-101/gis_tools.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:d6b500af509e1720059457013e554297cca119256cc2727750306711aad9f1ba"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Makefile\n",
      "Creating a *`Makefile`* allows others to recreate necessary components of your work flow while managing unecessary steps.\n",
      "\n",
      "###Why do this now?\n",
      "The *`Makefile`* process can be useful when working with shapefiles. Instead of carrying large .shp files around in a repo, we'll give anyone the opportunity to create the files.  \n",
      "\n",
      "###Makefile terminology and symbols\n",
      "Pseudo code for each section of a *`Makefile`* is as follows.  \n",
      "\n",
      "`Target: dependency`  \n",
      "&nbsp; `commands to run`\n",
      "\n",
      "The following variables come in handy when writing make files.  \n",
      "&nbsp; `$@` =  `target`   \n",
      "&nbsp; `$(dir $@)` = `target directory`  \n",
      "&nbsp; `$(notdir $@)` = `target filename`  \n",
      "&nbsp; `$<` = `dependency`   \n",
      "\n",
      "###Example\n",
      "For example, we could create a target named `somePath/someFile.txt` that depends on first building `stuff`.  \n",
      "\n",
      "`stuff`:  \n",
      "&nbsp; `echo $@`  \n",
      "&nbsp; `echo $(notdir $@)`  \n",
      "\n",
      "`somePath/someFile.txt: stuff`  \n",
      "&nbsp; `echo $@`  \n",
      "&nbsp; `echo $(dir $@)`  \n",
      "&nbsp; `echo $(notdir $@)`  \n",
      "&nbsp; `echo $<`  \n",
      " \n",
      " To see this example execute, run the following script from the command line within the rev_geo directory. \n",
      " <pre> make somePath/someFile.txt </pre> "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Shapefiles\n",
      "\n",
      "The [terminology](https://www.census.gov/geo/maps-data/data/pdfs/tiger/tgrshp2013/TGRSHP2013_TechDoc_Ch4.pdf) for shapefiles can be a little confusing, but the main ideas will be discussed below.\n",
      "\n",
      "First, let's get some data. To download a zip file containing county data and then upzip it in a new `data` directory, please run the following script from the command line within the rev_geo directory.  \n",
      "<pre> make data/tl_2013_us_county.shp </pre>  \n",
      "When it finishes, you can view the files in `geo/data`, which were downloaded from the US Census Bureau. \n",
      "\n",
      "* [Click](https://www.census.gov/geo/maps-data/data/tiger.html) for more information about the census data used in this example.   \n",
      "\n",
      "### What's inside a shapefile? \n",
      "We'll use python's `fiona` and `shapely` packages to view the contents of the shapefile and determine geometric properties of the corresponding features. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pprint\n",
      "import fiona\n",
      "\n",
      "\n",
      "with fiona.open(\"data/tl_2013_us_county.shp\") as fc: \n",
      "    print \"File type:\",fc.driver\n",
      "    print \"Schema:\",pprint.pprint(fc.schema)\n",
      "    print \"Number of records:\", len(fc)\n",
      "    print \"Bounds of all records:\", fc.bounds #visit http://boundingbox.klokantech.com/ to view these coords\n",
      "    print \"Additional info (coordinate reference system):\", fc.crs\n",
      "    print \"Field type map:\"\n",
      "    pprint.pprint(fiona.FIELD_TYPES_MAP)\n",
      "    records=list(fc)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Features\n",
      "* The collection (fc) contains 3,234 records or features. \n",
      "* Each feature is a Python dict structured exactly like a [GeoJSON](http://geojson.org/) Feature. \n",
      "* Each feature starts and ends with the same point (eg. complete polygon)\n",
      "* Each feature (in this example) contains information about a specific county.\n",
      "* [Click](http://toblerity.org/fiona/manual.html) for more information about `fiona`.   \n",
      "* [Click](http://www.census.gov/geo/reference/ansi.html) for codes related to the features' 'properites'.  \n",
      "* [Click](https://www.census.gov/geo/maps-data/data/pdfs/tiger/tgrshp2013/TGRSHP2013_TechDoc_A.pdf) for codes specific to county features (page A-80)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#import fiona\n",
      "#import pprint\n",
      "\n",
      "\n",
      "pprint.pprint(records[0])\n",
      "print \n",
      "print\n",
      "print \"###########################\"\n",
      "print \"alternative method below\"\n",
      "print \"###########################\"\n",
      "print \n",
      "print \n",
      "with fiona.open(\"data/tl_2013_us_county.shp\") as fc: \n",
      "    first=True\n",
      "    for record in fc:\n",
      "        if first:\n",
      "            print pprint.pprint(record)\n",
      "            break\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Analyzing features\n",
      "* use Python's 'shapely' package.\n",
      "* [Click](http://toblerity.org/shapely/manual.html#object.intersects) for documentation on evaluating if shapes intersect. \n",
      "* [Click](http://toblerity.org/shapely/manual.html#object.contains) for documentation on evaluating if points are contained within shapes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#import fiona\n",
      "from shapely.geometry import Point, shape, Polygon, box\n",
      "\n",
      "#print records[0]\n",
      "#print shape(records[0][\"geometry\"]) # one feature converted to shapley polygon object\n",
      "\n",
      "print \"{} Shapley object: {}\".format(records[0]['properties']['NAMELSAD'],box(*shape(records[0][\"geometry\"]).bounds))\n",
      "print\n",
      "shp1 = box(*shape(records[0][\"geometry\"]).bounds)\n",
      "print \"{} Shapley object: {}\".format(records[1]['properties']['NAMELSAD'],box(*shape(records[1][\"geometry\"]).bounds))\n",
      "print \n",
      "shp2 = box(*shape(records[1][\"geometry\"]).bounds)\n",
      "print \"Does {} intersect {}?\".format(records[0]['properties']['NAMELSAD'],records[1]['properties']['NAMELSAD'])\n",
      "print \"{}\".format(shp1.intersects(shp2)) \n",
      "print \n",
      "print \"Is (-123.4688,46.2674) in {}?\".format(records[0]['properties']['NAMELSAD'])\n",
      "p1=Point(-123.4688,46.2674)\n",
      "print \"{}\".format(shape(records[0][\"geometry\"]).contains(p1))\n",
      "print\n",
      "print \"Is (-123.4688,46.2674) in {}?\".format(records[1]['properties']['NAMELSAD'])\n",
      "print \"{}\".format(shape(records[1][\"geometry\"]).contains(p1))\n",
      "print \n",
      "print \"Create a list of centroids for every county.\"\n",
      "centroids=[]\n",
      "with fiona.open(\"data/tl_2013_us_county.shp\") as fc: \n",
      "    with open(\"data/centroids.txt\",\"wb\") as cntr:\n",
      "        for feature in fc:\n",
      "            pt=(float(feature[\"properties\"][\"INTPTLON\"]),float(feature[\"properties\"][\"INTPTLAT\"]))\n",
      "            centroids.append(pt)\n",
      "            cntr.write(str(pt)+'\\n')\n",
      "print \"{} centroids created\".format(len(centroids))\n",
      "    \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Utilities\n",
      "Two tools that can be handy when working with shapefiles:\n",
      "* ogr2ogr - part of [GDAL](http://www.gdal.org/) that merges multiple .shp files into a single [GeoJson](http://geojson.org/) file (among many other features).\n",
      "* [topojson](https://github.com/mbostock/topojson/wiki) - shrinks GeoJson files to a convenient format. \n",
      "\n",
      "[Example](http://bost.ocks.org/mike/map) that uses both tools."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Using rev_geo.py\n",
      "This utility takes coordinates and returns location information, but requires the user to point at the tl_2013_us_county.shp file. \n",
      "* currently only available for US state/county info\n",
      "\n",
      "###Standard options\n",
      "The options are rather limited at this point.\n",
      "* points outside of the US are currently excluded.  \n",
      "* considerations:   \n",
      "&nbsp; 1.  general solution for an efficiency grid indpendent of shapefile used.\n",
      "\n",
      "<pre>\n",
      "usage: rev_geo.py [-h] [-b GRID_BOUNDARIES] [-d DELTA] [-g]\n",
      "                  [-s SHAPE_FILE_PATH] [-t]\n",
      "                  [file_name]\n",
      "\n",
      "Reverse geo coder returns location info given a set of lon,lat\n",
      "\n",
      "positional arguments:\n",
      "  file_name             Input file name (optional).\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  -b GRID_BOUNDARIES, --bounding-box GRID_BOUNDARIES\n",
      "                        Set bounding box for region to include (default:\n",
      "                        [-185,15,-65,70])\n",
      "  -d DELTA, --delta DELTA\n",
      "                        Set the number of degrees between grid coords\n",
      "                        (default: 5)\n",
      "  -g, --use-saved-grid  Save grid or use previously saved version in\n",
      "                        data/grid.json\n",
      "  -s SHAPE_FILE_PATH, --shape-file-path SHAPE_FILE_PATH\n",
      "                        Set shapefile path (default:\n",
      "                        data/tl_2013_us_county.shp)\n",
      "  -t, --tweet-input     Set input as tweet payload instead of coordinates (in\n",
      "                        progress)\n",
      "</pre>\n",
      "\n",
      "###Standard output\n",
      "The general form of the output includes the following:\n",
      "* county: str [county name]  \n",
      "* centroid: (longitude, latitude)  [center of county]  \n",
      "* coords: (longitude, latitude) [specific coords passed to rev_geo.py]  \n",
      "* GEOID: int(5) [state_code + county_code]  \n",
      "<pre>\n",
      "{\"county\": \"Wahkiakum\", \n",
      " \"centroid\": [-123.4244583, 46.2946377], \n",
      " \"coords\": [-123.4244583, 46.2946377], \n",
      " \"GEOID\": \"53069\"}\n",
      "</pre>\n",
      "\n",
      "###Example script\n",
      "<pre>\n",
      "$head data/centroids.txt | ./rev_geo.py -g > info.json\n",
      "$cat info.json\n",
      "{\"county\": \"Cuming\", \"centroid\": [-96.7885168, 41.9158651], \"coords\": [-96.7885168, 41.9158651], \"GEOID\": \"31039\"}\n",
      "{\"county\": \"Wahkiakum\", \"centroid\": [-123.4244583, 46.2946377], \"coords\": [-123.4244583, 46.2946377], \"GEOID\": \"53069\"}\n",
      "{\"county\": \"De Baca\", \"centroid\": [-104.3686961, 34.3592729], \"coords\": [-104.3686961, 34.3592729], \"GEOID\": \"35011\"}\n",
      "{\"county\": \"Lancaster\", \"centroid\": [-96.6886584, 40.7835474], \"coords\": [-96.6886584, 40.7835474], \"GEOID\": \"31109\"}\n",
      "{\"county\": \"Nuckolls\", \"centroid\": [-98.0468422, 40.1764918], \"coords\": [-98.0468422, 40.1764918], \"GEOID\": \"31129\"}\n",
      "{\"county\": \"Las Piedras\", \"centroid\": [-65.871189, 18.1871483], \"coords\": [-65.871189, 18.1871483], \"GEOID\": \"72085\"}\n",
      "{\"county\": \"Minnehaha\", \"centroid\": [-96.7957261, 43.6674723], \"coords\": [-96.7957261, 43.6674723], \"GEOID\": \"46099\"}\n",
      "{\"county\": \"Menard\", \"centroid\": [-99.8539896, 30.8843655], \"coords\": [-99.8539896, 30.8843655], \"GEOID\": \"48327\"}\n",
      "{\"county\": \"Sierra\", \"centroid\": [-120.5219926, 39.5769252], \"coords\": [-120.5219926, 39.5769252], \"GEOID\": \"06091\"}\n",
      "{\"county\": \"Clinton\", \"centroid\": [-85.1534262, 36.7288647], \"coords\": [-85.1534262, 36.7288647], \"GEOID\": \"21053\"}\n",
      "</pre>\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}