From 51d44e561883ceaab1c30856a6579f55e3692b0e Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Wed, 5 Sep 2012 11:20:00 -0600 Subject: [PATCH 1/9] add cea projection. --- doc/users/cea.rst | 8 ++++++++ doc/users/figures/cea.py | 17 +++++++++++++++++ doc/users/mapsetup.rst | 1 + 3 files changed, 26 insertions(+) create mode 100644 doc/users/cea.rst create mode 100644 doc/users/figures/cea.py diff --git a/doc/users/cea.rst b/doc/users/cea.rst new file mode 100644 index 000000000..7271c7f64 --- /dev/null +++ b/doc/users/cea.rst @@ -0,0 +1,8 @@ +.. _cea: + +Cylindrial Equal-Area Projection +================================ + +It is what is says. + +.. plot:: users/figures/cea.py diff --git a/doc/users/figures/cea.py b/doc/users/figures/cea.py new file mode 100644 index 000000000..e15dbb775 --- /dev/null +++ b/doc/users/figures/cea.py @@ -0,0 +1,17 @@ +from mpl_toolkits.basemap import Basemap +import numpy as np +import matplotlib.pyplot as plt +# llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon +# are the lat/lon values of the lower left and upper right corners +# of the map. +# resolution = 'c' means use crude resolution coastlines. +m = Basemap(projection='cea',llcrnrlat=-90,urcrnrlat=90,\ + llcrnrlon=-180,urcrnrlon=180,resolution='c') +m.drawcoastlines() +m.fillcontinents(color='coral',lake_color='aqua') +# draw parallels and meridians. +m.drawparallels(np.arange(-90.,91.,30.)) +m.drawmeridians(np.arange(-180.,181.,60.)) +m.drawmapboundary(fill_color='aqua') +plt.title("Cylindrical Equal-Area Projection") +plt.show() diff --git a/doc/users/mapsetup.rst b/doc/users/mapsetup.rst index 18d73a850..8a796cd0a 100644 --- a/doc/users/mapsetup.rst +++ b/doc/users/mapsetup.rst @@ -53,6 +53,7 @@ time, many compromise between the two. poly.rst mill.rst gall.rst + cea.rst lcc.rst laea.rst stere.rst From 66e8bbb1e4c61af4f0647040db69c11616e45412 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Mon, 17 Sep 2012 13:23:47 -0600 Subject: [PATCH 2/9] have drawcoastlines use line segments instead of polygons to avoid 'thickening' of lines around edges of map. --- Changelog | 2 + lib/mpl_toolkits/basemap/__init__.py | 76 ++++++++++++++++------------ 2 files changed, 45 insertions(+), 33 deletions(-) diff --git a/Changelog b/Changelog index ab98c8cff..1f35bd109 100644 --- a/Changelog +++ b/Changelog @@ -1,5 +1,7 @@ version 1.0.6 (not yet released) -------------------------------- +* have drawcoastlines use line segments instead of coastline polygons, to + avoid 'thickening' of lines around edges of map. * added support for cylindrical equal area ('cea') projection. * add 'arcgisimage' method for displaying background image retrieved from an ArcGIS server using the REST API (requires using 'epsg' keyword to define projection). diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index b66a63dfc..936de046d 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -1034,28 +1034,9 @@ def __init__(self, llcrnrlon=None, llcrnrlat=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 = list(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 _cylproj: - # 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(list(range(iprev,i))): - # coastsegs.append(list(zip(x[iprev:i],y[iprev:i]))) - # iprev = i - #else: - # coastsegs.append(seg) - #self.coastsegs = coastsegs # create geos Polygon structures for land areas. # currently only used in is_land method. self.landpolygons=[] @@ -1124,7 +1105,7 @@ def makegrid(self,nx,ny,returnxy=False): """ return self.projtran.makegrid(nx,ny,returnxy=returnxy) - def _readboundarydata(self,name): + def _readboundarydata(self,name,as_polygons=True): """ read boundary data, clip to map projection region. """ @@ -1176,7 +1157,7 @@ def _readboundarydata(self,name): if np.abs(int(lat0)) == 90: lon0=0. maptran = pyproj.Proj(proj='stere',lon_0=lon0,lat_0=lat0,R=re) # boundary polygon for ortho/gnom/nsper projection - # in stereographic coorindates. + # in stereographic coordinates. b = self._boundarypolyll.boundary blons = b[:,0]; blats = b[:,1] b[:,0], b[:,1] = maptran(blons, blats) @@ -1271,7 +1252,12 @@ def _readboundarydata(self,name): # if polygon instersects map projection # region, process it. if poly.intersects(boundarypolyll): - geoms = poly.intersection(boundarypolyll) + if name != 'gshhs' or as_polygons: + geoms = poly.intersection(boundarypolyll) + else: + # convert polygons to line segments + poly = _geoslib.LineString(poly.boundary) + geoms = poly.intersection(boundarypolyll) # iterate over geometries in intersection. for psub in geoms: b = psub.boundary @@ -1295,7 +1281,12 @@ def _readboundarydata(self,name): # if polygon instersects map projection # region, process it. if poly.intersects(boundarypolyll): - geoms = poly.intersection(boundarypolyll) + if name != 'gshhs' or as_polygons: + geoms = poly.intersection(boundarypolyll) + else: + # convert polygons to line segments + poly = _geoslib.LineString(poly.boundary) + geoms = poly.intersection(boundarypolyll) # iterate over geometries in intersection. for psub in geoms: # only coastlines are polygons, @@ -1319,7 +1310,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 tostere: + if name == 'gshhs' and as_polygons and self.projection in tostere: b[:,0], b[:,1] = maptran(b[:,0], b[:,1]) else: b[:,0], b[:,1] = self(b[:,0], b[:,1]) @@ -1327,22 +1318,36 @@ def _readboundarydata(self,name): # if less than two points are valid in # map proj coords, skip this geometry. if np.sum(goodmask) <= 1: continue - if name != 'gshhs': + if name != 'gshhs' or (name == 'gshhs' and not as_polygons): # 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 ortho/gnom/nsper projection, all points - # outside map projection region are eliminated - # with the above step, so we're done. - if name != 'gshhs' or\ - (self.projection in tostere and len(bx) > 4): + # split coastline segments that jump across entire plot. + xd = (bx[1:]-bx[0:-1])**2 + yd = (by[1:]-by[0:-1])**2 + dist = np.sqrt(xd+yd) + split = dist > 5000000. + if np.sum(split) and self.projection not in _cylproj: + 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(list(range(iprev,i))): + polygons.append(list(zip(bx[iprev:i],by[iprev:i]))) + iprev = i + else: polygons.append(list(zip(bx,by))) - polygon_types.append(typ) + polygon_types.append(typ) continue # create a GEOS geometry object. - poly = Shape(b) + if name == 'gshhs' and not as_polygons: + # convert polygons to line segments + poly = _geoslib.LineString(poly.boundary) + else: + poly = Shape(b) # this is a workaround to avoid # "GEOS_ERROR: TopologyException: # found non-noded intersection between ..." @@ -1773,9 +1778,14 @@ def drawcoastlines(self,linewidth=1.,linestyle='solid',color='k',antialiased=1,a """ if self.resolution is None: raise AttributeError('there are no boundary datasets associated with this Basemap instance') + # read in coastlines as line segments, only keeping those that + # intersect map boundary polygon. + if not hasattr(self,'coastsegs2'): + self.coastsegs2, types =\ + self._readboundarydata('gshhs',as_polygons=False) # get current axes instance (if none specified). ax = ax or self._check_ax() - coastlines = LineCollection(self.coastsegs,antialiaseds=(antialiased,)) + coastlines = LineCollection(self.coastsegs2,antialiaseds=(antialiased,)) coastlines.set_color(color) coastlines.set_linestyle(linestyle) coastlines.set_linewidth(linewidth) From 21b5f5be34517ad59b1adc5624984c6f4ea398e1 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Mon, 17 Sep 2012 13:36:28 -0600 Subject: [PATCH 3/9] process coastsegs and coast polygons ahead of time --- lib/mpl_toolkits/basemap/__init__.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 936de046d..388b8cd14 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -1037,6 +1037,9 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, for seg in self.coastsegs: x, y = list(zip(*seg)) self.coastpolygons.append((x,y)) + # replace coastsegs with line segments (instead of polygons) + self.coastsegs, types =\ + self._readboundarydata('gshhs',as_polygons=False) # create geos Polygon structures for land areas. # currently only used in is_land method. self.landpolygons=[] @@ -1778,14 +1781,9 @@ def drawcoastlines(self,linewidth=1.,linestyle='solid',color='k',antialiased=1,a """ if self.resolution is None: raise AttributeError('there are no boundary datasets associated with this Basemap instance') - # read in coastlines as line segments, only keeping those that - # intersect map boundary polygon. - if not hasattr(self,'coastsegs2'): - self.coastsegs2, types =\ - self._readboundarydata('gshhs',as_polygons=False) # get current axes instance (if none specified). ax = ax or self._check_ax() - coastlines = LineCollection(self.coastsegs2,antialiaseds=(antialiased,)) + coastlines = LineCollection(self.coastsegs,antialiaseds=(antialiased,)) coastlines.set_color(color) coastlines.set_linestyle(linestyle) coastlines.set_linewidth(linewidth) From 37fcff7c02210d716a136cce6633ae4fb881055a Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 18 Sep 2012 08:19:34 -0600 Subject: [PATCH 4/9] code cleanup --- lib/mpl_toolkits/basemap/__init__.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 388b8cd14..58cc88b7a 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -1031,7 +1031,8 @@ def __init__(self, llcrnrlon=None, llcrnrlat=None, # 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') + self.coastsegs, self.coastpolygontypes =\ + self._readboundarydata('gshhs',as_polygons=True) # reformat for use in matplotlib.patches.Polygon. self.coastpolygons = [] for seg in self.coastsegs: @@ -1108,7 +1109,7 @@ def makegrid(self,nx,ny,returnxy=False): """ return self.projtran.makegrid(nx,ny,returnxy=returnxy) - def _readboundarydata(self,name,as_polygons=True): + def _readboundarydata(self,name,as_polygons=False): """ read boundary data, clip to map projection region. """ @@ -1279,7 +1280,6 @@ def _readboundarydata(self,name,as_polygons=True): polys = [poly1,poly,poly2] for poly in polys: # try to fix "non-noded intersection" errors. - #if not poly.is_valid(): poly=poly.simplify(1.e-10) if not poly.is_valid(): poly=poly.fix() # if polygon instersects map projection # region, process it. @@ -1292,16 +1292,12 @@ def _readboundarydata(self,name,as_polygons=True): geoms = poly.intersection(boundarypolyll) # 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) - if name != 'gshhs' or len(bx) > 4: + if not as_polygons or len(bx) > 4: polygons.append(list(zip(bx,by))) polygon_types.append(typ) # if map boundary polygon is not valid in lat/lon @@ -1354,7 +1350,6 @@ def _readboundarydata(self,name,as_polygons=True): # this is a workaround to avoid # "GEOS_ERROR: TopologyException: # found non-noded intersection between ..." - #if not poly.is_valid(): poly=poly.simplify(1.e-10) if not poly.is_valid(): poly=poly.fix() # if geometry instersects map projection # region, and doesn't have any invalid points, process it. @@ -1385,7 +1380,7 @@ def _readboundarydata(self,name,as_polygons=True): b[:,0], b[:,1] = maptran(b[:,0], b[:,1], inverse=True) # orthographic/gnomonic/nsper. b[:,0], b[:,1]= self(b[:,0], b[:,1]) - if name != 'gshhs' or b.shape[0] > 4: + if not as_polygons or len(bx) > 4: polygons.append(list(zip(b[:,0],b[:,1]))) polygon_types.append(typ) return polygons, polygon_types From 8faf74b9a3dc49b3d1418894970a46ddba109436 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 18 Sep 2012 08:33:55 -0600 Subject: [PATCH 5/9] fix for disappearing Eurasia line segments --- lib/mpl_toolkits/basemap/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 58cc88b7a..9c48c625a 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -1288,7 +1288,9 @@ def _readboundarydata(self,name,as_polygons=False): geoms = poly.intersection(boundarypolyll) else: # convert polygons to line segments - poly = _geoslib.LineString(poly.boundary) + # note: use fix method here or Eurasia + # line segments sometimes disappear. + poly = _geoslib.LineString(poly.fix().boundary) geoms = poly.intersection(boundarypolyll) # iterate over geometries in intersection. for psub in geoms: From b29ce6edcc6050bb09ad297602d28782d478c04d Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 18 Sep 2012 08:38:28 -0600 Subject: [PATCH 6/9] fix typo --- 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 9c48c625a..bd11a90cc 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -1382,7 +1382,7 @@ def _readboundarydata(self,name,as_polygons=False): b[:,0], b[:,1] = maptran(b[:,0], b[:,1], inverse=True) # orthographic/gnomonic/nsper. b[:,0], b[:,1]= self(b[:,0], b[:,1]) - if not as_polygons or len(bx) > 4: + if not as_polygons or len(b) > 4: polygons.append(list(zip(b[:,0],b[:,1]))) polygon_types.append(typ) return polygons, polygon_types From 53bc4d88f553f8034ec5da948c929c2fe1485a03 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 18 Sep 2012 08:44:42 -0600 Subject: [PATCH 7/9] only gshhs coastlines can be polygons --- lib/mpl_toolkits/basemap/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index bd11a90cc..cf04e78df 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -1119,6 +1119,8 @@ def _readboundarydata(self,name,as_polygons=False): If you are requesting a 'full' resolution dataset, you may need to download and install those files separately (see the basemap README for details).""") + # only gshhs coastlines can be polygons. + if name != 'gshhs': as_polygons=False 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') From 686fc5aed675f868031473541d3a4e26cc41d4bb Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 18 Sep 2012 09:06:43 -0600 Subject: [PATCH 8/9] fix streamplot --- lib/mpl_toolkits/basemap/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index cf04e78df..d89c1e98a 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -1119,7 +1119,7 @@ def _readboundarydata(self,name,as_polygons=False): If you are requesting a 'full' resolution dataset, you may need to download and install those files separately (see the basemap README for details).""") - # only gshhs coastlines can be polygons. + # only gshhs coastlines can be polygons. if name != 'gshhs': as_polygons=False try: bdatfile = open(os.path.join(basemap_datadir,name+'_'+self.resolution+'.dat'),'rb') @@ -3674,8 +3674,8 @@ def streamplot(self, x, y, u, v, *args, **kwargs): ax.hold(b) raise ax.hold(b) - if plt is not None and ret.get_array() is not None: - plt.sci(ret) + if plt is not None and ret.lines.get_array() is not None: + plt.sci(ret.lines) # clip for round polar plots. # streamplot arrows not returned in matplotlib 1.1.1, so clip all # FancyArrow patches attached to axes instance. From 836d5072c221121c0e6391449ec5b380bc4cc4d8 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 18 Sep 2012 09:18:13 -0600 Subject: [PATCH 9/9] fix for recent changes in tripcolor --- examples/ploticos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ploticos.py b/examples/ploticos.py index 85c9f68ca..dbc2c08f4 100644 --- a/examples/ploticos.py +++ b/examples/ploticos.py @@ -14,7 +14,7 @@ map.drawmapboundary() # tri=True forwards to axes.tripcolor #z = ma.masked_where(z < 1.e-5, z) # for testing masked arrays. -map.pcolor(x,y,z,tri=True,shading='faceted',vmin=0,vmax=3000) +map.pcolor(x,y,z,tri=True,shading='flat',edgecolors='k',cmap=plt.cm.hot_r,vmin=0,vmax=3000) #map.contourf(x,y,z,np.arange(0,3000,150),tri=True) #map.contour(x,y,z,np.arange(0,3000,150),tri=True) plt.title('pcolor plot on a global icosahedral mesh')