Menu
▾
▴
SF.net SVN: matplotlib:[6719] trunk/toolkits/basemap
|
From: <js...@us...> - 2008-12-31 14:37:48
|
Revision: 6719
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=6719&view=rev
Author: jswhit
Date: 2008-12-31 14:37:44 +0000 (Wed, 31 Dec 2008)
Log Message:
-----------
whole world aeqd added
Modified Paths:
--------------
trunk/toolkits/basemap/Changelog
trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py
Modified: trunk/toolkits/basemap/Changelog
===================================================================
--- trunk/toolkits/basemap/Changelog 2008-12-31 13:20:50 UTC (rev 6718)
+++ trunk/toolkits/basemap/Changelog 2008-12-31 14:37:44 UTC (rev 6719)
@@ -1,4 +1,7 @@
version 0.99.3 (not yet released)
+ * if upper right/lower left corners nor width/height given for
+ azimuthal equidistant ('aeqd') the whole world is drawn in
+ a circle (only works for perfect spheres, not ellipsoids).
* have setup.py check for already installed pyshapelib (just
like it does for httplib2 and pydap).
* Basemap will now look for it's data in BASEMAPDATA.
Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py 2008-12-31 13:20:50 UTC (rev 6718)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py 2008-12-31 14:37:44 UTC (rev 6719)
@@ -625,10 +625,17 @@
raise ValueError, 'must specify lat_0 and lon_0 for Azimuthal Equidistant basemap'
if not using_corners:
if width is None or height is None:
- raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
+ self._fulldisk = True
+ llcrnrlon = -180.
+ llcrnrlat = -90.
+ urcrnrlon = 180
+ urcrnrlat = 90.
+ else:
+ self._fulldisk = False
if lon_0 is None or lat_0 is None:
raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
+ if not self._fulldisk:
+ llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection in _cylproj:
@@ -1046,12 +1053,8 @@
if self.projection in ['ortho','geos']:
# circular region.
thetas = np.linspace(0.,2.*np.pi,2*nx*ny)[:-1]
- if self.projection == 'ortho':
- rminor = self.rmajor
- rmajor = self.rmajor
- else:
- rminor = self._height
- rmajor = self._width
+ rminor = self._height
+ rmajor = self._width
x = rmajor*np.cos(thetas) + rmajor
y = rminor*np.sin(thetas) + rminor
b = np.empty((len(x),2),np.float64)
@@ -1075,6 +1078,16 @@
projparms['x_0']=-llcrnrx
projparms['y_0']=-llcrnry
maptran = pyproj.Proj(projparms)
+ elif self.projection == 'aeqd' and self._fulldisk:
+ # circular region.
+ thetas = np.linspace(0.,2.*np.pi,2*nx*ny)[:-1]
+ rminor = self._height
+ rmajor = self._width
+ x = rmajor*np.cos(thetas) + rmajor
+ y = rminor*np.sin(thetas) + rminor
+ b = np.empty((len(x),2),np.float64)
+ b[:,0]=x; b[:,1]=y
+ boundaryxy = _geoslib.Polygon(b)
elif self.projection in _pseudocyl:
# quasi-elliptical region.
lon_0 = self.projparams['lon_0']
@@ -1188,11 +1201,10 @@
elif ax is None and self.ax is not None:
ax = self.ax
limb = None
- if self.projection == 'ortho':
- limb = Circle((self.rmajor,self.rmajor),self.rmajor)
- elif self.projection == 'geos':
+ if self.projection in ['ortho','geos'] or (self.projection=='aeqd' and\
+ self._fulldisk):
limb = Ellipse((self._width,self._height),2.*self._width,2.*self._height)
- if self.projection in ['ortho','geos'] and self._fulldisk:
+ if self.projection in ['ortho','geos','aeqd'] and self._fulldisk:
# elliptical region.
ax.add_patch(limb)
if fill_color is None:
@@ -1822,7 +1834,7 @@
linecolls[circ] = (lines,[])
# draw labels for parallels
# parallels not labelled for fulldisk orthographic or geostationary
- if self.projection in ['ortho','geos','vandg'] and max(labels):
+ if self.projection in ['ortho','geos','vandg','aeqd'] and max(labels):
if self.projection == 'vandg' or self._fulldisk:
print 'Warning: Cannot label parallels on %s basemap' % _projnames[self.projection]
labels = [0,0,0,0]
@@ -2068,9 +2080,12 @@
if self.projection in ['sinu','moll','vandg'] and max(labels):
print 'Warning: Cannot label meridians on %s basemap' % _projnames[self.projection]
labels = [0,0,0,0]
- if self.projection in ['ortho','geos'] and max(labels):
+ if self.projection in ['ortho','geos','aeqd'] and max(labels):
if self._fulldisk:
- print 'Warning: Cannot label meridians on full-disk Geostationary or Orthographic basemap'
+ print dedent(
+ """'Warning: Cannot label meridians on full-disk
+ Geostationary, Orthographic or Azimuthal equidistant basemap
+ """)
labels = [0,0,0,0]
# search along edges of map to see if parallels intersect.
# if so, find x,y location of intersection and draw a label there.
@@ -2535,7 +2550,7 @@
# turn off axes frame for non-rectangular projections.
if self.projection in _pseudocyl:
ax.set_frame_on(False)
- if self.projection in ['ortho','geos'] and self._fulldisk:
+ if self.projection in ['ortho','geos','aeqd'] and self._fulldisk:
ax.set_frame_on(False)
# make sure aspect ratio of map preserved.
# plot is re-centered in bounding rectangle.
Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py 2008-12-31 13:20:50 UTC (rev 6718)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/proj.py 2008-12-31 14:37:44 UTC (rev 6719)
@@ -1,6 +1,7 @@
import numpy as np
import pyproj
import math
+from matplotlib.cbook import dedent
__version__ = '1.2.2'
_dg2rad = math.radians(1.)
@@ -70,19 +71,18 @@
self.esq = (self.rmajor**2 - self.rminor**2)/self.rmajor**2
self.llcrnrlon = llcrnrlon
self.llcrnrlat = llcrnrlat
- if self.projection not in ['ortho','geos','cyl'] + _pseudocyl:
- self._proj4 = pyproj.Proj(projparams)
- llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
- elif self.projection == 'cyl':
+ if self.projection == 'cyl':
llcrnrx = llcrnrlon
llcrnry = llcrnrlat
- elif self.projection == 'ortho':
+ elif self.projection in 'ortho':
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
+ self._width = 0.5*(self.rmajor+self.rminor)
+ self._height = 0.5*(self.rmajor+self.rminor)
urcrnrx = -llcrnrx
urcrnry = -llcrnry
else:
@@ -91,6 +91,22 @@
llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
if llcrnrx > 1.e20 or llcrnry > 1.e20:
raise ValueError(_lower_left_out_of_bounds)
+ elif self.projection == 'aeqd' and\
+ (llcrnrlon == -180 and llcrnrlat == -90 and urcrnrlon == 180 and\
+ urcrnrlat == 90):
+ self._fulldisk = True
+ self._proj4 = pyproj.Proj(projparams)
+ if self.ellipsoid:
+ msg = dedent("""
+ full disk (whole world) Azimuthal Equidistant projection can
+ only be drawn for a perfect sphere""")
+ raise ValueError(msg)
+ llcrnrx = -0.5*(self.rmajor+self.rminor)
+ llcrnry = -0.5*(self.rmajor+self.rminor)
+ self._width = -llcrnrx
+ self._height = -llcrnry
+ urcrnrx = -llcrnrx
+ urcrnry = -llcrnry
elif self.projection == 'geos':
self._proj4 = pyproj.Proj(projparams)
# find major and minor axes of ellipse defining map proj region.
@@ -129,6 +145,10 @@
urcrnrx,xtmp = self(projparams['lon_0']+180.,0)
llcrnrx = -urcrnrx
llcrnry = -urcrnry
+ else:
+ self._proj4 = pyproj.Proj(projparams)
+ llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
+ if self.projection == 'aeqd': self._fulldisk=False
# compute x_0, y_0 so ll corner of domain is x=0,y=0.
# note that for 'cyl' x,y == lon,lat
self.projparams['x_0']=-llcrnrx
@@ -144,18 +164,10 @@
if urcrnrislatlon:
self.urcrnrlon = urcrnrlon
self.urcrnrlat = urcrnrlat
- if self.projection not in ['ortho','geos'] + _pseudocyl:
+ if self.projection not in ['ortho','geos','aeqd'] + _pseudocyl:
urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
- elif self.projection == 'ortho':
+ elif self.projection in ['ortho','geos','aeqd']:
if self._fulldisk:
- urcrnrx = 2.*self.rmajor
- urcrnry = 2.*self.rmajor
- else:
- urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
- if urcrnrx > 1.e20 or urcrnry > 1.e20:
- raise ValueError(_upper_right_out_of_bounds)
- elif self.projection == 'geos':
- if self._fulldisk:
urcrnrx = 2.*self._width
urcrnry = 2.*self._height
else:
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
