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.
-------------------------------------------------------------------------
Sponsored by: SourceForge.net Community Choice Awards: VOTE NOW!
Studies have shown that voting for your favorite open source project,
along with a healthy diet, reduces your potential for chronic lameness
and boredom. Vote Now at http://www.sourceforge.net/community/cca08
_______________________________________________
Matplotlib-checkins mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/matplotlib-checkins
