Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

avoid thickening of coastlines around edges of plot #78

Merged
merged 9 commits into from Sep 18, 2012
Merged

avoid thickening of coastlines around edges of plot #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
51d44e5
add cea projection.
jswhit Sep 5, 2012
66e8bbb
have drawcoastlines use line segments instead of polygons to
jswhit Sep 17, 2012
21b5f5b
process coastsegs and coast polygons ahead of time
jswhit Sep 17, 2012
37fcff7
code cleanup
jswhit Sep 18, 2012
8faf74b
fix for disappearing Eurasia line segments
jswhit Sep 18, 2012
b29ce6e
fix typo
jswhit Sep 18, 2012
53bc4d8
only gshhs coastlines can be polygons
jswhit Sep 18, 2012
686fc5a
fix streamplot
jswhit Sep 18, 2012
836d507
fix for recent changes in tripcolor
jswhit Sep 18, 2012
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 2 additions & 0 deletions Changelog
View file
Open in desktop
@@ -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).
Expand Down
8 changes: 8 additions & 0 deletions doc/users/cea.rst
View file
Open in desktop
@@ -0,0 +1,8 @@
.. _cea:

Cylindrial Equal-Area Projection
================================

It is what is says.

.. plot:: users/figures/cea.py
17 changes: 17 additions & 0 deletions doc/users/figures/cea.py
View file
Open in desktop
@@ -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()
1 change: 1 addition & 0 deletions doc/users/mapsetup.rst
View file
Open in desktop
Expand Up @@ -53,6 +53,7 @@ time, many compromise between the two.
poly.rst
mill.rst
gall.rst
cea.rst
lcc.rst
laea.rst
stere.rst
Expand Down
2 changes: 1 addition & 1 deletion examples/ploticos.py
View file
Open in desktop
Expand Up @@ -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')
Expand Down
93 changes: 50 additions & 43 deletions lib/mpl_toolkits/basemap/__init__.py
View file
Open in desktop
Expand Up @@ -1031,31 +1031,16 @@ 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 = []
# 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
# 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=[]
Expand Down Expand Up @@ -1124,7 +1109,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=False):
"""
read boundary data, clip to map projection region.
"""
Expand All @@ -1134,6 +1119,8 @@ def _readboundarydata(self,name):
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')
Expand Down Expand Up @@ -1176,7 +1163,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)
Expand Down Expand Up @@ -1271,7 +1258,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
Expand All @@ -1290,24 +1282,26 @@ def _readboundarydata(self,name):
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.
if poly.intersects(boundarypolyll):
geoms = poly.intersection(boundarypolyll)
if name != 'gshhs' or as_polygons:
geoms = poly.intersection(boundarypolyll)
else:
# convert polygons to line segments
# 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:
# 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
Expand All @@ -1319,34 +1313,47 @@ 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])
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 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 ..."
#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.
Expand Down Expand Up @@ -1377,7 +1384,7 @@ def _readboundarydata(self,name):
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(b) > 4:
polygons.append(list(zip(b[:,0],b[:,1])))
polygon_types.append(typ)
return polygons, polygon_types
Expand Down Expand Up @@ -3667,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.
Expand Down