posts/pyproj.mdwn

   1 [pyproj][] is a Python wrapper around [PROJ.4][].  Here's a quick
   2 walkthrough.
   3
   4 Initialize a [geodetic][] converter:
   5
   6     >>> from pyproj import Geod
   7     >>> g = Geod(ellps='clrk66')
   8
   9 where `ellps='clrk66'` selects [Clarke's 1866][clrk66] [reference
  10 ellipsoid][rell].  `help(Geod.__new__)` gives a list of possible
  11 ellipsoids.
  12
  13 Calculate the distance between two points, as well as the local
  14 heading, try
  15
  16     >>> lat1,lng1 = (40.7143528, -74.0059731)  # New York, NY
  17     >>> lat2,lng2 = (49.261226, -123.1139268)   # Vancouver, Canada
  18     >>> az12,az21,dist = g.inv(lng1,lat1,lng2,lat2)
  19     >>> az12,az21,dist
  20     (-59.10918706123901, 84.99453463527395, 3914198.2912370963)
  21
  22 which gives forward and back [azimuths][] as well as the geodesic
  23 distance in meters.  Not that longitude comes *before* latitude in the
  24 these pyproj argument lists.
  25
  26 Calculate the terminus of a geodesic from an initial point, azimuth,
  27 and distance with:
  28
  29     >>> lng3,lat3,az3 = g.fwd(lng1,lat1,az12, dist)
  30     >>> lat3,lng3,az3
  31     (49.26122600306212, -123.11392684861474, 84.99453467574762)
  32
  33 Plan your trip with:
  34
  35     >>> pts = g.npts(lng1,lat1,lng2,lat2,npts=5)
  36     >>> pts.insert(0, (lng1, lat1))
  37     >>> pts.append((lng2, lat2))
  38     >>> import numpy
  39     >>> npts = numpy.array(pts)
  40     >>> npts
  41     array([[ -74.0059731 ,   40.7143528 ],
  42            [ -80.93566289,   43.52686057],
  43            [ -88.48167748,   45.87969433],
  44            [ -96.61187851,   47.6930911 ],
  45            [-105.22271807,   48.89347605],
  46            [-114.13503215,   49.42510006],
  47            [-123.1139268 ,   49.261226  ]])
  48
  49 To plot the above New York to Vancouver route on a flat map, we need a
  50 `Proj` instance:
  51
  52     >>> from pyproj import Proj
  53     >>> awips221 = Proj(proj='lcc', R=6371200, lat_1=50, lat_2=50,
  54     ...     lon_0=-107, ellps='clrk66')
  55     >>> awips218 = Proj(proj='lcc', R=6371200, lat_1=25, lat_2=25,
  56     ...     lon_0=-95, ellps='clrk66')  #x_0=-llcrnrx,y_0=-llcrnry)
  57
  58     #llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower
  59     #    left hand corner of projection region.
  60
  61 where `proj='lcc'` selects the [Lambert conformal conic][lcc]
  62 projection for the x/y points, and `ellps='clrk66'` selects the
  63 reference ellipsoid for the lat/lng coordinates.  The other
  64 coordinates are LCC parameters that select the [AWIPS 221][awips221]
  65 and [AWIPS 226][awips226] projections respectively (`lat_1`
  66 corresponds to `Latin1`, `lat_2` corresponds to `Latin2`, and `lon_0`
  67 corresponds to `Lov`; see [this description][lcc-param] of the
  68 two-standard-parallel LCC and its PROJ.4 parameters).
  69
  70 Convert our lat/lng pairs into grid points:
  71
  72     >>> awips221(lng1, lat1)
  73     (2725283.842678774, 5823260.730665273)
  74     >>> x221,y221 = awips221(npts[:,0], npts[:,1])
  75     >>> # xy221 = numpy.concatenate((a1, a2, ...), axis=0)  # numpy-2.0
  76     >>> xy221 = numpy.ndarray(shape=npts.shape, dtype=npts.dtype)
  77     >>> xy221[:,0] = x221
  78     >>> xy221[:,1] = y221
  79     >>> xy221
  80     array([[ 2725283.84267877,  5823260.73066527],
  81            [ 2071820.3526011 ,  5892518.49630526],
  82            [ 1422529.71760395,  5967565.49899035],
  83            [  775650.03731228,  6046475.43928965],
  84            [  129946.46495299,  6127609.80532071],
  85            [ -515306.57275941,  6209785.69230076],
  86            [-1160447.80254759,  6292455.41884832]])
  87
  88 Finally, you can convert points from one projection to another.
  89
  90     >>> from pyproj import transform
  91     >>> x218,y218 = transform(awips221, awips218, x221, y221)
  92     >>> xy218 = numpy.ndarray(shape=npts.shape, dtype=npts.dtype)
  93     >>> xy218[:,0] = x218
  94     >>> xy218[:,1] = y218
  95     >>> xy218
  96     array([[ 1834251.59591561,  4780900.70184736],
  97            [ 1197541.13209718,  5028862.9881648 ],
  98            [  542391.04388716,  5258740.71523961],
  99            [ -131577.34942316,  5464828.45934687],
 100            [ -822685.42269932,  5641393.59760613],
 101            [-1527077.85176048,  5783597.16169582],
 102            [-2239159.34620498,  5888495.91009021]])
 103
 104 Another useful coordinate system is the [Universal Transverse
 105 Mercator][UTM] projection which slices the world into [zones][].
 106
 107     >>> p = Proj(proj='utm', zone=10, ellps='clrk66')
 108
 109 Putting everything together, here's a route map based on digital
 110 lat/lng pairs stored in a text file:
 111
 112     >>> from numpy import array
 113     >>> from pylab import plot, show
 114     >>> from pyproj import Geod, Proj
 115     >>> latlng = array([[float(x) for x in ln.split()]
 116     ...                for ln in open('coords', 'r')
 117     ...                if not ln.startswith('#')])
 118     >>> g = Geod(ellps='WGS84')
 119     >>> az12s,az21s,dists = g.inv(latlng[:-1,1], latlng[:-1,0],
 120     ...                           latlng[1:,1], latlng[1:,0])
 121     >>> print('total distance: %g m' % dists.sum())
 122     total distance: 2078.93 m
 123     >>> mlng = latlng[:,1].mean()
 124     >>> zone = int(round((mlng + 180) / 6.))
 125     >>> p = Proj(proj='utm', zone=zone, ellps='WGS84')
 126     >>> xs,ys = p(latlng[:,1], latlng[:,0])
 127     >>> lines = plot(xs, ys, 'r.-')
 128     >>> show()
 129
 130 I've written up a simple script using this approach:
 131 [[geoscript/maproute.py]].  I've also written up a simple script to
 132 draw a map with labeled points: [[geoscript/maplabel.py]].
 133
 134 Note that you can easily get lat/lng pairs using [geopy][] (ebuild in
 135 my [[Gentoo overlay]]):
 136
 137     >>> import geopy
 138     >>> g = geopy.geocoders.Google()
 139     >>> place1,(lat1,lng1) = g.geocode("New York, NY")
 140     >>> place2,(lat2,lng2) = g.geocode("Vancouver, Canada")
 141     >>> place1,(lat1,lng1)
 142     (u'New York, NY, USA', (40.7143528, -74.0059731))
 143     >>> place2,(lat2,lng2)
 144     (u'Vancouver, BC, Canada', (49.261226, -123.1139268))
 145
 146 If you're looking for a more compact [[C++]] package for geographic
 147 conversions, [GeographicLib][] looks promising.
 148
 149 [pyproj]: http://code.google.com/p/pyproj/
 150 [PROJ.4]: http://trac.osgeo.org/proj/
 151 [geodetic]: http://en.wikipedia.org/wiki/Geodesy
 152 [clrk66]: http://en.wikipedia.org/wiki/Alexander_Ross_Clarke
 153 [rell]: http://en.wikipedia.org/wiki/Reference_ellipsoid
 154 [azimuths]: http://en.wikipedia.org/wiki/Azimuth
 155 [LCC]: http://en.wikipedia.org/wiki/Lambert_conformal_conic_projection
 156 [awips221]: http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID221
 157 [awips226]: http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID218
 158 [lcc-param]: http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html
 159 [UTM]: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
 160 [zones]: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#UTM_zone
 161 [geopy]: http://code.google.com/p/geopy/
 162 [GeographicLib]: http://geographiclib.sourceforge.net/html/
 163
 164 [[!tag tags/fun]]
 165 [[!tag tags/linux]]
 166 [[!tag tags/python]]
 167 [[!tag tags/tools]]