From c49a29fc125cdd0169d06002431dae8fa4b61726 Mon Sep 17 00:00:00 2001 From: Jeffrey Whitaker Date: Mon, 16 Jan 2012 14:39:39 -0700 Subject: [PATCH 01/15] add round.py example (motivated by a question on the mailing list). --- MANIFEST.in | 1 + examples/README | 3 +++ examples/round.py | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 50 insertions(+) create mode 100644 examples/round.py diff --git a/MANIFEST.in b/MANIFEST.in index 57560b6e9..e6e8b0e00 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -11,6 +11,7 @@ include Changelog include setup.py include nad2bin.c include src/* +include examples/round.py include examples/allskymap.py include examples/allskymap_cr_example.py include examples/plothighsandlows.py diff --git a/examples/README b/examples/README index 3a2af55ce..c79cdf4d1 100644 --- a/examples/README +++ b/examples/README @@ -144,3 +144,6 @@ ploticos.py demonstrates plotting on unstructured grids. lic_demo.py shows how to use vectorplot scikit to visualize vector fields with Line Integral Convolutions (LIC). + +round.py shows how to make a polar plot round (by using a clip path to mask everything +outside the bounding latitude). diff --git a/examples/round.py b/examples/round.py new file mode 100644 index 000000000..9b9c0f887 --- /dev/null +++ b/examples/round.py @@ -0,0 +1,46 @@ +# hack to draw a round polar stereographic plot. +from mpl_toolkits.basemap import Basemap, cm +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.patches as patches + +# read in topo data (on a regular lat/lon grid) +# longitudes go from 20 to 380. +etopo = np.loadtxt('etopo20data.gz') +lons = np.loadtxt('etopo20lons.gz') +lats = np.loadtxt('etopo20lats.gz') + +print 'min/max etopo20 data:' +print etopo.min(),etopo.max() + +# create projection. +m = Basemap(boundinglat=20,lon_0=270,projection='npstere') +# 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,np.linspace(-7500,4500,41),cmap=cm.GMT_haxby,extend='both') +# colorbar +cb = m.colorbar() +# draw coastlines. +coasts = m.drawcoastlines() +# draw parallels and meridians. +parallels = m.drawparallels(np.arange(20.,90,20.)) +merids = m.drawmeridians(np.arange(0.,360.,60.)) +plt.box(on=False) # don't draw axes frame. +# create clip path. +clipit = patches.Circle((0.5*m.xmax,0.5*m.ymax),radius=0.5*m.xmax,fc='none') +ax = plt.gca() +ax.add_patch(clipit) +# clip coastlines. +coasts.set_clip_path(clipit) +# clip contours. +for cntr in cs.collections: + cntr.set_clip_path(clipit) +# clip meridian lines. +for merid in merids: + lines,labels = merids[merid] + for l in lines: + l.set_clip_path(clipit) +plt.title('a round polar plot') +plt.show() From 16787132d4f48196995f1ccb2fb1b996349b466e Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Mon, 16 Jan 2012 18:47:10 -0700 Subject: [PATCH 02/15] update --- examples/round.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/examples/round.py b/examples/round.py index 9b9c0f887..64e831e02 100644 --- a/examples/round.py +++ b/examples/round.py @@ -14,7 +14,9 @@ print etopo.min(),etopo.max() # create projection. -m = Basemap(boundinglat=20,lon_0=270,projection='npstere') +#m = Basemap(boundinglat=-20,lon_0=90,projection='spstere') +m = Basemap(boundinglat=20,lon_0=270,projection='npstere) +print m.xmin, m.xmax,m.ymin,m.ymax # compute native map projection coordinates for lat/lon grid. x,y = m(*np.meshgrid(lons,lats)) # make filled contour plot. @@ -29,9 +31,11 @@ merids = m.drawmeridians(np.arange(0.,360.,60.)) plt.box(on=False) # don't draw axes frame. # create clip path. -clipit = patches.Circle((0.5*m.xmax,0.5*m.ymax),radius=0.5*m.xmax,fc='none') +clipit =\ +patches.Circle((0.5*(m.xmax+m.xmin),0.5*(m.ymax+m.ymin)),radius=0.5*(m.xmax-m.xmin),fc='none') ax = plt.gca() -ax.add_patch(clipit) +p = ax.add_patch(clipit) +p.set_clip_on(False) # clip coastlines. coasts.set_clip_path(clipit) # clip contours. From 91e26d90751e01a745376f806d20f4c4f600a967 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Mon, 16 Jan 2012 19:40:57 -0700 Subject: [PATCH 03/15] add 'round' keyword to Basemap.__init__ for pole-centered projections. --- Changelog | 3 + lib/mpl_toolkits/basemap/__init__.py | 169 ++++++++++++++++++++++++++- 2 files changed, 170 insertions(+), 2 deletions(-) diff --git a/Changelog b/Changelog index 08bf71f07..eb3fd40e8 100644 --- a/Changelog +++ b/Changelog @@ -1,4 +1,7 @@ version 1.0.3 (not yet released) + * add 'round' keyword to Basemap.__init__ for pole-centered + projections to make them round (clipped at boundinglat) instead + of square. * fix broken daynight terminator function. * added kav7 (Kavrayskiy VII) and eck4 (Eckert IV) projections (https://github.com/matplotlib/basemap/issues/9). diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index c581503ac..6a62579ba 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -313,6 +313,10 @@ class methods such as drawcoastlines will raise an 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. + round cut off pole-centered projection at boundinglat + (so plot is a circle instead of a square). Only + relevant for npstere,spstere,nplaea,splaea,npaeqd + or spaeqd projections. satellite_height height of satellite (in m) above equator - only relevant for geostationary and near-sided perspective (``geos`` or ``nsper``) @@ -442,6 +446,7 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, fix_aspect=True, anchor='C', celestial=False, + round=False, ax=None): # docstring is added after __init__ method definition @@ -866,6 +871,14 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, if type == 2: self.lakepolygons.append(_geoslib.Polygon(b)) #if type == 3: self.islandinlakepolygons.append(_geoslib.Polygon(b)) #if type == 4: self.lakeinislandinlakepolygons.append(_geoslib.Polygon(b)) + # set clipping path for round polar plots. + self.round = False + if (self.projection.startswith('np') or + self.projection.startswith('sp')) and round: + self.clipcircle =\ + Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)), + radius=0.5*(self.xmax-self.xmin),fc='none') + self.round = True # set __init__'s docstring __init__.__doc__ = _Basemap_init_doc @@ -1343,6 +1356,17 @@ def drawmapboundary(self,color='k',linewidth=1.0,fill_color=None,\ limb.set_clip_on(False) if zorder is not None: limb.set_zorder(zorder) + elif self.round: + limb = self.clipcircle + 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. # use axesPatch for fill_color, frame for border line props. try: @@ -1468,6 +1492,13 @@ def fillcontinents(self,color='0.8',lake_color=None,ax=None,zorder=None): npoly = npoly + 1 # set axes limits to fit map region. self.set_axes_limits(ax=ax) + # clip continent polygons for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + for poly in polys: + poly.set_clip_path(self.clipcircle) return polys def drawcoastlines(self,linewidth=1.,color='k',antialiased=1,ax=None,zorder=None): @@ -1500,6 +1531,12 @@ def drawcoastlines(self,linewidth=1.,color='k',antialiased=1,ax=None,zorder=None coastlines.set_label('_nolabel_') if zorder is not None: coastlines.set_zorder(zorder) + # clip coastlines for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + coastlines.set_clip_path(self.clipcircle) ax.add_collection(coastlines) # set axes limits to fit map region. self.set_axes_limits(ax=ax) @@ -1541,6 +1578,12 @@ def drawcountries(self,linewidth=0.5,color='k',antialiased=1,ax=None,zorder=None if zorder is not None: countries.set_zorder(zorder) ax.add_collection(countries) + # clip countries for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + countries.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return countries @@ -1581,6 +1624,12 @@ def drawstates(self,linewidth=0.5,color='k',antialiased=1,ax=None,zorder=None): if zorder is not None: states.set_zorder(zorder) ax.add_collection(states) + # clip states for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + states.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return states @@ -1621,6 +1670,12 @@ def drawrivers(self,linewidth=0.5,color='k',antialiased=1,ax=None,zorder=None): if zorder is not None: rivers.set_zorder(zorder) ax.add_collection(rivers) + # clip rivers for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + rivers.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return rivers @@ -1783,6 +1838,12 @@ def readshapefile(self,shapefile,name,drawbounds=True,zorder=None, if zorder is not None: lines.set_zorder(zorder) ax.add_collection(lines) + # clip boundaries for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + lines.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) info = info + (lines,) @@ -2088,7 +2149,17 @@ def drawparallels(self,circles,color='k',linewidth=1.,zorder=None, \ else: linecolls[k] = _tup(linecolls[k]) # override __delitem__ in dict to call remove() on values. - return _dict(linecolls) + pardict = _dict(linecolls) + # clip meridian lines. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + for merid in pardict: + lines,labels = pardict[merid] + for l in lines: + l.set_clip_path(self.clipcircle) + return pardict def drawmeridians(self,meridians,color='k',linewidth=1., zorder=None,\ dashes=[1,1],labels=[0,0,0,0],labelstyle=None,\ @@ -2357,7 +2428,17 @@ def addlon(meridians,madd): # add a remove method to each tuple. linecolls[k] = _tup(linecolls[k]) # override __delitem__ in dict to call remove() on values. - return _dict(linecolls) + meridict = _dict(linecolls) + # clip meridian lines. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + for merid in meridict: + lines,labels = meridict[merid] + for l in lines: + l.set_clip_path(self.clipcircle) + return meridict def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs): """ @@ -2394,6 +2475,12 @@ def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs): seg.append((x,y)) poly = Polygon(seg,**kwargs) ax.add_patch(poly) + # clip polygons for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + poly.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return poly @@ -2723,6 +2810,9 @@ def set_axes_limits(self,ax=None): ax.set_frame_on(False) if self.projection in ['ortho','geos','nsper','aeqd'] and self._fulldisk: ax.set_frame_on(False) + # for round polar plots, turn off frame. + if self.round: + ax.set_frame_on(False) # make sure aspect ratio of map preserved. # plot is re-centered in bounding rectangle. # (anchor instance var determines where plot is placed) @@ -2763,6 +2853,12 @@ def scatter(self, *args, **kwargs): # reset current active image (only if pyplot is imported). if plt: plt.sci(ret) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + ret.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return ret @@ -2788,6 +2884,12 @@ def plot(self, *args, **kwargs): ax.hold(b) raise ax.hold(b) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + ret.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return ret @@ -2825,6 +2927,12 @@ def imshow(self, *args, **kwargs): # reset current active image (only if pyplot is imported). if plt: plt.sci(ret) + # clip image for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + ret.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return ret @@ -2891,6 +2999,12 @@ def pcolor(self,x,y,data,tri=False,**kwargs): # reset current active image (only if pyplot is imported). if plt: plt.sci(ret) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + ret.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return ret @@ -2919,6 +3033,12 @@ def pcolormesh(self,x,y,data,**kwargs): # reset current active image (only if pyplot is imported). if plt: plt.sci(ret) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + ret.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return ret @@ -3001,6 +3121,13 @@ def contour(self,x,y,data,*args,**kwargs): # reset current active image (only if pyplot is imported). if plt and CS.get_array() is not None: plt.sci(CS) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + for cntr in CS.collections: + cntr.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return CS @@ -3096,6 +3223,13 @@ def contourf(self,x,y,data,*args,**kwargs): plt.sci(CS) # set axes limits to fit map region. self.set_axes_limits(ax=ax) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + for cntr in CS.collections: + cntr.set_clip_path(self.clipcircle) return CS def quiver(self, x, y, u, v, *args, **kwargs): @@ -3121,6 +3255,12 @@ def quiver(self, x, y, u, v, *args, **kwargs): ax.hold(b) if plt is not None and ret.get_array() is not None: plt.sci(ret) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + ret.set_clip_path(self.clipcircle) # set axes limits to fit map region. self.set_axes_limits(ax=ax) return ret @@ -3166,6 +3306,12 @@ def barbs(self, x, y, u, v, *args, **kwargs): #if plt is not None and ret.get_array() is not None: # plt.sci(retnh) # set axes limits to fit map region. + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + ret.set_clip_path(self.clipcircle) self.set_axes_limits(ax=ax) return retnh,retsh @@ -3304,6 +3450,12 @@ def drawlsmask(self,land_color="0.8",ocean_color="w",lsmask=None, rgba[:,:,3] = np.where(mask==255,0,rgba[:,:,3]) # plot mask as rgba image. im = self.imshow(rgba,interpolation='nearest',ax=ax,**kwargs) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + im.set_clip_path(self.clipcircle) return im def bluemarble(self,ax=None,scale=None,**kwargs): @@ -3515,6 +3667,12 @@ def warpimage(self,image="bluemarble",scale=None,**kwargs): else: # bmproj True, no interpolation necessary. im = self.imshow(self._bm_rgba,ax=ax,**kwargs) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + im.set_clip_path(self.clipcircle) return im def drawmapscale(self,lon,lat,lon0,lat0,length,barstyle='simple',\ @@ -3789,6 +3947,13 @@ def nightshade(self,date,color="k",delta=0.25,alpha=0.5,ax=None,zorder=2): # is on top. for c in CS.collections: c.set_zorder(zorder) + # clip for round polar plots. + if self.round: + if self.clipcircle not in ax.patches: + p = ax.add_patch(self.clipcircle) + p.set_clip_on(False) + for cntr in CS.collections: + cntr.set_clip_path(self.clipcircle) return CS def _check_ax(self): From 7486bce72dc4b238c30a0dd093f0c22fc63faf22 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 17 Jan 2012 10:08:05 -0700 Subject: [PATCH 04/15] remove round.py example, modify polarmaps.py to produce 'round' maps. --- MANIFEST.in | 1 - examples/README | 3 -- examples/polarmaps.py | 6 +-- examples/round.py | 94 ++++++++++++++++++++++++++----------------- 4 files changed, 60 insertions(+), 44 deletions(-) diff --git a/MANIFEST.in b/MANIFEST.in index e6e8b0e00..57560b6e9 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -11,7 +11,6 @@ include Changelog include setup.py include nad2bin.c include src/* -include examples/round.py include examples/allskymap.py include examples/allskymap_cr_example.py include examples/plothighsandlows.py diff --git a/examples/README b/examples/README index c79cdf4d1..3a2af55ce 100644 --- a/examples/README +++ b/examples/README @@ -144,6 +144,3 @@ ploticos.py demonstrates plotting on unstructured grids. lic_demo.py shows how to use vectorplot scikit to visualize vector fields with Line Integral Convolutions (LIC). - -round.py shows how to make a polar plot round (by using a clip path to mask everything -outside the bounding latitude). diff --git a/examples/polarmaps.py b/examples/polarmaps.py index 1050f1157..47cb045bb 100644 --- a/examples/polarmaps.py +++ b/examples/polarmaps.py @@ -29,12 +29,12 @@ lon_0 = 130. lon_0_ortho = lon_0 - 180. lat_0 = -90. - bounding_lat = -20. + bounding_lat = -1. elif hem == 'North': lon_0 = -90. lon_0_ortho = lon_0 lat_0 = 90. - bounding_lat = 20. + bounding_lat = 1. # loop over projections, one for each panel of the figure. fig = plt.figure(figsize=(8,8)) npanel = 0 @@ -51,7 +51,7 @@ resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho) else: m = Basemap(boundinglat=bounding_lat,lon_0=lon_0,\ - resolution='c',area_thresh=10000.,projection=projection) + resolution='c',area_thresh=10000.,projection=projection,round=True) # compute native map projection coordinates for lat/lon grid. x,y = m(*np.meshgrid(lons,lats)) ax = fig.add_subplot(2,2,npanel) diff --git a/examples/round.py b/examples/round.py index 64e831e02..47cb045bb 100644 --- a/examples/round.py +++ b/examples/round.py @@ -1,8 +1,13 @@ -# hack to draw a round polar stereographic plot. -from mpl_toolkits.basemap import Basemap, cm +# make plots of etopo bathymetry/topography data on +# various map projections, drawing coastlines, state and +# country boundaries, filling continents and drawing +# parallels/meridians + +# illustrates special-case polar-centric projections. + +from mpl_toolkits.basemap import Basemap import numpy as np import matplotlib.pyplot as plt -import matplotlib.patches as patches # read in topo data (on a regular lat/lon grid) # longitudes go from 20 to 380. @@ -13,38 +18,53 @@ print 'min/max etopo20 data:' print etopo.min(),etopo.max() -# create projection. -#m = Basemap(boundinglat=-20,lon_0=90,projection='spstere') -m = Basemap(boundinglat=20,lon_0=270,projection='npstere) -print m.xmin, m.xmax,m.ymin,m.ymax -# 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,np.linspace(-7500,4500,41),cmap=cm.GMT_haxby,extend='both') -# colorbar -cb = m.colorbar() -# draw coastlines. -coasts = m.drawcoastlines() -# draw parallels and meridians. -parallels = m.drawparallels(np.arange(20.,90,20.)) -merids = m.drawmeridians(np.arange(0.,360.,60.)) -plt.box(on=False) # don't draw axes frame. -# create clip path. -clipit =\ -patches.Circle((0.5*(m.xmax+m.xmin),0.5*(m.ymax+m.ymin)),radius=0.5*(m.xmax-m.xmin),fc='none') -ax = plt.gca() -p = ax.add_patch(clipit) -p.set_clip_on(False) -# clip coastlines. -coasts.set_clip_path(clipit) -# clip contours. -for cntr in cs.collections: - cntr.set_clip_path(clipit) -# clip meridian lines. -for merid in merids: - lines,labels = merids[merid] - for l in lines: - l.set_clip_path(clipit) -plt.title('a round polar plot') +# these are the 4 polar projections +projs = ['laea','stere','aeqd','ortho'] # short names +# long names +projnames = ['Lambert Azimuthal Equal Area','Stereographic','Azimuthal Equidistant','Orthographic'] +# loop over hemispheres, make a 4-panel plot for each hemisphere +# showing all four polar projections. +for hem in ['North','South']: + if hem == 'South': + lon_0 = 130. + lon_0_ortho = lon_0 - 180. + lat_0 = -90. + bounding_lat = -1. + elif hem == 'North': + lon_0 = -90. + lon_0_ortho = lon_0 + lat_0 = 90. + bounding_lat = 1. + # loop over projections, one for each panel of the figure. + fig = plt.figure(figsize=(8,8)) + npanel = 0 + for proj,projname in zip(projs,projnames): + npanel = npanel + 1 + if hem == 'South': + projection = 'sp'+proj + elif hem == 'North': + projection = 'np'+proj + # setup map projection + # centered on Australia (for SH) or US (for NH). + if proj == 'ortho': + m = Basemap(projection='ortho', + resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho) + else: + m = Basemap(boundinglat=bounding_lat,lon_0=lon_0,\ + resolution='c',area_thresh=10000.,projection=projection,round=True) + # compute native map projection coordinates for lat/lon grid. + x,y = m(*np.meshgrid(lons,lats)) + ax = fig.add_subplot(2,2,npanel) + # make filled contour plot. + cs = m.contourf(x,y,etopo,20,cmap=plt.cm.jet) + # draw coastlines. + m.drawcoastlines() + # draw parallels and meridians. + m.drawparallels(np.arange(-80.,90,20.)) + m.drawmeridians(np.arange(0.,360.,60.)) + # draw boundary around map region. + m.drawmapboundary() + # draw title. + plt.title(hem+' Polar '+projname,y=1.05,fontsize=12) + print 'plotting '+hem+' Polar '+projname+' basemap ...' plt.show() From 9cc6607162019c74d8f710191486871cb7842535 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 17 Jan 2012 10:44:49 -0700 Subject: [PATCH 05/15] update docstring for round option. --- lib/mpl_toolkits/basemap/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 6a62579ba..36b9cee09 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -316,7 +316,7 @@ class methods such as drawcoastlines will raise an round cut off pole-centered projection at boundinglat (so plot is a circle instead of a square). Only relevant for npstere,spstere,nplaea,splaea,npaeqd - or spaeqd projections. + or spaeqd projections. Default False. satellite_height height of satellite (in m) above equator - only relevant for geostationary and near-sided perspective (``geos`` or ``nsper``) From a6c167ddd139a18126df4d33f4ef8b5f844a114d Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 17 Jan 2012 10:51:08 -0700 Subject: [PATCH 06/15] set boundinglat instance var. --- lib/mpl_toolkits/basemap/__init__.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 36b9cee09..fc888ca48 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -459,6 +459,8 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, self.celestial = celestial # map projection. self.projection = projection + # bounding lat (for pole-centered plots) + self.boundinglat = boundinglat # set up projection parameter dict. projparams = {} @@ -581,6 +583,9 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, 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 lat_0 == 90 or lat_0 == -90: + # for ortho plot centered on pole, set boundinglat to equator. + self.boundinglat = 0. if width is not None or height is not None: sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection]) if not using_corners: From 635e1d497fff8bb71fcc643f149cea338c711c95 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 17 Jan 2012 10:54:59 -0700 Subject: [PATCH 07/15] change boundinglat --- examples/polarmaps.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/polarmaps.py b/examples/polarmaps.py index 47cb045bb..10089b392 100644 --- a/examples/polarmaps.py +++ b/examples/polarmaps.py @@ -29,12 +29,14 @@ lon_0 = 130. lon_0_ortho = lon_0 - 180. lat_0 = -90. - bounding_lat = -1. + # Lambert Azimuth bounding lat must not extend into opposite hem. + bounding_lat = -0.01 elif hem == 'North': lon_0 = -90. lon_0_ortho = lon_0 lat_0 = 90. - bounding_lat = 1. + # Lambert Azimuth bounding lat must not extend into opposite hem. + bounding_lat = 0.01 # loop over projections, one for each panel of the figure. fig = plt.figure(figsize=(8,8)) npanel = 0 From d940cb65a306a0d6c275f25bfffc49a3b758ec37 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 17 Jan 2012 10:58:11 -0700 Subject: [PATCH 08/15] change colormap --- examples/polarmaps.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/polarmaps.py b/examples/polarmaps.py index 10089b392..ff5ed48bf 100644 --- a/examples/polarmaps.py +++ b/examples/polarmaps.py @@ -5,7 +5,7 @@ # illustrates special-case polar-centric projections. -from mpl_toolkits.basemap import Basemap +from mpl_toolkits.basemap import Basemap, cm import numpy as np import matplotlib.pyplot as plt @@ -58,7 +58,7 @@ x,y = m(*np.meshgrid(lons,lats)) ax = fig.add_subplot(2,2,npanel) # make filled contour plot. - cs = m.contourf(x,y,etopo,20,cmap=plt.cm.jet) + cs = m.contourf(x,y,etopo,np.linspace(-7500,4500,41),cmap=cm.GMT_haxby) # draw coastlines. m.drawcoastlines() # draw parallels and meridians. From befda9269f307aed044ddcdef27843d42070b9d4 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 17 Jan 2012 14:47:26 -0700 Subject: [PATCH 09/15] update docstring --- lib/mpl_toolkits/basemap/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index fc888ca48..bb904322d 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -585,7 +585,8 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, raise ValueError('must specify lat_0 and lon_0 for Orthographic basemap') if lat_0 == 90 or lat_0 == -90: # for ortho plot centered on pole, set boundinglat to equator. - self.boundinglat = 0. + # (so meridian labels can be drawn in this special case). + self.boundinglat = 0 if width is not None or height is not None: sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection]) if not using_corners: From 2824fa0e8dd09393091d5c949558105442fe8a32 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Wed, 18 Jan 2012 06:47:06 -0700 Subject: [PATCH 10/15] add meridians labels for pole-centered plots. --- examples/polarmaps.py | 8 +- lib/mpl_toolkits/basemap/__init__.py | 213 ++++++++++++++++++--------- 2 files changed, 145 insertions(+), 76 deletions(-) diff --git a/examples/polarmaps.py b/examples/polarmaps.py index ff5ed48bf..1310374d1 100644 --- a/examples/polarmaps.py +++ b/examples/polarmaps.py @@ -26,13 +26,13 @@ # showing all four polar projections. for hem in ['North','South']: if hem == 'South': - lon_0 = 130. + lon_0 = -130. lon_0_ortho = lon_0 - 180. lat_0 = -90. # Lambert Azimuth bounding lat must not extend into opposite hem. - bounding_lat = -0.01 + bounding_lat = -0.01 elif hem == 'North': - lon_0 = -90. + lon_0 = -130. lon_0_ortho = lon_0 lat_0 = 90. # Lambert Azimuth bounding lat must not extend into opposite hem. @@ -63,7 +63,7 @@ m.drawcoastlines() # draw parallels and meridians. m.drawparallels(np.arange(-80.,90,20.)) - m.drawmeridians(np.arange(0.,360.,60.)) + m.drawmeridians(np.arange(0.,340.,30.),labels=[1,0,0,0],fontsize=7) # draw boundary around map region. m.drawmapboundary() # draw title. diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index bb904322d..33082f437 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -461,6 +461,8 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, self.projection = projection # bounding lat (for pole-centered plots) self.boundinglat = boundinglat + # is a round pole-centered plot desired? + self.round = round # set up projection parameter dict. projparams = {} @@ -587,6 +589,7 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, # for ortho plot centered on pole, set boundinglat to equator. # (so meridian labels can be drawn in this special case). self.boundinglat = 0 + self.round = True if width is not None or height is not None: sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection]) if not using_corners: @@ -878,9 +881,9 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, #if type == 3: self.islandinlakepolygons.append(_geoslib.Polygon(b)) #if type == 4: self.lakeinislandinlakepolygons.append(_geoslib.Polygon(b)) # set clipping path for round polar plots. - self.round = False if (self.projection.startswith('np') or - self.projection.startswith('sp')) and round: + self.projection.startswith('sp') or + self.projection == 'ortho') and self.round: self.clipcircle =\ Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)), radius=0.5*(self.xmax-self.xmin),fc='none') @@ -2064,39 +2067,7 @@ def drawparallels(self,circles,color='k',linewidth=1.,zorder=None, \ nl = _searchlist(lats,lat) nr = _searchlist(lats[::-1],lat) if nr != -1: nr = len(lons)-nr-1 - try: # fmt is a function that returns a formatted string - latlab = fmt(lat) - except: # fmt is a format string. - if lat<0: - if rcParams['text.usetex']: - if labelstyle=='+/-': - latlabstr = r'${\/-%s\/^{\circ}}$'%fmt - else: - latlabstr = r'${%s\/^{\circ}\/S}$'%fmt - else: - if labelstyle=='+/-': - latlabstr = u'-%s\N{DEGREE SIGN}'%fmt - else: - latlabstr = u'%s\N{DEGREE SIGN}S'%fmt - latlab = latlabstr%np.fabs(lat) - elif lat>0: - if rcParams['text.usetex']: - if labelstyle=='+/-': - latlabstr = r'${\/+%s\/^{\circ}}$'%fmt - else: - latlabstr = r'${%s\/^{\circ}\/N}$'%fmt - else: - if labelstyle=='+/-': - latlabstr = u'+%s\N{DEGREE SIGN}'%fmt - else: - latlabstr = u'%s\N{DEGREE SIGN}N'%fmt - latlab = latlabstr%lat - else: - if rcParams['text.usetex']: - latlabstr = r'${%s\/^{\circ}}$'%fmt - else: - latlabstr = u'%s\N{DEGREE SIGN}'%fmt - latlab = latlabstr%lat + latlab = _setlatlab(fmt,lat,labelstyle) # parallels can intersect each map edge twice. for i,n in enumerate([nl,nr]): # don't bother if close to the first label. @@ -2193,6 +2164,9 @@ def drawmeridians(self,meridians,color='k',linewidth=1., zorder=None,\ example labels=[1,0,0,1] will cause meridians to be labelled where they intersect the left and and bottom of the plot, but not the right and top. + For round pole-centered plots, if element in labels + evaluates to True, then meridians are labelled all around + the bounding latitude. labelstyle if set to "+/-", east and west longitudes are labelled with "+" and "-", otherwise they are labelled with "E" and "W". @@ -2316,7 +2290,7 @@ def addlon(meridians,madd): sys.stdout.write('Warning: Cannot label meridians on %s basemap' % _projnames[self.projection]) labels = [0,0,0,0] if self.projection in ['ortho','geos','nsper','aeqd'] and max(labels): - if self._fulldisk: + if self._fulldisk and self.boundinglat is None: sys.stdout.write(dedent( """'Warning: Cannot label meridians on full-disk Geostationary, Orthographic or Azimuthal equidistant basemap @@ -2331,7 +2305,7 @@ def addlon(meridians,madd): xmin,ymin = self(lon_0-179.9,-90) xmax,ymax = self(lon_0+179.9,90) for dolab,side in zip(labels,['l','r','t','b']): - if not dolab: continue + if not dolab or self.round: continue # for cylindrical projections, don't draw meridians on left or right. if self.projection in _cylproj + _pseudocyl and side in ['l','r']: continue if side in ['l','r']: @@ -2371,39 +2345,7 @@ def addlon(meridians,madd): nl = _searchlist(lons,lon2) nr = _searchlist(lons[::-1],lon2) if nr != -1: nr = len(lons)-nr-1 - try: # fmt is a function that returns a formatted string - lonlab = fmt(lon) - except: # fmt is a format string. - if lon2>180: - if rcParams['text.usetex']: - if labelstyle=='+/-': - lonlabstr = r'${\/-%s\/^{\circ}}$'%fmt - else: - lonlabstr = r'${%s\/^{\circ}\/W}$'%fmt - else: - if labelstyle=='+/-': - lonlabstr = u'-%s\N{DEGREE SIGN}'%fmt - else: - lonlabstr = u'%s\N{DEGREE SIGN}W'%fmt - lonlab = lonlabstr%np.fabs(lon2-360) - elif lon2<180 and lon2 != 0: - if rcParams['text.usetex']: - if labelstyle=='+/-': - lonlabstr = r'${\/+%s\/^{\circ}}$'%fmt - else: - lonlabstr = r'${%s\/^{\circ}\/E}$'%fmt - else: - if labelstyle=='+/-': - lonlabstr = u'+%s\N{DEGREE SIGN}'%fmt - else: - lonlabstr = u'%s\N{DEGREE SIGN}E'%fmt - lonlab = lonlabstr%lon2 - else: - if rcParams['text.usetex']: - lonlabstr = r'${%s\/^{\circ}}$'%fmt - else: - lonlabstr = u'%s\N{DEGREE SIGN}'%fmt - lonlab = lonlabstr%lon2 + lonlab = _setlonlab(fmt,merid,labelstyle) # meridians can intersect each map edge twice. for i,n in enumerate([nl,nr]): lat = lats[n]/100. @@ -2422,7 +2364,7 @@ def addlon(meridians,madd): else: t = ax.text(xx[n],self.urcrnry+yoffset,lonlab,horizontalalignment='center',verticalalignment='bottom',**kwargs) - if t is not None: linecolls[lon][1].append(t) + if t is not None: linecolls[lon2][1].append(t) # set axes limits to fit map region. self.set_axes_limits(ax=ax) # remove empty values from linecolls dictionary @@ -2435,15 +2377,71 @@ def addlon(meridians,madd): linecolls[k] = _tup(linecolls[k]) # override __delitem__ in dict to call remove() on values. meridict = _dict(linecolls) - # clip meridian lines. + # clip meridian lines and label them. if self.round: if self.clipcircle not in ax.patches: p = ax.add_patch(self.clipcircle) p.set_clip_on(False) + # label desired? + label = False + for lab in labels: + if lab: label = True for merid in meridict: lines,labels = meridict[merid] + # clip lines. for l in lines: l.set_clip_path(self.clipcircle) + if not label: continue + # label + lonlab = _setlonlab(fmt,merid,labelstyle) + x,y = self(merid,self.boundinglat) + r = np.sqrt((x-0.5*(self.xmin+self.xmax))**2+ + (y-0.5*(self.ymin+self.ymax))**2) + r = r + np.sqrt(xoffset**2+yoffset**2) + if self.projection.startswith('np'): + pole = 1 + elif self.projection.startswith('sp'): + pole = -1 + elif self.projection == 'ortho' and self.round: + pole = 1 + if pole == 1: + if self.projection != 'ortho': + theta = (np.pi/180.)*(merid-self.projparams['lon_0']-90) + else: + theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90) + x = r*np.cos(theta)+0.5*(self.xmin+self.xmax) + y = r*np.sin(theta)+0.5*(self.ymin+self.ymax) + if x > 0.5*(self.xmin+self.xmax)+xoffset: + horizalign = 'left' + elif x < 0.5*(self.xmin+self.xmax)-xoffset: + horizalign = 'right' + else: + horizalign = 'center' + if y > 0.5*(self.ymin+self.ymax)+yoffset: + vertalign = 'bottom' + elif y < 0.5*(self.ymin+self.ymax)-yoffset: + vertalign = 'top' + else: + vertalign = 'center' + elif pole == -1: + theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90) + x = r*np.cos(theta)+0.5*(self.xmin+self.xmax) + y = r*np.sin(theta)+0.5*(self.ymin+self.ymax) + if x > 0.5*(self.xmin+self.xmax)-xoffset: + horizalign = 'right' + elif x < 0.5*(self.xmin+self.xmax)+xoffset: + horizalign = 'left' + else: + horizalign = 'center' + if y > 0.5*(self.ymin+self.ymax)-yoffset: + vertalign = 'top' + elif y < 0.5*(self.ymin+self.ymax)+yoffset: + vertalign = 'bottom' + else: + vertalign = 'center' + t =\ + ax.text(x,y,lonlab,horizontalalignment=horizalign,verticalalignment=vertalign,**kwargs) + meridict[merid][1].append(t) return meridict def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs): @@ -4337,3 +4335,74 @@ class _dict(dict): def __delitem__(self,key): self[key].remove() super(_dict, self).__delitem__(key) + +def _setlonlab(fmt,lon,labelstyle): + try: # fmt is a function that returns a formatted string + lonlab = fmt(lon) + except: # fmt is a format string. + if lon>180: + if rcParams['text.usetex']: + if labelstyle=='+/-': + lonlabstr = r'${\/-%s\/^{\circ}}$'%fmt + else: + lonlabstr = r'${%s\/^{\circ}\/W}$'%fmt + else: + if labelstyle=='+/-': + lonlabstr = u'-%s\N{DEGREE SIGN}'%fmt + else: + lonlabstr = u'%s\N{DEGREE SIGN}W'%fmt + lonlab = lonlabstr%np.fabs(lon-360) + elif lon<180 and lon != 0: + if rcParams['text.usetex']: + if labelstyle=='+/-': + lonlabstr = r'${\/+%s\/^{\circ}}$'%fmt + else: + lonlabstr = r'${%s\/^{\circ}\/E}$'%fmt + else: + if labelstyle=='+/-': + lonlabstr = u'+%s\N{DEGREE SIGN}'%fmt + else: + lonlabstr = u'%s\N{DEGREE SIGN}E'%fmt + lonlab = lonlabstr%lon + else: + if rcParams['text.usetex']: + lonlabstr = r'${%s\/^{\circ}}$'%fmt + else: + lonlabstr = u'%s\N{DEGREE SIGN}'%fmt + lonlab = lonlabstr%lon + return lonlab + +def _setlatlab(fmt,lat,labelstyle): + try: # fmt is a function that returns a formatted string + latlab = fmt(lat) + except: # fmt is a format string. + if lat<0: + if rcParams['text.usetex']: + if labelstyle=='+/-': + latlabstr = r'${\/-%s\/^{\circ}}$'%fmt + else: + latlabstr = r'${%s\/^{\circ}\/S}$'%fmt + else: + if labelstyle=='+/-': + latlabstr = u'-%s\N{DEGREE SIGN}'%fmt + else: + latlabstr = u'%s\N{DEGREE SIGN}S'%fmt + latlab = latlabstr%np.fabs(lat) + elif lat>0: + if rcParams['text.usetex']: + if labelstyle=='+/-': + latlabstr = r'${\/+%s\/^{\circ}}$'%fmt + else: + latlabstr = r'${%s\/^{\circ}\/N}$'%fmt + else: + if labelstyle=='+/-': + latlabstr = u'+%s\N{DEGREE SIGN}'%fmt + else: + latlabstr = u'%s\N{DEGREE SIGN}N'%fmt + latlab = latlabstr%lat + else: + if rcParams['text.usetex']: + latlabstr = r'${%s\/^{\circ}}$'%fmt + else: + latlabstr = u'%s\N{DEGREE SIGN}'%fmt + latlab = latlabstr%lat From 2b4c05ed7179a00d80ab7f8cc4df04a36fe20746 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Wed, 18 Jan 2012 09:02:41 -0700 Subject: [PATCH 11/15] update pole-centered meridian labelling --- lib/mpl_toolkits/basemap/__init__.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 33082f437..1bf0f5149 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -2164,9 +2164,6 @@ def drawmeridians(self,meridians,color='k',linewidth=1., zorder=None,\ example labels=[1,0,0,1] will cause meridians to be labelled where they intersect the left and and bottom of the plot, but not the right and top. - For round pole-centered plots, if element in labels - evaluates to True, then meridians are labelled all around - the bounding latitude. labelstyle if set to "+/-", east and west longitudes are labelled with "+" and "-", otherwise they are labelled with "E" and "W". @@ -2387,7 +2384,7 @@ def addlon(meridians,madd): for lab in labels: if lab: label = True for merid in meridict: - lines,labels = meridict[merid] + lines,labs = meridict[merid] # clip lines. for l in lines: l.set_clip_path(self.clipcircle) @@ -2423,6 +2420,11 @@ def addlon(meridians,madd): vertalign = 'top' else: vertalign = 'center' + # labels [l,r,t,b] + if labels[0] and x >= 0.5*(self.xmin+self.xmax)+xoffset: continue + if labels[1] and x <= 0.5*(self.xmin+self.xmax)-xoffset: continue + if labels[2] and y <= 0.5*(self.ymin+self.ymax)-yoffset: continue + if labels[3] and y >= 0.5*(self.ymin+self.ymax)+yoffset: continue elif pole == -1: theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90) x = r*np.cos(theta)+0.5*(self.xmin+self.xmax) @@ -2439,6 +2441,11 @@ def addlon(meridians,madd): vertalign = 'bottom' else: vertalign = 'center' + # labels [l,r,t,b] + if labels[0] and x <= 0.5*(self.xmin+self.xmax)+xoffset: continue + if labels[1] and x >= 0.5*(self.xmin+self.xmax)-xoffset: continue + if labels[2] and y >= 0.5*(self.ymin+self.ymax)-yoffset: continue + if labels[3] and y <= 0.5*(self.ymin+self.ymax)+yoffset: continue t =\ ax.text(x,y,lonlab,horizontalalignment=horizalign,verticalalignment=vertalign,**kwargs) meridict[merid][1].append(t) @@ -4337,6 +4344,7 @@ def __delitem__(self,key): super(_dict, self).__delitem__(key) def _setlonlab(fmt,lon,labelstyle): + # set lon label string (called by Basemap.drawmeridians) try: # fmt is a function that returns a formatted string lonlab = fmt(lon) except: # fmt is a format string. @@ -4373,6 +4381,7 @@ def _setlonlab(fmt,lon,labelstyle): return lonlab def _setlatlab(fmt,lat,labelstyle): + # set lat label string (called by Basemap.drawparallels) try: # fmt is a function that returns a formatted string latlab = fmt(lat) except: # fmt is a format string. From 072ba45ce9bdef728ce4a9361cdbdb2880af08c2 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Wed, 18 Jan 2012 09:16:47 -0700 Subject: [PATCH 12/15] more fixes for meridian labelling for round polar plots. --- examples/polarmaps.py | 5 +++-- lib/mpl_toolkits/basemap/__init__.py | 22 +++++++++++----------- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/examples/polarmaps.py b/examples/polarmaps.py index 1310374d1..5eedde276 100644 --- a/examples/polarmaps.py +++ b/examples/polarmaps.py @@ -26,7 +26,7 @@ # showing all four polar projections. for hem in ['North','South']: if hem == 'South': - lon_0 = -130. + lon_0 = 130. lon_0_ortho = lon_0 - 180. lat_0 = -90. # Lambert Azimuth bounding lat must not extend into opposite hem. @@ -63,7 +63,8 @@ m.drawcoastlines() # draw parallels and meridians. m.drawparallels(np.arange(-80.,90,20.)) - m.drawmeridians(np.arange(0.,340.,30.),labels=[1,0,0,0],fontsize=7) + #labels = [l,r,t,b] + m.drawmeridians(np.arange(0.,340.,30.),labels=[1,1,1,1],fontsize=7) # draw boundary around map region. m.drawmapboundary() # draw title. diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 1bf0f5149..663830d99 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -2402,9 +2402,9 @@ def addlon(meridians,madd): elif self.projection == 'ortho' and self.round: pole = 1 if pole == 1: - if self.projection != 'ortho': - theta = (np.pi/180.)*(merid-self.projparams['lon_0']-90) - else: + theta = (np.pi/180.)*(merid-self.projparams['lon_0']-90) + if self.projection == 'ortho' and\ + self.projparams['lat_0'] == -90: theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90) x = r*np.cos(theta)+0.5*(self.xmin+self.xmax) y = r*np.sin(theta)+0.5*(self.ymin+self.ymax) @@ -2421,10 +2421,10 @@ def addlon(meridians,madd): else: vertalign = 'center' # labels [l,r,t,b] - if labels[0] and x >= 0.5*(self.xmin+self.xmax)+xoffset: continue - if labels[1] and x <= 0.5*(self.xmin+self.xmax)-xoffset: continue - if labels[2] and y <= 0.5*(self.ymin+self.ymax)-yoffset: continue - if labels[3] and y >= 0.5*(self.ymin+self.ymax)+yoffset: continue + if labels[0] and not labels[1] and x >= 0.5*(self.xmin+self.xmax)+xoffset: continue + if labels[1] and not labels[0] and x <= 0.5*(self.xmin+self.xmax)-xoffset: continue + if labels[2] and not labels[3] and y <= 0.5*(self.ymin+self.ymax)-yoffset: continue + if labels[3] and not labels[2]and y >= 0.5*(self.ymin+self.ymax)+yoffset: continue elif pole == -1: theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90) x = r*np.cos(theta)+0.5*(self.xmin+self.xmax) @@ -2442,10 +2442,10 @@ def addlon(meridians,madd): else: vertalign = 'center' # labels [l,r,t,b] - if labels[0] and x <= 0.5*(self.xmin+self.xmax)+xoffset: continue - if labels[1] and x >= 0.5*(self.xmin+self.xmax)-xoffset: continue - if labels[2] and y >= 0.5*(self.ymin+self.ymax)-yoffset: continue - if labels[3] and y <= 0.5*(self.ymin+self.ymax)+yoffset: continue + if labels[0] and not labels[1] and x <= 0.5*(self.xmin+self.xmax)+xoffset: continue + if labels[1] and not labels[0] and x >= 0.5*(self.xmin+self.xmax)-xoffset: continue + if labels[2] and not labels[3] and y >= 0.5*(self.ymin+self.ymax)-yoffset: continue + if labels[3] and not labels[2] and y <= 0.5*(self.ymin+self.ymax)+yoffset: continue t =\ ax.text(x,y,lonlab,horizontalalignment=horizalign,verticalalignment=vertalign,**kwargs) meridict[merid][1].append(t) From ef9a980559c9233d34c1301d31df3147a82f6854 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Wed, 18 Jan 2012 09:24:22 -0700 Subject: [PATCH 13/15] do coastline clipping in stereographic coords for polar aeqd, laea. --- examples/polarmaps.py | 4 ++-- lib/mpl_toolkits/basemap/__init__.py | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/polarmaps.py b/examples/polarmaps.py index 5eedde276..213437f51 100644 --- a/examples/polarmaps.py +++ b/examples/polarmaps.py @@ -26,13 +26,13 @@ # showing all four polar projections. for hem in ['North','South']: if hem == 'South': - lon_0 = 130. + lon_0 = -130. lon_0_ortho = lon_0 - 180. lat_0 = -90. # Lambert Azimuth bounding lat must not extend into opposite hem. bounding_lat = -0.01 elif hem == 'North': - lon_0 = -130. + lon_0 = 130. lon_0_ortho = lon_0 lat_0 = 90. # Lambert Azimuth bounding lat must not extend into opposite hem. diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 663830d99..390249eb8 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -981,7 +981,8 @@ def _readboundarydata(self,name): # coordinates, then transform back. This is # because these projections are only defined on a hemisphere, and # some boundary features (like Eurasia) would be undefined otherwise. - if self.projection in ['ortho','gnom','nsper'] and name == 'gshhs': + tostere = ['ortho','gnom','nsper','nplaea','npaeqd','splaea','spaeqd'] + if self.projection in tostere and name == 'gshhs': containsPole = True lon_0=self.projparams['lon_0'] lat_0=self.projparams['lat_0'] @@ -1102,8 +1103,7 @@ def _readboundarydata(self,name): # to map projection coordinates. # special case for ortho/gnom/nsper, compute coastline polygon # vertices in stereographic coords. - if name == 'gshhs' and self.projection in\ - ['ortho','gnom','nsper']: + if name == 'gshhs' and self.projection in tostere: b[:,0], b[:,1] = maptran(b[:,0], b[:,1]) else: b[:,0], b[:,1] = self(b[:,0], b[:,1]) @@ -1120,7 +1120,7 @@ def _readboundarydata(self,name): # for ortho/gnom/nsper projection, all points # outside map projection region are eliminated # with the above step, so we're done. - if self.projection in ['ortho','gnom','nsper']: + if self.projection in tostere: polygons.append(list(zip(bx,by))) polygon_types.append(type) continue @@ -1147,7 +1147,7 @@ def _readboundarydata(self,name): # if projection in ['ortho','gnom','nsper'], # transform polygon from stereographic # to ortho/gnom/nsper coordinates. - if self.projection in ['ortho','gnom','nsper']: + if self.projection in tostere: # if coastline polygon covers more than 99% # of map region for fulldisk projection, # it's probably bogus, so skip it. From 7c478298eba39afdd615eae0ee16bb3fe04d7a05 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Wed, 18 Jan 2012 09:27:56 -0700 Subject: [PATCH 14/15] update --- Changelog | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Changelog b/Changelog index eb3fd40e8..c3d6f0018 100644 --- a/Changelog +++ b/Changelog @@ -1,4 +1,7 @@ version 1.0.3 (not yet released) + * clip coastlines for nplaea,npaeqd,splaea,spaeqd in stereographic + coordinates to avoid S. America disappearing in some south polar + plots. * add 'round' keyword to Basemap.__init__ for pole-centered projections to make them round (clipped at boundinglat) instead of square. From 4f3a30aeefb29501ec73444d6fddfb1d1b5af41d Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Wed, 18 Jan 2012 10:05:49 -0700 Subject: [PATCH 15/15] don't blow away labels list --- lib/mpl_toolkits/basemap/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 390249eb8..c58be0d9a 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -2133,7 +2133,7 @@ def drawparallels(self,circles,color='k',linewidth=1.,zorder=None, \ p = ax.add_patch(self.clipcircle) p.set_clip_on(False) for merid in pardict: - lines,labels = pardict[merid] + lines,labs = pardict[merid] for l in lines: l.set_clip_path(self.clipcircle) return pardict