{
 "metadata": {
  "name": "",
  "signature": "sha256:4d754bcea3d58808b0848f3e45708522b0bf3b135b99e028b2f67fdc390578aa"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Geospatial Data in Python: Database, Desktop, and the Web\n",
      "## Tutorial (Part 2)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from shapely.geometry import LineString\n",
      "\n",
      "# Dilating a line\n",
      "line = LineString([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])\n",
      "dilated = line.buffer(0.5)\n",
      "eroded = dilated.buffer(-0.3)\n",
      "\n",
      "# Demonstate Python Geo Interface\n",
      "print(line.__geo_interface__)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# Exploring the path of Hurican Katrina\n",
      "\n",
      "The data was originally sourced from the HURDAT2 dataset from [AOML/NOAA](http://www.aoml.noaa.gov/hrd/hurdat/newhurdat-all.html), and the Python lists are from the [cartopy documentation](http://scitools.org.uk/cartopy/docs/latest/examples/hurricane_katrina.html)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from shapely.geometry import LineString\n",
      "\n",
      "lons = [-75.1, -75.7, -76.2, -76.5, -76.9, -77.7, -78.4, -79.0,\n",
      "        -79.6, -80.1, -80.3, -81.3, -82.0, -82.6, -83.3, -84.0,\n",
      "        -84.7, -85.3, -85.9, -86.7, -87.7, -88.6, -89.2, -89.6,\n",
      "        -89.6, -89.6, -89.6, -89.6, -89.1, -88.6, -88.0, -87.0,\n",
      "        -85.3, -82.9]\n",
      "lats = [23.1, 23.4, 23.8, 24.5, 25.4, 26.0, 26.1, 26.2, 26.2, 26.0,\n",
      "        25.9, 25.4, 25.1, 24.9, 24.6, 24.4, 24.4, 24.5, 24.8, 25.2,\n",
      "        25.7, 26.3, 27.2, 28.2, 29.3, 29.5, 30.2, 31.1, 32.6, 34.1,\n",
      "        35.6, 37.0, 38.6, 40.1]\n",
      "\n",
      "# Turn the lons and lats into a shapely LineString\n",
      "katrina_track = LineString(zip(lons, lats))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "source": [
      "Buffer the linestring by two degrees (doesn't really make sense!). This *should* be about 200kms, but as we'll see, it's not quite accurate... **Why not**?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "katrina_buffer = katrina_track.buffer(2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "What if we reproject the lon/lats to a projection that preserves distances better?\n",
      "\n",
      "We *could* use `EPSG:32616` (UTM Zone 16), which covers where Katrina meets New Orleans, but we're probably better off using a custom `proj4` string based on a Lambert Conformal Conic projection. **Why**?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pyproj import Proj, transform\n",
      "from fiona.crs import from_string\n",
      "\n",
      "# Create custom proj4 string\n",
      "proj = from_string('+ellps=WGS84 +proj=lcc +lon_0=-96.0 +lat_0=39.0 '\n",
      "                   '+x_0=0.0 +y_0=0.0 +lat_1=33 +lat_2=45 +no_defs')\n",
      "\n",
      "# Create source and destination Proj objects (source is WGS84 lons/lats)\n",
      "src = Proj(init='epsg:4326')\n",
      "dst = Proj(proj)\n",
      "\n",
      "# Create a LineString from the transformed points\n",
      "proj_track = LineString(zip(*transform(src, dst, lons, lats)))\n",
      "# Buffer the LineString by 200 km (multiply by 1000 to work in meters)\n",
      "proj_buffer = proj_track.buffer(200*1000)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "## Aside: Coordinate tuples and x, y sequences\n",
      "\n",
      "`zip` is your friend! Use it with tuple unpacking to change between sequences of `(x, y)` pairs and seperate `x` and `y` sequences."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pts = [(0, 0), (1, 0), (1, 1), (2, 1), (2, 2)]\n",
      "x, y = zip(*pts)\n",
      "print x, y"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "source": [
      "Also, instead of calling `f(x, y)`, you can just use"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Skip this slide, not really needed for demo, just need Python function\n",
      "# Simple function to add 0.5 to each coordinate\n",
      "def f(x, y):\n",
      "    new_x = [i + 0.5 for i in x]\n",
      "    new_y = [j + 0.5 for j in y]\n",
      "    return new_x, new_y"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "skip"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(*zip(*pts)) # Function f adds 0.5 to each coordinate"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Plotting the path of Hurican Katrina"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(8, 8))\n",
      "ax = plt.axes(projection=ccrs.LambertConformal())\n",
      "ax.stock_img()  # Add background image (slow)\n",
      "ax.coastlines(resolution='110m')\n",
      "ax.add_geometries([katrina_buffer], ccrs.PlateCarree(), facecolor='blue', alpha=0.5)\n",
      "ax.add_geometries([katrina_track], ccrs.PlateCarree(), facecolor='none')\n",
      "\n",
      "# Let's add the projected buffer for comparison\n",
      "ax.add_geometries([proj_buffer], ccrs.LambertConformal(), facecolor='green', alpha=0.5)\n",
      "\n",
      "ax.set_extent([-125, -66.5, 20, 50], ccrs.PlateCarree())\n",
      "ax.gridlines()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Which `shapely` geometry method could we use to find where the tracks *differ*?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# Simple Transformation for Georeferencing a Raster\n",
      "\n",
      "What makes geospatial raster datasets different from other raster files is that their pixels map to regions of the Earth. In this case, we have a raster image which maps to Midtown Manhattan."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.image as mpimg\n",
      "import os\n",
      "\n",
      "# Read in a regular PNG image of Manhattan\n",
      "png_file = os.path.join('..', 'data', 'manhattan.png')\n",
      "img = mpimg.imread(png_file)\n",
      "\n",
      "# Take a quick look at the shape\n",
      "print(img.shape)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Specify the affine transformation\n",
      "from affine import Affine # You might not have this library installed\n",
      "A = Affine(0.999948245999997, 0.0, 583057.357, 0.0, -0.999948245999997, 4516255.36)\n",
      "\n",
      "# Compute the upper left and lower right corners\n",
      "ul = (0, 0)\n",
      "lr = img.shape[:2][::-1]\n",
      "\n",
      "print(\"Upper left:  \" + str(A * ul))\n",
      "print(\"Lower right: \" + str(A * lr))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here, the coordinate reference system is `EPSG:26918` (NAD83 / UTM Zone 18N), and the affine transformation matrix is given (in later examples, we'll get this information directly from the input raster datasets)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Upper left and bottom right corners in UTM coords\n",
      "left, top = A * ul\n",
      "right, bottom = A * lr\n",
      "\n",
      "# Plot showing original PNG image (axes correspond to rows and cols) on left\n",
      "# and 'transformed' PNG (axes correspond to UTM coords) on the right\n",
      "f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))\n",
      "ax1.imshow(img)\n",
      "ax1.set_title(\"PNG with row, col bounds\")\n",
      "ax2.imshow(img, extent=(left, right, bottom, top), aspect=\"equal\")\n",
      "ax2.set_title(\"PNG with correct bounds\")\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Time to work on Notebook:\n",
      "\n",
      "[`Plotting Great Circles in Python`](../exercises/Plotting Great Circles in Python.ipynb)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}