1 [pyproj][] is a Python wrapper around [PROJ.4][]. Here's a quick
4 Initialize a [geodetic][] converter:
6 >>> from pyproj import Geod
7 >>> g = Geod(ellps='clrk66')
9 where `ellps='clrk66'` selects [Clarke's 1866][clrk66] [reference
10 ellipsoid][rell]. `help(Geod.__new__)` gives a list of possible
13 Calculate the distance between two points, as well as the local
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)
20 (-59.10918706123901, 84.99453463527395, 3914198.2912370963)
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.
26 Calculate the terminus of a geodesic from an initial point, azimuth,
29 >>> lng3,lat3,az3 = g.fwd(lng1,lat1,az12, dist)
31 (49.26122600306212, -123.11392684861474, 84.99453467574762)
35 >>> pts = g.npts(lng1,lat1,lng2,lat2,npts=5)
36 >>> pts.insert(0, (lng1, lat1))
37 >>> pts.append((lng2, lat2))
39 >>> npts = numpy.array(pts)
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 ]])
49 To plot the above New York to Vancouver route on a flat map, we need a
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)
58 #llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower
59 # left hand corner of projection region.
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).
70 Convert our lat/lng pairs into grid points:
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)
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]])
88 Finally, you can convert points from one projection to another.
90 >>> from pyproj import transform
91 >>> x218,y218 = transform(awips221, awips218, x221, y221)
92 >>> xy218 = numpy.ndarray(shape=npts.shape, dtype=npts.dtype)
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]])
104 Another useful coordinate system is the [Universal Transverse
105 Mercator][UTM] projection which slices the world into [zones][].
107 >>> p = Proj(proj='utm', zone=10, ellps='clrk66')
109 Putting everything together, here's a route map based on digital
110 lat/lng pairs stored in a text file:
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.-')
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]].
134 Note that you can easily get lat/lng pairs using [geopy][] (ebuild in
135 my [[Gentoo overlay]]):
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))
146 If you're looking for a more compact [[C++]] package for geographic
147 conversions, [GeographicLib][] looks promising.
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/