Revision: 4071
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4071&view=rev
Author: efiring
Date: 2007-10-30 23:00:56 -0700 (Tue, 30 Oct 2007)
Log Message:
-----------
Numpification and some reformatting of pyproj and proj
Modified Paths:
--------------
trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/proj.py
trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/pyproj.py
Modified: trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/proj.py
===================================================================
--- trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/proj.py
2007-10-31 05:56:06 UTC (rev 4070)
+++ trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/proj.py
2007-10-31 06:00:56 UTC (rev 4071)
@@ -1,4 +1,4 @@
-import matplotlib.numerix as NX
+import numpy as npy
import pyproj
import math
@@ -6,36 +6,44 @@
_dg2rad = math.radians(1.)
_rad2dg = math.degrees(1.)
+_upper_right_out_of_bounds = (
+ 'the upper right corner of the plot is not in the map projection region')
+
+_lower_left_out_of_bounds = (
+ 'the lower left corner of the plot is not in the map projection region')
+
+
class Proj(object):
"""
- peforms cartographic transformations (converts from longitude,latitude
- to native map projection x,y coordinates and vice versa) using proj
- (http://proj.maptools.org/)
- Uses a pyrex generated C-interface to libproj.
-
- __init__ method sets up projection information.
- __call__ method compute transformations.
- See docstrings for __init__ and __call__ for details.
+ peforms cartographic transformations (converts from longitude,latitude
+ to native map projection x,y coordinates and vice versa) using proj
+ (http://proj.maptools.org/)
+ Uses a pyrex generated C-interface to libproj.
- Contact: Jeff Whitaker <[EMAIL PROTECTED]>
+ __init__ method sets up projection information.
+ __call__ method compute transformations.
+ See docstrings for __init__ and __call__ for details.
+
+ Contact: Jeff Whitaker <[EMAIL PROTECTED]>
"""
- def
__init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislatlon=True):
+ def __init__(self,projparams,llcrnrlon,llcrnrlat,
+ urcrnrlon,urcrnrlat,urcrnrislatlon=True):
"""
- initialize a Proj class instance.
+ initialize a Proj class instance.
- Input 'projparams' is a dictionary containing proj map
- projection control parameter key/value pairs.
- See the proj documentation (http://www.remotesensing.org/proj/)
- for details.
+ Input 'projparams' is a dictionary containing proj map
+ projection control parameter key/value pairs.
+ See the proj documentation (http://www.remotesensing.org/proj/)
+ for details.
- llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower left hand corner
- of projection region.
+ llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower
+ left hand corner of projection region.
- urcrnrlon,urcrnrlat are lon and lat (in degrees) of upper right hand corner
- of projection region if urcrnrislatlon=True (default). Otherwise,
- urcrnrlon,urcrnrlat are x,y in projection coordinates (units meters),
- assuming the lower left corner is x=0,y=0.
+ urcrnrlon,urcrnrlat are lon and lat (in degrees) of upper
+ right hand corner of projection region if urcrnrislatlon=True
+ (default). Otherwise, urcrnrlon,urcrnrlat are x,y in projection
+ coordinates (units meters), assuming the lower left corner is x=0,y=0.
"""
self.projparams = projparams
self.projection = projparams['proj']
@@ -63,8 +71,9 @@
llcrnrx = llcrnrlon
llcrnry = llcrnrlat
elif self.projection == 'ortho':
- if llcrnrlon == -180 and llcrnrlat == -90 and urcrnrlon == 180
and urcrnrlat == 90:
- self._fulldisk = True
+ if (llcrnrlon == -180 and llcrnrlat == -90 and
+ urcrnrlon == 180 and urcrnrlat == 90):
+ self._fulldisk = True
self._proj4 = pyproj.Proj(projparams)
llcrnrx = -self.rmajor
llcrnry = -self.rmajor
@@ -75,27 +84,28 @@
self._proj4 = pyproj.Proj(projparams)
llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
if llcrnrx > 1.e20 or llcrnry > 1.e20:
- raise ValueError('the lower left corner of the plot is not
in the map projection region')
+ raise ValueError(_lower_left_out_of_bounds)
elif self.projection == 'geos':
self._proj4 = pyproj.Proj(projparams)
# find major and minor axes of ellipse defining map proj region.
delta = 0.01
- lats = NX.arange(0,90,delta)
+ lats = npy.arange(0,90,delta)
lon_0 = projparams['lon_0']
- lons = lon_0*NX.ones(len(lats),'d')
+ lons = lon_0*npy.ones(len(lats),'d')
x, y = self._proj4(lons, lats)
yi = (y > 1.e20).tolist()
ny = yi.index(1)-1
height = y[ny]
- lons = NX.arange(lon_0,lon_0+90,delta)
- lats = NX.zeros(len(lons),'d')
+ lons = npy.arange(lon_0,lon_0+90,delta)
+ lats = npy.zeros(len(lons),'d')
x, y = self(lons, lats)
xi = (x > 1.e20).tolist()
nx = xi.index(1)-1
width = x[nx]
self._height = height
self._width = width
- if llcrnrlon == -180 and llcrnrlat == -90 and urcrnrlon == 180
and urcrnrlat == 90:
+ if (llcrnrlon == -180 and llcrnrlat == -90 and
+ urcrnrlon == 180 and urcrnrlat == 90):
self._fulldisk = True
llcrnrx = -width
llcrnry = -height
@@ -105,7 +115,7 @@
self._fulldisk = False
llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
if llcrnrx > 1.e20 or llcrnry > 1.e20:
- raise ValueError('the lower left corner of the plot is not
in the map projection region')
+ raise ValueError(_lower_left_out_of_bounds)
elif self.projection in ['moll','robin','sinu']:
self._proj4 = pyproj.Proj(projparams)
xtmp,urcrnry = self(projparams['lon_0'],90.)
@@ -116,7 +126,7 @@
# note that for 'cyl' x,y == lon,lat
self.projparams['x_0']=-llcrnrx
self.projparams['y_0']=-llcrnry
- # reset with x_0, y_0.
+ # reset with x_0, y_0.
if self.projection != 'cyl':
self._proj4 = pyproj.Proj(projparams)
llcrnry = 0.
@@ -136,7 +146,7 @@
else:
urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
if urcrnrx > 1.e20 or urcrnry > 1.e20:
- raise ValueError('the upper right corner of the plot
is not in the map projection region')
+ raise ValueError(_upper_right_out_of_bounds)
elif self.projection == 'geos':
if self._fulldisk:
urcrnrx = 2.*self._width
@@ -144,7 +154,7 @@
else:
urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
if urcrnrx > 1.e20 or urcrnry > 1.e20:
- raise ValueError('the upper right corner of the plot
is not in the map projection region')
+ raise ValueError(_upper_right_out_of_bounds)
elif self.projection in ['moll','robin','sinu']:
xtmp,urcrnry = self(projparams['lon_0'],90.)
urcrnrx,xtmp = self(projparams['lon_0']+180.,0)
@@ -172,24 +182,38 @@
self.ymax = llcrnry
self.ymin = urcrnry
- def __call__(self,x,y,inverse=False):
+ def __call__(self, *args, **kw):
+ # x,y,inverse=False):
"""
- Calling a Proj class instance with the arguments lon, lat will
- convert lon/lat (in degrees) to x/y native map projection
- coordinates (in meters). If optional keyword 'inverse' is
- True (default is False), the inverse transformation from x/y
- to lon/lat is performed.
+ Calling a Proj class instance with the arguments lon, lat will
+ convert lon/lat (in degrees) to x/y native map projection
+ coordinates (in meters). If optional keyword 'inverse' is
+ True (default is False), the inverse transformation from x/y
+ to lon/lat is performed.
- For cylindrical equidistant projection ('cyl'), this
- does nothing (i.e. x,y == lon,lat).
+ For cylindrical equidistant projection ('cyl'), this
+ does nothing (i.e. x,y == lon,lat).
- lon,lat can be either scalar floats or N arrays.
+ lon,lat can be either scalar floats or N arrays.
"""
+ if len(args) == 1:
+ xy = args[0]
+ onearray = True
+ else:
+ x,y = args
+ onearray = False
if self.projection == 'cyl': # for cyl x,y == lon,lat
- return x,y
+ if onearray:
+ return xy
+ else:
+ return x,y
+ inverse = kw.get('inverse', False)
+ if onearray:
+ outxy = self._proj4(xy, inverse=inverse)
+ else:
+ outx,outy = self._proj4(x, y, inverse=inverse)
if inverse:
- outx,outy = self._proj4(x,y,inverse=True)
- if self.projection in ['merc','mill']:
+ if self.projection in ['merc','mill']:
if self.projection == 'merc':
coslat = math.cos(math.radians(self.projparams['lat_ts']))
sinlat = math.sin(math.radians(self.projparams['lat_ts']))
@@ -199,12 +223,14 @@
# radius of curvature of the ellipse perpendicular to
# the plane of the meridian.
rcurv = self.rmajor*coslat/math.sqrt(1.-self.esq*sinlat**2)
- try: # x a scalar or an array
- outx = _rad2dg*(x/rcurv) + self.llcrnrlon
- except: # x a sequence
- outx = [_rad2dg*(xi/rcurv) + self.llcrnrlon for xi in x]
+ if onearray:
+ outxy[:,0] = _rad2dg*(xy[:,0]/rcurv) + self.llcrnrlon
+ else:
+ try: # x a scalar or an array
+ outx = _rad2dg*(x/rcurv) + self.llcrnrlon
+ except: # x a sequence
+ outx = [_rad2dg*(xi/rcurv) + self.llcrnrlon for xi in
x]
else:
- outx,outy = self._proj4(x,y)
if self.projection in ['merc','mill']:
if self.projection == 'merc':
coslat = math.cos(math.radians(self.projparams['lat_ts']))
@@ -215,28 +241,52 @@
# radius of curvature of the ellipse perpendicular to
# the plane of the meridian.
rcurv = self.rmajor*coslat/math.sqrt(1.-self.esq*sinlat**2)
- try: # x is a scalar or an array
- outx = rcurv*_dg2rad*(x-self.llcrnrlon)
- except: # x is a sequence.
- outx = [rcurv*_dg2rad*(xi-self.llcrnrlon) for xi in x]
- return outx,outy
+ if onearray:
+ outxy[:,0] = rcurv*_dg2rad*(xy[:,0]-self.llcrnrlon)
+ else:
+ try: # x is a scalar or an array
+ outx = rcurv*_dg2rad*(x-self.llcrnrlon)
+ except: # x is a sequence.
+ outx = [rcurv*_dg2rad*(xi-self.llcrnrlon) for xi in x]
+ if onearray:
+ return outxy
+ else:
+ return outx, outy
def makegrid(self,nx,ny,returnxy=False):
"""
- return arrays of shape (ny,nx) containing lon,lat coordinates of
- an equally spaced native projection grid.
- if returnxy=True, the x,y values of the grid are returned also.
+ return arrays of shape (ny,nx) containing lon,lat coordinates of
+ an equally spaced native projection grid.
+ if returnxy=True, the x,y values of the grid are returned also.
"""
dx = (self.urcrnrx-self.llcrnrx)/(nx-1)
dy = (self.urcrnry-self.llcrnry)/(ny-1)
- x = self.llcrnrx+dx*NX.indices((ny,nx),NX.Float32)[1,:,:]
- y = self.llcrnry+dy*NX.indices((ny,nx),NX.Float32)[0,:,:]
+ x = self.llcrnrx+dx*npy.indices((ny,nx),npy.float32)[1,:,:]
+ y = self.llcrnry+dy*npy.indices((ny,nx),npy.float32)[0,:,:]
lons, lats = self(x, y, inverse=True)
if returnxy:
return lons, lats, x, y
else:
return lons, lats
+ def makegrid3d(self,nx,ny,returnxy=False):
+ """
+ return array of shape (ny,nx, 2) containing lon,lat coordinates of
+ an equally spaced native projection grid.
+ if returnxy=True, the x,y values of the grid are returned also.
+ """
+ dx = (self.urcrnrx-self.llcrnrx)/(nx-1)
+ dy = (self.urcrnry-self.llcrnry)/(ny-1)
+ xy = npy.empty((ny,nx,2), npy.float64)
+ xy[...,0] = self.llcrnrx+dx*npy.indices((ny,nx),npy.float32)[1,:,:]
+ xy[...,1] = self.llcrnry+dy*npy.indices((ny,nx),npy.float32)[0,:,:]
+ lonlat = self(xy, inverse=True)
+ if returnxy:
+ return lonlat, xy
+ else:
+ return lonlat
+
+
if __name__ == "__main__":
params = {}
@@ -247,10 +297,10 @@
params['lon_0'] = -107
nx = 349; ny = 277; dx = 32463.41; dy = dx
awips221 = Proj(params,-145.5,1.0,(nx-1)*dx,(ny-1)*dy,urcrnrislatlon=False)
-# AWIPS grid 221 parameters
-# (from http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html)
+ # AWIPS grid 221 parameters
+ # (from http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html)
llcornerx, llcornery = awips221(-145.5,1.)
-# find 4 lon/lat corners of AWIPS grid 221.
+ # find 4 lon/lat corners of AWIPS grid 221.
llcornerx = 0.; llcornery = 0.
lrcornerx = dx*(nx-1); lrcornery = 0.
ulcornerx = 0.; ulcornery = dy*(ny-1)
@@ -270,13 +320,22 @@
print ' -68.318 0.897'
print ' -2.566 46.352'
print ' 148.639 46.635'
-# compute lons and lats for the whole AWIPS grid 221 (377x249).
+ # compute lons and lats for the whole AWIPS grid 221 (377x249).
import time; t1 = time.clock()
lons, lats = awips221.makegrid(nx,ny)
t2 = time.clock()
print 'compute lats/lons for all points on AWIPS 221 grid (%sx%s)' %(nx,ny)
print 'max/min lons'
- print min(NX.ravel(lons)),max(NX.ravel(lons))
+ print min(npy.ravel(lons)),max(npy.ravel(lons))
print 'max/min lats'
- print min(NX.ravel(lats)),max(NX.ravel(lats))
+ print min(npy.ravel(lats)),max(npy.ravel(lats))
print 'took',t2-t1,'secs'
+ print 'Same thing but with a single 3-D array'
+ t1 = time.clock()
+ lonlat, xy = awips221.makegrid3d(nx,ny, returnxy=True)
+ t2 = time.clock()
+ print 'took',t2-t1,'secs'
+
+ assert (lons==lonlat[...,0]).all(), "The longitudes are different"
+ assert (lats==lonlat[...,1]).all(), "The latitudes are different"
+
Modified: trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/pyproj.py
===================================================================
--- trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/pyproj.py
2007-10-31 05:56:06 UTC (rev 4070)
+++ trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/pyproj.py
2007-10-31 06:00:56 UTC (rev 4071)
@@ -1,5 +1,5 @@
"""
-Pyrex wrapper to provide python interfaces to
+Pyrex wrapper to provide python interfaces to
PROJ.4 (http://proj.maptools.org) functions.
Performs cartographic transformations and geodetic computations.
@@ -53,68 +53,69 @@
from array import array
from types import TupleType, ListType, NoneType
import os
+import numpy as npy
pyproj_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
set_datapath(pyproj_datadir)
class Proj(_Proj):
"""
-performs cartographic transformations (converts from
-longitude,latitude to native map projection x,y coordinates and
-vice versa) using proj (http://proj.maptools.org/)
+ performs cartographic transformations (converts from
+ longitude,latitude to native map projection x,y coordinates and
+ vice versa) using proj (http://proj.maptools.org/)
-A Proj class instance is initialized with proj map projection
-control parameter key/value pairs. The key/value pairs can
-either be passed in a dictionary, or as keyword arguments. See
-http://www.remotesensing.org/geotiff/proj_list for examples of
-key/value pairs defining different map projections.
+ A Proj class instance is initialized with proj map projection
+ control parameter key/value pairs. The key/value pairs can
+ either be passed in a dictionary, or as keyword arguments. See
+ http://www.remotesensing.org/geotiff/proj_list for examples of
+ key/value pairs defining different map projections.
-Calling a Proj class instance with the arguments lon, lat will
-convert lon/lat (in degrees) to x/y native map projection
-coordinates (in meters). If optional keyword 'inverse' is True
-(default is False), the inverse transformation from x/y to
-lon/lat is performed. If optional keyword 'radians' is True
-(default is False) lon/lat are interpreted as radians instead of
-degrees. If optional keyword 'errcheck' is True (default is
-False) an exception is raised if the transformation is invalid.
-If errcheck=False and the transformation is invalid, no
-exception is raised and 1.e30 is returned.
+ Calling a Proj class instance with the arguments lon, lat will
+ convert lon/lat (in degrees) to x/y native map projection
+ coordinates (in meters). If optional keyword 'inverse' is True
+ (default is False), the inverse transformation from x/y to
+ lon/lat is performed. If optional keyword 'radians' is True
+ (default is False) lon/lat are interpreted as radians instead of
+ degrees. If optional keyword 'errcheck' is True (default is
+ False) an exception is raised if the transformation is invalid.
+ If errcheck=False and the transformation is invalid, no
+ exception is raised and 1.e30 is returned.
-Works with numpy and regular python array objects, python
-sequences and scalars.
+ Works with numpy and regular python array objects, python
+ sequences and scalars.
"""
def __new__(self, projparams=None, **kwargs):
"""
-initialize a Proj class instance.
+ initialize a Proj class instance.
-Proj4 projection control parameters must either be given in a
-dictionary 'projparams' or as keyword arguments. See the proj
-documentation (http://proj.maptools.org) for more information
-about specifying projection parameters.
+ Proj4 projection control parameters must either be given in a
+ dictionary 'projparams' or as keyword arguments. See the proj
+ documentation (http://proj.maptools.org) for more information
+ about specifying projection parameters.
-Example usage:
+ Example usage:
->>> from pyproj import Proj
->>> p = Proj(proj='utm',zone=10,ellps='WGS84')
->>> x,y = p(-120.108, 34.36116666)
->>> print 'x=%9.3f y=%11.3f' % (x,y)
-x=765975.641 y=3805993.134
->>> print 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True)
-lon=-120.108 lat=34.361
->>> # do 3 cities at a time in a tuple (Fresno, LA, SF)
->>> lons = (-119.72,-118.40,-122.38)
->>> lats = (36.77, 33.93, 37.62 )
->>> x,y = p(lons, lats)
->>> print 'x: %9.3f %9.3f %9.3f' % x
-x: 792763.863 925321.537 554714.301
->>> print 'y: %9.3f %9.3f %9.3f' % y
-y: 4074377.617 3763936.941 4163835.303
->>> lons, lats = p(x, y, inverse=True) # inverse transform
->>> print 'lons: %8.3f %8.3f %8.3f' % lons
-lons: -119.720 -118.400 -122.380
->>> print 'lats: %8.3f %8.3f %8.3f' % lats
-lats: 36.770 33.930 37.620
+ >>> from pyproj import Proj
+ >>> p = Proj(proj='utm',zone=10,ellps='WGS84')
+ >>> x,y = p(-120.108, 34.36116666)
+ >>> print 'x=%9.3f y=%11.3f' % (x,y)
+ x=765975.641 y=3805993.134
+ >>> print 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True)
+ lon=-120.108 lat=34.361
+ >>> # do 3 cities at a time in a tuple (Fresno, LA, SF)
+ >>> lons = (-119.72,-118.40,-122.38)
+ >>> lats = (36.77, 33.93, 37.62 )
+ >>> x,y = p(lons, lats)
+ >>> print 'x: %9.3f %9.3f %9.3f' % x
+ x: 792763.863 925321.537 554714.301
+ >>> print 'y: %9.3f %9.3f %9.3f' % y
+ y: 4074377.617 3763936.941 4163835.303
+ >>> lons, lats = p(x, y, inverse=True) # inverse transform
+ >>> print 'lons: %8.3f %8.3f %8.3f' % lons
+ lons: -119.720 -118.400 -122.380
+ >>> print 'lats: %8.3f %8.3f %8.3f' % lats
+ lats: 36.770 33.930 37.620
"""
# if projparams is None, use kwargs.
if projparams is None:
@@ -130,25 +131,42 @@
projparams['units']='m'
return _Proj.__new__(self, projparams)
- def __call__(self,lon,lat,inverse=False,radians=False,errcheck=False):
+ def __call__(self, *args, **kw):
+ #,lon,lat,inverse=False,radians=False,errcheck=False):
"""
-Calling a Proj class instance with the arguments lon, lat will
-convert lon/lat (in degrees) to x/y native map projection
-coordinates (in meters). If optional keyword 'inverse' is True
-(default is False), the inverse transformation from x/y to
-lon/lat is performed. If optional keyword 'radians' is True
-(default is False) the units of lon/lat are radians instead of
-degrees. If optional keyword 'errcheck' is True (default is
-False) an exception is raised if the transformation is invalid.
-If errcheck=False and the transformation is invalid, no
-execption is raised and 1.e30 is returned.
+ Calling a Proj class instance with the arguments lon, lat will
+ convert lon/lat (in degrees) to x/y native map projection
+ coordinates (in meters). If optional keyword 'inverse' is True
+ (default is False), the inverse transformation from x/y to
+ lon/lat is performed. If optional keyword 'radians' is True
+ (default is False) the units of lon/lat are radians instead of
+ degrees. If optional keyword 'errcheck' is True (default is
+ False) an exception is raised if the transformation is invalid.
+ If errcheck=False and the transformation is invalid, no
+ exception is raised and 1.e30 is returned.
-Inputs should be doubles (they will be cast to doubles if they
-are not, causing a slight performance hit).
+ Instead of calling with lon, lat, a single ndarray of
+ shape n,2 may be used, and one of the same shape will
+ be returned; this is more efficient.
-Works with numpy and regular python array objects, python
-sequences and scalars, but is fastest for array objects.
+ Inputs should be doubles (they will be cast to doubles if they
+ are not, causing a slight performance hit).
+
+ Works with numpy and regular python array objects, python
+ sequences and scalars, but is fastest for array objects.
"""
+ inverse = kw.get('inverse', False)
+ radians = kw.get('radians', False)
+ errcheck = kw.get('errcheck', False)
+ if len(args) == 1:
+ latlon = npy.array(args[0], copy=True,
+ order='C', dtype=float, ndmin=2)
+ if inverse:
+ _Proj._invn(self, latlon, radians=radians, errcheck=errcheck)
+ else:
+ _Proj._fwdn(self, latlon, radians=radians, errcheck=errcheck)
+ return latlon
+ lon, lat = args
# process inputs, making copies that support buffer API.
inx, xisfloat, xislist, xistuple = _copytobuffer(lon)
iny, yisfloat, yislist, yistuple = _copytobuffer(lat)
@@ -172,68 +190,68 @@
def transform(p1, p2, x, y, z=None, radians=False):
"""
-x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False)
+ x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False)
-Transform points between two coordinate systems defined by the
-Proj instances p1 and p2.
+ Transform points between two coordinate systems defined by the
+ Proj instances p1 and p2.
-The points x1,y1,z1 in the coordinate system defined by p1 are
-transformed to x2,y2,z2 in the coordinate system defined by p2.
+ The points x1,y1,z1 in the coordinate system defined by p1 are
+ transformed to x2,y2,z2 in the coordinate system defined by p2.
-z1 is optional, if it is not set it is assumed to be zero (and
-only x2 and y2 are returned).
+ z1 is optional, if it is not set it is assumed to be zero (and
+ only x2 and y2 are returned).
-In addition to converting between cartographic and geographic
-projection coordinates, this function can take care of datum
-shifts (which cannot be done using the __call__ method of the
-Proj instances). It also allows for one of the coordinate
-systems to be geographic (proj = 'latlong').
+ In addition to converting between cartographic and geographic
+ projection coordinates, this function can take care of datum
+ shifts (which cannot be done using the __call__ method of the
+ Proj instances). It also allows for one of the coordinate
+ systems to be geographic (proj = 'latlong').
-If optional keyword 'radians' is True (default is False) and p1
-is defined in geographic coordinate (pj.is_latlong() is True),
-x1,y1 is interpreted as radians instead of the default degrees.
-Similarly, if p2 is defined in geographic coordinates and
-radians=True, x2, y2 are returned in radians instead of degrees.
-if p1.is_latlong() and p2.is_latlong() both are False, the
-radians keyword has no effect.
+ If optional keyword 'radians' is True (default is False) and p1
+ is defined in geographic coordinate (pj.is_latlong() is True),
+ x1,y1 is interpreted as radians instead of the default degrees.
+ Similarly, if p2 is defined in geographic coordinates and
+ radians=True, x2, y2 are returned in radians instead of degrees.
+ if p1.is_latlong() and p2.is_latlong() both are False, the
+ radians keyword has no effect.
-x,y and z can be numpy or regular python arrays, python
-lists/tuples or scalars. Arrays are fastest. For projections in
-geocentric coordinates, values of x and y are given in meters.
-z is always meters.
+ x,y and z can be numpy or regular python arrays, python
+ lists/tuples or scalars. Arrays are fastest. For projections in
+ geocentric coordinates, values of x and y are given in meters.
+ z is always meters.
-Example usage:
+ Example usage:
->>> # projection 1: UTM zone 15, grs80 ellipse, NAD83 datum
->>> # (defined by epsg code 26915)
->>> p1 = Proj(init='epsg:26915')
->>> # projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum
->>> p2 = Proj(init='epsg:26715')
->>> # find x,y of Jefferson City, MO.
->>> x1, y1 = p1(-92.199881,38.56694)
->>> # transform this point to projection 2 coordinates.
->>> x2, y2 = transform(p1,p2,x1,y1)
->>> print '%9.3f %11.3f' % (x1,y1)
-569704.566 4269024.671
->>> print '%9.3f %11.3f' % (x2,y2)
-569706.333 4268817.680
->>> print '%8.3f %5.3f' % p2(x2,y2,inverse=True)
- -92.200 38.567
->>> # process 3 points at a time in a tuple
->>> lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri
->>> lons = (-92.22,-94.72,-90.37)
->>> x1, y1 = p1(lons,lats)
->>> x2, y2 = transform(p1,p2,x1,y1)
->>> xy = x1+y1
->>> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
-567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005
->>> xy = x2+y2
->>> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
-567705.072 351727.113 728558.917 4297993.157 4353490.111 4292111.678
->>> lons, lats = p2(x2,y2,inverse=True)
->>> xy = lons+lats
->>> print '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy
- -92.220 -94.720 -90.370 38.830 39.320 38.750
+ >>> # projection 1: UTM zone 15, grs80 ellipse, NAD83 datum
+ >>> # (defined by epsg code 26915)
+ >>> p1 = Proj(init='epsg:26915')
+ >>> # projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum
+ >>> p2 = Proj(init='epsg:26715')
+ >>> # find x,y of Jefferson City, MO.
+ >>> x1, y1 = p1(-92.199881,38.56694)
+ >>> # transform this point to projection 2 coordinates.
+ >>> x2, y2 = transform(p1,p2,x1,y1)
+ >>> print '%9.3f %11.3f' % (x1,y1)
+ 569704.566 4269024.671
+ >>> print '%9.3f %11.3f' % (x2,y2)
+ 569706.333 4268817.680
+ >>> print '%8.3f %5.3f' % p2(x2,y2,inverse=True)
+ -92.200 38.567
+ >>> # process 3 points at a time in a tuple
+ >>> lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri
+ >>> lons = (-92.22,-94.72,-90.37)
+ >>> x1, y1 = p1(lons,lats)
+ >>> x2, y2 = transform(p1,p2,x1,y1)
+ >>> xy = x1+y1
+ >>> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
+ 567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005
+ >>> xy = x2+y2
+ >>> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
+ 567705.072 351727.113 728558.917 4297993.157 4353490.111 4292111.678
+ >>> lons, lats = p2(x2,y2,inverse=True)
+ >>> xy = lons+lats
+ >>> print '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy
+ -92.220 -94.720 -90.370 38.830 39.320 38.750
"""
# process inputs, making copies that support buffer API.
inx, xisfloat, xislist, xistuple = _copytobuffer(x)
@@ -254,13 +272,13 @@
return outx, outy
def _copytobuffer(x):
- """
-return a copy of x as an object that supports the python Buffer
-API (python array if input is float, list or tuple, numpy array
-if input is a numpy array). returns copyofx, isfloat, islist,
-istuple (islist is True if input is a list, istuple is true if
-input is a tuple, isfloat is true if input is a float).
"""
+ return a copy of x as an object that supports the python Buffer
+ API (python array if input is float, list or tuple, numpy array
+ if input is a numpy array). returns copyofx, isfloat, islist,
+ istuple (islist is True if input is a list, istuple is true if
+ input is a tuple, isfloat is true if input is a float).
+ """
# make sure x supports Buffer API and contains doubles.
isfloat = False; islist = False; istuple = False
# first, if it's a numpy array scalar convert to float
@@ -281,7 +299,7 @@
try:
x.typecode
inx = array('d',x)
- except:
+ except:
# try to convert to python array
# a list.
if type(x) is ListType:
@@ -315,107 +333,107 @@
class Geod(_Geod):
"""
-performs forward and inverse geodetic, or Great Circle,
-computations. The forward computation (using the 'fwd' method)
-involves determining latitude, longitude and back azimuth of a
-terminus point given the latitude and longitude of an initial
-point, plus azimuth and distance. The inverse computation (using
-the 'inv' method) involves determining the forward and back
-azimuths and distance given the latitudes and longitudes of an
-initial and terminus point.
+ performs forward and inverse geodetic, or Great Circle,
+ computations. The forward computation (using the 'fwd' method)
+ involves determining latitude, longitude and back azimuth of a
+ terminus point given the latitude and longitude of an initial
+ point, plus azimuth and distance. The inverse computation (using
+ the 'inv' method) involves determining the forward and back
+ azimuths and distance given the latitudes and longitudes of an
+ initial and terminus point.
"""
def __new__(self, initparams=None, **kwargs):
"""
-initialize a Geod class instance.
+ initialize a Geod class instance.
-Geodetic parameters for specifying the ellipsoid or sphere to
-use must either be given in a dictionary 'initparams' or as
-keyword arguments. Following is a list of the ellipsoids that
-may be defined using the 'ellps' keyword:
+ Geodetic parameters for specifying the ellipsoid or sphere to
+ use must either be given in a dictionary 'initparams' or as
+ keyword arguments. Following is a list of the ellipsoids that
+ may be defined using the 'ellps' keyword:
- MERIT a=6378137.0 rf=298.257 MERIT 1983
- SGS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85
- GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980)
- IAU76 a=6378140.0 rf=298.257 IAU 1976
- airy a=6377563.396 b=6356256.910 Airy 1830
- APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
- NWL9D a=6378145.0. rf=298.25 Naval Weapons Lab., 1965
-mod_airy a=6377340.189 b=6356034.446 Modified Airy
- andrae a=6377104.43 rf=300.0 Andrae 1876 (Den., Iclnd.)
- aust_SA a=6378160.0 rf=298.25 Australian Natl & S. Amer. 1969
- GRS67 a=6378160.0 rf=298.2471674270 GRS 67(IUGG 1967)
- bessel a=6377397.155 rf=299.1528128 Bessel 1841
-bess_nam a=6377483.865 rf=299.1528128 Bessel 1841 (Namibia)
- clrk66 a=6378206.4 b=6356583.8 Clarke 1866
- clrk80 a=6378249.145 rf=293.4663 Clarke 1880 mod.
- CPM a=6375738.7 rf=334.29 Comm. des Poids et Mesures 1799
- delmbr a=6376428. rf=311.5 Delambre 1810 (Belgium)
- engelis a=6378136.05 rf=298.2566 Engelis 1985
- evrst30 a=6377276.345 rf=300.8017 Everest 1830
- evrst48 a=6377304.063 rf=300.8017 Everest 1948
- evrst56 a=6377301.243 rf=300.8017 Everest 1956
- evrst69 a=6377295.664 rf=300.8017 Everest 1969
- evrstSS a=6377298.556 rf=300.8017 Everest (Sabah & Sarawak)
- fschr60 a=6378166. rf=298.3 Fischer (Mercury Datum) 1960
-fschr60m a=6378155. rf=298.3 Modified Fischer 1960
- fschr68 a=6378150. rf=298.3 Fischer 1968
- helmert a=6378200. rf=298.3 Helmert 1906
- hough a=6378270.0 rf=297. Hough
- intl a=6378388.0 rf=297. International 1909 (Hayford)
- krass a=6378245.0 rf=298.3 Krassovsky, 1942
- kaula a=6378163. rf=298.24 Kaula 1961
- lerch a=6378139. rf=298.257 Lerch 1979
- mprts a=6397300. rf=191. Maupertius 1738
-new_intl a=6378157.5 b=6356772.2 New International 1967
- plessis a=6376523. b=6355863. Plessis 1817 (France)
- SEasia a=6378155.0 b=6356773.3205 Southeast Asia
- walbeck a=6376896.0 b=6355834.8467 Walbeck
- WGS60 a=6378165.0 rf=298.3 WGS 60
- WGS66 a=6378145.0 rf=298.25 WGS 66
- WGS72 a=6378135.0 rf=298.26 WGS 72
- WGS84 a=6378137.0 rf=298.257223563 WGS 84
- sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)
+ MERIT a=6378137.0 rf=298.257 MERIT 1983
+ SGS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85
+ GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980)
+ IAU76 a=6378140.0 rf=298.257 IAU 1976
+ airy a=6377563.396 b=6356256.910 Airy 1830
+ APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
+ NWL9D a=6378145.0. rf=298.25 Naval Weapons Lab., 1965
+ mod_airy a=6377340.189 b=6356034.446 Modified Airy
+ andrae a=6377104.43 rf=300.0 Andrae 1876 (Den., Iclnd.)
+ aust_SA a=6378160.0 rf=298.25 Australian Natl & S. Amer.
1969
+ GRS67 a=6378160.0 rf=298.2471674270 GRS 67(IUGG 1967)
+ bessel a=6377397.155 rf=299.1528128 Bessel 1841
+ bess_nam a=6377483.865 rf=299.1528128 Bessel 1841 (Namibia)
+ clrk66 a=6378206.4 b=6356583.8 Clarke 1866
+ clrk80 a=6378249.145 rf=293.4663 Clarke 1880 mod.
+ CPM a=6375738.7 rf=334.29 Comm. des Poids et Mesures
1799
+ delmbr a=6376428. rf=311.5 Delambre 1810 (Belgium)
+ engelis a=6378136.05 rf=298.2566 Engelis 1985
+ evrst30 a=6377276.345 rf=300.8017 Everest 1830
+ evrst48 a=6377304.063 rf=300.8017 Everest 1948
+ evrst56 a=6377301.243 rf=300.8017 Everest 1956
+ evrst69 a=6377295.664 rf=300.8017 Everest 1969
+ evrstSS a=6377298.556 rf=300.8017 Everest (Sabah & Sarawak)
+ fschr60 a=6378166. rf=298.3 Fischer (Mercury Datum) 1960
+ fschr60m a=6378155. rf=298.3 Modified Fischer 1960
+ fschr68 a=6378150. rf=298.3 Fischer 1968
+ helmert a=6378200. rf=298.3 Helmert 1906
+ hough a=6378270.0 rf=297. Hough
+ intl a=6378388.0 rf=297. International 1909 (Hayford)
+ krass a=6378245.0 rf=298.3 Krassovsky, 1942
+ kaula a=6378163. rf=298.24 Kaula 1961
+ lerch a=6378139. rf=298.257 Lerch 1979
+ mprts a=6397300. rf=191. Maupertius 1738
+ new_intl a=6378157.5 b=6356772.2 New International 1967
+ plessis a=6376523. b=6355863. Plessis 1817 (France)
+ SEasia a=6378155.0 b=6356773.3205 Southeast Asia
+ walbeck a=6376896.0 b=6355834.8467 Walbeck
+ WGS60 a=6378165.0 rf=298.3 WGS 60
+ WGS66 a=6378145.0 rf=298.25 WGS 66
+ WGS72 a=6378135.0 rf=298.26 WGS 72
+ WGS84 a=6378137.0 rf=298.257223563 WGS 84
+ sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)
-The parameters of the ellipsoid may also be set directly using
-the 'a' (semi-major or equatorial axis radius) keyword, and
-any one of the following keywords: 'b' (semi-minor,
-or polar axis radius), 'e' (eccentricity), 'es' (eccentricity
-squared), 'f' (flattening), or 'rf' (reciprocal flattening).
+ The parameters of the ellipsoid may also be set directly using
+ the 'a' (semi-major or equatorial axis radius) keyword, and
+ any one of the following keywords: 'b' (semi-minor,
+ or polar axis radius), 'e' (eccentricity), 'es' (eccentricity
+ squared), 'f' (flattening), or 'rf' (reciprocal flattening).
-See the proj documentation (http://proj.maptools.org) for more
-information about specifying ellipsoid parameters (specifically,
-the chapter 'Specifying the Earth's figure' in the main Proj
-users manual).
+ See the proj documentation (http://proj.maptools.org) for more
+ information about specifying ellipsoid parameters (specifically,
+ the chapter 'Specifying the Earth's figure' in the main Proj
+ users manual).
-Example usage:
+ Example usage:
->>> from pyproj import Geod
->>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
->>> # specify the lat/lons of some cities.
->>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
->>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
->>> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
->>> london_lat = 51.+(32./60.); london_lon = -(5./60.)
->>> # compute forward and back azimuths, plus distance
->>> # between Boston and Portland.
->>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
->>> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
--66.531 75.654 4164192.708
->>> # compute latitude, longitude and back azimuth of Portland,
->>> # given Boston lat/lon, forward azimuth and distance to Portland.
->>> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
->>> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
-45.517 -123.683 75.654
->>> # compute the azimuths, distances from New York to several
->>> # cities (pass a list)
->>> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
->>> lons2 = [boston_lon, portland_lon, london_lon]
->>> lats2 = [boston_lat, portland_lat, london_lat]
->>> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
->>> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" %
(faz,baz,d)
- 54.663 -123.448 288303.720
--65.463 79.342 4013037.318
- 51.254 -71.576 5579916.649
+ >>> from pyproj import Geod
+ >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
+ >>> # specify the lat/lons of some cities.
+ >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
+ >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
+ >>> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
+ >>> london_lat = 51.+(32./60.); london_lon = -(5./60.)
+ >>> # compute forward and back azimuths, plus distance
+ >>> # between Boston and Portland.
+ >>> az12,az21,dist =
g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
+ >>> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
+ -66.531 75.654 4164192.708
+ >>> # compute latitude, longitude and back azimuth of Portland,
+ >>> # given Boston lat/lon, forward azimuth and distance to Portland.
+ >>> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
+ >>> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
+ 45.517 -123.683 75.654
+ >>> # compute the azimuths, distances from New York to several
+ >>> # cities (pass a list)
+ >>> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
+ >>> lons2 = [boston_lon, portland_lon, london_lon]
+ >>> lats2 = [boston_lat, portland_lat, london_lat]
+ >>> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
+ >>> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" %
(faz,baz,d)
+ 54.663 -123.448 288303.720
+ -65.463 79.342 4013037.318
+ 51.254 -71.576 5579916.649
"""
# if projparams is None, use kwargs.
if initparams is None:
@@ -433,16 +451,16 @@
def fwd(self, lons, lats, az, dist, radians=False):
"""
-forward transformation - Returns longitudes, latitudes and back
-azimuths of terminus points given longitudes (lons) and
-latitudes (lats) of initial points, plus forward azimuths (az)
-and distances (dist).
+ forward transformation - Returns longitudes, latitudes and back
+ azimuths of terminus points given longitudes (lons) and
+ latitudes (lats) of initial points, plus forward azimuths (az)
+ and distances (dist).
-Works with numpy and regular python array objects, python
-sequences and scalars.
+ Works with numpy and regular python array objects, python
+ sequences and scalars.
-if radians=True, lons/lats and azimuths are radians instead of
-degrees. Distances are in meters.
+ if radians=True, lons/lats and azimuths are radians instead of
+ degrees. Distances are in meters.
"""
# process inputs, making copies that support buffer API.
inx, xisfloat, xislist, xistuple = _copytobuffer(lons)
@@ -459,15 +477,15 @@
def inv(self, lons1, lats1, lons2, lats2, radians=False):
"""
-inverse transformation - Returns forward and back azimuths, plus
-distances between initial points (specified by lons1, lats1) and
-terminus points (specified by lons2, lats2).
+ inverse transformation - Returns forward and back azimuths, plus
+ distances between initial points (specified by lons1, lats1) and
+ terminus points (specified by lons2, lats2).
-Works with numpy and regular python array objects, python
-sequences and scalars.
+ Works with numpy and regular python array objects, python
+ sequences and scalars.
-if radians=True, lons/lats and azimuths are radians instead of
-degrees. Distances are in meters.
+ if radians=True, lons/lats and azimuths are radians instead of
+ degrees. Distances are in meters.
"""
# process inputs, making copies that support buffer API.
inx, xisfloat, xislist, xistuple = _copytobuffer(lons1)
@@ -484,35 +502,35 @@
def npts(self, lon1, lat1, lon2, lat2, npts, radians=False):
"""
-Given a single initial point and terminus point (specified by
-python floats lon1,lat1 and lon2,lat2), returns a list of
-longitude/latitude pairs describing npts equally spaced
-intermediate points along the geodesic between the initial and
-terminus points.
+ Given a single initial point and terminus point (specified by
+ python floats lon1,lat1 and lon2,lat2), returns a list of
+ longitude/latitude pairs describing npts equally spaced
+ intermediate points along the geodesic between the initial and
+ terminus points.
-if radians=True, lons/lats are radians instead of degrees.
+ if radians=True, lons/lats are radians instead of degrees.
-Example usage:
+ Example usage:
->>> from pyproj import Geod
->>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
->>> # specify the lat/lons of Boston and Portland.
->>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
->>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
->>> # find ten equally spaced points between Boston and Portland.
->>> lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
->>> for lon,lat in lonlats: print '%6.3f %7.3f' % (lat, lon)
-43.528 -75.414
-44.637 -79.883
-45.565 -84.512
-46.299 -89.279
-46.830 -94.156
-47.149 -99.112
-47.251 -104.106
-47.136 -109.100
-46.805 -114.051
-46.262 -118.924
- """
+ >>> from pyproj import Geod
+ >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
+ >>> # specify the lat/lons of Boston and Portland.
+ >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
+ >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
+ >>> # find ten equally spaced points between Boston and Portland.
+ >>> lonlats =
g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
+ >>> for lon,lat in lonlats: print '%6.3f %7.3f' % (lat, lon)
+ 43.528 -75.414
+ 44.637 -79.883
+ 45.565 -84.512
+ 46.299 -89.279
+ 46.830 -94.156
+ 47.149 -99.112
+ 47.251 -104.106
+ 47.136 -109.100
+ 46.805 -114.051
+ 46.262 -118.924
+ """
lons, lats = _Geod._npts(self,lon1,lat1,lon2,lat2,npts,radians=radians)
return zip(lons, lats)
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/
_______________________________________________
Matplotlib-checkins mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/matplotlib-checkins
