Pierre - That’s great. As soon as I get the zorder problem figured out (can’t have land AND ocean masks on!) , I’ll give your solution a try.
···
On Nov 2, 2007, at 10:06 AM, Pierre GM wrote:
On Friday 02 November 2007 11:51:55 Michael Hearne wrote:
- Has anyone figured out a way to make an ocean mask? I need my
map to look like this
Assuming that you have gdal installed, you can use this set of functions.
Basically, you use OGR to compute the differences between the background and
the land polygons… Works fine for my applications (I interpolate some
rainfall fields over AL/GA/FL, and need to hide the interpolated values in
the Gulf of Mexico), but might still be a tad buggy. Please don’t hesittate
to contact me if you need help with that.
#####--------------------------------------------------------------------------
#---- — OGR conversion —
#####--------------------------------------------------------------------------
import ogr
def polygon_to_geometry(polygon):
"""Creates a new ogr.Geometry object from a matplolib.Polygon."""
if not isinstance(polygon, Polygon):
raise TypeError, "The input data should be a valid Polygon object!"
listvert = ["%s %s" % (x,y) for (x,y) in polygon.get_verts()]
listvert.append(listvert[0])
return ogr.CreateGeometryFromWkt("POLYGON ((%s))" % ','.join(listvert))
#-----------------------------------------------------------------------------
#---- OGR to matplotlib —
#-----------------------------------------------------------------------------
def _geometry_to_vertices(geometry):
"""Creates a list of vertices (x,y) from the current geometry."""
verts = []
for nl in range(geometry.GetGeometryCount()):
line = geometry.GetGeometryRef(nl)
points = [(line.GetX(i), line.GetY(i)) for i in
range(line.GetPointCount())]
verts.append(points)
return verts
def geometry_to_vertices(geometry):
"""Creates lists of vertices (x,y) from the current geometry."""
if not isinstance(geometry, ogr.Geometry):
raise TypeError, "The input data should be a valid ogr.Geometry
object!"
vertices = []
#
if geometry.GetGeometryType() == ogr.wkbPolygon:
vertices.extend(_geometry_to_vertices(geometry))
elif geometry.GetGeometryType() == ogr.wkbMultiPolygon:
for np in range(geometry.GetGeometryCount()):
poly = geometry.GetGeometryRef(np)
vertices.extend(_geometry_to_vertices(poly))
return vertices
#####--------------------------------------------------------------------------
#---- — Add-ons
#####--------------------------------------------------------------------------
def fillwaterbodies(basemap, color=‘blue’, inlands=True, ax=None,
zorder=None):
"""Fills the water bodies with color.
If inlands is True, inland water bodies are also filled.
:Inputs:
basemap : Basemap
The basemap on which to fill.
color : string/RGBA tuple *['blue']*
Filling color
inlands : boolean *[True]*
Whether inland water bodies should be filled.
ax : Axes instance *[None]*
Axe on which to plot. If None, the current axe is selected.
zorder : integer *[None]*
zorder of the water bodies polygons.
"""
if not isinstance(basemap, Basemap):
raise TypeError, "The input basemap should be a valid Basemap object!"
#
if ax is None and basemap.ax is None:
try:
ax = pylab.gca()
except:
import pylab
ax = pylab.gca()
#
coastpolygons = basemap.coastpolygons
coastpolygontypes = basemap.coastpolygontypes
(llx, lly) = (basemap.llcrnrx, basemap.llcrnry)
(urx, ury) = (basemap.urcrnrx,basemap.urcrnry)
background = Polygon([(llx, lly), (urx, lly), (urx, ury), (llx, ury)])
#
ogr_background = polygon_to_geometry(background)
inland_polygons = []
#
for (poly, polytype) in zip(coastpolygons, coastpolygontypes):
if polytype != 2:
verts = ["%s %s" % (x,y) for (x,y) in zip(*poly)]
ogr_poly = ogr.CreateGeometryFromWkt("POLYGON
((%s))" % ‘,’.join(verts))
ogr_background = ogr_background.Difference(ogr_poly)
else:
inland_polygons.append(Polygon(zip(*poly),
facecolor=color,
edgecolor=color,
linewidth=0))
#
background = geometry_to_vertices(ogr_background)
for xy in background:
patch = Polygon(xy, facecolor=color, edgecolor=color, linewidth=0)
if zorder is not None:
patch.set_zorder(zorder)
ax.add_patch(patch)
#
if inlands:
for poly in inland_polygons:
if zorder is not None:
poly.set_zorder(zorder)
ax.add_patch(poly)
basemap.set_axes_limits(ax=ax)
return
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
Michael Hearne
mhearne@…924…
(303) 273-8620
USGS National Earthquake Information Center
1711 Illinois St. Golden CO 80401
Senior Software Engineer
Synergetics, Inc.
