Package pyproj

Source Code for Package pyproj

   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 # python 2.6+ 
  62  except ImportError: 
  63      pass # python 3.x 
  64  #import numpy as np 
  65   
  66  # Python 2/3 compatibility 
  67  if sys.version_info[0] == 2:            # Python 2 
  68     string_types = (basestring,) 
  69  else:                                   # Python 3 
  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  #if not os.path.isdir(pyproj_datadir): 
 259  #    msg="proj data directory not found. Expecting it at: %s"%pyproj_datadir 
 260  #    raise IOError(msg) 
 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 # if projparams is None, use kwargs. 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 # if projparams is a string or a unicode string, interpret as a proj4 init string. 344 projstring = projparams 345 else: # projparams a dict 346 projstring = _dict2string(projparams) 347 # make sure units are meters if preserve_units is False. 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 # look for EPSG, replace with epsg (EPSG only works 360 # on case-insensitive filesystems). 361 projstring = projstring.replace('EPSG','epsg') 362 return _proj.Proj.__new__(self, projstring)
363
364 - def __call__(self, *args, **kw):
365 #,lon,lat,inverse=False,radians=False,errcheck=False): 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 #if len(args) == 1: 388 # latlon = np.array(args[0], copy=True, 389 # order='C', dtype=float, ndmin=2) 390 # if inverse: 391 # _proj.Proj._invn(self, latlon, radians=radians, errcheck=errcheck) 392 # else: 393 # _proj.Proj._fwdn(self, latlon, radians=radians, errcheck=errcheck) 394 # return latlon 395 lon, lat = args 396 # process inputs, making copies that support buffer API. 397 inx, xisfloat, xislist, xistuple = _copytobuffer(lon) 398 iny, yisfloat, yislist, yistuple = _copytobuffer(lat) 399 # call proj4 functions. inx and iny modified in place. 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 # if inputs were lists, tuples or floats, convert back. 405 outx = _convertback(xisfloat,xislist,xistuple,inx) 406 outy = _convertback(yisfloat,yislist,xistuple,iny) 407 return outx, outy
408
409 - def to_latlong(self):
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 # Decode bytes --> Unicode String (for Python3) 414 if not isinstance(string_def, str): 415 string_def = string_def.decode() 416 return Proj(string_def)
417
418 - def is_latlong(self):
419 """returns True if projection in geographic (lon/lat) coordinates""" 420 return _proj.Proj.is_latlong(self)
421
422 - def is_geocent(self):
423 """returns True if projection in geocentric (x/y) coordinates""" 424 return _proj.Proj.is_geocent(self)
425
426 - def definition_string(self):
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
441 -def transform(p1, p2, x, y, z=None, radians=False):
442 """ 443 x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False) 444 445 Transform points between two coordinate systems defined by the 446 Proj instances p1 and p2. 447 448 The points x1,y1,z1 in the coordinate system defined by p1 are 449 transformed to x2,y2,z2 in the coordinate system defined by p2. 450 451 z1 is optional, if it is not set it is assumed to be zero (and 452 only x2 and y2 are returned). 453 454 In addition to converting between cartographic and geographic 455 projection coordinates, this function can take care of datum 456 shifts (which cannot be done using the __call__ method of the 457 Proj instances). It also allows for one of the coordinate 458 systems to be geographic (proj = 'latlong'). 459 460 If optional keyword 'radians' is True (default is False) and p1 461 is defined in geographic coordinate (pj.is_latlong() is True), 462 x1,y1 is interpreted as radians instead of the default degrees. 463 Similarly, if p2 is defined in geographic coordinates and 464 radians=True, x2, y2 are returned in radians instead of degrees. 465 if p1.is_latlong() and p2.is_latlong() both are False, the 466 radians keyword has no effect. 467 468 x,y and z can be numpy or regular python arrays, python 469 lists/tuples or scalars. Arrays are fastest. For projections in 470 geocentric coordinates, values of x and y are given in meters. 471 z is always meters. 472 473 Example usage: 474 475 >>> # projection 1: UTM zone 15, grs80 ellipse, NAD83 datum 476 >>> # (defined by epsg code 26915) 477 >>> p1 = Proj(init='epsg:26915') 478 >>> # projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum 479 >>> p2 = Proj(init='epsg:26715') 480 >>> # find x,y of Jefferson City, MO. 481 >>> x1, y1 = p1(-92.199881,38.56694) 482 >>> # transform this point to projection 2 coordinates. 483 >>> x2, y2 = transform(p1,p2,x1,y1) 484 >>> '%9.3f %11.3f' % (x1,y1) 485 '569704.566 4269024.671' 486 >>> '%9.3f %11.3f' % (x2,y2) 487 '569722.342 4268814.027' 488 >>> '%8.3f %5.3f' % p2(x2,y2,inverse=True) 489 ' -92.200 38.567' 490 >>> # process 3 points at a time in a tuple 491 >>> lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri 492 >>> lons = (-92.22,-94.72,-90.37) 493 >>> x1, y1 = p1(lons,lats) 494 >>> x2, y2 = transform(p1,p2,x1,y1) 495 >>> xy = x1+y1 496 >>> '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy 497 '567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005' 498 >>> xy = x2+y2 499 >>> '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy 500 '567721.149 351747.558 728569.133 4297989.112 4353489.644 4292106.305' 501 >>> lons, lats = p2(x2,y2,inverse=True) 502 >>> xy = lons+lats 503 >>> '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy 504 ' -92.220 -94.720 -90.370 38.830 39.320 38.750' 505 >>> # test datum shifting, installation of extra datum grid files. 506 >>> p1 = Proj(proj='latlong',datum='WGS84') 507 >>> x1 = -111.5; y1 = 45.25919444444 508 >>> p2 = Proj(proj="utm",zone=10,datum='NAD27') 509 >>> x2, y2 = transform(p1, p2, x1, y1) 510 >>> "%s %s" % (str(x2)[:9],str(y2)[:9]) 511 '1402285.9 5076292.4' 512 """ 513 # check that p1 and p2 are from the Proj class 514 if not isinstance(p1, Proj): 515 raise TypeError("p1 must be a Proj class") 516 if not isinstance(p2, Proj): 517 raise TypeError("p2 must be a Proj class") 518 519 # process inputs, making copies that support buffer API. 520 inx, xisfloat, xislist, xistuple = _copytobuffer(x) 521 iny, yisfloat, yislist, yistuple = _copytobuffer(y) 522 if z is not None: 523 inz, zisfloat, zislist, zistuple = _copytobuffer(z) 524 else: 525 inz = None 526 # call pj_transform. inx,iny,inz buffers modified in place. 527 _proj._transform(p1,p2,inx,iny,inz,radians) 528 # if inputs were lists, tuples or floats, convert back. 529 outx = _convertback(xisfloat,xislist,xistuple,inx) 530 outy = _convertback(yisfloat,yislist,xistuple,iny) 531 if inz is not None: 532 outz = _convertback(zisfloat,zislist,zistuple,inz) 533 return outx, outy, outz 534 else: 535 return outx, outy
536
537 -def itransform(p1, p2, points, radians=False, switch=False):
538 """ 539 points2 = transform(p1, p2, points1, radians=False) 540 Iterator/generator version of the function pyproj.transform. 541 Transform points between two coordinate systems defined by the 542 Proj instances p1 and p2. This function can be used as an alternative 543 to pyproj.transform when there is a need to transform a big number of 544 coordinates lazily, for example when reading and processing from a file. 545 Points1 is an iterable/generator of coordinates x1,y1(,z1) or lon1,lat1(,z1) 546 in the coordinate system defined by p1. Points2 is an iterator that returns tuples 547 of x2,y2(,z2) or lon2,lat2(,z2) coordinates in the coordinate system defined by p2. 548 z are provided optionally. 549 550 Points1 can be: 551 - a tuple/list of tuples/lists i.e. for 2d points: [(xi,yi),(xj,yj),....(xn,yn)] 552 - a Nx3 or Nx2 2d numpy array where N is the point number 553 - a generator of coordinates (xi,yi) for 2d points or (xi,yi,zi) for 3d 554 555 If optional keyword 'switch' is True (default is False) then x, y or lon,lat coordinates 556 of points are switched to y, x or lat, lon. 557 If optional keyword 'radians' is True (default is False) and p1 558 is defined in geographic coordinate (pj.is_latlong() is True), 559 x1,y1 is interpreted as radians instead of the default degrees. 560 Similarly, if p2 is defined in geographic coordinates and 561 radians=True, x2, y2 are returned in radians instead of degrees. 562 if p1.is_latlong() and p2.is_latlong() both are False, the 563 radians keyword has no effect. 564 Example usage: 565 >>> # projection 1: WGS84 566 >>> # (defined by epsg code 4326) 567 >>> p1 = Proj(init='epsg:4326') 568 >>> # projection 2: GGRS87 / Greek Grid 569 >>> p2 = Proj(init='epsg:2100') 570 >>> # Three points with coordinates lon, lat in p1 571 >>> points = [(22.95, 40.63), (22.81, 40.53), (23.51, 40.86)] 572 >>> # transform this point to projection 2 coordinates. 573 >>> for pt in itransform(p1,p2,points): '%6.3f %7.3f' % pt 574 '411050.470 4497928.574' 575 '399060.236 4486978.710' 576 '458553.243 4523045.485' 577 """ 578 if not isinstance(p1, Proj): 579 raise TypeError("p1 must be a Proj class") 580 if not isinstance(p2, Proj): 581 raise TypeError("p2 must be a Proj class") 582 583 it = iter(points) # point iterator 584 # get first point to check stride 585 try: 586 fst_pt = next(it) 587 except StopIteration: 588 raise ValueError("iterable must contain at least one point") 589 590 stride = len(fst_pt) 591 if stride not in (2,3): 592 raise ValueError("points can contain up to 3 coordinates") 593 594 # create a coordinate sequence generator etc. x1,y1,z1,x2,y2,z2,.... 595 # chain so the generator returns the first point that was already acquired 596 coord_gen = chain(fst_pt, (coords[c] for coords in it for c in range(stride))) 597 598 while True: 599 # create a temporary buffer storage for the next 64 points (64*stride*8 bytes) 600 buff = array('d', islice(coord_gen, 0, 64*stride)) 601 if len(buff) == 0: 602 break; 603 604 _proj._transform_sequence(p1, p2, stride, buff, radians, switch) 605 606 for pt in zip(*([iter(buff)]*stride)): 607 yield pt
608
609 -def _copytobuffer_return_scalar(x):
610 try: 611 # inx,isfloat,islist,istuple 612 return array('d',(float(x),)),True,False,False 613 except: 614 raise TypeError('input must be an array, list, tuple or scalar')
615
616 -def _copytobuffer(x):
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 # make sure x supports Buffer API and contains doubles. 625 isfloat = False; islist = False; istuple = False 626 # first, if it's a numpy array scalar convert to float 627 # (array scalars don't support buffer API) 628 if hasattr(x,'shape'): 629 if x.shape == (): 630 return _copytobuffer_return_scalar(x) 631 else: 632 try: 633 # typecast numpy arrays to double. 634 # (this makes a copy - which is crucial 635 # since buffer is modified in place) 636 x.dtype.char 637 # Basemap issue 638 # https://github.com/matplotlib/basemap/pull/223/files 639 # (deal with input array in fortran order) 640 inx = x.copy(order="C").astype('d') 641 # inx,isfloat,islist,istuple 642 return inx,False,False,False 643 except: 644 try: # perhaps they are Numeric/numarrays? 645 # sorry, not tested yet. 646 # i don't know Numeric/numarrays has `shape'. 647 x.typecode() 648 inx = x.astype('d') 649 # inx,isfloat,islist,istuple 650 return inx,False,False,False 651 except: 652 raise TypeError('input must be an array, list, tuple or scalar') 653 else: 654 # perhaps they are regular python arrays? 655 if hasattr(x, 'typecode'): 656 #x.typecode 657 inx = array('d',x) 658 # try to convert to python array 659 # a list. 660 elif type(x) == list: 661 inx = array('d',x) 662 islist = True 663 # a tuple. 664 elif type(x) == tuple: 665 inx = array('d',x) 666 istuple = True 667 # a scalar? 668 else: 669 return _copytobuffer_return_scalar(x) 670 return inx,isfloat,islist,istuple
671
672 -def _convertback(isfloat,islist,istuple,inx):
673 # if inputs were lists, tuples or floats, convert back to original type. 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
683 -def _dict2string(projparams):
684 # convert a dict to a proj4 string. 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 # if initparams is a proj-type init string, 806 # convert to dict. 807 ellpsd = {} 808 if initstring is not None: 809 for kvpair in initstring.split(): 810 # Actually only +a and +b are needed 811 # We can ignore safely any parameter that doesn't have a value 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 # merge this dict with kwargs dict. 820 kwargs = dict(list(kwargs.items()) + list(ellpsd.items())) 821 sphere = False 822 if 'ellps' in kwargs: 823 # ellipse name given, look up in pj_ellps dict 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 # a (semi-major axis) and one of 838 # b the semi-minor axis 839 # rf the reciprocal flattening 840 # f flattening 841 # es eccentricity squared 842 # must be given. 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 #msg='ellipse name or a, plus one of f,es,b must be given' 869 #raise ValueError(msg) 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 # process inputs, making copies that support buffer API. 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 # if inputs were lists, tuples or floats, convert back. 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 # process inputs, making copies that support buffer API. 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 # if inputs were lists, tuples or floats, convert back. 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
975 - def __repr__(self):
976 # search for ellipse name 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 # self.sphere is True when self.f is zero or very close to 982 # zero (0), so prevent divide by zero. 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 # no ellipse name found, call super class 990 return _proj.Geod.__repr__(self)
991
992 - def __eq__(self, other):
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