1 """
2 Cython wrapper to provide python interfaces to
3 PROJ.4 (https://github.com/OSGeo/proj.4/wiki) 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://python.org/pypi/pyproj
27
28 Requirements: Python 2.6, 2.7, 3.2 or higher version.
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 import sys
51 from pyproj import _proj
52 from pyproj.datadir import pyproj_datadir
53 __version__ = _proj.__version__
54 proj_version_str = _proj.proj_version_str
55 geodesic_version_str = _proj.geodesic_version_str
56 set_datapath = _proj.set_datapath
57 from array import array
58 import os, math
59 from itertools import islice, chain
60 try:
61 from future_builtins import zip
62 except ImportError:
63 pass
64
65
66
67 if sys.version_info[0] == 2:
68 string_types = (basestring,)
69 else:
70 string_types = (str,)
71
72 pj_list={
73 'aea': "Albers Equal Area",
74 'aeqd': "Azimuthal Equidistant",
75 'airy': "Airy",
76 'aitoff': "Aitoff",
77 'alsk': "Mod. Stererographics of Alaska",
78 'apian': "Apian Globular I",
79 'august': "August Epicycloidal",
80 'bacon': "Bacon Globular",
81 'bipc': "Bipolar conic of western hemisphere",
82 'boggs': "Boggs Eumorphic",
83 'bonne': "Bonne (Werner lat_1=90)",
84 'cass': "Cassini",
85 'cc': "Central Cylindrical",
86 'cea': "Equal Area Cylindrical",
87 'chamb': "Chamberlin Trimetric",
88 'collg': "Collignon",
89 'comill': "Compact Miller",
90 'crast': "Craster Parabolic (Putnins P4)",
91 'denoy': "Denoyer Semi-Elliptical",
92 'eck1': "Eckert I",
93 'eck2': "Eckert II",
94 'eck3': "Eckert III",
95 'eck4': "Eckert IV",
96 'eck5': "Eckert V",
97 'eck6': "Eckert VI",
98 'eqc': "Equidistant Cylindrical (Plate Caree)",
99 'eqdc': "Equidistant Conic",
100 'etmerc': "Extended Transverse Mercator" ,
101 'euler': "Euler",
102 'fahey': "Fahey",
103 'fouc': "Foucaut",
104 'fouc_s': "Foucaut Sinusoidal",
105 'gall': "Gall (Gall Stereographic)",
106 'geocent': "Geocentric",
107 'geos': "Geostationary Satellite View",
108 'gins8': "Ginsburg VIII (TsNIIGAiK)",
109 'gn_sinu': "General Sinusoidal Series",
110 'gnom': "Gnomonic",
111 'goode': "Goode Homolosine",
112 'gs48': "Mod. Stererographics of 48 U.S.",
113 'gs50': "Mod. Stererographics of 50 U.S.",
114 'hammer': "Hammer & Eckert-Greifendorff",
115 'hatano': "Hatano Asymmetrical Equal Area",
116 'healpix': "HEALPix",
117 'rhealpix': "rHEALPix",
118 'igh': "Interrupted Goode Homolosine",
119 'imw_p': "Internation Map of the World Polyconic",
120 'isea': "Icosahedral Snyder Equal Area",
121 'kav5': "Kavraisky V",
122 'kav7': "Kavraisky VII",
123 'krovak': "Krovak",
124 'labrd': "Laborde",
125 'laea': "Lambert Azimuthal Equal Area",
126 'lagrng': "Lagrange",
127 'larr': "Larrivee",
128 'lask': "Laskowski",
129 'lonlat': "Lat/long (Geodetic)",
130 'latlon': "Lat/long (Geodetic alias)",
131 'latlong': "Lat/long (Geodetic alias)",
132 'longlat': "Lat/long (Geodetic alias)",
133 'lcc': "Lambert Conformal Conic",
134 'lcca': "Lambert Conformal Conic Alternative",
135 'leac': "Lambert Equal Area Conic",
136 'lee_os': "Lee Oblated Stereographic",
137 'loxim': "Loximuthal",
138 'lsat': "Space oblique for LANDSAT",
139 'mbt_s': "McBryde-Thomas Flat-Polar Sine",
140 'mbt_fps': "McBryde-Thomas Flat-Pole Sine (No. 2)",
141 'mbtfpp': "McBride-Thomas Flat-Polar Parabolic",
142 'mbtfpq': "McBryde-Thomas Flat-Polar Quartic",
143 'mbtfps': "McBryde-Thomas Flat-Polar Sinusoidal",
144 'merc': "Mercator",
145 'mil_os': "Miller Oblated Stereographic",
146 'mill': "Miller Cylindrical",
147 'moll': "Mollweide",
148 'murd1': "Murdoch I",
149 'murd2': "Murdoch II",
150 'murd3': "Murdoch III",
151 'natearth': "Natural Earth",
152 'natearth2': "Natural Earth 2",
153 'nell': "Nell",
154 'nell_h': "Nell-Hammer",
155 'nicol': "Nicolosi Globular",
156 'nsper': "Near-sided perspective",
157 'nzmg': "New Zealand Map Grid",
158 'ob_tran': "General Oblique Transformation",
159 'ocea': "Oblique Cylindrical Equal Area",
160 'oea': "Oblated Equal Area",
161 'omerc': "Oblique Mercator",
162 'ortel': "Ortelius Oval",
163 'ortho': "Orthographic",
164 'patterson': "Patterson Cylindrical",
165 'pconic': "Perspective Conic",
166 'poly': "Polyconic (American)",
167 'putp1': "Putnins P1",
168 'putp2': "Putnins P2",
169 'putp3': "Putnins P3",
170 'putp3p': "Putnins P3'",
171 'putp4p': "Putnins P4'",
172 'putp5': "Putnins P5",
173 'putp5p': "Putnins P5'",
174 'putp6': "Putnins P6",
175 'putp6p': "Putnins P6'",
176 'qua_aut': "Quartic Authalic",
177 'robin': "Robinson",
178 'rouss': "Roussilhe Stereographic",
179 'rpoly': "Rectangular Polyconic",
180 'sinu': "Sinusoidal (Sanson-Flamsteed)",
181 'somerc': "Swiss. Obl. Mercator",
182 'stere': "Stereographic",
183 'sterea': "Oblique Stereographic Alternative",
184 'gstmerc': "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)",
185 'tcc': "Transverse Central Cylindrical",
186 'tcea': "Transverse Cylindrical Equal Area",
187 'times': "Times",
188 'tissot': "Tissot Conic",
189 'tmerc': "Transverse Mercator",
190 'tpeqd': "Two Point Equidistant",
191 'tpers': "Tilted perspective",
192 'ups': "Universal Polar Stereographic",
193 'urm5': "Urmaev V",
194 'urmfps': "Urmaev Flat-Polar Sinusoidal",
195 'utm': "Universal Transverse Mercator (UTM)",
196 'vandg': "van der Grinten (I)",
197 'vandg2': "van der Grinten II",
198 'vandg3': "van der Grinten III",
199 'vandg4': "van der Grinten IV",
200 'vitk1': "Vitkovsky I",
201 'wag1': "Wagner I (Kavraisky VI)",
202 'wag2': "Wagner II",
203 'wag3': "Wagner III",
204 'wag4': "Wagner IV",
205 'wag5': "Wagner V",
206 'wag6': "Wagner VI",
207 'wag7': "Wagner VII",
208 'weren': "Werenskiold I",
209 'wink1': "Winkel I",
210 'wink2': "Winkel II",
211 'wintri': "Winkel Tripel"}
212
213 pj_ellps={
214 "MERIT": {'a':6378137.0,'rf':298.257,'description':"MERIT 1983"},
215 "SGS85": {'a':6378136.0,'rf':298.257,'description':"Soviet Geodetic System 85"},
216 "GRS80": {'a':6378137.0,'rf':298.257222101,'description':"GRS 1980(IUGG, 1980)"},
217 "IAU76": {'a':6378140.0,'rf':298.257,'description':"IAU 1976"},
218 "airy": {'a':6377563.396,'b':6356256.910,'description':"Airy 1830"},
219 "APL4.9": {'a':6378137.0,'rf':298.25,'description':"Appl. Physics. 1965"},
220 "NWL9D": {'a':6378145.0,'rf':298.25,'description':" Naval Weapons Lab., 1965"},
221 "mod_airy": {'a':6377340.189,'b':6356034.446,'description':"Modified Airy"},
222 "andrae": {'a':6377104.43,'rf':300.0,'description':"Andrae 1876 (Den., Iclnd.)"},
223 "aust_SA": {'a':6378160.0,'rf':298.25,'description':"Australian Natl & S. Amer. 1969"},
224 "GRS67": {'a':6378160.0,'rf':298.2471674270,'description':"GRS 67(IUGG 1967)"},
225 "bessel": {'a':6377397.155,'rf':299.1528128,'description':"Bessel 1841"},
226 "bess_nam": {'a':6377483.865,'rf':299.1528128,'description':"Bessel 1841 (Namibia)"},
227 "clrk66": {'a':6378206.4,'b':6356583.8,'description':"Clarke 1866"},
228 "clrk80": {'a':6378249.145,'rf':293.4663,'description':"Clarke 1880 mod."},
229 "CPM": {'a':6375738.7,'rf':334.29,'description':"Comm. des Poids et Mesures 1799"},
230 "delmbr": {'a':6376428.,'rf':311.5,'description':"Delambre 1810 (Belgium)"},
231 "engelis": {'a':6378136.05,'rf':298.2566,'description':"Engelis 1985"},
232 "evrst30": {'a':6377276.345,'rf':300.8017,'description':"Everest 1830"},
233 "evrst48": {'a':6377304.063,'rf':300.8017,'description':"Everest 1948"},
234 "evrst56": {'a':6377301.243,'rf':300.8017,'description':"Everest 1956"},
235 "evrst69": {'a':6377295.664,'rf':300.8017,'description':"Everest 1969"},
236 "evrstSS": {'a':6377298.556,'rf':300.8017,'description':"Everest (Sabah & Sarawak)"},
237 "fschr60": {'a':6378166.,'rf':298.3,'description':"Fischer (Mercury Datum) 1960"},
238 "fschr60m": {'a':6378155.,'rf':298.3,'description':"Modified Fischer 1960"},
239 "fschr68": {'a':6378150.,'rf':298.3,'description':"Fischer 1968"},
240 "helmert": {'a':6378200.,'rf':298.3,'description':"Helmert 1906"},
241 "hough": {'a':6378270.0,'rf':297.,'description':"Hough"},
242 "intl": {'a':6378388.0,'rf':297.,'description':"International 1909 (Hayford)"},
243 "krass": {'a':6378245.0,'rf':298.3,'description':"Krassovsky, 1942"},
244 "kaula": {'a':6378163.,'rf':298.24,'description':"Kaula 1961"},
245 "lerch": {'a':6378139.,'rf':298.257,'description':"Lerch 1979"},
246 "mprts": {'a':6397300.,'rf':191.,'description':"Maupertius 1738"},
247 "new_intl": {'a':6378157.5,'b':6356772.2,'description':"New International 1967"},
248 "plessis": {'a':6376523.,'b':6355863.,'description':"Plessis 1817 (France)"},
249 "SEasia": {'a':6378155.0,'b':6356773.3205,'description':"Southeast Asia"},
250 "walbeck": {'a':6376896.0,'b':6355834.8467,'description':"Walbeck"},
251 "WGS60": {'a':6378165.0,'rf':298.3,'description':"WGS 60"},
252 "WGS66": {'a':6378145.0,'rf':298.25,'description':"WGS 66"},
253 "WGS72": {'a':6378135.0,'rf':298.26,'description':"WGS 72"},
254 "WGS84": {'a':6378137.0,'rf':298.257223563,'description':"WGS 84"},
255 "sphere": {'a':6370997.0,'b':6370997.0,'description':"Normal Sphere"},
256 }
257
258
259
260
261
262 set_datapath(pyproj_datadir)
263
264 -class Proj(_proj.Proj):
265 """
266 performs cartographic transformations (converts from
267 longitude,latitude to native map projection x,y coordinates and
268 vice versa) using proj (https://github.com/OSGeo/proj.4/wiki).
269
270 A Proj class instance is initialized with proj map projection
271 control parameter key/value pairs. The key/value pairs can
272 either be passed in a dictionary, or as keyword arguments,
273 or as a proj4 string (compatible with the proj command). See
274 http://www.remotesensing.org/geotiff/proj_list for examples of
275 key/value pairs defining different map projections.
276
277 Calling a Proj class instance with the arguments lon, lat will
278 convert lon/lat (in degrees) to x/y native map projection
279 coordinates (in meters). If optional keyword 'inverse' is True
280 (default is False), the inverse transformation from x/y to
281 lon/lat is performed. If optional keyword 'radians' is True
282 (default is False) lon/lat are interpreted as radians instead of
283 degrees. If optional keyword 'errcheck' is True (default is
284 False) an exception is raised if the transformation is invalid.
285 If errcheck=False and the transformation is invalid, no
286 exception is raised and 1.e30 is returned. If the optional keyword
287 'preserve_units' is True, the units in map projection coordinates
288 are not forced to be meters.
289
290 Works with numpy and regular python array objects, python
291 sequences and scalars.
292 """
293
294 - def __new__(self, projparams=None, preserve_units=False, **kwargs):
295 """
296 initialize a Proj class instance.
297
298 Proj4 projection control parameters must either be given in a
299 dictionary 'projparams' or as keyword arguments. See the proj
300 documentation (https://github.com/OSGeo/proj.4/wiki) for more information
301 about specifying projection parameters.
302
303 Example usage:
304
305 >>> from pyproj import Proj
306 >>> p = Proj(proj='utm',zone=10,ellps='WGS84') # use kwargs
307 >>> x,y = p(-120.108, 34.36116666)
308 >>> 'x=%9.3f y=%11.3f' % (x,y)
309 'x=765975.641 y=3805993.134'
310 >>> 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True)
311 'lon=-120.108 lat=34.361'
312 >>> # do 3 cities at a time in a tuple (Fresno, LA, SF)
313 >>> lons = (-119.72,-118.40,-122.38)
314 >>> lats = (36.77, 33.93, 37.62 )
315 >>> x,y = p(lons, lats)
316 >>> 'x: %9.3f %9.3f %9.3f' % x
317 'x: 792763.863 925321.537 554714.301'
318 >>> 'y: %9.3f %9.3f %9.3f' % y
319 'y: 4074377.617 3763936.941 4163835.303'
320 >>> lons, lats = p(x, y, inverse=True) # inverse transform
321 >>> 'lons: %8.3f %8.3f %8.3f' % lons
322 'lons: -119.720 -118.400 -122.380'
323 >>> 'lats: %8.3f %8.3f %8.3f' % lats
324 'lats: 36.770 33.930 37.620'
325 >>> p2 = Proj('+proj=utm +zone=10 +ellps=WGS84') # use proj4 string
326 >>> x,y = p2(-120.108, 34.36116666)
327 >>> 'x=%9.3f y=%11.3f' % (x,y)
328 'x=765975.641 y=3805993.134'
329 >>> p = Proj(init="epsg:32667")
330 >>> 'x=%12.3f y=%12.3f (meters)' % p(-114.057222, 51.045)
331 'x=-1783486.760 y= 6193833.196 (meters)'
332 >>> p = Proj("+init=epsg:32667",preserve_units=True)
333 >>> 'x=%12.3f y=%12.3f (feet)' % p(-114.057222, 51.045)
334 'x=-5851322.810 y=20320934.409 (feet)'
335 """
336
337 if projparams is None:
338 if len(kwargs) == 0:
339 raise RuntimeError('no projection control parameters specified')
340 else:
341 projstring = _dict2string(kwargs)
342 elif isinstance(projparams, string_types):
343
344 projstring = projparams
345 else:
346 projstring = _dict2string(projparams)
347
348 if not projstring.count('+units=') and not preserve_units:
349 projstring = '+units=m '+projstring
350 else:
351 kvpairs = []
352 for kvpair in projstring.split():
353 if kvpair.startswith('+units') and not preserve_units:
354 k,v = kvpair.split('=')
355 kvpairs.append(k+'=m ')
356 else:
357 kvpairs.append(kvpair+' ')
358 projstring = ''.join(kvpairs)
359
360
361 projstring = projstring.replace('EPSG','epsg')
362 return _proj.Proj.__new__(self, projstring)
363
365
366 """
367 Calling a Proj class instance with the arguments lon, lat will
368 convert lon/lat (in degrees) to x/y native map projection
369 coordinates (in meters). If optional keyword 'inverse' is True
370 (default is False), the inverse transformation from x/y to
371 lon/lat is performed. If optional keyword 'radians' is True
372 (default is False) the units of lon/lat are radians instead of
373 degrees. If optional keyword 'errcheck' is True (default is
374 False) an exception is raised if the transformation is invalid.
375 If errcheck=False and the transformation is invalid, no
376 exception is raised and 1.e30 is returned.
377
378 Inputs should be doubles (they will be cast to doubles if they
379 are not, causing a slight performance hit).
380
381 Works with numpy and regular python array objects, python
382 sequences and scalars, but is fastest for array objects.
383 """
384 inverse = kw.get('inverse', False)
385 radians = kw.get('radians', False)
386 errcheck = kw.get('errcheck', False)
387
388
389
390
391
392
393
394
395 lon, lat = args
396
397 inx, xisfloat, xislist, xistuple = _copytobuffer(lon)
398 iny, yisfloat, yislist, yistuple = _copytobuffer(lat)
399
400 if inverse:
401 _proj.Proj._inv(self, inx, iny, radians=radians, errcheck=errcheck)
402 else:
403 _proj.Proj._fwd(self, inx, iny, radians=radians, errcheck=errcheck)
404
405 outx = _convertback(xisfloat,xislist,xistuple,inx)
406 outy = _convertback(yisfloat,yislist,xistuple,iny)
407 return outx, outy
408
410 """returns an equivalent Proj in the corresponding lon/lat
411 coordinates. (see pj_latlong_from_proj() in the Proj.4 C API)"""
412 string_def = _proj.Proj.to_latlong_def(self)
413
414 if not isinstance(string_def, str):
415 string_def = string_def.decode()
416 return Proj(string_def)
417
419 """returns True if projection in geographic (lon/lat) coordinates"""
420 return _proj.Proj.is_latlong(self)
421
423 """returns True if projection in geocentric (x/y) coordinates"""
424 return _proj.Proj.is_geocent(self)
425
427 """Returns formal definition string for projection
428
429 >>> Proj('+init=epsg:4326').definition_string()
430 ' +units=m +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'
431 >>>
432 """
433
434 string_def = _proj.Proj.definition_string(self)
435
436 if not isinstance(string_def, str):
437 string_def = string_def.decode()
438
439 return string_def
440
536
608
610 try:
611
612 return array('d',(float(x),)),True,False,False
613 except:
614 raise TypeError('input must be an array, list, tuple or scalar')
615
617 """
618 return a copy of x as an object that supports the python Buffer
619 API (python array if input is float, list or tuple, numpy array
620 if input is a numpy array). returns copyofx, isfloat, islist,
621 istuple (islist is True if input is a list, istuple is true if
622 input is a tuple, isfloat is true if input is a float).
623 """
624
625 isfloat = False; islist = False; istuple = False
626
627
628 if hasattr(x,'shape'):
629 if x.shape == ():
630 return _copytobuffer_return_scalar(x)
631 else:
632 try:
633
634
635
636 x.dtype.char
637
638
639
640 inx = x.copy(order="C").astype('d')
641
642 return inx,False,False,False
643 except:
644 try:
645
646
647 x.typecode()
648 inx = x.astype('d')
649
650 return inx,False,False,False
651 except:
652 raise TypeError('input must be an array, list, tuple or scalar')
653 else:
654
655 if hasattr(x, 'typecode'):
656
657 inx = array('d',x)
658
659
660 elif type(x) == list:
661 inx = array('d',x)
662 islist = True
663
664 elif type(x) == tuple:
665 inx = array('d',x)
666 istuple = True
667
668 else:
669 return _copytobuffer_return_scalar(x)
670 return inx,isfloat,islist,istuple
671
673
674 if isfloat:
675 return inx[0]
676 elif islist:
677 return inx.tolist()
678 elif istuple:
679 return tuple(inx)
680 else:
681 return inx
682
684
685 pjargs = []
686 for key,value in projparams.items():
687 pjargs.append('+'+key+"="+str(value)+' ')
688 return ''.join(pjargs)
689
690 -class Geod(_proj.Geod):
691 """
692 performs forward and inverse geodetic, or Great Circle,
693 computations. The forward computation (using the 'fwd' method)
694 involves determining latitude, longitude and back azimuth of a
695 computations. The forward computation (using the 'fwd' method)
696 involves determining latitude, longitude and back azimuth of a
697 terminus point given the latitude and longitude of an initial
698 point, plus azimuth and distance. The inverse computation (using
699 the 'inv' method) involves determining the forward and back
700 azimuths and distance given the latitudes and longitudes of an
701 initial and terminus point.
702 """
703 - def __new__(cls, initstring=None, **kwargs):
704 """
705 initialize a Geod class instance.
706
707 Geodetic parameters for specifying the ellipsoid
708 can be given in a dictionary 'initparams', as keyword arguments,
709 or as as proj4 geod initialization string.
710 Following is a list of the ellipsoids that may be defined using the
711 'ellps' keyword (these are stored in the model variable pj_ellps)::
712
713 MERIT a=6378137.0 rf=298.257 MERIT 1983
714 SGS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85
715 GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980)
716 IAU76 a=6378140.0 rf=298.257 IAU 1976
717 airy a=6377563.396 b=6356256.910 Airy 1830
718 APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
719 airy a=6377563.396 b=6356256.910 Airy 1830
720 APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
721 NWL9D a=6378145.0. rf=298.25 Naval Weapons Lab., 1965
722 mod_airy a=6377340.189 b=6356034.446 Modified Airy
723 andrae a=6377104.43 rf=300.0 Andrae 1876 (Den., Iclnd.)
724 aust_SA a=6378160.0 rf=298.25 Australian Natl & S. Amer. 1969
725 GRS67 a=6378160.0 rf=298.247167427 GRS 67(IUGG 1967)
726 bessel a=6377397.155 rf=299.1528128 Bessel 1841
727 bess_nam a=6377483.865 rf=299.1528128 Bessel 1841 (Namibia)
728 clrk66 a=6378206.4 b=6356583.8 Clarke 1866
729 clrk80 a=6378249.145 rf=293.4663 Clarke 1880 mod.
730 CPM a=6375738.7 rf=334.29 Comm. des Poids et Mesures 1799
731 delmbr a=6376428. rf=311.5 Delambre 1810 (Belgium)
732 engelis a=6378136.05 rf=298.2566 Engelis 1985
733 evrst30 a=6377276.345 rf=300.8017 Everest 1830
734 evrst48 a=6377304.063 rf=300.8017 Everest 1948
735 evrst56 a=6377301.243 rf=300.8017 Everest 1956
736 evrst69 a=6377295.664 rf=300.8017 Everest 1969
737 evrstSS a=6377298.556 rf=300.8017 Everest (Sabah & Sarawak)
738 fschr60 a=6378166. rf=298.3 Fischer (Mercury Datum) 1960
739 fschr60m a=6378155. rf=298.3 Modified Fischer 1960
740 fschr68 a=6378150. rf=298.3 Fischer 1968
741 helmert a=6378200. rf=298.3 Helmert 1906
742 hough a=6378270.0 rf=297. Hough
743 helmert a=6378200. rf=298.3 Helmert 1906
744 hough a=6378270.0 rf=297. Hough
745 intl a=6378388.0 rf=297. International 1909 (Hayford)
746 krass a=6378245.0 rf=298.3 Krassovsky, 1942
747 kaula a=6378163. rf=298.24 Kaula 1961
748 lerch a=6378139. rf=298.257 Lerch 1979
749 mprts a=6397300. rf=191. Maupertius 1738
750 new_intl a=6378157.5 b=6356772.2 New International 1967
751 plessis a=6376523. b=6355863. Plessis 1817 (France)
752 SEasia a=6378155.0 b=6356773.3205 Southeast Asia
753 walbeck a=6376896.0 b=6355834.8467 Walbeck
754 WGS60 a=6378165.0 rf=298.3 WGS 60
755 WGS66 a=6378145.0 rf=298.25 WGS 66
756 WGS72 a=6378135.0 rf=298.26 WGS 72
757 WGS84 a=6378137.0 rf=298.257223563 WGS 84
758 sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)
759
760 The parameters of the ellipsoid may also be set directly using
761 the 'a' (semi-major or equatorial axis radius) keyword, and
762 any one of the following keywords: 'b' (semi-minor,
763 or polar axis radius), 'e' (eccentricity), 'es' (eccentricity
764 squared), 'f' (flattening), or 'rf' (reciprocal flattening).
765
766 See the proj documentation (https://github.com/OSGeo/proj.4/wiki) for more
767 information about specifying ellipsoid parameters (specifically,
768 the chapter 'Specifying the Earth's figure' in the main Proj
769 users manual).
770
771 Example usage:
772
773 >>> from pyproj import Geod
774 >>> g = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid.
775 >>> # specify the lat/lons of some cities.
776 >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
777 >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
778 >>> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
779 >>> london_lat = 51.+(32./60.); london_lon = -(5./60.)
780 >>> # compute forward and back azimuths, plus distance
781 >>> # between Boston and Portland.
782 >>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
783 >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist)
784 '-66.531 75.654 4164192.708'
785 >>> # compute latitude, longitude and back azimuth of Portland,
786 >>> # given Boston lat/lon, forward azimuth and distance to Portland.
787 >>> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
788 >>> "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
789 '45.517 -123.683 75.654'
790 >>> # compute the azimuths, distances from New York to several
791 >>> # cities (pass a list)
792 >>> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
793 >>> lons2 = [boston_lon, portland_lon, london_lon]
794 >>> lats2 = [boston_lat, portland_lat, london_lat]
795 >>> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
796 >>> for faz,baz,d in list(zip(az12,az21,dist)): "%7.3f %7.3f %9.3f" % (faz,baz,d)
797 ' 54.663 -123.448 288303.720'
798 '-65.463 79.342 4013037.318'
799 ' 51.254 -71.576 5579916.651'
800 >>> g2 = Geod('+ellps=clrk66') # use proj4 style initialization string
801 >>> az12,az21,dist = g2.inv(boston_lon,boston_lat,portland_lon,portland_lat)
802 >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist)
803 '-66.531 75.654 4164192.708'
804 """
805
806
807 ellpsd = {}
808 if initstring is not None:
809 for kvpair in initstring.split():
810
811
812 if kvpair.find('=') == -1:
813 continue
814 k,v = kvpair.split('=')
815 k = k.lstrip('+')
816 if k in ['a','b','rf','f','es','e']:
817 v = float(v)
818 ellpsd[k] = v
819
820 kwargs = dict(list(kwargs.items()) + list(ellpsd.items()))
821 sphere = False
822 if 'ellps' in kwargs:
823
824 ellps_dict = pj_ellps[kwargs['ellps']]
825 a = ellps_dict['a']
826 if ellps_dict['description']=='Normal Sphere':
827 sphere = True
828 if 'b' in ellps_dict:
829 b = ellps_dict['b']
830 es = 1. - (b * b) / (a * a)
831 f = (a - b)/a
832 elif 'rf' in ellps_dict:
833 f = 1./ellps_dict['rf']
834 b = a*(1. - f)
835 es = 1. - (b * b) / (a * a)
836 else:
837
838
839
840
841
842
843 a = kwargs['a']
844 if 'b' in kwargs:
845 b = kwargs['b']
846 es = 1. - (b * b) / (a * a)
847 f = (a - b)/a
848 elif 'rf' in kwargs:
849 f = 1./kwargs['rf']
850 b = a*(1. - f)
851 es = 1. - (b * b) / (a * a)
852 elif 'f' in kwargs:
853 f = kwargs['f']
854 b = a*(1. - f)
855 es = 1. - (b/a)**2
856 elif 'es' in kwargs:
857 es = kwargs['es']
858 b = math.sqrt(a**2 - es*a**2)
859 f = (a - b)/a
860 elif 'e' in kwargs:
861 es = kwargs['e']**2
862 b = math.sqrt(a**2 - es*a**2)
863 f = (a - b)/a
864 else:
865 b = a
866 f = 0.
867 es = 0.
868
869
870 if math.fabs(f) < 1.e-8: sphere = True
871
872 return _proj.Geod.__new__(cls, a, f, sphere, b, es)
873
874 - def fwd(self, lons, lats, az, dist, radians=False):
875 """
876 forward transformation - Returns longitudes, latitudes and back
877 azimuths of terminus points given longitudes (lons) and
878 latitudes (lats) of initial points, plus forward azimuths (az)
879 and distances (dist).
880 latitudes (lats) of initial points, plus forward azimuths (az)
881 and distances (dist).
882
883 Works with numpy and regular python array objects, python
884 sequences and scalars.
885
886 if radians=True, lons/lats and azimuths are radians instead of
887 degrees. Distances are in meters.
888 """
889
890 inx, xisfloat, xislist, xistuple = _copytobuffer(lons)
891 iny, yisfloat, yislist, yistuple = _copytobuffer(lats)
892 inz, zisfloat, zislist, zistuple = _copytobuffer(az)
893 ind, disfloat, dislist, distuple = _copytobuffer(dist)
894 _proj.Geod._fwd(self, inx, iny, inz, ind, radians=radians)
895
896 outx = _convertback(xisfloat,xislist,xistuple,inx)
897 outy = _convertback(yisfloat,yislist,xistuple,iny)
898 outz = _convertback(zisfloat,zislist,zistuple,inz)
899 return outx, outy, outz
900
901 - def inv(self,lons1,lats1,lons2,lats2,radians=False):
902 """
903 inverse transformation - Returns forward and back azimuths, plus
904 distances between initial points (specified by lons1, lats1) and
905 terminus points (specified by lons2, lats2).
906
907 Works with numpy and regular python array objects, python
908 sequences and scalars.
909
910 if radians=True, lons/lats and azimuths are radians instead of
911 degrees. Distances are in meters.
912 """
913
914 inx, xisfloat, xislist, xistuple = _copytobuffer(lons1)
915 iny, yisfloat, yislist, yistuple = _copytobuffer(lats1)
916 inz, zisfloat, zislist, zistuple = _copytobuffer(lons2)
917 ind, disfloat, dislist, distuple = _copytobuffer(lats2)
918 _proj.Geod._inv(self,inx,iny,inz,ind,radians=radians)
919
920 outx = _convertback(xisfloat,xislist,xistuple,inx)
921 outy = _convertback(yisfloat,yislist,xistuple,iny)
922 outz = _convertback(zisfloat,zislist,zistuple,inz)
923 return outx, outy, outz
924
925 - def npts(self, lon1, lat1, lon2, lat2, npts, radians=False):
926 """
927 Given a single initial point and terminus point (specified by
928 python floats lon1,lat1 and lon2,lat2), returns a list of
929 longitude/latitude pairs describing npts equally spaced
930 intermediate points along the geodesic between the initial and
931 terminus points.
932
933 if radians=True, lons/lats are radians instead of degrees.
934
935 Example usage:
936
937 >>> from pyproj import Geod
938 >>> g = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid.
939 >>> # specify the lat/lons of Boston and Portland.
940 >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
941 >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
942 >>> # find ten equally spaced points between Boston and Portland.
943 >>> lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
944 >>> for lon,lat in lonlats: '%6.3f %7.3f' % (lat, lon)
945 '43.528 -75.414'
946 '44.637 -79.883'
947 '45.565 -84.512'
948 '46.299 -89.279'
949 '46.830 -94.156'
950 '47.149 -99.112'
951 '47.251 -104.106'
952 '47.136 -109.100'
953 '46.805 -114.051'
954 '46.262 -118.924'
955 >>> # test with radians=True (inputs/outputs in radians, not degrees)
956 >>> import math
957 >>> dg2rad = math.radians(1.)
958 >>> rad2dg = math.degrees(1.)
959 >>> lonlats = g.npts(dg2rad*boston_lon,dg2rad*boston_lat,dg2rad*portland_lon,dg2rad*portland_lat,10,radians=True)
960 >>> for lon,lat in lonlats: '%6.3f %7.3f' % (rad2dg*lat, rad2dg*lon)
961 '43.528 -75.414'
962 '44.637 -79.883'
963 '45.565 -84.512'
964 '46.299 -89.279'
965 '46.830 -94.156'
966 '47.149 -99.112'
967 '47.251 -104.106'
968 '47.136 -109.100'
969 '46.805 -114.051'
970 '46.262 -118.924'
971 """
972 lons, lats = _proj.Geod._npts(self, lon1, lat1, lon2, lat2, npts, radians=radians)
973 return list(zip(lons, lats))
974
976
977 for (ellps, vals) in pj_ellps.items():
978 if self.a == vals['a']:
979 b = vals.get('b', None)
980 rf = vals.get('rf', None)
981
982
983 if self.b == b or (not self.sphere and (1.0/self.f) == rf):
984 return "{modname}.{classname}(ellps={ellps!r})" \
985 "".format(modname=self.__module__,
986 classname=self.__class__.__name__,
987 ellps=ellps)
988
989
990 return _proj.Geod.__repr__(self)
991
993 """
994 equality operator == for Geod objects
995
996 Example usage:
997
998 >>> from pyproj import Geod
999 >>> gclrk1 = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid.
1000 >>> gclrk2 = Geod(a=6378206.4, b=6356583.8) # Define Clarke 1866 using parameters
1001 >>> gclrk1 == gclrk2
1002 True
1003 >>> gwgs66 = Geod('+ellps=WGS66') # WGS 66 ellipsoid, Proj.4 style
1004 >>> gnwl9d = Geod('+ellps=NWL9D') # Naval Weapons Lab., 1965 ellipsoid
1005 >>> # these ellipsoids are the same
1006 >>> gnwl9d == gwgs66
1007 True
1008 >>> gclrk1 != gnwl9d # Clarke 1866 is unlike NWL9D
1009 True
1010 """
1011 if not isinstance(other, _proj.Geod):
1012 return False
1013
1014 return self.__repr__() == other.__repr__()
1015
1016 -def test(**kwargs):
1017 """run the examples in the docstrings using the doctest module"""
1018 import doctest
1019 import pyproj
1020
1021 verbose = kwargs.get('verbose')
1022 failure_count, test_count = doctest.testmod(pyproj, verbose=verbose)
1023
1024 return failure_count
1025
1026 if __name__ == "__main__":
1027 sys.exit(test(verbose=True))
1028