{
 "metadata": {
  "name": "",
  "signature": "sha256:e997bcc65d4ed8d4e17e808f973e7759a9bfabc34f78eeedbdd31c81671f79ac"
 },
 "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 1)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "# Converting coordinates with PyProj"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pyproj import Proj\n",
      "\n",
      "# Create projection transformation object\n",
      "p = Proj(init='epsg:3857') # EPSG code for Web Mercator\n",
      "\n",
      "# Convert from long/lat to Mercator and back\n",
      "print(p(-97.740372, 30.282642))\n",
      "print(p(-10880408.440985134, 3539932.8204972977, inverse=True))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Fiona (which we will be using shortly) has several \n",
      "# helper functions for working with proj4 strings\n",
      "from fiona.crs import to_string, from_epsg, from_string\n",
      "\n",
      "# Create a crs dict from a proj4 string\n",
      "crs = from_string('+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 '\n",
      "                  '+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 '\n",
      "                  '+y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs')\n",
      "\n",
      "# Using a proj4 string\n",
      "nyc_proj = Proj(crs, preserve_units=True)\n",
      "\n",
      "# Using an EPSG code\n",
      "nyc_epsg = Proj(init='epsg:2263', preserve_units=True)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Here's about where my office in NYC is located (in long/lat)\n",
      "office = (-73.9637, 40.7684)\n",
      "\n",
      "# Are they close?\n",
      "print(nyc_proj(*office))\n",
      "print(nyc_epsg(*office))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "source": [
      "## Plotting Eyjafjallaj\u00f6kull volcano with Cartopy"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "import cartopy.crs as ccrs\n",
      "import matplotlib.pyplot as plt\n",
      "import cartopy.io.img_tiles as cimgt\n",
      "\n",
      "# Create a MapQuest open aerial instance\n",
      "map_quest_aerial = cimgt.MapQuestOpenAerial()\n",
      "\n",
      "# What is the projection?\n",
      "print(map_quest_aerial.crs.proj4_init)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Specify the lon/lat for the volcano\n",
      "volcano = (-19.613333, 63.62)\n",
      "\n",
      "# Define the plotting extent of the map\n",
      "extent = [-22, -15, 63, 65]\n",
      "\n",
      "# Specify the transform to use when plotting\n",
      "transform=ccrs.Geodetic()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(10,8))\n",
      "# Create a GeoAxes in the tile's projection\n",
      "ax = plt.axes(projection=map_quest_aerial.crs)\n",
      "ax.set_extent(extent)\n",
      "# Add the MapQuest data at zoom level 8\n",
      "ax.add_image(map_quest_aerial, 8)\n",
      "ax.plot(*volcano, marker='o', color='yellow', markersize=12,\n",
      "        alpha=0.7, transform=transform)\n",
      "ax.set_title(u'Eyjafjallaj\u00f6kull Volcano')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "subslide"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Time to work on Notebook:\n",
      "\n",
      "[`Working with Projections in Python`](../exercises/Working with Projections in Python.ipynb)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}