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.