Re: [Matplotlib-users] Basemap questions

[Michael H. <mh...@us...>]

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.

One question regarding - I found this page:
http://trac.osgeo.org/gdal/wiki/GdalOgrInPython

but it wasn't immediately apparent how to _get_ it.  Do I have to  
download/build ogr first, then install the python bindings?

--Mike
On Nov 2, 2007, at 10:06 AM, Pierre GM wrote:

> On Friday 02 November 2007 11:51:55 Michael Hearne wrote:
>> 2) 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
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users




------------------------------------------------------
Michael Hearne
mh...@us...
(303) 273-8620
USGS National Earthquake Information Center
1711 Illinois St. Golden CO 80401
Senior Software Engineer
Synergetics, Inc.
------------------------------------------------------