1 """
2 Cython wrapper to provide python interfaces to
3 PROJ.4 (http://trac.osgeo.org/proj/) functions.
4
5 Performs cartographic transformations and geodetic computations.
6
7 The Proj class can convert from geographic (longitude,latitude)
8 to native map projection (x,y) coordinates and vice versa, or
9 from one map projection coordinate system directly to another.
10 The module variable pj_list is a dictionary containing all the
11 available projections and their descriptions.
12
13 The Geod class can perform forward and inverse geodetic, or
14 Great Circle, computations. The forward computation involves
15 determining latitude, longitude and back azimuth of a terminus
16 point given the latitude and longitude of an initial point, plus
17 azimuth and distance. The inverse computation involves
18 determining the forward and back azimuths and distance given the
19 latitudes and longitudes of an initial and terminus point.
20
21 Input coordinates can be given as python arrays, lists/tuples,
22 scalars or numpy/Numeric/numarray arrays. Optimized for objects
23 that support the Python buffer protocol (regular python and
24 numpy array objects).
25
26 Download: http://code.google.com/p/pyproj/downloads/list
27
28 Requirements: python 2.4 or higher.
29
30 Example scripts are in 'test' subdirectory of source distribution.
31 The 'test()' function will run the examples in the docstrings.
32
33 Contact: Jeffrey Whitaker <jeffrey.s.whitaker@noaa.gov
34
35 copyright (c) 2006 by Jeffrey Whitaker.
36
37 Permission to use, copy, modify, and distribute this software
38 and its documentation for any purpose and without fee is hereby
39 granted, provided that the above copyright notice appear in all
40 copies and that both the copyright notice and this permission
41 notice appear in supporting documentation. THE AUTHOR DISCLAIMS
42 ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL
43 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT
44 SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT OR
45 CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
46 LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
47 NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
48 CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. """
49
50 from pyproj import _proj
51 from pyproj.geodesic import Geodesic
52 __version__ = _proj.__version__
53 set_datapath = _proj.set_datapath
54 from array import array
55 import os, math
56
57 pj_list={
58 'aea': "Albers Equal Area",
59 'aeqd': "Azimuthal Equidistant",
60 'airy': "Airy",
61 'aitoff': "Aitoff",
62 'alsk': "Mod. Stererographics of Alaska",
63 'apian': "Apian Globular I",
64 'august': "August Epicycloidal",
65 'bacon': "Bacon Globular",
66 'bipc': "Bipolar conic of western hemisphere",
67 'boggs': "Boggs Eumorphic",
68 'bonne': "Bonne (Werner lat_1=90)",
69 'cass': "Cassini",
70 'cc': "Central Cylindrical",
71 'cea': "Equal Area Cylindrical",
72 'chamb': "Chamberlin Trimetric",
73 'collg': "Collignon",
74 'crast': "Craster Parabolic (Putnins P4)",
75 'denoy': "Denoyer Semi-Elliptical",
76 'eck1': "Eckert I",
77 'eck2': "Eckert II",
78 'eck3': "Eckert III",
79 'eck4': "Eckert IV",
80 'eck5': "Eckert V",
81 'eck6': "Eckert VI",
82 'eqc': "Equidistant Cylindrical (Plate Caree)",
83 'eqdc': "Equidistant Conic",
84 'etmerc': "Extended Transverse Mercator" ,
85 'euler': "Euler",
86 'fahey': "Fahey",
87 'fouc': "Foucaut",
88 'fouc_s': "Foucaut Sinusoidal",
89 'gall': "Gall (Gall Stereographic)",
90 'geocent': "Geocentric",
91 'geos': "Geostationary Satellite View",
92 'gins8': "Ginsburg VIII (TsNIIGAiK)",
93 'gn_sinu': "General Sinusoidal Series",
94 'gnom': "Gnomonic",
95 'goode': "Goode Homolosine",
96 'gs48': "Mod. Stererographics of 48 U.S.",
97 'gs50': "Mod. Stererographics of 50 U.S.",
98 'hammer': "Hammer & Eckert-Greifendorff",
99 'hatano': "Hatano Asymmetrical Equal Area",
100 'healpix': "HEALPix",
101 'rhealpix': "rHEALPix",
102 'igh': "Interrupted Goode Homolosine",
103 'imw_p': "Internation Map of the World Polyconic",
104 'kav5': "Kavraisky V",
105 'kav7': "Kavraisky VII",
106 'krovak': "Krovak",
107 'labrd': "Laborde",
108 'laea': "Lambert Azimuthal Equal Area",
109 'lagrng': "Lagrange",
110 'larr': "Larrivee",
111 'lask': "Laskowski",
112 'lonlat': "Lat/long (Geodetic)",
113 'latlon': "Lat/long (Geodetic alias)",
114 'latlong': "Lat/long (Geodetic alias)",
115 'longlat': "Lat/long (Geodetic alias)",
116 'lcc': "Lambert Conformal Conic",
117 'lcca': "Lambert Conformal Conic Alternative",
118 'leac': "Lambert Equal Area Conic",
119 'lee_os': "Lee Oblated Stereographic",
120 'loxim': "Loximuthal",
121 'lsat': "Space oblique for LANDSAT",
122 'mbt_s': "McBryde-Thomas Flat-Polar Sine",
123 'mbt_fps': "McBryde-Thomas Flat-Pole Sine (No. 2)",
124 'mbtfpp': "McBride-Thomas Flat-Polar Parabolic",
125 'mbtfpq': "McBryde-Thomas Flat-Polar Quartic",
126 'mbtfps': "McBryde-Thomas Flat-Polar Sinusoidal",
127 'merc': "Mercator",
128 'mil_os': "Miller Oblated Stereographic",
129 'mill': "Miller Cylindrical",
130 'moll': "Mollweide",
131 'murd1': "Murdoch I",
132 'murd2': "Murdoch II",
133 'murd3': "Murdoch III",
134 'nell': "Nell",
135 'nell_h': "Nell-Hammer",
136 'nicol': "Nicolosi Globular",
137 'nsper': "Near-sided perspective",
138 'nzmg': "New Zealand Map Grid",
139 'ob_tran': "General Oblique Transformation",
140 'ocea': "Oblique Cylindrical Equal Area",
141 'oea': "Oblated Equal Area",
142 'omerc': "Oblique Mercator",
143 'ortel': "Ortelius Oval",
144 'ortho': "Orthographic",
145 'pconic': "Perspective Conic",
146 'poly': "Polyconic (American)",
147 'putp1': "Putnins P1",
148 'putp2': "Putnins P2",
149 'putp3': "Putnins P3",
150 'putp3p': "Putnins P3'",
151 'putp4p': "Putnins P4'",
152 'putp5': "Putnins P5",
153 'putp5p': "Putnins P5'",
154 'putp6': "Putnins P6",
155 'putp6p': "Putnins P6'",
156 'qua_aut': "Quartic Authalic",
157 'robin': "Robinson",
158 'rouss': "Roussilhe Stereographic",
159 'rpoly': "Rectangular Polyconic",
160 'sinu': "Sinusoidal (Sanson-Flamsteed)",
161 'somerc': "Swiss. Obl. Mercator",
162 'stere': "Stereographic",
163 'sterea': "Oblique Stereographic Alternative",
164 'gstmerc': "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)",
165 'tcc': "Transverse Central Cylindrical",
166 'tcea': "Transverse Cylindrical Equal Area",
167 'tissot': "Tissot Conic",
168 'tmerc': "Transverse Mercator",
169 'tpeqd': "Two Point Equidistant",
170 'tpers': "Tilted perspective",
171 'ups': "Universal Polar Stereographic",
172 'urm5': "Urmaev V",
173 'urmfps': "Urmaev Flat-Polar Sinusoidal",
174 'utm': "Universal Transverse Mercator (UTM)",
175 'vandg': "van der Grinten (I)",
176 'vandg2': "van der Grinten II",
177 'vandg3': "van der Grinten III",
178 'vandg4': "van der Grinten IV",
179 'vitk1': "Vitkovsky I",
180 'wag1': "Wagner I (Kavraisky VI)",
181 'wag2': "Wagner II",
182 'wag3': "Wagner III",
183 'wag4': "Wagner IV",
184 'wag5': "Wagner V",
185 'wag6': "Wagner VI",
186 'wag7': "Wagner VII",
187 'weren': "Werenskiold I",
188 'wink1': "Winkel I",
189 'wink2': "Winkel II",
190 'wintri': "Winkel Tripel"}
191
192 pj_ellps={
193 "MERIT": {'a':6378137.0,'rf':298.257,'description':"MERIT 1983"},
194 "SGS85": {'a':6378136.0,'rf':298.257,'description':"Soviet Geodetic System 85"},
195 "GRS80": {'a':6378137.0,'rf':298.257222101,'description':"GRS 1980(IUGG, 1980)"},
196 "IAU76": {'a':6378140.0,'rf':298.257,'description':"IAU 1976"},
197 "airy": {'a':6377563.396,'b':6356256.910,'description':"Airy 1830"},
198 "APL4.9": {'a':6378137.0,'rf':298.25,'description':"Appl. Physics. 1965"},
199 "NWL9D": {'a':6378145.0,'rf':298.25,'description':" Naval Weapons Lab., 1965"},
200 "mod_airy": {'a':6377340.189,'b':6356034.446,'description':"Modified Airy"},
201 "andrae": {'a':6377104.43,'rf':300.0,'description':"Andrae 1876 (Den., Iclnd.)"},
202 "aust_SA": {'a':6378160.0,'rf':298.25,'description':"Australian Natl & S. Amer. 1969"},
203 "GRS67": {'a':6378160.0,'rf':298.2471674270,'description':"GRS 67(IUGG 1967)"},
204 "bessel": {'a':6377397.155,'rf':299.1528128,'description':"Bessel 1841"},
205 "bess_nam": {'a':6377483.865,'rf':299.1528128,'description':"Bessel 1841 (Namibia)"},
206 "clrk66": {'a':6378206.4,'b':6356583.8,'description':"Clarke 1866"},
207 "clrk80": {'a':6378249.145,'rf':293.4663,'description':"Clarke 1880 mod."},
208 "CPM": {'a':6375738.7,'rf':334.29,'description':"Comm. des Poids et Mesures 1799"},
209 "delmbr": {'a':6376428.,'rf':311.5,'description':"Delambre 1810 (Belgium)"},
210 "engelis": {'a':6378136.05,'rf':298.2566,'description':"Engelis 1985"},
211 "evrst30": {'a':6377276.345,'rf':300.8017,'description':"Everest 1830"},
212 "evrst48": {'a':6377304.063,'rf':300.8017,'description':"Everest 1948"},
213 "evrst56": {'a':6377301.243,'rf':300.8017,'description':"Everest 1956"},
214 "evrst69": {'a':6377295.664,'rf':300.8017,'description':"Everest 1969"},
215 "evrstSS": {'a':6377298.556,'rf':300.8017,'description':"Everest (Sabah & Sarawak)"},
216 "fschr60": {'a':6378166.,'rf':298.3,'description':"Fischer (Mercury Datum) 1960"},
217 "fschr60m": {'a':6378155.,'rf':298.3,'description':"Modified Fischer 1960"},
218 "fschr68": {'a':6378150.,'rf':298.3,'description':"Fischer 1968"},
219 "helmert": {'a':6378200.,'rf':298.3,'description':"Helmert 1906"},
220 "hough": {'a':6378270.0,'rf':297.,'description':"Hough"},
221 "intl": {'a':6378388.0,'rf':297.,'description':"International 1909 (Hayford)"},
222 "krass": {'a':6378245.0,'rf':298.3,'description':"Krassovsky, 1942"},
223 "kaula": {'a':6378163.,'rf':298.24,'description':"Kaula 1961"},
224 "lerch": {'a':6378139.,'rf':298.257,'description':"Lerch 1979"},
225 "mprts": {'a':6397300.,'rf':191.,'description':"Maupertius 1738"},
226 "new_intl": {'a':6378157.5,'b':6356772.2,'description':"New International 1967"},
227 "plessis": {'a':6376523.,'b':6355863.,'description':"Plessis 1817 (France)"},
228 "SEasia": {'a':6378155.0,'b':6356773.3205,'description':"Southeast Asia"},
229 "walbeck": {'a':6376896.0,'b':6355834.8467,'description':"Walbeck"},
230 "WGS60": {'a':6378165.0,'rf':298.3,'description':"WGS 60"},
231 "WGS66": {'a':6378145.0,'rf':298.25,'description':"WGS 66"},
232 "WGS72": {'a':6378135.0,'rf':298.26,'description':"WGS 72"},
233 "WGS84": {'a':6378137.0,'rf':298.257223563,'description':"WGS 84"},
234 "sphere": {'a':6370997.0,'b':6370997.0,'description':"Normal Sphere"},
235 }
236
237 pyproj_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
238 if not os.path.isdir(pyproj_datadir):
239 msg="proj data directory not found. Expecting it at: %s"%pyproj_datadir
240 raise IOError(msg)
241
242 set_datapath(pyproj_datadir)
243
244 -class Proj(_proj.Proj):
245 """
246 performs cartographic transformations (converts from
247 longitude,latitude to native map projection x,y coordinates and
248 vice versa) using proj (http://trac.osgeo.org/proj/).
249
250 A Proj class instance is initialized with proj map projection
251 control parameter key/value pairs. The key/value pairs can
252 either be passed in a dictionary, or as keyword arguments,
253 or as a proj4 string (compatible with the proj command). See
254 http://www.remotesensing.org/geotiff/proj_list for examples of
255 key/value pairs defining different map projections.
256
257 Calling a Proj class instance with the arguments lon, lat will
258 convert lon/lat (in degrees) to x/y native map projection
259 coordinates (in meters). If optional keyword 'inverse' is True
260 (default is False), the inverse transformation from x/y to
261 lon/lat is performed. If optional keyword 'radians' is True
262 (default is False) lon/lat are interpreted as radians instead of
263 degrees. If optional keyword 'errcheck' is True (default is
264 False) an exception is raised if the transformation is invalid.
265 If errcheck=False and the transformation is invalid, no
266 exception is raised and 1.e30 is returned. If the optional keyword
267 'preserve_units' is True, the units in map projection coordinates
268 are not forced to be meters.
269
270 Works with numpy and regular python array objects, python
271 sequences and scalars.
272 """
273
274 - def __new__(self, projparams=None, preserve_units=False, **kwargs):
275 """
276 initialize a Proj class instance.
277
278 Proj4 projection control parameters must either be given in a
279 dictionary 'projparams' or as keyword arguments. See the proj
280 documentation (http://trac.osgeo.org/proj/) for more information
281 about specifying projection parameters.
282
283 Example usage:
284
285 >>> from pyproj import Proj
286 >>> p = Proj(proj='utm',zone=10,ellps='WGS84') # use kwargs
287 >>> x,y = p(-120.108, 34.36116666)
288 >>> 'x=%9.3f y=%11.3f' % (x,y)
289 'x=765975.641 y=3805993.134'
290 >>> 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True)
291 'lon=-120.108 lat=34.361'
292 >>> # do 3 cities at a time in a tuple (Fresno, LA, SF)
293 >>> lons = (-119.72,-118.40,-122.38)
294 >>> lats = (36.77, 33.93, 37.62 )
295 >>> x,y = p(lons, lats)
296 >>> 'x: %9.3f %9.3f %9.3f' % x
297 'x: 792763.863 925321.537 554714.301'
298 >>> 'y: %9.3f %9.3f %9.3f' % y
299 'y: 4074377.617 3763936.941 4163835.303'
300 >>> lons, lats = p(x, y, inverse=True) # inverse transform
301 >>> 'lons: %8.3f %8.3f %8.3f' % lons
302 'lons: -119.720 -118.400 -122.380'
303 >>> 'lats: %8.3f %8.3f %8.3f' % lats
304 'lats: 36.770 33.930 37.620'
305 >>> p2 = Proj('+proj=utm +zone=10 +ellps=WGS84') # use proj4 string
306 >>> x,y = p2(-120.108, 34.36116666)
307 >>> 'x=%9.3f y=%11.3f' % (x,y)
308 'x=765975.641 y=3805993.134'
309 >>> p = Proj(init="epsg:32667")
310 >>> 'x=%12.3f y=%12.3f (meters)' % p(-114.057222, 51.045)
311 'x=-1783486.760 y= 6193833.196 (meters)'
312 >>> p = Proj("+init=epsg:32667",preserve_units=True)
313 >>> 'x=%12.3f y=%12.3f (feet)' % p(-114.057222, 51.045)
314 'x=-5851322.810 y=20320934.409 (feet)'
315 """
316
317 if projparams is None:
318 if len(kwargs) == 0:
319 raise RuntimeError('no projection control parameters specified')
320 else:
321 projstring = _dict2string(kwargs)
322 elif type(projparams) == str:
323
324 projstring = projparams
325 else:
326 projstring = _dict2string(projparams)
327
328 if not projstring.count('+units=') and not preserve_units:
329 projstring = '+units=m '+projstring
330 else:
331 kvpairs = []
332 for kvpair in projstring.split():
333 if kvpair.startswith('+units') and not preserve_units:
334 k,v = kvpair.split('=')
335 kvpairs.append(k+'=m ')
336 else:
337 kvpairs.append(kvpair+' ')
338 projstring = ''.join(kvpairs)
339
340
341 projstring = projstring.replace('EPSG','epsg')
342 return _proj.Proj.__new__(self, projstring)
343
345
346 """
347 Calling a Proj class instance with the arguments lon, lat will
348 convert lon/lat (in degrees) to x/y native map projection
349 coordinates (in meters). If optional keyword 'inverse' is True
350 (default is False), the inverse transformation from x/y to
351 lon/lat is performed. If optional keyword 'radians' is True
352 (default is False) the units of lon/lat are radians instead of
353 degrees. If optional keyword 'errcheck' is True (default is
354 False) an exception is raised if the transformation is invalid.
355 If errcheck=False and the transformation is invalid, no
356 exception is raised and 1.e30 is returned.
357
358 Instead of calling with lon, lat, a single ndarray of
359 shape n,2 may be used, and one of the same shape will
360 be returned; this is more efficient.
361
362 Inputs should be doubles (they will be cast to doubles if they
363 are not, causing a slight performance hit).
364
365 Works with numpy and regular python array objects, python
366 sequences and scalars, but is fastest for array objects.
367 """
368 inverse = kw.get('inverse', False)
369 radians = kw.get('radians', False)
370 errcheck = kw.get('errcheck', False)
371
372
373
374
375
376
377
378
379 lon, lat = args
380
381 inx, xisfloat, xislist, xistuple = _copytobuffer(lon)
382 iny, yisfloat, yislist, yistuple = _copytobuffer(lat)
383
384 if inverse:
385 _proj.Proj._inv(self, inx, iny, radians=radians, errcheck=errcheck)
386 else:
387 _proj.Proj._fwd(self, inx, iny, radians=radians, errcheck=errcheck)
388
389 outx = _convertback(xisfloat,xislist,xistuple,inx)
390 outy = _convertback(yisfloat,yislist,xistuple,iny)
391 return outx, outy
392
394 """returns True if projection in geographic (lon/lat) coordinates"""
395 return _proj.Proj.is_latlong(self)
396
398 """returns True if projection in geocentric (x/y) coordinates"""
399 return _proj.Proj.is_geocent(self)
400
490
492 """
493 return a copy of x as an object that supports the python Buffer
494 API (python array if input is float, list or tuple, numpy array
495 if input is a numpy array). returns copyofx, isfloat, islist,
496 istuple (islist is True if input is a list, istuple is true if
497 input is a tuple, isfloat is true if input is a float).
498 """
499
500 isfloat = False; islist = False; istuple = False
501
502
503 if hasattr(x,'shape') and x.shape == (): x = float(x)
504 try:
505
506
507
508 x.dtype.char
509 inx = x.astype('d')
510 except:
511 try:
512 x.typecode()
513 inx = x.astype('d')
514 except:
515
516 try:
517 x.typecode
518 inx = array('d',x)
519 except:
520
521
522 if type(x) == list:
523 inx = array('d',x)
524 islist = True
525
526 elif type(x) == tuple:
527 inx = array('d',x)
528 istuple = True
529
530 else:
531 try:
532 x = float(x)
533 inx = array('d',(x,))
534 isfloat = True
535 except:
536 raise TypeError('input must be an array, list, tuple or scalar')
537 return inx,isfloat,islist,istuple
538
540
541 if isfloat:
542 return inx[0]
543 elif islist:
544 return inx.tolist()
545 elif istuple:
546 return tuple(inx)
547 else:
548 return inx
549
551
552 pjargs = []
553 for key,value in projparams.items():
554 pjargs.append('+'+key+"="+str(value)+' ')
555 return ''.join(pjargs)
556
558 """
559 performs forward and inverse geodetic, or Great Circle,
560 computations. The forward computation (using the 'fwd' method)
561 involves determining latitude, longitude and back azimuth of a
562 computations. The forward computation (using the 'fwd' method)
563 involves determining latitude, longitude and back azimuth of a
564 terminus point given the latitude and longitude of an initial
565 point, plus azimuth and distance. The inverse computation (using
566 the 'inv' method) involves determining the forward and back
567 azimuths and distance given the latitudes and longitudes of an
568 initial and terminus point.
569 """
570 - def __init__(self, initstring=None, **kwargs):
571 """
572 initialize a Geod class instance.
573
574 Geodetic parameters for specifying the ellipsoid
575 can be given in a dictionary 'initparams', as keyword arguments,
576 or as as proj4 geod initialization string.
577 Following is a list of the ellipsoids that may be defined using the
578 'ellps' keyword (these are stored in the model variable pj_ellps)::
579
580 MERIT a=6378137.0 rf=298.257 MERIT 1983
581 SGS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85
582 GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980)
583 IAU76 a=6378140.0 rf=298.257 IAU 1976
584 airy a=6377563.396 b=6356256.910 Airy 1830
585 APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
586 airy a=6377563.396 b=6356256.910 Airy 1830
587 APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
588 NWL9D a=6378145.0. rf=298.25 Naval Weapons Lab., 1965
589 mod_airy a=6377340.189 b=6356034.446 Modified Airy
590 andrae a=6377104.43 rf=300.0 Andrae 1876 (Den., Iclnd.)
591 aust_SA a=6378160.0 rf=298.25 Australian Natl & S. Amer. 1969
592 GRS67 a=6378160.0 rf=298.247167427 GRS 67(IUGG 1967)
593 bessel a=6377397.155 rf=299.1528128 Bessel 1841
594 bess_nam a=6377483.865 rf=299.1528128 Bessel 1841 (Namibia)
595 clrk66 a=6378206.4 b=6356583.8 Clarke 1866
596 clrk80 a=6378249.145 rf=293.4663 Clarke 1880 mod.
597 CPM a=6375738.7 rf=334.29 Comm. des Poids et Mesures 1799
598 delmbr a=6376428. rf=311.5 Delambre 1810 (Belgium)
599 engelis a=6378136.05 rf=298.2566 Engelis 1985
600 evrst30 a=6377276.345 rf=300.8017 Everest 1830
601 evrst48 a=6377304.063 rf=300.8017 Everest 1948
602 evrst56 a=6377301.243 rf=300.8017 Everest 1956
603 evrst69 a=6377295.664 rf=300.8017 Everest 1969
604 evrstSS a=6377298.556 rf=300.8017 Everest (Sabah & Sarawak)
605 fschr60 a=6378166. rf=298.3 Fischer (Mercury Datum) 1960
606 fschr60m a=6378155. rf=298.3 Modified Fischer 1960
607 fschr68 a=6378150. rf=298.3 Fischer 1968
608 helmert a=6378200. rf=298.3 Helmert 1906
609 hough a=6378270.0 rf=297. Hough
610 helmert a=6378200. rf=298.3 Helmert 1906
611 hough a=6378270.0 rf=297. Hough
612 intl a=6378388.0 rf=297. International 1909 (Hayford)
613 krass a=6378245.0 rf=298.3 Krassovsky, 1942
614 kaula a=6378163. rf=298.24 Kaula 1961
615 lerch a=6378139. rf=298.257 Lerch 1979
616 mprts a=6397300. rf=191. Maupertius 1738
617 new_intl a=6378157.5 b=6356772.2 New International 1967
618 plessis a=6376523. b=6355863. Plessis 1817 (France)
619 SEasia a=6378155.0 b=6356773.3205 Southeast Asia
620 walbeck a=6376896.0 b=6355834.8467 Walbeck
621 WGS60 a=6378165.0 rf=298.3 WGS 60
622 WGS66 a=6378145.0 rf=298.25 WGS 66
623 WGS72 a=6378135.0 rf=298.26 WGS 72
624 WGS84 a=6378137.0 rf=298.257223563 WGS 84
625 sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)
626
627 The parameters of the ellipsoid may also be set directly using
628 the 'a' (semi-major or equatorial axis radius) keyword, and
629 any one of the following keywords: 'b' (semi-minor,
630 or polar axis radius), 'e' (eccentricity), 'es' (eccentricity
631 squared), 'f' (flattening), or 'rf' (reciprocal flattening).
632
633 See the proj documentation (http://trac.osgeo.org/proj/) for more
634
635 See the proj documentation (http://trac.osgeo.org/proj/) for more
636 information about specifying ellipsoid parameters (specifically,
637 the chapter 'Specifying the Earth's figure' in the main Proj
638 users manual).
639
640 Example usage:
641
642 >>> from pyproj import Geod
643 >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
644 >>> # specify the lat/lons of some cities.
645 >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
646 >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
647 >>> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
648 >>> london_lat = 51.+(32./60.); london_lon = -(5./60.)
649 >>> # compute forward and back azimuths, plus distance
650 >>> # between Boston and Portland.
651 >>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
652 >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist)
653 '-66.531 75.654 4164192.708'
654 >>> # compute latitude, longitude and back azimuth of Portland,
655 >>> # given Boston lat/lon, forward azimuth and distance to Portland.
656 >>> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
657 >>> "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
658 '45.517 -123.683 75.654'
659 >>> # compute the azimuths, distances from New York to several
660 >>> # cities (pass a list)
661 >>> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
662 >>> lons2 = [boston_lon, portland_lon, london_lon]
663 >>> lats2 = [boston_lat, portland_lat, london_lat]
664 >>> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
665 >>> for faz,baz,d in list(zip(az12,az21,dist)): "%7.3f %7.3f %9.3f" % (faz,baz,d)
666 ' 54.663 -123.448 288303.720'
667 '-65.463 79.342 4013037.318'
668 ' 51.254 -71.576 5579916.651'
669 >>> g2 = Geod('+ellps=clrk66') # use proj4 style initialization string
670 >>> az12,az21,dist = g2.inv(boston_lon,boston_lat,portland_lon,portland_lat)
671 >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist)
672 '-66.531 75.654 4164192.708'
673 """
674
675
676 ellpsd = {}
677 if initstring is not None:
678 for kvpair in initstring.split():
679 k,v = kvpair.split('=')
680 k = k.lstrip('+')
681 if k in ['a','b','rf','f','es','e']:
682 v = float(v)
683 ellpsd[k] = v
684
685 kwargs = dict(list(kwargs.items()) + list(ellpsd.items()))
686 self.sphere = False
687 if 'ellps' in kwargs:
688
689 ellps_dict = pj_ellps[kwargs['ellps']]
690 a = ellps_dict['a']
691 if ellps_dict['description']=='Normal Sphere':
692 self.sphere = True
693 if 'b' in ellps_dict:
694 b = ellps_dict['b']
695 es = 1. - (b * b) / (a * a)
696 f = (a - b)/a
697 elif 'rf' in ellps_dict:
698 f = 1./ellps_dict['rf']
699 b = a*(1. - f)
700 es = 1. - (b * b) / (a * a)
701 else:
702
703
704
705
706
707
708 a = kwargs['a']
709 if 'b' in kwargs:
710 b = kwargs['b']
711 es = 1. - (b * b) / (a * a)
712 f = (a - b)/a
713 elif 'rf' in kwargs:
714 f = 1./kwargs['rf']
715 b = a*(1. - f)
716 es = 1. - (b * b) / (a * a)
717 elif 'f' in kwargs:
718 f = kwargs['f']
719 b = a*(1. - f)
720 es = 1. - (b/a)**2
721 elif 'es' in kwargs:
722 es = kwargs['es']
723 b = math.sqrt(a**2 - es*a**2)
724 f = (a - b)/a
725 elif 'e' in kwargs:
726 es = kwargs['e']**2
727 b = math.sqrt(a**2 - es*a**2)
728 f = (a - b)/a
729 else:
730 b = a
731 f = 0.
732 es = 0.
733
734
735 if math.fabs(f) < 1.e-8: self.sphere = True
736 self.a = a
737 self.b = b
738 self.f = f
739 self.es = es
740 self.G = Geodesic(self.a, self.f)
741
742 - def fwd(self, lons, lats, az, dist, radians=False):
743 """
744 forward transformation - Returns longitudes, latitudes and back
745 azimuths of terminus points given longitudes (lons) and
746 latitudes (lats) of initial points, plus forward azimuths (az)
747 and distances (dist).
748 latitudes (lats) of initial points, plus forward azimuths (az)
749 and distances (dist).
750
751 Works with numpy and regular python array objects, python
752 sequences and scalars.
753
754 if radians=True, lons/lats and azimuths are radians instead of
755 degrees. Distances are in meters.
756 """
757
758 inx, xisfloat, xislist, xistuple = _copytobuffer(lons)
759 iny, yisfloat, yislist, yistuple = _copytobuffer(lats)
760 inz, zisfloat, zislist, zistuple = _copytobuffer(az)
761 ind, disfloat, dislist, distuple = _copytobuffer(dist)
762 n = 0
763 zipin = zip(inx,iny,inz,ind)
764 for lon,lat,az,dist in zipin:
765 result = self.G.Direct(lat, lon, az, dist)
766 inx[n] = result['lon2']
767 iny[n] = result['lat2']
768 inz[n] = result['azi2']
769 az = result['azi2']
770 if az > 0:
771 inz[n] = az-180.
772 elif az < 0:
773 inz[n] = az+180.
774 else:
775 inz[n] = az
776 n = n + 1
777
778 outx = _convertback(xisfloat,xislist,xistuple,inx)
779 outy = _convertback(yisfloat,yislist,xistuple,iny)
780 outz = _convertback(zisfloat,zislist,zistuple,inz)
781 outy = _convertback(yisfloat,yislist,xistuple,iny)
782 outz = _convertback(zisfloat,zislist,zistuple,inz)
783 return outx, outy, outz
784
785 - def inv(self,lons1,lats1,lons2,lats2,radians=False):
786 """
787 inverse transformation - Returns forward and back azimuths, plus
788 distances between initial points (specified by lons1, lats1) and
789 terminus points (specified by lons2, lats2).
790
791 Works with numpy and regular python array objects, python
792 sequences and scalars.
793
794 if radians=True, lons/lats and azimuths are radians instead of
795 degrees. Distances are in meters.
796 """
797
798 inx, xisfloat, xislist, xistuple = _copytobuffer(lons1)
799 iny, yisfloat, yislist, yistuple = _copytobuffer(lats1)
800 inz, zisfloat, zislist, zistuple = _copytobuffer(lons2)
801 ind, disfloat, dislist, distuple = _copytobuffer(lats2)
802 n = 0
803 zipin = zip(inx,iny,inz,ind)
804 for lon1,lat1,lon2,lat2 in zipin:
805 result = self.G.Inverse(lat1, lon1, lat2, lon2)
806 inx[n] = result['azi1']
807 az = result['azi2']
808 if az > 0:
809 iny[n] = az-180.
810 elif az < 0:
811 iny[n] = az+180.
812 else:
813 iny[n] = az
814 inz[n] = result['s12']
815 n = n + 1
816
817 outx = _convertback(xisfloat,xislist,xistuple,inx)
818 outy = _convertback(yisfloat,yislist,xistuple,iny)
819 outz = _convertback(zisfloat,zislist,zistuple,inz)
820 return outx, outy, outz
821
822 - def npts(self, lon1, lat1, lon2, lat2, npts, radians=False):
823 """
824 Given a single initial point and terminus point (specified by
825 python floats lon1,lat1 and lon2,lat2), returns a list of
826 longitude/latitude pairs describing npts equally spaced
827 intermediate points along the geodesic between the initial and
828 terminus points.
829
830 if radians=True, lons/lats are radians instead of degrees.
831
832 Example usage:
833
834 >>> from pyproj import Geod
835 >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
836 >>> # specify the lat/lons of Boston and Portland.
837 >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
838 >>> # specify the lat/lons of Boston and Portland.
839 >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
840 >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
841 >>> # find ten equally spaced points between Boston and Portland.
842 >>> lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
843 >>> for lon,lat in lonlats: '%6.3f %7.3f' % (lat, lon)
844 '43.528 -75.414'
845 '44.637 -79.883'
846 '45.565 -84.512'
847 '46.299 -89.279'
848 '46.830 -94.156'
849 '47.149 -99.112'
850 '47.251 -104.106'
851 '47.136 -109.100'
852 '46.805 -114.051'
853 '46.262 -118.924'
854 """
855 result = self.G.Inverse(lat1, lon1, lat2, lon2)
856 dist = result['s12']
857 az = result['azi1']
858
859 del_s = dist/(npts+1)
860
861 del_s = dist/(npts+1)
862
863 lats = ()
864 lons = ()
865
866 for i in range(1,npts+1):
867 S = i*del_s
868 result = self.G.Direct(lat1, lon1, az, S)
869 if radians:
870 lats = lats + (_dg2rad*result['lat2'],)
871 lons = lons + (_dg2rad*result['lon2'],)
872 else:
873 lats = lats + (result['lat2'],)
874 lons = lons + (result['lon2'],)
875 return list(zip(lons, lats))
876
878 """run the examples in the docstrings using the doctest module"""
879 import doctest, pyproj
880 doctest.testmod(pyproj,verbose=True)
881
882 if __name__ == "__main__": test()
883