examples/pyproj example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transformations with pyproj\n",
    "\n",
    "`pyproj` is a python interface to the venerable PROJ.4 library. Transformations to and from geographic coordinates can be done by a `Proj` instance, which can be defined in a few different ways.\n",
    "\n",
    "Transforming from longitude, latitude to map coordinates is a call on a `Proj` instance that we have initiallized by its EPSG code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(586912.6635852784, 4514845.724019826)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "from pyproj import Proj, transform\n",
    "epsg32618 = Proj(init='epsg:32618')\n",
    "lat, lon = 40.78, -73.97\n",
    "x, y = epsg32618(lon, lat)\n",
    "x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse transformation is done by the same object with the keyword parameter `inverse=True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-73.96999999999984, 40.77999999999991)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epsg32618(x, y, inverse=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can initialize a UTM projection by its zone. Zone 18 is equivalent to EPSG above, so we should ge the same result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(586912.6635846694, 4514845.7241347665)\n"
     ]
    }
   ],
   "source": [
    "utm18 = Proj(proj=\"utm\", zone=\"18\")\n",
    "print(utm18(lon, lat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A string like `\"+proj=utm +zone=18\"` can be used as well (this is compatible with `PROJ.4` command line arguments)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(586912.6635846694, 4514845.7241347665)\n"
     ]
    }
   ],
   "source": [
    "utm18 = Proj(\"+proj=utm +zone=18\")\n",
    "print(utm18(lon, lat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do datum transformations as well as map projections. To create a projection in the NAD27 datum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(586914.9904254354, 4514634.228909769)\n"
     ]
    }
   ],
   "source": [
    "nad27 = Proj(proj=\"utm\", zone=\"18\", ellps=\"clrk66\", datum=\"NAD27\")\n",
    "print(nad27(lon, lat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could also do this transfromation between the two reference frames directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(586914.9904260582, 4514634.228794821)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "old_x, old_y = transform(utm18, nad27, x, y)\n",
    "old_x, old_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Libraries including Fiona, rasterio, and GeoPandas use a lightweight python dictionary of parameters that can be passed to the `Proj` constructor when it is needed for a transformation. Here we use the explicit parameters equivalent to EPSG:2263 (NY state plane Long Island, US feet):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(302532.4183897138, 68108.70199755549)\n"
     ]
    }
   ],
   "source": [
    "crs = {'lon_0': -74, 'datum': 'NAD83', 'y_0': 0, 'no_defs': True, 'proj': 'lcc',\n",
    "       'x_0': 300000, 'units': 'us-ft', 'lat_2': 41.03333333333333,\n",
    "       'lat_1': 40.66666666666666, 'lat_0': 40.16666666666666}\n",
    "nyc = Proj(**crs)\n",
    "print(nyc(lon, lat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mini-exercise:\n",
    "\n",
    "Transform the locations of the Texas state capitol (97.74035° W, 30.27467° N) and the AT&T Center (97.74034° W, 30.28240° N) to UTM zone 14 and calculate the distance between them in meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "name": "python",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}