Revision: 5737
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5737&view=rev
Author: jswhit
Date: 2008-07-11 08:03:53 -0700 (Fri, 11 Jul 2008)
Log Message:
-----------
don't use shapefiles.
Modified Paths:
--------------
trunk/toolkits/basemap/examples/plot_tissot.py
Modified: trunk/toolkits/basemap/examples/plot_tissot.py
===================================================================
--- trunk/toolkits/basemap/examples/plot_tissot.py 2008-07-11 14:36:49 UTC
(rev 5736)
+++ trunk/toolkits/basemap/examples/plot_tissot.py 2008-07-11 15:03:53 UTC
(rev 5737)
@@ -1,6 +1,7 @@
import numpy as np
import matplotlib.pyplot as plt
-from mpl_toolkits.basemap import Basemap as Basemap
+from mpl_toolkits.basemap import Basemap
+from mpl_toolkits.basemap import pyproj
from matplotlib.patches import Polygon
# Tissot's Indicatrix (http://en.wikipedia.org/wiki/Tissot's_Indicatrix).
@@ -11,90 +12,62 @@
# Tissot's indicatrix have all unit area, although their shapes and
# orientations vary with location.
-# adapted from http://www.perrygeo.net/wordpress/?p=4
+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 new figure
-fig=plt.figure()
-m = Basemap(llcrnrlon=-180,llcrnrlat=-80,urcrnrlon=180,urcrnrlat=80,
- projection='cyl')
-shp_info = m.readshapefile('tissot','tissot',drawbounds=True)
-ax = plt.gca()
-for nshape,seg in enumerate(m.tissot):
- poly = Polygon(seg,facecolor='green',zorder=10)
- ax.add_patch(poly)
-# draw meridians and parallels.
-m.drawparallels(np.arange(-90,91,30),labels=[1,0,0,0])
-m.drawmeridians(np.arange(-180,180,60),labels=[0,0,0,1])
-m.drawcoastlines()
-m.fillcontinents()
-plt.title('Tissot Diagram - Cylindrical Equal Area')
-print 'plot Cylindrical Equidistant Equal Area Tissot diagram ...'
+# create Basemap instances with several different projections
+m1 = Basemap2(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,
+ 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')
-# create new figure
-fig=plt.figure()
-m = Basemap(llcrnrlon=-180,llcrnrlat=-70,urcrnrlon=180,urcrnrlat=70,
- projection='merc',lat_ts=20)
-shp_info = m.readshapefile('tissot','tissot',drawbounds=True)
-ax = plt.gca()
-for nshape,seg in enumerate(m.tissot):
- poly = Polygon(seg,facecolor='green',zorder=10)
- ax.add_patch(poly)
-# draw meridians and parallels.
-m.drawparallels(np.arange(-90,91,30),labels=[1,0,0,0])
-m.drawmeridians(np.arange(-180,180,60),labels=[0,0,0,1])
-m.drawcoastlines()
-m.fillcontinents()
-plt.title('Tissot Diagram - Mercator Conformal')
-print 'plot Mercator Conformal Tissot diagram ...'
+for m in [m1,m2,m3,m4,m5]:
+ # make a new figure.
+ fig = plt.figure()
+ ax = plt.gca()
+ # draw "circles" at specified longitudes and latitudes.
+ for parallel in range(-60,61,30):
+ for meridian in range(-165,166,30):
+ seg = m.tissot(meridian,parallel,6,100)
+ poly = Polygon(seg,facecolor='green',zorder=10)
+ ax.add_patch(poly)
+ # draw meridians and parallels.
+ m.drawparallels(np.arange(-60,61,30))
+ m.drawmeridians(np.arange(-180,180,60))
+ # draw coastlines, fill continents, plot title.
+ m.drawcoastlines()
+ m.fillcontinents()
+ title = 'Tissot Diagram: projection = %s' % m.projection
+ print title
+ plt.title(title)
-# create new figure
-fig=plt.figure()
-m = Basemap(lon_0=-60,lat_0=45,projection='ortho')
-shp_info = m.readshapefile('tissot','tissot',drawbounds=False)
-ax = plt.gca()
-for nshape,seg in enumerate(m.tissot):
- xx,yy = zip(*seg)
- if max(xx) < 1.e20 and max(yy) < 1.e20:
- poly = Polygon(seg,facecolor='green',zorder=10)
- ax.add_patch(poly)
-m.drawcoastlines()
-m.fillcontinents()
-m.drawparallels(np.arange(-90,91,30))
-m.drawmeridians(np.arange(-180,180,30))
-plt.title('Tissot Diagram - Orthographic')
-m.drawmapboundary()
-plt.gca().set_frame_on(True)
-print 'plot Orthographic Tissot diagram ...'
-
-# create new figure
-fig=plt.figure()
-m = Basemap(lon_0=270,lat_0=90,boundinglat=10,projection='npstere')
-shp_info = m.readshapefile('tissot','tissot',drawbounds=True)
-ax = plt.gca()
-for nshape,seg in enumerate(m.tissot):
- poly = Polygon(seg,facecolor='green',zorder=10)
- ax.add_patch(poly)
-# draw meridians and parallels.
-m.drawparallels(np.arange(-90,91,30),labels=[1,0,0,0])
-m.drawmeridians(np.arange(-180,180,30),labels=[0,0,0,1])
-m.drawcoastlines()
-m.fillcontinents()
-plt.title('Tissot Diagram - North Polar Stereographic Conformal')
-print 'plot North Polar Stereographic Conformal Tissot diagram ...'
-
-# create new figure
-fig=plt.figure()
-m = Basemap(lon_0=270,lat_0=90,boundinglat=10,projection='nplaea')
-shp_info = m.readshapefile('tissot','tissot',drawbounds=True)
-ax = plt.gca()
-for nshape,seg in enumerate(m.tissot):
- poly = Polygon(seg,facecolor='green',zorder=10)
- ax.add_patch(poly)
-# draw meridians and parallels.
-m.drawparallels(np.arange(-90,91,30),labels=[1,0,0,0])
-m.drawmeridians(np.arange(-180,180,30),labels=[0,0,0,1])
-m.drawcoastlines()
-m.fillcontinents()
-plt.title('Tissot Diagram - North Polar Lambert Azimuthal Equal Area')
-print 'plot North Polar Lambert Azimuthal Equal Area Tissot diagram ...'
plt.show()
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
