From: W. Trevor King Date: Thu, 19 May 2011 12:27:32 +0000 (-0400) Subject: Add pyproj post. X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=8bdee63103ad2be9c56c41cf60f9a242b41110c3;p=blog.git Add pyproj post. --- diff --git a/posts/pyproj.mdwn b/posts/pyproj.mdwn new file mode 100644 index 0000000..9a10e40 --- /dev/null +++ b/posts/pyproj.mdwn @@ -0,0 +1,160 @@ +[pyproj][] is a Python wrapper around [PROJ.4][]. Here's a quick +walkthrough. + +Initialize a [geodetic][] converter: + + >>> from pyproj import Geod + >>> g = Geod(ellps='clrk66') + +where `ellps='clrk66'` selects [Clarke's 1866][clrk66] [reference +ellipsoid][rell]. `help(Geod.__new__)` gives a list of possible +ellipsoids. + +Calculate the distance between two points, as well as the local +heading, try + + >>> lat1,lng1 = (40.7143528, -74.0059731) # New York, NY + >>> lat2,lng2 = (49.261226, -123.1139268) # Vancouver, Canada + >>> az12,az21,dist = g.inv(lng1,lat1,lng2,lat2) + >>> az12,az21,dist + (-59.10918706123901, 84.99453463527395, 3914198.2912370963) + +which gives forward and back [azimuths][] as well as the geodesic +distance in meters. Not that longitude comes *before* latitude in the +these pyproj argument lists. + +Calculate the terminus of a geodesic from an initial point, azimuth, +and distance with: + + >>> lng3,lat3,az3 = g.fwd(lng1,lat1,az12, dist) + >>> lat3,lng3,az3 + (49.26122600306212, -123.11392684861474, 84.99453467574762) + +Plan your trip with: + + >>> pts = g.npts(lng1,lat1,lng2,lat2,npts=5) + >>> pts.insert(0, (lng1, lat1)) + >>> pts.append((lng2, lat2)) + >>> import numpy + >>> npts = numpy.array(pts) + >>> npts + array([[ -74.0059731 , 40.7143528 ], + [ -80.93566289, 43.52686057], + [ -88.48167748, 45.87969433], + [ -96.61187851, 47.6930911 ], + [-105.22271807, 48.89347605], + [-114.13503215, 49.42510006], + [-123.1139268 , 49.261226 ]]) + +To plot the above New York to Vancouver route on a flat map, we need a +`Proj` instance: + + >>> from pyproj import Proj + >>> awips221 = Proj(proj='lcc', R=6371200, lat_1=50, lat_2=50, + ... lon_0=-107, ellps='clrk66') + >>> awips218 = Proj(proj='lcc', R=6371200, lat_1=25, lat_2=25, + ... lon_0=-95, ellps='clrk66') #x_0=-llcrnrx,y_0=-llcrnry) + + #llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower + # left hand corner of projection region. + +where `proj='lcc` selects the [Lambert conformal conic][lcc] +projection for the x/y points, and `ellps='clrk66'` selects the +reference ellipsoid for the lat/lng coordinates. The other +coordinates are LCC parameters that select the [AWIPS 221][awips221] +and [AWIPS 226][awips226] projections respectively (`lat_1` +corresponds to `Latin1`, `lat_2` corresponds to `Latin2`, and `lon_0` +corresponds to `Lov`). + +Convert our lat/lng pairs into grid points: + + >>> awips221(lng1, lat1) + (2725283.842678774, 5823260.730665273) + >>> x221,y221 = awips221(npts[:,0], npts[:,1]) + >>> # xy221 = numpy.concatenate((a1, a2, ...), axis=0) # numpy-2.0 + >>> xy221 = numpy.ndarray(shape=npts.shape, dtype=npts.dtype) + >>> xy221[:,0] = x221 + >>> xy221[:,1] = y221 + >>> xy221 + array([[ 2725283.84267877, 5823260.73066527], + [ 2071820.3526011 , 5892518.49630526], + [ 1422529.71760395, 5967565.49899035], + [ 775650.03731228, 6046475.43928965], + [ 129946.46495299, 6127609.80532071], + [ -515306.57275941, 6209785.69230076], + [-1160447.80254759, 6292455.41884832]]) + +Finally, you can convert points from one projection to another. + + >>> from pyproj import transform + >>> x218,y218 = transform(awips221, awips218, x221, y221) + >>> xy218 = numpy.ndarray(shape=npts.shape, dtype=npts.dtype) + >>> xy218[:,0] = x218 + >>> xy218[:,1] = y218 + >>> xy218 + array([[ 1834251.59591561, 4780900.70184736], + [ 1197541.13209718, 5028862.9881648 ], + [ 542391.04388716, 5258740.71523961], + [ -131577.34942316, 5464828.45934687], + [ -822685.42269932, 5641393.59760613], + [-1527077.85176048, 5783597.16169582], + [-2239159.34620498, 5888495.91009021]]) + +Another useful coordinate system is the [Universal Transverse +Mercator][UTM] projection which slices the world into [zones][]. + + >>> p = Proj(proj='utm', zone=10, ellps='clrk66') + +Putting everything together, here's a route map based on digital +lat/lng pairs stored in a text file: + + >>> from numpy import array + >>> from pylab import plot, show + >>> from pyproj import Geod, Proj + >>> latlng = array([[float(x) for x in ln.split()] + ... for ln in open('coords', 'r') + ... if not ln.startswith('#')]) + >>> g = Geod(ellps='WGS84') + >>> az12s,az21s,dists = g.inv(latlng[:-1,1], latlng[:-1,0], + ... latlng[1:,1], latlng[1:,0]) + >>> print('total distance: %g m' % dists.sum()) + total distance: 2078.93 m + >>> mlng = latlng[:,1].mean() + >>> zone = int(round((mlng + 180) / 6.)) + >>> p = Proj(proj='utm', zone=zone, ellps='WGS84') + >>> xs,ys = p(latlng[:,1], latlng[:,0]) + >>> lines = plot(xs, ys, 'r.-') + >>> show() + +Note that you can easily get lat/lng pairs using [geopy][]: + + >>> import geopy + >>> g = geopy.geocoders.Google() + >>> place1,(lat1,lng1) = g.geocode("New York, NY") + >>> place2,(lat2,lng2) = g.geocode("Vancouver, Canada") + >>> place1,(lat1,lng1) + (u'New York, NY, USA', (40.7143528, -74.0059731)) + >>> place2,(lat2,lng2) + (u'Vancouver, BC, Canada', (49.261226, -123.1139268)) + +If you're looking for a more compact [[C++]] package for geographic +conversions, [GeographicLib][] looks promising. + +[pyproj]: http://code.google.com/p/pyproj/ +[PROJ.4]: http://trac.osgeo.org/proj/ +[geodetic]: http://en.wikipedia.org/wiki/Geodesy +[clrk66]: http://en.wikipedia.org/wiki/Alexander_Ross_Clarke +[rell]: http://en.wikipedia.org/wiki/Reference_ellipsoid +[azimuths]: http://en.wikipedia.org/wiki/Azimuth +[LCC]: http://en.wikipedia.org/wiki/Lambert_conformal_conic_projection +[awips221]: http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID221 +[awips226]: http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID218 +[UTM]: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system +[zone]: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#UTM_zone +[geopy]: http://code.google.com/p/geopy/ +[GeographicLib]: http://geographiclib.sourceforge.net/html/ + +[[!tag tags/fun]] +[[!tag tags/linux]] +[[!tag tags/python]] +[[!tag tags/tools]]