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

[<js...@us...>]

Revision: 5741
          http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5741&view=rev
Author:   jswhit
Date:     2008-07-11 08:39:08 -0700 (Fri, 11 Jul 2008)

Log Message:
-----------
added tissot method.

Modified Paths:
--------------
    trunk/toolkits/basemap/Changelog
    trunk/toolkits/basemap/examples/plot_tissot.py
    trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py

Modified: trunk/toolkits/basemap/Changelog
===================================================================
--- trunk/toolkits/basemap/Changelog	2008-07-11 15:05:09 UTC (rev 5740)
+++ trunk/toolkits/basemap/Changelog	2008-07-11 15:39:08 UTC (rev 5741)
@@ -1,3 +1,5 @@
+           * added "tissot" method for generating Tissot's indicatrix
+             (see example plot_tissot.py).
            * fixed processing of coastlines for gnomonic projection.
            * don't try to use PyNIO in NetCDFFile (it was causing
              too many suprises).

Modified: trunk/toolkits/basemap/examples/plot_tissot.py
===================================================================
--- trunk/toolkits/basemap/examples/plot_tissot.py	2008-07-11 15:05:09 UTC (rev 5740)
+++ trunk/toolkits/basemap/examples/plot_tissot.py	2008-07-11 15:39:08 UTC (rev 5741)
@@ -12,43 +12,14 @@
 # Tissot's indicatrix have all unit area, although their shapes and 
 # orientations vary with location.
 
-class Basemap2(Basemap):
-    def tissot(self,lon_0,lat_0,radius_deg,npts):
-        # create list of npts lon,lat pairs that are equidistant on the
-        # surface of the earth from central point lon_0,lat_0
-        # and has radius along lon_0 of radius_deg degrees of latitude.
-        # points are transformed to map projection coordinates.
-        # The ellipse that list of points represents is a
-        # Tissot's indicatrix
-        # (http://en.wikipedia.org/wiki/Tissot%27s_Indicatrix),
-        # which when drawn on a map shows the distortion
-        # inherent in the map projection.
-        g = pyproj.Geod(a=self.rmajor,b=self.rminor)
-        az12,az21,dist = g.inv(lon_0,lat_0,lon_0,lat_0+radius_deg)
-        seg = [self(lon_0,lat_0+radius_deg)]
-        delaz = 360./npts
-        az = az12
-        for n in range(npts):
-            az = az+delaz
-            # skip segments along equator (Geod can't handel equatorial arcs)
-            if np.allclose(0.,lat_0) and (np.allclose(90.,az) or np.allclose(270.,az)):
-                continue
-            else:
-                lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
-            x,y = self(lon,lat)
-            # add segment if it is in the map projection region.
-            if x < 1.e20 and y < 1.e20:
-                seg.append((x,y))
-        return seg
-
 # create Basemap instances with several different projections
-m1 = Basemap2(llcrnrlon=-180,llcrnrlat=-80,urcrnrlon=180,urcrnrlat=80,
+m1 = Basemap(llcrnrlon=-180,llcrnrlat=-80,urcrnrlon=180,urcrnrlat=80,
               projection='cyl')
-m2 = Basemap2(lon_0=-60,lat_0=45,projection='ortho')
-m3 = Basemap2(llcrnrlon=-180,llcrnrlat=-70,urcrnrlon=180,urcrnrlat=70,
+m2 = Basemap(lon_0=-60,lat_0=45,projection='ortho')
+m3 = Basemap(llcrnrlon=-180,llcrnrlat=-70,urcrnrlon=180,urcrnrlat=70,
               projection='merc',lat_ts=20,rsphere=(6378137.0,6356752.3142))
-m4 = Basemap2(lon_0=270,lat_0=90,boundinglat=10,projection='npstere')
-m5 = Basemap2(lon_0=270,lat_0=90,boundinglat=10,projection='nplaea')
+m4 = Basemap(lon_0=270,lat_0=90,boundinglat=10,projection='npstere')
+m5 = Basemap(lon_0=270,lat_0=90,boundinglat=10,projection='nplaea')
 
 for m in [m1,m2,m3,m4,m5]:
     # make a new figure.

Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2008-07-11 15:05:09 UTC (rev 5740)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2008-07-11 15:39:08 UTC (rev 5741)
@@ -2146,6 +2146,34 @@
             if v == ([], []): del linecolls[k]
         return linecolls
 
+    def tissot(self,lon_0,lat_0,radius_deg,npts):
+        """
+        create list of ``npts`` x,y pairs that are equidistant on the
+        surface of the earth from central point ``lon_0,lat_0`` and form
+        an ellipse with radius of ``radius_deg`` degrees of latitude along
+        longitude ``lon_0``.
+        The ellipse represents a Tissot's indicatrix
+        (http://en.wikipedia.org/wiki/Tissot%27s_Indicatrix),
+        which when drawn on a map shows the distortion
+        inherent in the map projection."""
+        g = pyproj.Geod(a=self.rmajor,b=self.rminor)
+        az12,az21,dist = g.inv(lon_0,lat_0,lon_0,lat_0+radius_deg)
+        seg = [self(lon_0,lat_0+radius_deg)]
+        delaz = 360./npts
+        az = az12
+        for n in range(npts):
+            az = az+delaz
+            # skip segments along equator (Geod can't handel equatorial arcs)
+            if np.allclose(0.,lat_0) and (np.allclose(90.,az) or np.allclose(270.,az)):
+                continue
+            else:
+                lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
+            x,y = self(lon,lat)
+            # add segment if it is in the map projection region.
+            if x < 1.e20 and y < 1.e20:
+                seg.append((x,y))
+        return seg
+
     def gcpoints(self,lon1,lat1,lon2,lat2,npoints):
         """
         compute ``points`` points along a great circle with endpoints


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.