Revision: 4886 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4886&view=rev Author: jswhit Date: 2008年01月23日 05:20:03 -0800 (2008年1月23日) Log Message: ----------- make a copy of projparams dict, so original is not modified. Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年01月22日 19:43:18 UTC (rev 4885) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年01月23日 13:20:03 UTC (rev 4886) @@ -930,7 +930,7 @@ boundaryxy = _geos.Polygon(b) # compute proj instance for full disk, if necessary. if not self._fulldisk: - projparms = self.projparams + projparms = self.projparams.copy() del projparms['x_0'] del projparms['y_0'] if self.projection == 'ortho': This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 4897 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4897&view=rev Author: jswhit Date: 2008年01月26日 05:06:47 -0800 (2008年1月26日) Log Message: ----------- add llcrnrx,llcrnry,urcrnrx,urcrnry keyword for Basemap.__init__, only relevant for 'ortho' and 'geos' (to enable plotting of geostationary sector images). Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年01月26日 00:11:36 UTC (rev 4896) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年01月26日 13:06:47 UTC (rev 4897) @@ -30,6 +30,7 @@ from matplotlib.collections import LineCollection from matplotlib.patches import Ellipse, Circle, Polygon from matplotlib.lines import Line2D +from matplotlib.transforms import Bbox import pyproj, sys, os, math, dbflib from proj import Proj import numpy as npy @@ -94,8 +95,8 @@ 'spstere' : 'bounding_lat,lon_0,lat_0,no corners or width/height', 'cass' : 'lon_0,lat_0', 'poly' : 'lon_0,lat_0', - 'ortho' : 'lon_0,lat_0', - 'geos' : 'lon_0,lat_0,satellite_height', + 'ortho' : 'lon_0,lat_0,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height', + 'geos' : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height', 'sinu' : 'lon_0,lat_0,no corners or width/height', 'moll' : 'lon_0,lat_0,no corners or width/height', 'robin' : 'lon_0,lat_0,no corners or width/height', @@ -143,7 +144,9 @@ projections except 'ortho' and 'geos', either the lat/lon values of the corners or width and height must be specified by the user. For 'ortho' and 'geos', the lat/lon values of the corners may be specified, - but if they are not, the entire globe is plotted. + or the x/y values of the corners (llcrnrx,llcrnry,urcrnrx,urcrnry) in the + coordinate system of the global projection. If the corners are not + specified, the entire globe is plotted. resolution - resolution of boundary database to use. Can be 'c' (crude), 'l' (low), 'i' (intermediate), 'h' (high), 'f' (full) or None. @@ -323,6 +326,8 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, urcrnrlon=None, urcrnrlat=None, + llcrnrx=None, llcrnry=None, + urcrnrx=None, urcrnry=None, width=None, height=None, projection='cyl', resolution='c', area_thresh=None, rsphere=6370997.0, @@ -605,11 +610,18 @@ self.aspect = (self.urcrnrlat-self.llcrnrlat)/(self.urcrnrlon-self.llcrnrlon) else: self.aspect = (proj.ymax-proj.ymin)/(proj.xmax-proj.xmin) - self.llcrnrx = proj.llcrnrx - self.llcrnry = proj.llcrnry - self.urcrnrx = proj.urcrnrx - self.urcrnry = proj.urcrnry - + if projection in ['geos','ortho'] and \ + None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]: + self.llcrnrx = llcrnrx + self.llcrnry = llcrnry + self.urcrnrx = urcrnrx + self.urcrnry = urcrnry + self._fulldisk = False + else: + self.llcrnrx = proj.llcrnrx + self.llcrnry = proj.llcrnry + self.urcrnrx = proj.urcrnrx + self.urcrnry = proj.urcrnry # set min/max lats for projection domain. if projection in ['mill','cyl','merc']: self.latmin = self.llcrnrlat @@ -1049,34 +1061,23 @@ ax = pylab.gca() elif ax is None and self.ax is not None: ax = self.ax - if self.projection == 'ortho' and self._fulldisk: # circular region. - # define a circle patch, add it to axes instance. - circle = Circle((self.rmajor,self.rmajor),self.rmajor) - ax.add_patch(circle) + if self.projection == 'ortho': + limb = Circle((self.rmajor,self.rmajor),self.rmajor) + elif self.projection == 'geos': + limb = Ellipse((self._width,self._height),2.*self._width,2.*self._height) + if self.projection in ['ortho','geos'] and self._fulldisk: + # elliptical region. + ax.add_patch(limb) if fill_color is None: - circle.set_fill(False) + limb.set_fill(False) else: - circle.set_facecolor(fill_color) - circle.set_zorder(0) - circle.set_edgecolor(color) - circle.set_linewidth(linewidth) - circle.set_clip_on(False) + limb.set_facecolor(fill_color) + limb.set_zorder(0) + limb.set_edgecolor(color) + limb.set_linewidth(linewidth) + limb.set_clip_on(False) if zorder is not None: - circle.set_zorder(zorder) - elif self.projection == 'geos' and self._fulldisk: # elliptical region - # define an Ellipse patch, add it to axes instance. - ellps = Ellipse((self._width,self._height),2.*self._width,2.*self._height) - ax.add_patch(ellps) - if fill_color is None: - ellps.set_fill(False) - else: - ellps.set_facecolor(fill_color) - ellps.set_zorder(0) - ellps.set_edgecolor(color) - ellps.set_linewidth(linewidth) - ellps.set_clip_on(False) - if zorder is not None: - ellps.set_zorder(0) + limb.set_zorder(zorder) elif self.projection in ['moll','robin','sinu']: # elliptical region. nx = 100; ny = 100 # quasi-elliptical region. @@ -1109,15 +1110,37 @@ poly.set_zorder(zorder) else: # all other projections are rectangular. ax.axesPatch.set_linewidth(linewidth) - if fill_color is None: + if self.projection not in ['geos','ortho']: + if fill_color is None: + ax.axesPatch.set_facecolor(ax.get_axis_bgcolor()) + else: + ax.axesPatch.set_facecolor(fill_color) + ax.axesPatch.set_zorder(0) ax.axesPatch.set_facecolor(ax.get_axis_bgcolor()) + ax.axesPatch.set_edgecolor(color) + ax.set_frame_on(True) + if zorder is not None: + ax.axesPatch.set_zorder(zorder) else: - ax.axesPatch.set_facecolor(fill_color) - ax.axesPatch.set_zorder(0) - ax.axesPatch.set_edgecolor(color) - ax.set_frame_on(True) - if zorder is not None: - ax.axesPatch.set_zorder(zorder) + ax.axesPatch.set_facecolor(ax.get_axis_bgcolor()) + ax.axesPatch.set_edgecolor(color) + ax.set_frame_on(True) + if zorder is not None: + ax.axesPatch.set_zorder(zorder) + # for geos or ortho projections, also + # draw and fill map projection limb, clipped + # to rectangular region. + ax.add_patch(limb) + if fill_color is None: + limb.set_fill(False) + else: + limb.set_facecolor(fill_color) + limb.set_zorder(0) + limb.set_edgecolor(color) + limb.set_linewidth(linewidth) + if zorder is not None: + limb.set_zorder(zorder) + limb.set_clip_on(True) # set axes limits to fit map region. self.set_axes_limits(ax=ax) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 4961 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4961&view=rev Author: jswhit Date: 2008年02月13日 13:24:21 -0800 (2008年2月13日) Log Message: ----------- added no_rot parameter (only relevant for omerc projection) Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年02月13日 17:15:11 UTC (rev 4960) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年02月13日 21:24:21 UTC (rev 4961) @@ -79,7 +79,7 @@ projection_params = {'cyl' : 'corners only (no width/height)', 'merc' : 'corners plus lat_ts (no width/height)', 'tmerc' : 'lon_0,lat_0', - 'omerc' : 'lat_0,lat_1,lat_2,lon_1,lon_2', + 'omerc' : 'lon_0,lat_0,lat_1,lat_2,lon_1,lon_2,no_rot', 'mill' : 'corners only (no width/height)', 'lcc' : 'lon_0,lat_0,lat_1,lat_2', 'laea' : 'lon_0,lat_0', @@ -215,6 +215,9 @@ for oblique mercator. lon_2 - longitude of one of the two points on the projection centerline for oblique mercator. + no_rot - only used by oblique mercator. If set to True, the map projection + coordinates will not be rotated to true North. Default is False (projection + coordinates are automatically rotated). lat_0 - central latitude (y-axis origin) - used by all projections, Must be equator for mercator projection. lon_0 - central meridian (x-axis origin) - used by all projections, @@ -336,6 +339,7 @@ lat_1=None, lat_2=None, lat_0=None, lon_0=None, lon_1=None, lon_2=None, + no_rot=False, suppress_ticks=True, satellite_height=35786000, boundinglat=None, @@ -527,13 +531,15 @@ self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat elif projection == 'omerc': - if lat_1 is None or lon_1 is None or lat_2 is None or lon_2 is None or lat_0 is None: - raise ValueError, 'must specify lat_0 and lat_1,lon_1 and lat_2,lon_2 for Oblique Mercator basemap' + if lat_1 is None or lon_1 is None or lat_2 is None or lon_2 is None: + raise ValueError, 'must specify lat_1,lon_1 and lat_2,lon_2 for Oblique Mercator basemap' projparams['lat_1'] = lat_1 projparams['lon_1'] = lon_1 projparams['lat_2'] = lat_2 projparams['lon_2'] = lon_2 projparams['lat_0'] = lat_0 + if no_rot: + projparams['no_rot']='' #if not using_corners: # raise ValueError, 'cannot specify map region with width and height keywords for this projection, please specify lat/lon values of corners' if not using_corners: This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 4967 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4967&view=rev Author: jswhit Date: 2008年02月14日 15:18:09 -0800 (2008年2月14日) Log Message: ----------- fix fill_color in drawmapboundary Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年02月14日 22:22:13 UTC (rev 4966) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年02月14日 23:18:09 UTC (rev 4967) @@ -1133,7 +1133,6 @@ else: ax.axesPatch.set_facecolor(fill_color) ax.axesPatch.set_zorder(0) - ax.axesPatch.set_facecolor(ax.get_axis_bgcolor()) ax.axesPatch.set_edgecolor(color) ax.set_frame_on(True) if zorder is not None: This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 4971 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4971&view=rev Author: jswhit Date: 2008年02月15日 05:45:41 -0800 (2008年2月15日) Log Message: ----------- fix typo in warpimage method Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年02月15日 13:43:19 UTC (rev 4970) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年02月15日 13:45:41 UTC (rev 4971) @@ -2686,7 +2686,7 @@ try: from PIL import Image except ImportError: - raise ImportError('bluemarble method requires PIL (http://www.pythonware.com/products/pil/)') + raise ImportError('warpimage method requires PIL (http://www.pythonware.com/products/pil)') from matplotlib.image import pil_to_array if not kwargs.has_key('ax') and self.ax is None: try: This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 4994 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4994&view=rev Author: jswhit Date: 2008年03月04日 14:58:49 -0800 (2008年3月04日) Log Message: ----------- update docstring for drawmapscale Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年03月04日 22:11:10 UTC (rev 4993) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年03月04日 22:58:49 UTC (rev 4994) @@ -2796,6 +2796,8 @@ fillcolor1 and fillcolor2 are only relevant for the 'fancy' scale bar. They are the colors of the alternating filled regions (default white and black). + + extra keyword 'ax' can be used to override the default axis instance. """ # get current axes instance (if none specified). if ax is None and self.ax is None: This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 5029 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5029&view=rev Author: jswhit Date: 2008年04月04日 08:02:48 -0700 (2008年4月04日) Log Message: ----------- workaround for problem reported by Rob Hetland (meridians don't reach edge of very tiny maps) Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年04月03日 15:24:20 UTC (rev 5028) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年04月04日 15:02:48 UTC (rev 5029) @@ -1587,10 +1587,13 @@ lats = circ*npy.ones(len(lons),npy.float32) x,y = self(lons,lats) # remove points outside domain. - testx = npy.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta) + # leave a little slop around edges (3*xdelta) + # don't really know why, but this appears to be needed to + # or lines sometimes don't reach edge of plot. + testx = npy.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta) x = npy.compress(testx, x) y = npy.compress(testx, y) - testy = npy.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta) + testy = npy.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta) x = npy.compress(testy, x) y = npy.compress(testy, y) lines = [] @@ -1812,10 +1815,13 @@ lons = merid*npy.ones(len(lats),npy.float32) x,y = self(lons,lats) # remove points outside domain. - testx = npy.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta) + # leave a little slop around edges (3*xdelta) + # don't really know why, but this appears to be needed to + # or lines sometimes don't reach edge of plot. + testx = npy.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta) x = npy.compress(testx, x) y = npy.compress(testx, y) - testy = npy.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta) + testy = npy.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta) x = npy.compress(testy, x) y = npy.compress(testy, y) lines = [] This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 5199 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5199&view=rev Author: jswhit Date: 2008年05月20日 08:01:46 -0700 (2008年5月20日) Log Message: ----------- convert to np/plt (from numpy, pylab) Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月20日 12:24:54 UTC (rev 5198) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月20日 15:01:46 UTC (rev 5199) @@ -33,8 +33,8 @@ from matplotlib.transforms import Bbox import pyproj, sys, os, math, dbflib from proj import Proj -import numpy as npy -from numpy import linspace, squeeze, ma +import numpy as np +from numpy import ma from shapelib import ShapeFile import _geoslib, pupynere, netcdftime @@ -504,7 +504,7 @@ self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat # FIXME: won't work for points exactly on equator?? - if npy.abs(lat_0) < 1.e-2: lat_0 = 1.e-2 + if np.abs(lat_0) < 1.e-2: lat_0 = 1.e-2 projparams['lat_0'] = lat_0 elif projection == 'geos': if lon_0 is None: @@ -652,7 +652,7 @@ self.latmin = lats.min() self.latmax = lats.max() - # if ax == None, pylab.gca may be used. + # if ax == None, pyplot.gca may be used. self.ax = ax self.lsmask = None @@ -685,13 +685,13 @@ for seg in self.coastsegs: x, y = zip(*seg) self.coastpolygons.append((x,y)) - x = npy.array(x,npy.float64); y = npy.array(y,npy.float64) + x = np.array(x,np.float64); y = np.array(y,np.float64) xd = (x[1:]-x[0:-1])**2 yd = (y[1:]-y[0:-1])**2 - dist = npy.sqrt(xd+yd) + dist = np.sqrt(xd+yd) split = dist > 5000000. - if npy.sum(split) and self.projection not in ['merc','cyl','mill']: - ind = (npy.compress(split,squeeze(split*npy.indices(xd.shape)))+1).tolist() + if np.sum(split) and self.projection not in ['merc','cyl','mill']: + ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist() iprev = 0 ind.append(len(xd)) for i in ind: @@ -779,9 +779,9 @@ re = self.projparams['R'] # center of stereographic projection restricted to be # nearest one of 6 points on the sphere (every 90 deg lat/lon). - lon0 = 90.*(npy.around(lon_0/90.)) - lat0 = 90.*(npy.around(lat_0/90.)) - if npy.abs(int(lat0)) == 90: lon0=0. + lon0 = 90.*(np.around(lon_0/90.)) + lat0 = 90.*(np.around(lat_0/90.)) + if np.abs(int(lat0)) == 90: lon0=0. maptran = pyproj.Proj(proj='stere',lon_0=lon0,lat_0=lat0,R=re) # boundary polygon for orthographic projection # in stereographic coorindates. @@ -806,7 +806,7 @@ # numpy array (first column is lons, second is lats). polystring = bdatfile.read(bytecount) # binary data is little endian. - b = npy.array(npy.fromstring(polystring,dtype='<f4'),'f8') + b = np.array(np.fromstring(polystring,dtype='<f4'),'f8') b.shape = (npts,2) b2 = b.copy() # if map boundary polygon is a valid one in lat/lon @@ -833,7 +833,7 @@ lats.insert(0,-90.) lons.append(lonend) lats.append(-90.) - b = npy.empty((len(lons),2),npy.float64) + b = np.empty((len(lons),2),np.float64) b[:,0] = lons; b[:,1] = lats poly = _geoslib.Polygon(b) antart = True @@ -887,16 +887,16 @@ b[:,0], b[:,1] = maptran(b[:,0], b[:,1]) else: b[:,0], b[:,1] = self(b[:,0], b[:,1]) - goodmask = npy.logical_and(b[:,0]<1.e20,b[:,1]<1.e20) + goodmask = np.logical_and(b[:,0]<1.e20,b[:,1]<1.e20) # if less than two points are valid in # map proj coords, skip this geometry. - if npy.sum(goodmask) <= 1: continue + if np.sum(goodmask) <= 1: continue if name != 'gshhs': # if not a polygon, # just remove parts of geometry that are undefined # in this map projection. - bx = npy.compress(goodmask, b[:,0]) - by = npy.compress(goodmask, b[:,1]) + bx = np.compress(goodmask, b[:,0]) + by = np.compress(goodmask, b[:,1]) # for orthographic projection, all points # outside map projection region are eliminated # with the above step, so we're done. @@ -947,16 +947,16 @@ maptran = self if self.projection in ['ortho','geos']: # circular region. - thetas = linspace(0.,2.*npy.pi,2*nx*ny)[:-1] + 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 - x = rmajor*npy.cos(thetas) + rmajor - y = rminor*npy.sin(thetas) + rminor - b = npy.empty((len(x),2),npy.float64) + 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) # compute proj instance for full disk, if necessary. @@ -981,41 +981,41 @@ # quasi-elliptical region. lon_0 = self.projparams['lon_0'] # left side - lats1 = linspace(-89.9,89.9,ny).tolist() + lats1 = np.linspace(-89.9,89.9,ny).tolist() lons1 = len(lats1)*[lon_0-179.9] # top. - lons2 = linspace(lon_0-179.9,lon_0+179.9,nx).tolist() + lons2 = np.linspace(lon_0-179.9,lon_0+179.9,nx).tolist() lats2 = len(lons2)*[89.9] # right side - lats3 = linspace(89.9,-89.9,ny).tolist() + lats3 = np.linspace(89.9,-89.9,ny).tolist() lons3 = len(lats3)*[lon_0+179.9] # bottom. - lons4 = linspace(lon_0+179.9,lon_0-179.9,nx).tolist() + lons4 = np.linspace(lon_0+179.9,lon_0-179.9,nx).tolist() lats4 = len(lons4)*[-89.9] - lons = npy.array(lons1+lons2+lons3+lons4,npy.float64) - lats = npy.array(lats1+lats2+lats3+lats4,npy.float64) + lons = np.array(lons1+lons2+lons3+lons4,np.float64) + lats = np.array(lats1+lats2+lats3+lats4,np.float64) x, y = maptran(lons,lats) - b = npy.empty((len(x),2),npy.float64) + b = np.empty((len(x),2),np.float64) b[:,0]=x; b[:,1]=y boundaryxy = _geoslib.Polygon(b) else: # all other projections are rectangular. # left side (x = xmin, ymin <= y <= ymax) - yy = linspace(self.ymin, self.ymax, ny)[:-1] + yy = np.linspace(self.ymin, self.ymax, ny)[:-1] x = len(yy)*[self.xmin]; y = yy.tolist() # top (y = ymax, xmin <= x <= xmax) - xx = npy.linspace(self.xmin, self.xmax, nx)[:-1] + xx = np.np.linspace(self.xmin, self.xmax, nx)[:-1] x = x + xx.tolist() y = y + len(xx)*[self.ymax] # right side (x = xmax, ymin <= y <= ymax) - yy = linspace(self.ymax, self.ymin, ny)[:-1] + yy = np.linspace(self.ymax, self.ymin, ny)[:-1] x = x + len(yy)*[self.xmax]; y = y + yy.tolist() # bottom (y = ymin, xmin <= x <= xmax) - xx = linspace(self.xmax, self.xmin, nx)[:-1] + xx = np.linspace(self.xmax, self.xmin, nx)[:-1] x = x + xx.tolist() y = y + len(xx)*[self.ymin] - x = npy.array(x,npy.float64) - y = npy.array(y,npy.float64) - b = npy.empty((4,2),npy.float64) + x = np.array(x,np.float64) + y = np.array(y,np.float64) + b = np.empty((4,2),np.float64) b[:,0]=[self.xmin,self.xmin,self.xmax,self.xmax] b[:,1]=[self.ymin,self.ymax,self.ymax,self.ymin] boundaryxy = _geoslib.Polygon(b) @@ -1032,7 +1032,7 @@ lons = [self.llcrnrlon, self.llcrnrlon, self.urcrnrlon, self.urcrnrlon] lats = [llcrnrlat, urcrnrlat, urcrnrlat, llcrnrlat] x, y = self(lons, lats) - b = npy.empty((len(x),2),npy.float64) + b = np.empty((len(x),2),np.float64) b[:,0]=x; b[:,1]=y boundaryxy = _geoslib.Polygon(b) else: @@ -1042,7 +1042,7 @@ n = 1 lonprev = lons[0] for lon,lat in zip(lons[1:],lats[1:]): - if npy.abs(lon-lonprev) > 90.: + if np.abs(lon-lonprev) > 90.: if lonprev < 0: lon = lon - 360. else: @@ -1050,7 +1050,7 @@ lons[n] = lon lonprev = lon n = n + 1 - b = npy.empty((len(lons),2),npy.float64) + b = np.empty((len(lons),2),np.float64) b[:,0]=lons; b[:,1]=lats boundaryll = _geoslib.Polygon(b) return boundaryll, boundaryxy @@ -1075,10 +1075,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax limb = None @@ -1104,19 +1104,19 @@ # quasi-elliptical region. lon_0 = self.projparams['lon_0'] # left side - lats1 = linspace(-89.9,89.9,ny).tolist() + lats1 = np.linspace(-89.9,89.9,ny).tolist() lons1 = len(lats1)*[lon_0-179.9] # top. - lons2 = linspace(lon_0-179.9,lon_0+179.9,nx).tolist() + lons2 = np.linspace(lon_0-179.9,lon_0+179.9,nx).tolist() lats2 = len(lons2)*[89.9] # right side - lats3 = linspace(89.9,-89.9,ny).tolist() + lats3 = np.linspace(89.9,-89.9,ny).tolist() lons3 = len(lats3)*[lon_0+179.9] # bottom. - lons4 = linspace(lon_0+179.9,lon_0-179.9,nx).tolist() + lons4 = np.linspace(lon_0+179.9,lon_0-179.9,nx).tolist() lats4 = len(lons4)*[-89.9] - lons = npy.array(lons1+lons2+lons3+lons4,npy.float64) - lats = npy.array(lats1+lats2+lats3+lats4,npy.float64) + lons = np.array(lons1+lons2+lons3+lons4,np.float64) + lats = np.array(lats1+lats2+lats3+lats4,np.float64) x, y = self(lons,lats) xy = zip(x,y) limb = Polygon(xy,edgecolor=color,linewidth=linewidth) @@ -1185,32 +1185,32 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax # get axis background color. axisbgc = ax.get_axis_bgcolor() np = 0 for x,y in self.coastpolygons: - xa = npy.array(x,npy.float32) - ya = npy.array(y,npy.float32) + xa = np.array(x,np.float32) + ya = np.array(y,np.float32) # check to see if all four corners of domain in polygon (if so, # don't draw since it will just fill in the whole map). delx = 10; dely = 10 if self.projection in ['cyl']: delx = 0.1 dely = 0.1 - test1 = npy.fabs(xa-self.urcrnrx) < delx - test2 = npy.fabs(xa-self.llcrnrx) < delx - test3 = npy.fabs(ya-self.urcrnry) < dely - test4 = npy.fabs(ya-self.llcrnry) < dely - hasp1 = npy.sum(test1*test3) - hasp2 = npy.sum(test2*test3) - hasp4 = npy.sum(test2*test4) - hasp3 = npy.sum(test1*test4) + test1 = np.fabs(xa-self.urcrnrx) < delx + test2 = np.fabs(xa-self.llcrnrx) < delx + test3 = np.fabs(ya-self.urcrnry) < dely + test4 = np.fabs(ya-self.llcrnry) < dely + hasp1 = np.sum(test1*test3) + hasp2 = np.sum(test2*test3) + hasp4 = np.sum(test2*test4) + hasp3 = np.sum(test1*test4) if not hasp1 or not hasp2 or not hasp3 or not hasp4: xy = zip(xa.tolist(),ya.tolist()) if self.coastpolygontypes[np] not in [2,4]: @@ -1245,10 +1245,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax coastlines = LineCollection(self.coastsegs,antialiaseds=(antialiased,)) @@ -1283,10 +1283,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax countries = LineCollection(self.cntrysegs,antialiaseds=(antialiased,)) @@ -1321,10 +1321,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax states = LineCollection(self.statesegs,antialiaseds=(antialiased,)) @@ -1359,10 +1359,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax rivers = LineCollection(self.riversegs,antialiaseds=(antialiased,)) @@ -1484,10 +1484,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax # make LineCollections for each polygon. @@ -1537,7 +1537,7 @@ ax - axes instance (overrides default axes instance) additional keyword arguments control text properties for labels (see - pylab.text documentation) + plt.text documentation) returns a dictionary whose keys are the parallels, and whose values are tuples containing lists of the Line2D and Text instances @@ -1546,10 +1546,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax # don't draw meridians past latmax, always draw parallel at latmax. @@ -1565,13 +1565,13 @@ xoffset = (self.urcrnrx-self.llcrnrx)/100. if self.projection in ['merc','cyl','mill','moll','robin','sinu']: - lons = npy.arange(self.llcrnrlon,self.urcrnrlon+0.01,0.01) + lons = np.arange(self.llcrnrlon,self.urcrnrlon+0.01,0.01) elif self.projection in ['tmerc']: lon_0 = self.projparams['lon_0'] # tmerc only defined within +/- 90 degrees of lon_0 - lons = npy.arange(lon_0-90,lon_0+90.01,0.01) + lons = np.arange(lon_0-90,lon_0+90.01,0.01) else: - lons = npy.arange(0,360.01,0.01) + lons = np.arange(0,360.01,0.01) # make sure latmax degree parallel is drawn if projection not merc or cyl or miller try: circlesl = circles.tolist() @@ -1586,28 +1586,28 @@ ydelta = 0.01*(self.ymax-self.ymin) linecolls = {} for circ in circlesl: - lats = circ*npy.ones(len(lons),npy.float32) + lats = circ*np.ones(len(lons),np.float32) x,y = self(lons,lats) # remove points outside domain. # leave a little slop around edges (3*xdelta) # don't really know why, but this appears to be needed to # or lines sometimes don't reach edge of plot. - testx = npy.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta) - x = npy.compress(testx, x) - y = npy.compress(testx, y) - testy = npy.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta) - x = npy.compress(testy, x) - y = npy.compress(testy, y) + testx = np.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta) + x = np.compress(testx, x) + y = np.compress(testx, y) + testy = np.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta) + x = np.compress(testy, x) + y = np.compress(testy, y) lines = [] if len(x) > 1 and len(y) > 1: # split into separate line segments if necessary. # (not necessary for mercator or cylindrical or miller). xd = (x[1:]-x[0:-1])**2 yd = (y[1:]-y[0:-1])**2 - dist = npy.sqrt(xd+yd) + dist = np.sqrt(xd+yd) split = dist > 500000. - if npy.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']: - ind = (npy.compress(split,squeeze(split*npy.indices(xd.shape)))+1).tolist() + if np.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']: + ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist() xl = [] yl = [] iprev = 0 @@ -1650,15 +1650,15 @@ if self.projection in ['cyl','merc','mill','moll','robin','sinu'] and side in ['t','b']: continue if side in ['l','r']: nmax = int((self.ymax-self.ymin)/dy+1) - yy = linspace(self.llcrnry,self.urcrnry,nmax) + yy = np.linspace(self.llcrnry,self.urcrnry,nmax) # mollweide inverse transform undefined at South Pole if self.projection == 'moll' and yy[0] < 1.e-4: yy[0] = 1.e-4 if side == 'l': - lons,lats = self(self.llcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True) + lons,lats = self(self.llcrnrx*np.ones(yy.shape,np.float32),yy,inverse=True) lons = lons.tolist(); lats = lats.tolist() else: - lons,lats = self(self.urcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True) + lons,lats = self(self.urcrnrx*np.ones(yy.shape,np.float32),yy,inverse=True) lons = lons.tolist(); lats = lats.tolist() if max(lons) > 1.e20 or max(lats) > 1.e20: raise ValueError,'inverse transformation undefined - please adjust the map projection region' @@ -1666,12 +1666,12 @@ lons = [(lon+360) % 360 for lon in lons] else: nmax = int((self.xmax-self.xmin)/dx+1) - xx = linspace(self.llcrnrx,self.urcrnrx,nmax) + xx = np.linspace(self.llcrnrx,self.urcrnrx,nmax) if side == 'b': - lons,lats = self(xx,self.llcrnry*npy.ones(xx.shape,npy.float32),inverse=True) + lons,lats = self(xx,self.llcrnry*np.ones(xx.shape,np.float32),inverse=True) lons = lons.tolist(); lats = lats.tolist() else: - lons,lats = self(xx,self.urcrnry*npy.ones(xx.shape,npy.float32),inverse=True) + lons,lats = self(xx,self.urcrnry*np.ones(xx.shape,np.float32),inverse=True) lons = lons.tolist(); lats = lats.tolist() if max(lons) > 1.e20 or max(lats) > 1.e20: raise ValueError,'inverse transformation undefined - please adjust the map projection region' @@ -1697,7 +1697,7 @@ latlabstr = u'-%s\N{DEGREE SIGN}'%fmt else: latlabstr = u'%s\N{DEGREE SIGN}S'%fmt - latlab = latlabstr%npy.fabs(lat) + latlab = latlabstr%np.fabs(lat) elif lat>0: if rcParams['text.usetex']: if labelstyle=='+/-': @@ -1779,7 +1779,7 @@ ax - axes instance (overrides default axes instance) additional keyword arguments control text properties for labels (see - pylab.text documentation) + plt.text documentation) returns a dictionary whose keys are the meridians, and whose values are tuples containing lists of the Line2D and Text instances @@ -1788,10 +1788,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax # don't draw meridians past latmax, always draw parallel at latmax. @@ -1807,35 +1807,35 @@ xoffset = (self.urcrnrx-self.llcrnrx)/100. if self.projection not in ['merc','cyl','mill','moll','robin','sinu']: - lats = npy.arange(-latmax,latmax+0.01,0.01) + lats = np.arange(-latmax,latmax+0.01,0.01) else: - lats = npy.arange(-90,90.01,0.01) + lats = np.arange(-90,90.01,0.01) xdelta = 0.01*(self.xmax-self.xmin) ydelta = 0.01*(self.ymax-self.ymin) linecolls = {} for merid in meridians: - lons = merid*npy.ones(len(lats),npy.float32) + lons = merid*np.ones(len(lats),np.float32) x,y = self(lons,lats) # remove points outside domain. # leave a little slop around edges (3*xdelta) # don't really know why, but this appears to be needed to # or lines sometimes don't reach edge of plot. - testx = npy.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta) - x = npy.compress(testx, x) - y = npy.compress(testx, y) - testy = npy.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta) - x = npy.compress(testy, x) - y = npy.compress(testy, y) + testx = np.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta) + x = np.compress(testx, x) + y = np.compress(testx, y) + testy = np.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta) + x = np.compress(testy, x) + y = np.compress(testy, y) lines = [] if len(x) > 1 and len(y) > 1: # split into separate line segments if necessary. # (not necessary for mercator or cylindrical or miller). xd = (x[1:]-x[0:-1])**2 yd = (y[1:]-y[0:-1])**2 - dist = npy.sqrt(xd+yd) + dist = np.sqrt(xd+yd) split = dist > 500000. - if npy.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']: - ind = (npy.compress(split,squeeze(split*npy.indices(xd.shape)))+1).tolist() + if np.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']: + ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist() xl = [] yl = [] iprev = 0 @@ -1884,12 +1884,12 @@ if self.projection in ['cyl','merc','mill','sinu','robin','moll'] and side in ['l','r']: continue if side in ['l','r']: nmax = int((self.ymax-self.ymin)/dy+1) - yy = linspace(self.llcrnry,self.urcrnry,nmax) + yy = np.linspace(self.llcrnry,self.urcrnry,nmax) if side == 'l': - lons,lats = self(self.llcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True) + lons,lats = self(self.llcrnrx*np.ones(yy.shape,np.float32),yy,inverse=True) lons = lons.tolist(); lats = lats.tolist() else: - lons,lats = self(self.urcrnrx*npy.ones(yy.shape,npy.float32),yy,inverse=True) + lons,lats = self(self.urcrnrx*np.ones(yy.shape,np.float32),yy,inverse=True) lons = lons.tolist(); lats = lats.tolist() if max(lons) > 1.e20 or max(lats) > 1.e20: raise ValueError,'inverse transformation undefined - please adjust the map projection region' @@ -1897,12 +1897,12 @@ lons = [(lon+360) % 360 for lon in lons] else: nmax = int((self.xmax-self.xmin)/dx+1) - xx = linspace(self.llcrnrx,self.urcrnrx,nmax) + xx = np.linspace(self.llcrnrx,self.urcrnrx,nmax) if side == 'b': - lons,lats = self(xx,self.llcrnry*npy.ones(xx.shape,npy.float32),inverse=True) + lons,lats = self(xx,self.llcrnry*np.ones(xx.shape,np.float32),inverse=True) lons = lons.tolist(); lats = lats.tolist() else: - lons,lats = self(xx,self.urcrnry*npy.ones(xx.shape,npy.float32),inverse=True) + lons,lats = self(xx,self.urcrnry*np.ones(xx.shape,np.float32),inverse=True) lons = lons.tolist(); lats = lats.tolist() if max(lons) > 1.e20 or max(lats) > 1.e20: raise ValueError,'inverse transformation undefined - please adjust the map projection region' @@ -1930,7 +1930,7 @@ lonlabstr = u'-%s\N{DEGREE SIGN}'%fmt else: lonlabstr = u'%s\N{DEGREE SIGN}W'%fmt - lonlab = lonlabstr%npy.fabs(lon2-360) + lonlab = lonlabstr%np.fabs(lon2-360) elif lon2<180 and lon2 != 0: if rcParams['text.usetex']: if labelstyle=='+/-': @@ -2002,7 +2002,7 @@ del_s - points on great circle computed every delta kilometers (default 100). Other keyword arguments (**kwargs) control plotting of great circle line, - see pylab.plot documentation for details. + see plt.plot documentation for details. Note: cannot handle situations in which the great circle intersects the edge of the map projection domain, and then re-enters the domain. @@ -2055,8 +2055,8 @@ raise ValueError, 'lons and lats must be increasing!' # check that lons in -180,180 for non-cylindrical projections. if self.projection not in ['cyl','merc','mill']: - lonsa = npy.array(lons) - count = npy.sum(lonsa < -180.00001) + npy.sum(lonsa > 180.00001) + lonsa = np.array(lons) + count = np.sum(lonsa < -180.00001) + np.sum(lonsa > 180.00001) if count > 1: raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)' # allow for wraparound point to be outside. @@ -2109,8 +2109,8 @@ raise ValueError, 'lons and lats must be increasing!' # check that lons in -180,180 for non-cylindrical projections. if self.projection not in ['cyl','merc','mill']: - lonsa = npy.array(lons) - count = npy.sum(lonsa < -180.00001) + npy.sum(lonsa > 180.00001) + lonsa = np.array(lons) + count = np.sum(lonsa < -180.00001) + np.sum(lonsa > 180.00001) if count > 1: raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)' # allow for wraparound point to be outside. @@ -2127,10 +2127,10 @@ vin = vin.filled(1) masked = True # override kwarg with reality uvc = uin + 1j*vin - uvmag = npy.abs(uvc) + uvmag = np.abs(uvc) delta = 0.1 # increment in longitude dlon = delta*uin/uvmag - dlat = delta*(vin/uvmag)*npy.cos(latsout*npy.pi/180.0) + dlat = delta*(vin/uvmag)*np.cos(latsout*np.pi/180.0) farnorth = latsout+dlat >= 90.0 somenorth = farnorth.any() if somenorth: @@ -2139,10 +2139,10 @@ lon1 = lonsout + dlon lat1 = latsout + dlat xn, yn = self(lon1, lat1) - vecangle = npy.arctan2(yn-y, xn-x) + vecangle = np.arctan2(yn-y, xn-x) if somenorth: - vecangle[farnorth] += npy.pi - uvcout = uvmag * npy.exp(1j*vecangle) + vecangle[farnorth] += np.pi + uvcout = uvmag * np.exp(1j*vecangle) uout = uvcout.real vout = uvcout.imag if masked: @@ -2184,10 +2184,10 @@ else: masked = False uvc = uin + 1j*vin - uvmag = npy.abs(uvc) + uvmag = np.abs(uvc) delta = 0.1 # increment in longitude dlon = delta*uin/uvmag - dlat = delta*(vin/uvmag)*npy.cos(lats*npy.pi/180.0) + dlat = delta*(vin/uvmag)*np.cos(lats*np.pi/180.0) farnorth = lats+dlat >= 90.0 somenorth = farnorth.any() if somenorth: @@ -2196,10 +2196,10 @@ lon1 = lons + dlon lat1 = lats + dlat xn, yn = self(lon1, lat1) - vecangle = npy.arctan2(yn-y, xn-x) + vecangle = np.arctan2(yn-y, xn-x) if somenorth: - vecangle[farnorth] += npy.pi - uvcout = uvmag * npy.exp(1j*vecangle) + vecangle[farnorth] += np.pi + uvcout = uvmag * np.exp(1j*vecangle) uout = uvcout.real vout = uvcout.imag if masked: @@ -2218,10 +2218,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax # update data limits for map domain. @@ -2250,15 +2250,15 @@ def scatter(self, *args, **kwargs): """ - Plot points with markers on the map (see pylab.scatter documentation). + Plot points with markers on the map (see plt.scatter documentation). extra keyword 'ax' can be used to override the default axes instance. """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2271,7 +2271,7 @@ try: ret = ax.scatter(*args, **kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: @@ -2284,15 +2284,15 @@ def plot(self, *args, **kwargs): """ - Draw lines and/or markers on the map (see pylab.plot documentation). + Draw lines and/or markers on the map (see plt.plot documentation). extra keyword 'ax' can be used to override the default axis instance. """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2305,7 +2305,7 @@ try: ret = ax.plot(*args, **kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: @@ -2318,17 +2318,17 @@ def imshow(self, *args, **kwargs): """ - Display an image over the map (see pylab.imshow documentation). + Display an image over the map (see plt.imshow documentation). extent and origin keywords set automatically so image will be drawn over map region. extra keyword 'ax' can be used to override the default axis instance. """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2345,16 +2345,16 @@ try: ret = ax.imshow(*args, **kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: ax.hold(b) raise ax.hold(b) - # reset current active image (only if pylab is imported). + # reset current active image (only if pyplot is imported). try: - pylab.gci._current = ret + plt.gci._current = ret except: pass # set axes limits to fit map region. @@ -2364,7 +2364,7 @@ def pcolor(self,x,y,data,**kwargs): """ Make a pseudo-color plot over the map. - see pylab.pcolor documentation for definition of **kwargs + see plt.pcolor documentation for definition of **kwargs If x or y are outside projection limb (i.e. they have values > 1.e20) they will be convert to masked arrays with those values masked. As a result, those values will not be plotted. @@ -2372,18 +2372,18 @@ """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: ax = kwargs.pop('ax') # make x,y masked arrays # (masked where data is outside of projection limb) - x = ma.masked_values(npy.where(x > 1.e20,1.e20,x), 1.e20) - y = ma.masked_values(npy.where(y > 1.e20,1.e20,y), 1.e20) + x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20) + y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20) # allow callers to override the hold state by passing hold=True|False b = ax.ishold() h = kwargs.pop('hold',None) @@ -2392,16 +2392,16 @@ try: ret = ax.pcolor(x,y,data,**kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: ax.hold(b) raise ax.hold(b) - # reset current active image (only if pylab is imported). + # reset current active image (only if pyplot is imported). try: - pylab.gci._current = ret + plt.gci._current = ret except: pass # set axes limits to fit map region. @@ -2411,15 +2411,15 @@ def pcolormesh(self,x,y,data,**kwargs): """ Make a pseudo-color plot over the map. - see pylab.pcolormesh documentation for definition of **kwargs + see plt.pcolormesh documentation for definition of **kwargs extra keyword 'ax' can be used to override the default axis instance. """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2432,16 +2432,16 @@ try: ret = ax.pcolormesh(x,y,data,**kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: ax.hold(b) raise ax.hold(b) - # reset current active image (only if pylab is imported). + # reset current active image (only if pyplot is imported). try: - pylab.gci._current = ret + plt.gci._current = ret except: pass # set axes limits to fit map region. @@ -2450,15 +2450,15 @@ def contour(self,x,y,data,*args,**kwargs): """ - Make a contour plot over the map (see pylab.contour documentation). + Make a contour plot over the map (see plt.contour documentation). extra keyword 'ax' can be used to override the default axis instance. """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2482,10 +2482,10 @@ function to adjust the data to be consistent with the map projection region (see examples/contour_demo.py).""") # mask for points outside projection limb. - xymask = npy.logical_or(npy.greater(x,1.e20),npy.greater(y,1.e20)) + xymask = np.logical_or(np.greater(x,1.e20),np.greater(y,1.e20)) data = ma.asarray(data) # combine with data mask. - mask = npy.logical_or(ma.getmaskarray(data),xymask) + mask = np.logical_or(ma.getmaskarray(data),xymask) data = ma.masked_array(data,mask=mask) # allow callers to override the hold state by passing hold=True|False b = ax.ishold() @@ -2495,7 +2495,7 @@ try: CS = ax.contour(x,y,data,*args,**kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: @@ -2504,29 +2504,29 @@ ax.hold(b) # set axes limits to fit map region. self.set_axes_limits(ax=ax) - # reset current active image (only if pylab is imported). + # reset current active image (only if pyplot is imported). try: try: # new contour. - if CS._A is not None: pylab.gci._current = CS + if CS._A is not None: plt.gci._current = CS except: # old contour. - if CS[1].mappable is not None: pylab.gci._current = CS[1].mappable + if CS[1].mappable is not None: plt.gci._current = CS[1].mappable except: pass return CS def contourf(self,x,y,data,*args,**kwargs): """ - Make a filled contour plot over the map (see pylab.contourf documentation). + Make a filled contour plot over the map (see plt.contourf documentation). If x or y are outside projection limb (i.e. they have values > 1.e20), the corresponing data elements will be masked. Extra keyword 'ax' can be used to override the default axis instance. """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2550,10 +2550,10 @@ function to adjust the data to be consistent with the map projection region (see examples/contour_demo.py).""") # mask for points outside projection limb. - xymask = npy.logical_or(npy.greater(x,1.e20),npy.greater(y,1.e20)) + xymask = np.logical_or(np.greater(x,1.e20),np.greater(y,1.e20)) data = ma.asarray(data) # combine with data mask. - mask = npy.logical_or(ma.getmaskarray(data),xymask) + mask = np.logical_or(ma.getmaskarray(data),xymask) data = ma.masked_array(data,mask=mask) # allow callers to override the hold state by passing hold=True|False b = ax.ishold() @@ -2563,7 +2563,7 @@ try: CS = ax.contourf(x,y,data,*args,**kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: @@ -2572,12 +2572,12 @@ ax.hold(b) # set axes limits to fit map region. self.set_axes_limits(ax=ax) - # reset current active image (only if pylab is imported). + # reset current active image (only if pyplot is imported). try: try: # new contour. - if CS._A is not None: pylab.gci._current = CS + if CS._A is not None: plt.gci._current = CS except: # old contour. - if CS[1].mappable is not None: pylab.gci._current = CS[1].mappable + if CS[1].mappable is not None: plt.gci._current = CS[1].mappable except: pass return CS @@ -2587,15 +2587,15 @@ Make a vector plot (u, v) with arrows on the map. Extra arguments (*args and **kwargs) passed to quiver Axes method (see - pylab.quiver documentation for details). + plt.quiver documentation for details). extra keyword 'ax' can be used to override the default axis instance. """ if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2608,7 +2608,7 @@ try: ret = ax.quiver(x,y,u,v,*args,**kwargs) try: - pylab.draw_if_interactive() + plt.draw_if_interactive() except: pass except: @@ -2656,13 +2656,13 @@ extra keyword 'ax' can be used to override the default axis instance. """ # look for axes instance (as keyword, an instance variable - # or from pylab.gca(). + # or from plt.gca(). if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2677,9 +2677,9 @@ lsmaskf = open(os.path.join(basemap_datadir,'5minmask.bin'),'rb') nlons = 4320; nlats = nlons/2 delta = 360./float(nlons) - lsmask_lons = npy.arange(-180+0.5*delta,180.,delta) - lsmask_lats = npy.arange(-90.+0.5*delta,90.,delta) - lsmask = npy.reshape(npy.fromstring(lsmaskf.read(),npy.uint8),(nlats,nlons)) + lsmask_lons = np.arange(-180+0.5*delta,180.,delta) + lsmask_lats = np.arange(-90.+0.5*delta,90.,delta) + lsmask = np.reshape(np.fromstring(lsmaskf.read(),np.uint8),(nlats,nlons)) lsmaskf.close() # instance variable lsmask is set on first invocation, # it contains the land-sea mask interpolated to the native @@ -2708,28 +2708,28 @@ lons, lats = self(x, y, inverse=True) lon_0 = self.projparams['lon_0'] lats = lats[:,nx/2] - lons1 = (lon_0+180.)*npy.ones(lons.shape[0],npy.float64) - lons2 = (lon_0-180.)*npy.ones(lons.shape[0],npy.float64) + lons1 = (lon_0+180.)*np.ones(lons.shape[0],np.float64) + lons2 = (lon_0-180.)*np.ones(lons.shape[0],np.float64) xmax,ytmp = self(lons1,lats) xmin,ytmp = self(lons2,lats) for j in range(lats.shape[0]): xx = x[j,:] - mask[j,:]=npy.where(npy.logical_or(xx<xmin[j],xx>xmax[j]),\ + mask[j,:]=np.where(np.logical_or(xx<xmin[j],xx>xmax[j]),\ 255,mask[j,:]) self.lsmask = mask # optionally, set lakes to ocean color. if lakes: - mask = npy.where(self.lsmask==2,0,self.lsmask) + mask = np.where(self.lsmask==2,0,self.lsmask) else: mask = self.lsmask ny, nx = mask.shape - rgba = npy.ones((ny,nx,4),npy.uint8) - rgba_land = npy.array(rgba_land,npy.uint8) - rgba_ocean = npy.array(rgba_ocean,npy.uint8) + rgba = np.ones((ny,nx,4),np.uint8) + rgba_land = np.array(rgba_land,np.uint8) + rgba_ocean = np.array(rgba_ocean,np.uint8) for k in range(4): - rgba[:,:,k] = npy.where(mask,rgba_land[k],rgba_ocean[k]) + rgba[:,:,k] = np.where(mask,rgba_land[k],rgba_ocean[k]) # make points outside projection limb transparent. - rgba[:,:,3] = npy.where(mask==255,0,rgba[:,:,3]) + rgba[:,:,3] = np.where(mask==255,0,rgba[:,:,3]) # plot mask as rgba image. im = self.imshow(rgba,interpolation='nearest',ax=ax,**kwargs) return im @@ -2760,10 +2760,10 @@ from matplotlib.image import pil_to_array if not kwargs.has_key('ax') and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif not kwargs.has_key('ax') and self.ax is not None: ax = self.ax else: @@ -2786,12 +2786,12 @@ pilImage = Image.open(self._bm_file) self._bm_rgba = pil_to_array(pilImage) # convert to normalized floats. - self._bm_rgba = self._bm_rgba.astype(npy.float32)/255. + self._bm_rgba = self._bm_rgba.astype(np.float32)/255. # define lat/lon grid that image spans. nlons = self._bm_rgba.shape[1]; nlats = self._bm_rgba.shape[0] delta = 360./float(nlons) - self._bm_lons = npy.arange(-180.+0.5*delta,180.,delta) - self._bm_lats = npy.arange(-90.+0.5*delta,90.,delta) + self._bm_lons = np.arange(-180.+0.5*delta,180.,delta) + self._bm_lats = np.arange(-90.+0.5*delta,90.,delta) if self.projection != 'cyl': if newfile or not hasattr(self,'_bm_rgba_warped'): @@ -2799,10 +2799,10 @@ # projection grid. # nx and ny chosen to have roughly the # same horizontal res as original image. - dx = 2.*npy.pi*self.rmajor/float(nlons) + dx = 2.*np.pi*self.rmajor/float(nlons) nx = int((self.xmax-self.xmin)/dx)+1 ny = int((self.ymax-self.ymin)/dx)+1 - self._bm_rgba_warped = npy.zeros((ny,nx,4),npy.float64) + self._bm_rgba_warped = np.zeros((ny,nx,4),np.float64) # interpolate rgba values from geographic coords (proj='cyl') # to map projection coords. # if masked=True, values outside of @@ -2814,8 +2814,8 @@ # for ortho,geos mask pixels outside projection limb. if self.projection in ['geos','ortho']: lonsr,latsr = self(x,y,inverse=True) - mask = ma.zeros((nx,ny,4),npy.int8) - mask[:,:,0] = npy.logical_or(lonsr>1.e20,latsr>1.e30) + mask = ma.zeros((nx,ny,4),np.int8) + mask[:,:,0] = np.logical_or(lonsr>1.e20,latsr>1.e30) for k in range(1,4): mask[:,:,k] = mask[:,:,0] self._bm_rgba_warped = \ @@ -2862,10 +2862,10 @@ # get current axes instance (if none specified). if ax is None and self.ax is None: try: - ax = pylab.gca() + ax = plt.gca() except: - import pylab - ax = pylab.gca() + import matplotlib.pyplot as plt + ax = plt.gca() elif ax is None and self.ax is not None: ax = self.ax # not valid for cylindrical projection @@ -3091,8 +3091,8 @@ else: # irregular (but still rectilinear) input grid. xoutflat = xout.flatten(); youtflat = yout.flatten() - ix = (npy.searchsorted(xin,xoutflat)-1).tolist() - iy = (npy.searchsorted(yin,youtflat)-1).tolist() + ix = (np.searchsorted(xin,xoutflat)-1).tolist() + iy = (np.searchsorted(yin,youtflat)-1).tolist() xoutflat = xoutflat.tolist(); xin = xin.tolist() youtflat = youtflat.tolist(); yin = yin.tolist() xcoords = []; ycoords = [] @@ -3110,33 +3110,33 @@ ycoords.append(len(yin)) # outside range on upper end else: ycoords.append(float(j)+(youtflat[m]-yin[j])/(yin[j+1]-yin[j])) - xcoords = npy.reshape(xcoords,xout.shape) - ycoords = npy.reshape(ycoords,yout.shape) + xcoords = np.reshape(xcoords,xout.shape) + ycoords = np.reshape(ycoords,yout.shape) # data outside range xin,yin will be clipped to # values on boundary. if masked: - xmask = npy.logical_or(npy.less(xcoords,0),npy.greater(xcoords,len(xin)-1)) - ymask = npy.logical_or(npy.less(ycoords,0),npy.greater(ycoords,len(yin)-1)) - xymask = npy.logical_or(xmask,ymask) - xcoords = npy.clip(xcoords,0,len(xin)-1) - ycoords = npy.clip(ycoords,0,len(yin)-1) + xmask = np.logical_or(np.less(xcoords,0),np.greater(xcoords,len(xin)-1)) + ymask = np.logical_or(np.less(ycoords,0),np.greater(ycoords,len(yin)-1)) + xymask = np.logical_or(xmask,ymask) + xcoords = np.clip(xcoords,0,len(xin)-1) + ycoords = np.clip(ycoords,0,len(yin)-1) # interpolate to output grid using bilinear interpolation. if order == 1: - xi = xcoords.astype(npy.int32) - yi = ycoords.astype(npy.int32) + xi = xcoords.astype(np.int32) + yi = ycoords.astype(np.int32) xip1 = xi+1 yip1 = yi+1 - xip1 = npy.clip(xip1,0,len(xin)-1) - yip1 = npy.clip(yip1,0,len(yin)-1) - delx = xcoords-xi.astype(npy.float32) - dely = ycoords-yi.astype(npy.float32) + xip1 = np.clip(xip1,0,len(xin)-1) + yip1 = np.clip(yip1,0,len(yin)-1) + delx = xcoords-xi.astype(np.float32) + dely = ycoords-yi.astype(np.float32) dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \ delx*dely*datain[yip1,xip1] + \ (1.-delx)*dely*datain[yip1,xi] + \ delx*(1.-dely)*datain[yi,xip1] elif order == 0: - xcoordsi = npy.around(xcoords).astype(npy.int32) - ycoordsi = npy.around(ycoords).astype(npy.int32) + xcoordsi = np.around(xcoords).astype(np.int32) + ycoordsi = np.around(ycoords).astype(np.int32) dataout = datain[ycoordsi,xcoordsi] else: raise ValueError,'order keyword must be 0 or 1' @@ -3145,7 +3145,7 @@ newmask = ma.mask_or(ma.getmask(dataout), xymask) dataout = ma.masked_array(dataout,mask=newmask) elif masked and is_scalar(masked): - dataout = npy.where(xymask,masked,dataout) + dataout = np.where(xymask,masked,dataout) return dataout def shiftgrid(lon0,datain,lonsin,start=True): @@ -3163,13 +3163,13 @@ returns dataout,lonsout (data and longitudes on shifted grid). """ - if npy.fabs(lonsin[-1]-lonsin[0]-360.) > 1.e-4: + if np.fabs(lonsin[-1]-lonsin[0]-360.) > 1.e-4: raise ValueError, 'cyclic point not included' if lon0 < lonsin[0] or lon0 > lonsin[-1]: raise ValueError, 'lon0 outside of range of lonsin' - i0 = npy.argmin(npy.fabs(lonsin-lon0)) - dataout = npy.zeros(datain.shape,datain.dtype) - lonsout = npy.zeros(lonsin.shape,lonsin.dtype) + i0 = np.argmin(np.fabs(lonsin-lon0)) + dataout = np.zeros(datain.shape,datain.dtype) + lonsout = np.zeros(lonsin.shape,lonsin.dtype) if start: lonsout[0:len(lonsin)-i0] = lonsin[i0:] else: @@ -3193,13 +3193,13 @@ if hasattr(arrin,'mask'): arrout = ma.zeros((nlats,nlons+1),arrin.dtype) else: - arrout = npy.zeros((nlats,nlons+1),arrin.dtype) + arrout = np.zeros((nlats,nlons+1),arrin.dtype) arrout[:,0:nlons] = arrin[:,:] arrout[:,nlons] = arrin[:,0] if hasattr(lonsin,'mask'): lonsout = ma.zeros(nlons+1,lonsin.dtype) else: - lonsout = npy.zeros(nlons+1,lonsin.dtype) + lonsout = np.zeros(nlons+1,lonsin.dtype) lonsout[0:nlons] = lonsin[:] lonsout[nlons] = lonsin[-1] + lonsin[1]-lonsin[0] return arrout,lonsout This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 5200 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5200&view=rev Author: jswhit Date: 2008年05月20日 08:08:21 -0700 (2008年5月20日) Log Message: ----------- fix problems caused by change to "import numpy as np" Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月20日 15:01:46 UTC (rev 5199) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月20日 15:08:21 UTC (rev 5200) @@ -1003,7 +1003,7 @@ yy = np.linspace(self.ymin, self.ymax, ny)[:-1] x = len(yy)*[self.xmin]; y = yy.tolist() # top (y = ymax, xmin <= x <= xmax) - xx = np.np.linspace(self.xmin, self.xmax, nx)[:-1] + xx = np.linspace(self.xmin, self.xmax, nx)[:-1] x = x + xx.tolist() y = y + len(xx)*[self.ymax] # right side (x = xmax, ymin <= y <= ymax) @@ -1193,7 +1193,7 @@ ax = self.ax # get axis background color. axisbgc = ax.get_axis_bgcolor() - np = 0 + npt = 0 for x,y in self.coastpolygons: xa = np.array(x,np.float32) ya = np.array(y,np.float32) @@ -1213,7 +1213,7 @@ hasp3 = np.sum(test1*test4) if not hasp1 or not hasp2 or not hasp3 or not hasp4: xy = zip(xa.tolist(),ya.tolist()) - if self.coastpolygontypes[np] not in [2,4]: + if self.coastpolygontypes[npt] not in [2,4]: poly = Polygon(xy,facecolor=color,edgecolor=color,linewidth=0) else: # lakes filled with background color by default if lake_color is None: @@ -1223,7 +1223,7 @@ if zorder is not None: poly.set_zorder(zorder) ax.add_patch(poly) - np = np + 1 + npt = npt + 1 # set axes limits to fit map region. self.set_axes_limits(ax=ax) return poly This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 5201 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5201&view=rev Author: jswhit Date: 2008年05月20日 08:14:39 -0700 (2008年5月20日) Log Message: ----------- change variable name Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月20日 15:08:21 UTC (rev 5200) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月20日 15:14:39 UTC (rev 5201) @@ -1193,7 +1193,7 @@ ax = self.ax # get axis background color. axisbgc = ax.get_axis_bgcolor() - npt = 0 + npoly = 0 for x,y in self.coastpolygons: xa = np.array(x,np.float32) ya = np.array(y,np.float32) @@ -1213,7 +1213,7 @@ hasp3 = np.sum(test1*test4) if not hasp1 or not hasp2 or not hasp3 or not hasp4: xy = zip(xa.tolist(),ya.tolist()) - if self.coastpolygontypes[npt] not in [2,4]: + if self.coastpolygontypes[npoly] not in [2,4]: poly = Polygon(xy,facecolor=color,edgecolor=color,linewidth=0) else: # lakes filled with background color by default if lake_color is None: @@ -1223,7 +1223,7 @@ if zorder is not None: poly.set_zorder(zorder) ax.add_patch(poly) - npt = npt + 1 + npoly = npoly + 1 # set axes limits to fit map region. self.set_axes_limits(ax=ax) return poly This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 5341 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5341&view=rev Author: jswhit Date: 2008年05月31日 05:57:10 -0700 (2008年5月31日) Log Message: ----------- fixed warpimage method to work with recent matplotlib change "2008-05-29 matplotlib.image.imread now no longer always returns RGBA" Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月31日 12:26:21 UTC (rev 5340) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月31日 12:57:10 UTC (rev 5341) @@ -2802,12 +2802,12 @@ dx = 2.*np.pi*self.rmajor/float(nlons) nx = int((self.xmax-self.xmin)/dx)+1 ny = int((self.ymax-self.ymin)/dx)+1 - self._bm_rgba_warped = np.zeros((ny,nx,4),np.float64) + self._bm_rgba_warped = np.ones((ny,nx,4),np.float64) # interpolate rgba values from geographic coords (proj='cyl') # to map projection coords. # if masked=True, values outside of # projection limb will be masked. - for k in range(4): + for k in range(3): self._bm_rgba_warped[:,:,k],x,y = \ self.transform_scalar(self._bm_rgba[:,:,k],\ self._bm_lons,self._bm_lats,nx,ny,returnxy=True) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 5343 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5343&view=rev Author: jswhit Date: 2008年05月31日 06:05:26 -0700 (2008年5月31日) Log Message: ----------- use "import numpy.ma as ma" instead of "from numpy import ma" Modified Paths: -------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月31日 13:03:28 UTC (rev 5342) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年05月31日 13:05:26 UTC (rev 5343) @@ -34,7 +34,7 @@ import pyproj, sys, os, math, dbflib from proj import Proj import numpy as np -from numpy import ma +import numpy.ma as ma from shapelib import ShapeFile import _geoslib, pupynere, netcdftime This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
Revision: 5547 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5547&view=rev Author: jswhit Date: 2008年06月15日 05:08:13 -0700 (2008年6月15日) Log Message: ----------- moved to __init__.py so Sphinx can find docstrings. Removed Paths: ------------- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py Deleted: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py =================================================================== --- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年06月15日 11:57:51 UTC (rev 5546) +++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/basemap.py 2008年06月15日 12:08:13 UTC (rev 5547) @@ -1,3344 +0,0 @@ -""" -Module for plotting data on maps with matplotlib. - -Contains the Basemap class (which does most of the -heavy lifting), and the following functions: - -NetCDFFile: Read local and remote NetCDF datasets. - -interp: bilinear interpolation between rectilinear grids. - -shiftgrid: shifts global lat/lon grids east or west. - -addcyclic: Add cyclic (wraparound) point in longitude. - -num2date: convert from a numeric time value to a datetime object. - -date2num: convert from a datetime object to a numeric time value. -""" -from matplotlib import __version__ as _matplotlib_version -from matplotlib.cbook import is_scalar, dedent -# check to make sure matplotlib is not too old. -_mpl_required_version = '0.98' -if _matplotlib_version < _mpl_required_version: - msg = dedent(""" - your matplotlib is too old - basemap requires version %s or - higher, you have version %s""" % - (_mpl_required_version,_matplotlib_version)) - raise ImportError(msg) -from matplotlib import rcParams, is_interactive, _pylab_helpers -from matplotlib.collections import LineCollection -from matplotlib.patches import Ellipse, Circle, Polygon -from matplotlib.lines import Line2D -from matplotlib.transforms import Bbox -import pyproj, sys, os, math, dbflib -from proj import Proj -import numpy as np -import numpy.ma as ma -from shapelib import ShapeFile -import _geoslib, pupynere, netcdftime - -# basemap data files now installed in lib/matplotlib/toolkits/basemap/data -basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data']) - -__version__ = '0.99' - -# supported map projections. -_projnames = {'cyl' : 'Cylindrical Equidistant', - 'merc' : 'Mercator', - 'tmerc' : 'Transverse Mercator', - 'omerc' : 'Oblique Mercator', - 'mill' : 'Miller Cylindrical', - 'lcc' : 'Lambert Conformal', - 'laea' : 'Lambert Azimuthal Equal Area', - 'nplaea' : 'North-Polar Lambert Azimuthal', - 'splaea' : 'South-Polar Lambert Azimuthal', - 'eqdc' : 'Equidistant Conic', - 'aeqd' : 'Azimuthal Equidistant', - 'npaeqd' : 'North-Polar Azimuthal Equidistant', - 'spaeqd' : 'South-Polar Azimuthal Equidistant', - 'aea' : 'Albers Equal Area', - 'stere' : 'Stereographic', - 'npstere' : 'North-Polar Stereographic', - 'spstere' : 'South-Polar Stereographic', - 'cass' : 'Cassini-Soldner', - 'poly' : 'Polyconic', - 'ortho' : 'Orthographic', - 'geos' : 'Geostationary', - 'sinu' : 'Sinusoidal', - 'moll' : 'Mollweide', - 'robin' : 'Robinson', - 'gnom' : 'Gnomonic', - } -supported_projections = [] -for _items in _projnames.iteritems(): - supported_projections.append("'%s' = %s\n" % (_items)) -supported_projections = ''.join(supported_projections) - -# projection specific parameters. -projection_params = {'cyl' : 'corners only (no width/height)', - 'merc' : 'corners plus lat_ts (no width/height)', - 'tmerc' : 'lon_0,lat_0', - 'omerc' : 'lon_0,lat_0,lat_1,lat_2,lon_1,lon_2,no_rot', - 'mill' : 'corners only (no width/height)', - 'lcc' : 'lon_0,lat_0,lat_1,lat_2', - 'laea' : 'lon_0,lat_0', - 'nplaea' : 'bounding_lat,lon_0,lat_0,no corners or width/height', - 'splaea' : 'bounding_lat,lon_0,lat_0,no corners or width/height', - 'eqdc' : 'lon_0,lat_0,lat_1,lat_2', - 'aeqd' : 'lon_0,lat_0', - 'npaeqd' : 'bounding_lat,lon_0,lat_0,no corners or width/height', - 'spaeqd' : 'bounding_lat,lon_0,lat_0,no corners or width/height', - 'aea' : 'lon_0,lat_0,lat_1', - 'stere' : 'lon_0,lat_0,lat_ts', - 'npstere' : 'bounding_lat,lon_0,lat_0,no corners or width/height', - 'spstere' : 'bounding_lat,lon_0,lat_0,no corners or width/height', - 'cass' : 'lon_0,lat_0', - 'poly' : 'lon_0,lat_0', - 'ortho' : 'lon_0,lat_0,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height', - 'geos' : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height', - 'sinu' : 'lon_0,lat_0,no corners or width/height', - 'moll' : 'lon_0,lat_0,no corners or width/height', - 'robin' : 'lon_0,lat_0,no corners or width/height', - 'gnom' : 'lon_0,lat_0', - } - -# The __init__ docstring is pulled out here because it is so long; -# Having it in the usual place makes it hard to get from the -# __init__ argument list to the code that uses the arguments. -_Basemap_init_doc = """ - create a Basemap instance. - - This sets up a basemap with specified map projection. - and creates the coastline data structures in native map projection - coordinates. - - arguments: - - projection - map projection. Supported projections are: -%(supported_projections)s - Default is 'cyl'. - - For most map projections, the map projection region can either be - specified by setting these keywords: - - llcrnrlon - longitude of lower left hand corner of the desired map domain (degrees). - llcrnrlat - latitude of lower left hand corner of the desired map domain (degrees). - urcrnrlon - longitude of upper right hand corner of the desired map domain (degrees). - urcrnrlat - latitude of upper right hand corner of the desired map domain (degrees). - - or these keywords: - - width - width of desired map domain in projection coordinates (meters). - height - height of desired map domain in projection coordinates (meters). - lon_0 - center of desired map domain (in degrees). - lat_0 - center of desired map domain (in degrees). - - For 'sinu', 'moll', 'npstere', 'spstere', 'nplaea', 'splaea', 'nplaea', - 'splaea', 'npaeqd', 'spaeqd' or 'robin', the values of - llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,width and height are ignored - (because either they are computed internally, or entire globe is - always plotted). For the cylindrical projections - ('cyl','merc' and 'mill'), the default is to use - llcrnrlon=-180,llcrnrlat=-90, urcrnrlon=180 and urcrnrlat=90). For all other - projections except 'ortho' and 'geos', either the lat/lon values of the - corners or width and height must be specified by the user. - For 'ortho' and 'geos', the lat/lon values of the corners may be specified, - or the x/y values of the corners (llcrnrx,llcrnry,urcrnrx,urcrnry) in the - coordinate system of the global projection (with x=0,y=0 at the center - of the global projection). If the corners are not specified, - the entire globe is plotted. - - resolution - resolution of boundary database to use. Can be 'c' (crude), - 'l' (low), 'i' (intermediate), 'h' (high), 'f' (full) or None. - If None, no boundary data will be read in (and class methods - such as drawcoastlines will raise an exception if invoked). - Resolution drops off by roughly 80%% - between datasets. Higher res datasets are much slower to draw. - Default 'c'. Coastline data is from the GSHHS - (http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html). - State, country and river datasets from the Generic Mapping - Tools (http://gmt.soest.hawaii.edu). - - area_thresh - coastline or lake with an area smaller than area_thresh - in km^2 will not be plotted. Default 10000,1000,100,10,1 for resolution - 'c','l','i','h','f'. - - rsphere - radius of the sphere used to define map projection (default - 6370997 meters, close to the arithmetic mean radius of the earth). If - given as a sequence, the first two elements are interpreted as - the the radii of the major and minor axes of an ellipsoid. Note: sometimes - an ellipsoid is specified by the major axis and an 'inverse flattening - parameter' (if). The minor axis (b) can be computed from the major axis (a) - and the inverse flattening parameter using the formula if = a/(a-b). - - suppress_ticks - suppress automatic drawing of axis ticks and labels - in map projection coordinates. Default False, so parallels and meridians - can be labelled instead. If parallel or meridian labelling is requested - (using drawparallels and drawmeridians methods), automatic tick labelling - will be supressed even is suppress_ticks=False. suppress_ticks=False - is useful if you want to use your own custom tick formatter, or - if you want to let matplotlib label the axes in meters - using native map projection coordinates - - anchor - determines how map is placed in axes rectangle (passed to - axes.set_aspect). Default is 'C', which means map is centered. - Allowed values are ['C', 'SW', 'S', 'SE', 'E', 'NE', 'N', 'NW', 'W']. - - ax - set default axes instance (default None - pylab.gca() may be used - to get the current axes instance). If you don't want pylab to be imported, - you can either set this to a pre-defined axes instance, or use the 'ax' - keyword in each Basemap method call that does drawing. In the first case, - all Basemap method calls will draw to the same axes instance. In the - second case, you can draw to different axes with the same Basemap instance. - You can also use the 'ax' keyword in individual method calls to - selectively override the default axes instance. - - The following parameters are map projection parameters which all default to - None. Not all parameters are used by all projections, some are ignored. - The module variable 'projection_params' is a dictionary which - lists which parameters apply to which projections. - - lat_ts - latitude of true scale for mercator projection, optional - for stereographic projection. - lat_1 - first standard parallel for lambert conformal, albers - equal area projection and equidistant conic projections. Latitude of one - of the two points on the projection centerline for oblique mercator. - If lat_1 is not given, but lat_0 is, lat_1 is set to lat_0 for - lambert conformal, albers equal area and equidistant conic. - lat_2 - second standard parallel for lambert conformal, albers - equal area projection and equidistant conic projections. Latitude of one - of the two points on the projection centerline for oblique mercator. - If lat_2 is not given, it is set to lat_1 for - lambert conformal, albers equal area and equidistant conic. - lon_1 - longitude of one of the two points on the projection centerline - for oblique mercator. - lon_2 - longitude of one of the two points on the projection centerline - for oblique mercator. - no_rot - only used by oblique mercator. If set to True, the map projection - coordinates will not be rotated to true North. Default is False (projection - coordinates are automatically rotated). - lat_0 - central latitude (y-axis origin) - used by all projections, - Must be equator for mercator projection. - lon_0 - central meridian (x-axis origin) - used by all projections, - boundinglat - bounding latitude for pole-centered projections (npstere,spstere, - nplaea,splaea,npaeqd,spaeqd). These projections are square regions centered - on the north or south pole. The longitude lon_0 is at 6-o'clock, and the - latitude circle boundinglat is tangent to the edge of the map at lon_0. - satellite_height - height of satellite (in m) above equator - - only relevant for geostationary projections ('geos'). Default 35,786 km. - - Here are the most commonly used class methods (see the docstring - for each for more details): - - To draw a graticule grid (labelled latitude and longitude lines) - use the drawparallels and drawmeridians methods. - - To draw coastline, rivers and political boundaries, use the - the drawcontinents, drawrivers, drawcountries and drawstates methods. - - To fill the continents and inland lakes, use the fillcontinents method. - - To draw the boundary of the map projection region, and fill the - interior a certain color, use the drawmapboundary method. - - The contour, contourf, pcolor, pcolormesh, plot, scatter - quiver and imshow methods use the corresponding matplotlib axes - methods to draw on the map. - - The transform_scalar method can be used to interpolate regular - lat/lon grids of scalar data to native map projection grids. - - The transform_vector method can be used to interpolate and rotate - regular lat/lon grids of vector data to native map projection grids. - - The rotate_vector method rotates a vector field from lat/lon - coordinates into map projections coordinates, without doing any - interpolation. - - The readshapefile method can be used to read in data from ESRI - shapefiles. - - The drawgreatcircle method draws great circles on the map. -""" % locals() - -# unsupported projection error message. -_unsupported_projection = ["'%s' is an unsupported projection.\n"] -_unsupported_projection.append("The supported projections are:\n") -_unsupported_projection.append(supported_projections) -_unsupported_projection = ''.join(_unsupported_projection) - -def _validated_ll(param, name, minval, maxval): - param = float(param) - if param > maxval or param < minval: - raise ValueError('%s must be between %f and %f degrees' % - (name, minval, maxval)) - return param - -def _insert_validated(d, param, name, minval, maxval): - if param is not None: - d[name] = _validated_ll(param, name, minval, maxval) - -class Basemap(object): - """ - Class for plotting data on map projections with matplotlib. - See __init__ docstring for details on how to create a class - instance for a given map projection. - - Useful instance variables: - - projection - map projection. Print the module variable - "supported_projections" to see a list. - aspect - map aspect ratio (size of y dimension / size of x dimension). - llcrnrlon - longitude of lower left hand corner of the desired map domain. - llcrnrlon - latitude of lower left hand corner of the desired map domain. - urcrnrlon - longitude of upper right hand corner of the desired map domain. - urcrnrlon - latitude of upper right hand corner of the desired map domain. - llcrnrx,llcrnry,urcrnrx,urcrnry - corners of map domain in projection coordinates. - rmajor,rminor - equatorial and polar radii of ellipsoid used (in meters). - resolution - resolution of boundary dataset being used ('c' for crude, - 'l' for low, etc.). If None, no boundary dataset is associated with the - Basemap instance. - srs - a string representing the 'spatial reference system' for the map - projection as defined by PROJ.4. - - Example Usage: - - >>> from mpl_toolkits.basemap import Basemap - >>> import numpy as np - >>> import matplotlib.pyplot as plt - >>> import matplotlib.mlab as mlab - >>> # read in topo data (on a regular lat/lon grid) - >>> etopo = mlab.load('etopo20data.gz') - >>> lons = mlab.load('etopo20lons.gz') - >>> lats = mlab.load('etopo20lats.gz') - >>> # create Basemap instance for Robinson projection. - >>> m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1])) - >>> # compute native map projection coordinates for lat/lon grid. - >>> x, y = m(*np.meshgrid(lons,lats)) - >>> # make filled contour plot. - >>> cs = m.contourf(x,y,etopo,30,cmap=plt.cm.jet) - >>> m.drawcoastlines() # draw coastlines - >>> m.drawmapboundary() # draw a line around the map region - >>> m.drawparallels(np.arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels - >>> m.drawmeridians(np.arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians - >>> plt.title('Robinson Projection') # add a title - >>> plt.show() - - [this example (simpletest.py) plus many others can be found in the - examples directory of source distribution. The "OO" version of this - example (which does not use pylab) is called "simpletest_oo.py".] - """ - - def __init__(self, llcrnrlon=None, llcrnrlat=None, - urcrnrlon=None, urcrnrlat=None, - llcrnrx=None, llcrnry=None, - urcrnrx=None, urcrnry=None, - width=None, height=None, - projection='cyl', resolution='c', - area_thresh=None, rsphere=6370997.0, - lat_ts=None, - lat_1=None, lat_2=None, - lat_0=None, lon_0=None, - lon_1=None, lon_2=None, - no_rot=False, - suppress_ticks=True, - satellite_height=35786000, - boundinglat=None, - anchor='C', - ax=None): - # docstring is added after __init__ method definition - - # where to put plot in figure (default is 'C' or center) - self.anchor = anchor - # map projection. - self.projection = projection - - # set up projection parameter dict. - projparams = {} - projparams['proj'] = projection - try: - if rsphere[0] > rsphere[1]: - projparams['a'] = rsphere[0] - projparams['b'] = rsphere[1] - else: - projparams['a'] = rsphere[1] - projparams['b'] = rsphere[0] - except: - if projection == 'tmerc': - # use bR_a instead of R because of obscure bug - # in proj4 for tmerc projection. - projparams['bR_a'] = rsphere - else: - projparams['R'] = rsphere - # set units to meters. - projparams['units']='m' - # check for sane values of lon_0, lat_0, lat_ts, lat_1, lat_2 - _insert_validated(projparams, lat_0, 'lat_0', -90, 90) - _insert_validated(projparams, lat_1, 'lat_1', -90, 90) - _insert_validated(projparams, lat_2, 'lat_2', -90, 90) - _insert_validated(projparams, lat_ts, 'lat_ts', -90, 90) - _insert_validated(projparams, lon_0, 'lon_0', -360, 720) - _insert_validated(projparams, lon_1, 'lon_1', -360, 720) - _insert_validated(projparams, lon_2, 'lon_2', -360, 720) - if projection == 'geos': - projparams['h'] = satellite_height - # check for sane values of projection corners. - using_corners = (None not in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]) - if using_corners: - self.llcrnrlon = _validated_ll(llcrnrlon, 'llcrnrlon', -360, 720) - self.urcrnrlon = _validated_ll(urcrnrlon, 'urcrnrlon', -360, 720) - self.llcrnrlat = _validated_ll(llcrnrlat, 'llcrnrlat', -90, 90) - self.urcrnrlat = _validated_ll(urcrnrlat, 'urcrnrlat', -90, 90) - - # for each of the supported projections, - # compute lat/lon of domain corners - # and set values in projparams dict as needed. - - if projection in ['lcc', 'eqdc', 'aea']: - # if lat_0 is given, but not lat_1, - # set lat_1=lat_0 - if lat_1 is None and lat_0 is not None: - lat_1 = lat_0 - projparams['lat_1'] = lat_1 - if lat_1 is None or lon_0 is None: - raise ValueError('must specify lat_1 or lat_0 and lon_0 for %s basemap (lat_2 is optional)' % _projnames[projection]) - if lat_2 is None: - projparams['lat_2'] = lat_1 - 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' - 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) - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection == 'stere': - if lat_0 is None or lon_0 is None: - raise ValueError, 'must specify lat_0 and lon_0 for Stereographic basemap (lat_ts is optional)' - 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' - 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) - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection in ['spstere', 'npstere', - 'splaea', 'nplaea', - 'spaeqd', 'npaeqd']: - if boundinglat is None or lon_0 is None: - raise ValueError('must specify boundinglat and lon_0 for %s basemap' % _projnames[projection]) - if projection[0] == 's': - sgn = -1 - else: - sgn = 1 - rootproj = projection[2:] - projparams['proj'] = rootproj - if rootproj == 'stere': - projparams['lat_ts'] = sgn * 90. - projparams['lat_0'] = sgn * 90. - self.llcrnrlon = lon_0 - sgn*45. - self.urcrnrlon = lon_0 + sgn*135. - proj = pyproj.Proj(projparams) - x,y = proj(lon_0,boundinglat) - lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True) - self.urcrnrlat = self.llcrnrlat - if width is not None or height is not None: - print 'warning: width and height keywords ignored for %s projection' % _projnames[projection] - elif projection == 'laea': - if lat_0 is None or lon_0 is None: - raise ValueError, 'must specify lat_0 and lon_0 for Lambert Azimuthal 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' - 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) - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection == 'merc': - if lat_ts is None: - raise ValueError, 'must specify lat_ts for Mercator basemap' - # clip plot region to be within -89.99S to 89.99N - # (mercator is singular at poles) - if not using_corners: - llcrnrlon = -180. - llcrnrlat = -90. - urcrnrlon = 180 - urcrnrlat = 90. - if llcrnrlat < -89.99: llcrnrlat = -89.99 - if llcrnrlat > 89.99: llcrnrlat = 89.99 - if urcrnrlat < -89.99: urcrnrlat = -89.99 - if urcrnrlat > 89.99: urcrnrlat = 89.99 - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - if width is not None or height is not None: - print 'warning: width and height keywords ignored for %s projection' % self.projection - elif projection in ['tmerc','gnom','cass','poly'] : - if lat_0 is None or lon_0 is None: - raise ValueError, 'must specify lat_0 and lon_0 for Transverse Mercator, Gnomonic, Cassini-Soldnerr Polyconic 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' - 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) - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection == 'ortho': - if not projparams.has_key('R'): - raise ValueError, 'orthographic projection only works for perfect spheres - not ellipsoids' - if lat_0 is None or lon_0 is None: - raise ValueError, 'must specify lat_0 and lon_0 for Orthographic basemap' - if width is not None or height is not None: - print 'warning: width and height keywords ignored for %s projection' % self.projection - if not using_corners: - llcrnrlon = -180. - llcrnrlat = -90. - urcrnrlon = 180 - urcrnrlat = 90. - self._fulldisk = True - else: - self._fulldisk = False - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - # FIXME: won't work for points exactly on equator?? - if np.abs(lat_0) < 1.e-2: lat_0 = 1.e-2 - projparams['lat_0'] = lat_0 - elif projection == 'geos': - if lon_0 is None: - raise ValueError, 'must specify lon_0 for Geostationary basemap' - if width is not None or height is not None: - print 'warning: width and height keywords ignored for %s projection' % self.projection - if not using_corners: - llcrnrlon = -180. - llcrnrlat = -90. - urcrnrlon = 180 - urcrnrlat = 90. - self._fulldisk = True - else: - self._fulldisk = False - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection in ['moll','robin','sinu']: - if lon_0 is None: - raise ValueError, 'must specify lon_0 for Robinson, Mollweide, or Sinusoidal basemap' - if width is not None or height is not None: - print 'warning: width and height keywords ignored for %s projection' % self.projection - llcrnrlon = -180. - llcrnrlat = -90. - urcrnrlon = 180 - urcrnrlat = 90. - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection == 'omerc': - if lat_1 is None or lon_1 is None or lat_2 is None or lon_2 is None: - raise ValueError, 'must specify lat_1,lon_1 and lat_2,lon_2 for Oblique Mercator basemap' - projparams['lat_1'] = lat_1 - projparams['lon_1'] = lon_1 - projparams['lat_2'] = lat_2 - projparams['lon_2'] = lon_2 - projparams['lat_0'] = lat_0 - if no_rot: - projparams['no_rot']='' - #if not using_corners: - # raise ValueError, 'cannot specify map region with width and height keywords for this projection, please specify lat/lon values of corners' - 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' - 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) - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection == 'aeqd': - if lat_0 is None or lon_0 is None: - 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' - 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) - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - elif projection == 'mill': - if not using_corners: - llcrnrlon = -180. - llcrnrlat = -90. - urcrnrlon = 180 - urcrnrlat = 90. - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - if width is not None or height is not None: - print 'warning: width and height keywords ignored for %s projection' % self.projection - elif projection == 'cyl': - if not using_corners: - llcrnrlon = -180. - llcrnrlat = -90. - urcrnrlon = 180 - urcrnrlat = 90. - self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat - self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat - if width is not None or height is not None: - print 'warning: width and height keywords ignored for %s projection' % self.projection - else: - raise ValueError(_unsupported_projection % projection) - - # initialize proj4 - proj = Proj(projparams,self.llcrnrlon,self.llcrnrlat,self.urcrnrlon,self.urcrnrlat) - - # make sure axis ticks are suppressed. - self.noticks = suppress_ticks - - # make Proj instance a Basemap instance variable. - self.projtran = proj - # copy some Proj attributes. - atts = ['rmajor','rminor','esq','flattening','ellipsoid','projparams'] - for att in atts: - self.__dict__[att] = proj.__dict__[att] - # these only exist for geostationary projection. - if hasattr(proj,'_width'): - self.__dict__['_width'] = proj.__dict__['_width'] - if hasattr(proj,'_height'): - self.__dict__['_height'] = proj.__dict__['_height'] - # spatial reference string (useful for georeferencing output - # images with gdal_translate). - if hasattr(self,'_proj4'): - self.srs = proj._proj4.srs - else: - pjargs = [] - for key,value in self.projparams.iteritems(): - # 'cyl' projection translates to 'eqc' in PROJ.4 - if projection == 'cyl' and key == 'proj': - value = 'eqc' - # ignore x_0 and y_0 settings for 'cyl' projection - # (they are not consistent with what PROJ.4 uses) - elif projection == 'cyl' and key in ['x_0','y_0']: - continue - pjargs.append('+'+key+"="+str(value)+' ') - self.srs = ''.join(pjargs) - # set instance variables defining map region. - self.xmin = proj.xmin - self.xmax = proj.xmax - self.ymin = proj.ymin - self.ymax = proj.ymax - if projection == 'cyl': - self.aspect = (self.urcrnrlat-self.llcrnrlat)/(self.urcrnrlon-self.llcrnrlon) - else: - self.aspect = (proj.ymax-proj.ymin)/(proj.xmax-proj.xmin) - if projection in ['geos','ortho'] and \ - None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]: - self.llcrnrx = llcrnrx+0.5*proj.xmax - self.llcrnry = llcrnry+0.5*proj.ymax - self.urcrnrx = urcrnrx+0.5*proj.xmax - self.urcrnry = urcrnry+0.5*proj.ymax - self._fulldisk = False - else: - self.llcrnrx = proj.llcrnrx - self.llcrnry = proj.llcrnry - self.urcrnrx = proj.urcrnrx - self.urcrnry = proj.urcrnry - # set min/max lats for projection domain. - if projection in ['mill','cyl','merc']: - self.latmin = self.llcrnrlat - self.latmax = self.urcrnrlat - elif projection in ['ortho','geos','moll','robin','sinu']: - self.latmin = -90. - self.latmax = 90. - else: - lons, lats = self.makegrid(101,101) - self.latmin = lats.min() - self.latmax = lats.max() - - # if ax == None, pyplot.gca may be used. - self.ax = ax - self.lsmask = None - - # set defaults for area_thresh. - self.resolution = resolution - if area_thresh is None and resolution is not None: - if resolution == 'c': - area_thresh = 10000. - elif resolution == 'l': - area_thresh = 1000. - elif resolution == 'i': - area_thresh = 100. - elif resolution == 'h': - area_thresh = 10. - elif resolution == 'f': - area_thresh = 1. - else: - raise ValueError, "boundary resolution must be one of 'c','l','i','h' or 'f'" - self.area_thresh = area_thresh - # define map boundary polygon (in lat/lon coordinates) - self._boundarypolyll, self._boundarypolyxy = self._getmapboundary() - # read in coastline polygons, only keeping those that - # intersect map boundary polygon. - if self.resolution is not None: - self.coastsegs, self.coastpolygontypes = self._readboundarydata('gshhs') - # reformat for use in matplotlib.patches.Polygon. - self.coastpolygons = [] - # also, split coastline segments that jump across entire plot. - coastsegs = [] - for seg in self.coastsegs: - x, y = zip(*seg) - self.coastpolygons.append((x,y)) - x = np.array(x,np.float64); y = np.array(y,np.float64) - xd = (x[1:]-x[0:-1])**2 - yd = (y[1:]-y[0:-1])**2 - dist = np.sqrt(xd+yd) - split = dist > 5000000. - if np.sum(split) and self.projection not in ['merc','cyl','mill']: - ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist() - iprev = 0 - ind.append(len(xd)) - for i in ind: - # don't add empty lists. - if len(range(iprev,i)): - coastsegs.append(zip(x[iprev:i],y[iprev:i])) - iprev = i - else: - coastsegs.append(seg) - self.coastsegs = coastsegs - # set __init__'s docstring - __init__.__doc__ = _Basemap_init_doc - - def __call__(self,x,y,inverse=False): - """ - Calling a Basemap 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 non-cylindrical projections, the inverse transformation - always returns longitudes between -180 and 180 degrees. For - cylindrical projections (self.projection == 'cyl','mill' or 'merc') - the inverse transformation will return longitudes between - self.llcrnrlon and self.llcrnrlat. - - input arguments lon, lat can be either scalar floats or N arrays. - """ - return self.projtran(x,y,inverse=inverse) - - 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 self.projtran.makegrid(nx,ny,returnxy=returnxy) - - def _readboundarydata(self,name): - """ - read boundary data, clip to map projection region. - """ - msg = dedent(""" - Unable to open boundary dataset file. Only the 'crude', 'low', - 'intermediate' and 'high' resolution datasets are installed by default. - If you are requesting a 'full' resolution dataset, you may need to - download and install those files separately - (see the basemap README for details).""") - try: - bdatfile = open(os.path.join(basemap_datadir,name+'_'+self.resolution+'.dat'),'rb') - bdatmetafile = open(os.path.join(basemap_datadir,name+'meta_'+self.resolution+'.dat'),'r') - except: - raise IOError, msg - polygons = [] - polygon_types = [] - # coastlines are polygons, other boundaries are line segments. - if name == 'gshhs': - Shape = _geoslib.Polygon - else: - Shape = _geoslib.LineString - # see if map projection region polygon contains a pole. - NPole = _geoslib.Point(self(0.,90.)) - SPole = _geoslib.Point(self(0.,-90.)) - boundarypolyxy = self._boundarypolyxy - boundarypolyll = self._boundarypolyll - hasNP = NPole.within(boundarypolyxy) - hasSP = SPole.within(boundarypolyxy) - containsPole = hasNP or hasSP - # these projections cannot cross pole. - if containsPole and\ - self.projection in ['merc','mill','cyl','robin','moll','sinu','geos']: - #self.projection in ['tmerc','omerc','cass','merc','mill','cyl','robin','moll','sinu','geos']: - raise ValueError('%s projection cannot cross pole'%(self.projection)) - # make sure orthographic projection has containsPole=True - # we will compute the intersections in stereographic - # coordinates, then transform to orthographic. - if self.projection == 'ortho' and name == 'gshhs': - containsPole = True - lon_0=self.projparams['lon_0'] - lat_0=self.projparams['lat_0'] - re = self.projparams['R'] - # center of stereographic projection restricted to be - # nearest one of 6 points on the sphere (every 90 deg lat/lon). - lon0 = 90.*(np.around(lon_0/90.)) - lat0 = 90.*(np.around(lat_0/90.)) - if np.abs(int(lat0)) == 90: lon0=0. - maptran = pyproj.Proj(proj='stere',lon_0=lon0,lat_0=lat0,R=re) - # boundary polygon for orthographic projection - # in stereographic coorindates. - b = self._boundarypolyll.boundary - blons = b[:,0]; blats = b[:,1] - b[:,0], b[:,1] = maptran(blons, blats) - boundarypolyxy = _geoslib.Polygon(b) - for line in bdatmetafile: - linesplit = line.split() - area = float(linesplit[1]) - south = float(linesplit[3]) - north = float(linesplit[4]) - if area < 0.: area = 1.e30 - useit = self.latmax>=south and self.latmin<=north and area>self.area_thresh - if useit: - type = int(linesplit[0]) - npts = int(linesplit[2]) - offsetbytes = int(linesplit[5]) - bytecount = int(linesplit[6]) - bdatfile.seek(offsetbytes,0) - # read in binary string convert into an npts by 2 - # numpy array (first column is lons, second is lats). - polystring = bdatfile.read(bytecount) - # binary data is little endian. - b = np.array(np.fromstring(polystring,dtype='<f4'),'f8') - b.shape = (npts,2) - b2 = b.copy() - # if map boundary polygon is a valid one in lat/lon - # coordinates (i.e. it does not contain either pole), - # the intersections of the boundary geometries - # and the map projection region can be computed before - # transforming the boundary geometry to map projection - # coordinates (this saves time, especially for small map - # regions and high-resolution boundary geometries). - if not containsPole: - # close Antarctica. - if name == 'gshhs' and south < -68 and area > 10000: - lons = b[:,0] - lats = b[:,1] - lons2 = lons[:-2][::-1] - lats2 = lats[:-2][::-1] - lons1 = lons2 - 360. - lons3 = lons2 + 360. - lons = lons1.tolist()+lons2.tolist()+lons3.tolist() - lats = lats2.tolist()+lats2.tolist()+lats2.tolist() - lonstart,latstart = lons[0], lats[0] - lonend,latend = lons[-1], lats[-1] - lons.insert(0,lonstart) - lats.insert(0,-90.) - lons.append(lonend) - lats.append(-90.) - b = np.empty((len(lons),2),np.float64) - b[:,0] = lons; b[:,1] = lats - poly = _geoslib.Polygon(b) - antart = True - else: - poly = Shape(b) - antart = False - # create duplicate polygons shifted by -360 and +360 - # (so as to properly treat polygons that cross - # Greenwich meridian). - if not antart: - b2[:,0] = b[:,0]-360 - poly1 = Shape(b2) - b2[:,0] = b[:,0]+360 - poly2 = Shape(b2) - polys = [poly1,poly,poly2] - else: # Antartica already extends from -360 to +720. - polys = [poly] - for poly in polys: - # if polygon instersects map projection - # region, process it. - if poly.intersects(boundarypolyll): - # if geometry intersection calculation fails, - # just move on. - try: - geoms = poly.intersection(boundarypolyll) - except: - continue - # iterate over geometries in intersection. - for psub in geoms: - # only coastlines are polygons, - # which have a 'boundary' attribute. - # otherwise, use 'coords' attribute - # to extract coordinates. - b = psub.boundary - blons = b[:,0]; blats = b[:,1] - # transformation from lat/lon to - # map projection coordinates. - bx, by = self(blons, blats) - polygons.append(zip(bx,by)) - polygon_types.append(type) - # if map boundary polygon is not valid in lat/lon - # coordinates, compute intersection between map - # projection region and boundary geometries in map - # projection coordinates. - else: - # transform coordinates from lat/lon - # to map projection coordinates. - # special case for ortho, compute coastline polygon - # vertices in stereographic coords. - if name == 'gshhs' and self.projection == 'ortho': - b[:,0], b[:,1] = maptran(b[:,0], b[:,1]) - else: - b[:,0], b[:,1] = self(b[:,0], b[:,1]) - goodmask = np.logical_and(b[:,0]<1.e20,b[:,1]<1.e20) - # if less than two points are valid in - # map proj coords, skip this geometry. - if np.sum(goodmask) <= 1: continue - if name != 'gshhs': - # if not a polygon, - # just remove parts of geometry that are undefined - # in this map projection. - bx = np.compress(goodmask, b[:,0]) - by = np.compress(goodmask, b[:,1]) - # for orthographic projection, all points - # outside map projection region are eliminated - # with the above step, so we're done. - if self.projection == 'ortho': - polygons.append(zip(bx,by)) - polygon_types.append(type) - continue - # create a GEOS geometry object. - poly = Shape(b) - # if geometry instersects map projection - # region, and doesn't have any invalid points, process it. - if goodmask.all() and poly.intersects(boundarypolyxy): - # if geometry intersection calculation fails, - # just move on. - try: - geoms = poly.intersection(boundarypolyxy) - except: - continue - # iterate over geometries in intersection. - for psub in geoms: - b = psub.boundary - # if projection == 'ortho', - # transform polygon from stereographic - # to orthographic coordinates. - if self.projection == 'ortho': - # if coastline polygon covers more than 99% - # of map region for fulldisk projection, - # it's probably bogus, so skip it. - areafrac = psub.area()/boundarypolyxy.area() - if name == 'gshhs' and\ - self._fulldisk and\ - areafrac > 0.99: continue - # inverse transform from stereographic - # to lat/lon. - b[:,0], b[:,1] = maptran(b[:,0], b[:,1], inverse=True) - # orthographic. - b[:,0], b[:,1]= self(b[:,0], b[:,1]) - polygons.append(zip(b[:,0],b[:,1])) - polygon_types.append(type) - return polygons, polygon_types - - def _getmapboundary(self): - """ - create map boundary polygon (in lat/lon and x/y coordinates) - """ - nx = 100 - ny = 100 - maptran = self - 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 - 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) - # compute proj instance for full disk, if necessary. - if not self._fulldisk: - projparms = self.projparams.copy() - del projparms['x_0'] - del projparms['y_0'] - if self.projection == 'ortho': - llcrnrx = -self.rmajor - llcrnry = -self.rmajor - urcrnrx = -llcrnrx - urcrnry = -llcrnry - else: - llcrnrx = -self._width - llcrnry = -self._height - urcrnrx = -llcrnrx - urcrnry = -llcrnry - projparms['x_0']=-llcrnrx - projparms['y_0']=-llcrnry - maptran = pyproj.Proj(projparms) - elif self.projection in ['moll','robin','sinu']: - # quasi-elliptical region. - lon_0 = self.projparams['lon_0'] - # left side - lats1 = np.linspace(-89.9,89.9,ny).tolist() - lons1 = len(lats1)*[lon_0-179.9] - # top. - lons2 = np.linspace(lon_0-179.9,lon_0+179.9,nx).tolist() - lats2 = len(lons2)*[89.9] - # right side - lats3 = np.linspace(89.9,-89.9,ny).tolist() - lons3 = len(lats3)*[lon_0+179.9] - # bottom. - lons4 = np.linspace(lon_0+179.9,lon_0-179.9,nx).tolist() - lats4 = len(lons4)*[-89.9] - lons = np.array(lons1+lons2+lons3+lons4,np.float64) - lats = np.array(lats1+lats2+lats3+lats4,np.float64) - x, y = maptran(lons,lats) - b = np.empty((len(x),2),np.float64) - b[:,0]=x; b[:,1]=y - boundaryxy = _geoslib.Polygon(b) - else: # all other projections are rectangular. - # left side (x = xmin, ymin <= y <= ymax) - yy = np.linspace(self.ymin, self.ymax, ny)[:-1] - x = len(yy)*[self.xmin]; y = yy.tolist() - # top (y = ymax, xmin <= x <= xmax) - xx = np.linspace(self.xmin, self.xmax, nx)[:-1] - x = x + xx.tolist() - y = y + len(xx)*[self.ymax] - # right side (x = xmax, ymin <= y <= ymax) - yy = np.linspace(self.ymax, self.ymin, ny)[:-1] - x = x + len(yy)*[self.xmax]; y = y + yy.tolist() - # bottom (y = ymin, xmin <= x <= xmax) - xx = np.linspace(self.xmax, self.xmin, nx)[:-1] - x = x + xx.tolist() - y = y + len(xx)*[self.ymin] - x = np.array(x,np.float64) - y = np.array(y,np.float64) - b = np.empty((4,2),np.float64) - b[:,0]=[self.xmin,self.xmin,self.xmax,self.xmax] - b[:,1]=[self.ymin,self.ymax,self.ymax,self.ymin] - boundaryxy = _geoslib.Polygon(b) - if self.projection in ['mill','merc','cyl']: - # make sure map boundary doesn't quite include pole. - if self.urcrnrlat > 89.9999: - urcrnrlat = 89.9999 - else: - urcrnrlat = self.urcrnrlat - if self.llcrnrlat < -89.9999: - llcrnrlat = -89.9999 - else: - llcrnrlat = self.llcrnrlat - lons = [self.llcrnrlon, self.llcrnrlon, self.urcrnrlon, self.urcrnrlon] - lats = [llcrnrlat, urcrnrlat, urcrnrlat, llcrnrlat] - x, y = self(lons, lats) - b = np.empty((len(x),2),np.float64) - b[:,0]=x; b[:,1]=y - boundaryxy = _geoslib.Polygon(b) - else: - if self.projection not in ['moll','robin','sinu']: - lons, lats = maptran(x,y,inverse=True) - # fix lons so there are no jumps. - n = 1 - lonprev = lons[0] - for lon,lat in zip(lons[1:],lats[1:]): - if np.abs(lon-lonprev) > 90.: - if lonprev < 0: - lon = lon - 360. - else: - lon = lon + 360 - lons[n] = lon - lonprev = lon - n = n + 1 - b = np.empty((len(lons),2),np.float64) - b[:,0]=lons; b[:,1]=lats - boundaryll = _geoslib.Polygon(b) - return boundaryll, boundaryxy - - - def drawmapboundary(self,color='k',linewidth=1.0,fill_color=None,\ - zorder=None,ax=None): - """ - draw boundary around map projection region, optionally - filling interior of region. - - linewidth - line width for boundary (default 1.) - color - color of boundary line (default black) - fill_color - fill the map region background with this - color (default is no fill or fill with axis background color). - zorder - sets the zorder for filling map background - (default 0). - ax - axes instance to use (default None, use default axes - instance). - returns PatchCollection representing map boundary. - """ - # get current axes instance (if none specified). - if ax is None and self.ax is None: - try: - ax = plt.gca() - except: - import matplotlib.pyplot as plt - ax = plt.gca() - 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': - limb = Ellipse((self._width,self._height),2.*self._width,2.*self._height) - if self.projection in ['ortho','geos'] and self._fulldisk: - # elliptical region. - ax.add_patch(limb) - if fill_color is None: - limb.set_fill(False) - else: - limb.set_facecolor(fill_color) - limb.set_zorder(0) - limb.set_edgecolor(color) - limb.set_linewidth(linewidth) - limb.set_clip_on(False) - if zorder is not None: - limb.set_zorder(zorder) - elif self.projection in ['moll','robin','sinu']: # elliptical region. - nx = 100; ny = 100 - # quasi-elliptical region. - lon_0 = self.projparams['lon_0'] - # left side - lats1 = np.linspace(-89.9,89.9,ny).tolist() - lons1 = len(lats1)*[lon_0-179.9] - # top. - lons2 = np.linspace(lon_0-179.9,lon_0+179.9,nx).tolist() - lats2 = len(lons2)*[89.9] - # right side - lats3 = np.linspace(89.9,-89.9,ny).tolist() - lons3 = len(lats3)*[lon_0+179.9] - # bottom. - lons4 = np.linspace(lon_0+179.9,lon_0-179.9,nx).tolist() - lats4 = len(lons4)*[-89.9] - lons = np.array(lons1+lons2+lons3+lons4,np.float64) - lats = np.array(lats1+lats2+lats3+lats4,np.float64) - x, y = self(lons,lats) - xy = zip(x,y) - limb = Polygon(xy,edgecolor=color,linewidth=linewidth) - ax.add_patch(limb) - if fill_color is None: - limb.set_fill(False) - else: - limb.set_facecolor(fill_color) - limb.set_zorder(0) - limb.set_clip_on(False) - if zorder is not None: - limb.set_zorder(zorder) - else: # all other projections are rectangular. - ax.axesPatch.set_linewidth(linewidth) - if self.projection not in ['geos','ortho']: - if fill_color is None: - ax.axesPatch.set_facecolor(ax.get_axis_bgcolor()) - else: - ax.axesPatch.set_facecolor(fill_color) - ax.axesPatch.set_zorder(0) - ax.axesPatch.set_edgecolor(color) - ax.set_frame_on(True) - if zorder is not None: - ax.axesPatch.set_zorder(zorder) - else: - ax.axesPatch.set_facecolor(ax.get_axis_bgcolor()) - ax.axesPatch.set_edgecolor(color) - ax.set_frame_on(True) - if zorder is not None: - ax.axesPatch.set_zorder(zorder) - # for geos or ortho projections, also - # draw and fill map projection limb, clipped - # to rectangular region. - ax.add_patch(limb) - if fill_color is None: - limb.set_fill(False) - else: - limb.set_facecolor(fill_color) - limb.set_zorder(0) - limb.set_edgecolor(color) - limb.set_linewidth(linewidth) - if zorder is not None: - limb.set_zorder(zorder) - limb.set_clip_on(True) - # set axes limits to fit map region. - self.set_axes_limits(ax=ax) - return limb - - def fillcontinents(self,color='0.8',lake_color=None,ax=None,zorder=None): - """ - Fill continents. - - color - color to fill continents (default gray). - lake_color - color to fill inland lakes (default axes background). - ax - axes instance (overrides default axes instance). - zorder - sets the zorder for the continent polygons (if not specified, - uses default zorder for a Polygon patch). Set to zero if you want to paint - over the filled continents). - - After filling continents, lakes are re-filled with axis background color. - - returns Polygon object. - """ - if self.resolution is None: - raise AttributeError, 'there are no boundary datasets associated with this Basemap instance' - # get current axes instance (if none specified). - if ax is None and self.ax is None: - try: - ax = plt.gca() - except: - import matplotlib.pyplot as plt - ax = plt.gca() - elif ax is None and self.ax is not None: - ax = self.ax - # get axis background color. - axisbgc = ax.get_axis_bgcolor() - npoly = 0 - for x,y in self.coastpolygons: - xa = np.array(x,np.float32) - ya = np.array(y,np.float32) - # check to see if all four corners of domain in polygon (if so, - # don't draw since it will just fill in the whole map). - delx = 10; dely = 10 - if self.projection in ['cyl']: - delx = 0.1 - dely = 0.1 - test1 = np.fabs(xa-self.urcrnrx) < delx - test2 = np.fabs(xa-self.llcrnrx) < delx - test3 = np.fabs(ya-self.urcrnry) < dely - test4 = np.fabs(ya-self.llcrnry) < dely - hasp1 = np.sum(test1*test3) - hasp2 = np.sum(test2*test3) - hasp4 = np.sum(test2*test4) - hasp3 = np.sum(test1*test4) - if not hasp1 or not hasp2 or not ... [truncated message content]