Package pyproj

Source Code for Package pyproj

  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  #import numpy as np 
 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 # if projparams is None, use kwargs. 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 # if projparams is a string, interpret as a proj4 init string. 324 projstring = projparams 325 else: # projparams a dict 326 projstring = _dict2string(projparams) 327 # make sure units are meters if preserve_units is False. 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 # look for EPSG, replace with epsg (EPSG only works 340 # on case-insensitive filesystems). 341 projstring = projstring.replace('EPSG','epsg') 342 return _proj.Proj.__new__(self, projstring)
343
344 - def __call__(self, *args, **kw):
345 #,lon,lat,inverse=False,radians=False,errcheck=False): 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 #if len(args) == 1: 372 # latlon = np.array(args[0], copy=True, 373 # order='C', dtype=float, ndmin=2) 374 # if inverse: 375 # _proj.Proj._invn(self, latlon, radians=radians, errcheck=errcheck) 376 # else: 377 # _proj.Proj._fwdn(self, latlon, radians=radians, errcheck=errcheck) 378 # return latlon 379 lon, lat = args 380 # process inputs, making copies that support buffer API. 381 inx, xisfloat, xislist, xistuple = _copytobuffer(lon) 382 iny, yisfloat, yislist, yistuple = _copytobuffer(lat) 383 # call proj4 functions. inx and iny modified in place. 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 # if inputs were lists, tuples or floats, convert back. 389 outx = _convertback(xisfloat,xislist,xistuple,inx) 390 outy = _convertback(yisfloat,yislist,xistuple,iny) 391 return outx, outy
392
393 - def is_latlong(self):
394 """returns True if projection in geographic (lon/lat) coordinates""" 395 return _proj.Proj.is_latlong(self)
396
397 - def is_geocent(self):
398 """returns True if projection in geocentric (x/y) coordinates""" 399 return _proj.Proj.is_geocent(self)
400
401 -def transform(p1, p2, x, y, z=None, radians=False):
402 """ 403 x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False) 404 405 Transform points between two coordinate systems defined by the 406 Proj instances p1 and p2. 407 408 The points x1,y1,z1 in the coordinate system defined by p1 are 409 transformed to x2,y2,z2 in the coordinate system defined by p2. 410 411 z1 is optional, if it is not set it is assumed to be zero (and 412 only x2 and y2 are returned). 413 414 In addition to converting between cartographic and geographic 415 projection coordinates, this function can take care of datum 416 shifts (which cannot be done using the __call__ method of the 417 Proj instances). It also allows for one of the coordinate 418 systems to be geographic (proj = 'latlong'). 419 420 If optional keyword 'radians' is True (default is False) and p1 421 is defined in geographic coordinate (pj.is_latlong() is True), 422 x1,y1 is interpreted as radians instead of the default degrees. 423 Similarly, if p2 is defined in geographic coordinates and 424 radians=True, x2, y2 are returned in radians instead of degrees. 425 if p1.is_latlong() and p2.is_latlong() both are False, the 426 radians keyword has no effect. 427 428 x,y and z can be numpy or regular python arrays, python 429 lists/tuples or scalars. Arrays are fastest. For projections in 430 geocentric coordinates, values of x and y are given in meters. 431 z is always meters. 432 433 Example usage: 434 435 >>> # projection 1: UTM zone 15, grs80 ellipse, NAD83 datum 436 >>> # (defined by epsg code 26915) 437 >>> p1 = Proj(init='epsg:26915') 438 >>> # projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum 439 >>> p2 = Proj(init='epsg:26715') 440 >>> # find x,y of Jefferson City, MO. 441 >>> x1, y1 = p1(-92.199881,38.56694) 442 >>> # transform this point to projection 2 coordinates. 443 >>> x2, y2 = transform(p1,p2,x1,y1) 444 >>> '%9.3f %11.3f' % (x1,y1) 445 '569704.566 4269024.671' 446 >>> '%9.3f %11.3f' % (x2,y2) 447 '569722.342 4268814.027' 448 >>> '%8.3f %5.3f' % p2(x2,y2,inverse=True) 449 ' -92.200 38.567' 450 >>> # process 3 points at a time in a tuple 451 >>> lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri 452 >>> lons = (-92.22,-94.72,-90.37) 453 >>> x1, y1 = p1(lons,lats) 454 >>> x2, y2 = transform(p1,p2,x1,y1) 455 >>> xy = x1+y1 456 >>> '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy 457 '567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005' 458 >>> xy = x2+y2 459 >>> '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy 460 '567721.149 351747.558 728569.133 4297989.112 4353489.644 4292106.305' 461 >>> lons, lats = p2(x2,y2,inverse=True) 462 >>> xy = lons+lats 463 >>> '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy 464 ' -92.220 -94.720 -90.370 38.830 39.320 38.750' 465 >>> # test datum shifting, installation of extra datum grid files. 466 >>> p1 = Proj(proj='latlong',datum='WGS84') 467 >>> x1 = -111.5; y1 = 45.25919444444 468 >>> p2 = Proj(proj="utm",zone=10,datum='NAD27') 469 >>> x2, y2 = transform(p1, p2, x1, y1) 470 >>> "%12.3f %12.3f" % (x2,y2) 471 ' 1402285.991 5076292.423' 472 """ 473 # process inputs, making copies that support buffer API. 474 inx, xisfloat, xislist, xistuple = _copytobuffer(x) 475 iny, yisfloat, yislist, yistuple = _copytobuffer(y) 476 if z is not None: 477 inz, zisfloat, zislist, zistuple = _copytobuffer(z) 478 else: 479 inz = None 480 # call pj_transform. inx,iny,inz buffers modified in place. 481 _proj._transform(p1,p2,inx,iny,inz,radians) 482 # if inputs were lists, tuples or floats, convert back. 483 outx = _convertback(xisfloat,xislist,xistuple,inx) 484 outy = _convertback(yisfloat,yislist,xistuple,iny) 485 if inz is not None: 486 outz = _convertback(zisfloat,zislist,zistuple,inz) 487 return outx, outy, outz 488 else: 489 return outx, outy
490
491 -def _copytobuffer(x):
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 # make sure x supports Buffer API and contains doubles. 500 isfloat = False; islist = False; istuple = False 501 # first, if it's a numpy array scalar convert to float 502 # (array scalars don't support buffer API) 503 if hasattr(x,'shape') and x.shape == (): x = float(x) 504 try: 505 # typecast numpy arrays to double. 506 # (this makes a copy - which is crucial 507 # since buffer is modified in place) 508 x.dtype.char 509 inx = x.astype('d') 510 except: 511 try: # perhaps they are Numeric/numarrays? 512 x.typecode() 513 inx = x.astype('d') 514 except: 515 # perhaps they are regular python arrays? 516 try: 517 x.typecode 518 inx = array('d',x) 519 except: 520 # try to convert to python array 521 # a list. 522 if type(x) == list: 523 inx = array('d',x) 524 islist = True 525 # a tuple. 526 elif type(x) == tuple: 527 inx = array('d',x) 528 istuple = True 529 # a scalar? 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
539 -def _convertback(isfloat,islist,istuple,inx):
540 # if inputs were lists, tuples or floats, convert back to original type. 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
550 -def _dict2string(projparams):
551 # convert a dict to a proj4 string. 552 pjargs = [] 553 for key,value in projparams.items(): 554 pjargs.append('+'+key+"="+str(value)+' ') 555 return ''.join(pjargs)
556
557 -class Geod(object):
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 # if initparams is a proj-type init string, 675 # convert to dict. 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 # merge this dict with kwargs dict. 685 kwargs = dict(list(kwargs.items()) + list(ellpsd.items())) 686 self.sphere = False 687 if 'ellps' in kwargs: 688 # ellipse name given, look up in pj_ellps dict 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 # a (semi-major axis) and one of 703 # b the semi-minor axis 704 # rf the reciprocal flattening 705 # f flattening 706 # es eccentricity squared 707 # must be given. 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 #msg='ellipse name or a, plus one of f,es,b must be given' 734 #raise ValueError(msg) 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 # process inputs, making copies that support buffer API. 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 # if inputs were lists, tuples or floats, convert back. 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 # process inputs, making copies that support buffer API. 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 # if inputs were lists, tuples or floats, convert back. 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 # distance increment. 859 del_s = dist/(npts+1) 860 # initialize output tuples. 861 del_s = dist/(npts+1) 862 # initialize output tuples. 863 lats = () 864 lons = () 865 # loop over intermediate points, compute lat/lons. 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
877 -def test():
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