SF.net SVN: matplotlib: [4071] trunk/toolkits/basemap/lib/matplotlib/ toolkits/basemap

[<ef...@us...>]

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 <jef...@no...>
+    __init__ method sets up projection information.
+    __call__ method compute transformations.
+    See docstrings for __init__ and __call__ for details.
+
+    Contact: Jeff Whitaker <jef...@no...>
     """
 
-    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.