Basemap questions

I have two questions:

  1. The fillcontinents() method has a zorder keyword parameter. Is this supposed to work with imshow()? I have the latest tarball from the website, and I can’t get my image to paint on top of the continents.

  2. Has anyone figured out a way to make an ocean mask? I need my map to look like this

(http://earthquake.usgs.gov/eqcenter/pager/us/2007itah/us/5/shakemap.600.jpg)

, where the image data (the yellow and orange bits) that extends out beyond the land boundary is masked by a blue ocean.

Thanks,

Mike

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.

···

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

#####--------------------------------------------------------------------------
#---- --- 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

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:

  1. 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.


Michael Hearne wrote:

I have two questions:

1) The fillcontinents() method has a zorder keyword parameter. Is this supposed to work with imshow()? I have the latest tarball from the website, and I can't get my image to paint on top of the continents.

2) Has anyone figured out a way to make an _ocean_ mask? I need my map to look like this

( http://earthquake.usgs.gov/eqcenter/pager/us/2007itah/us/5/shakemap.600.jpg)

, where the image data (the yellow and orange bits) that extends out beyond the land boundary is masked by a blue ocean.

Thanks,

Mike

Mike: You can set the background color for the axes instance to blue, and then draw the filled continents over it.

You'll also have to make your image a masked array, and set the missing values to be transparent.

The basemap examples plotmap_masked.py shows how to do both of these things.

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Jeff - I looked at that example file, and I think there’s a big difference - your etopo base data set is global, and you can plot over the data in the oceans by setting the mask on all pixels less than zero.

My dataset (a map of earthquake shaking) is not global, and actually has NO missing data. I think I need a way to “clip” the data by the land mask - that is, find all of the pixels that are NOT on land, and then mask them off.

Is there an easy way to do this with matplotlib/basemap tools?

Regarding my other issue - I used my script to test x/y offset values: [0.05,0.1,0.5,1.0,10] and couldn’t see any difference. I’d be more than happy to provide test output, or debugging information…

Just to be clear - these offsets are supposed to move the meridian and/or parallel labels around with respect to the map edge? My actual goal is to get the labels inside the edge of the map (I tried negative numbers to accomplish this, to no effect.)

On a positive note, I can make solid lines!

Thanks for all of your help,

Mike

···

On Nov 2, 2007, at 11:47 AM, Jeff Whitaker wrote:

Michael Hearne wrote:

I have two questions:

  1. The fillcontinents() method has a zorder keyword parameter. Is this supposed to work with imshow()? I have the latest tarball from the website, and I can’t get my image to paint on top of the continents.
  1. Has anyone figured out a way to make an ocean mask? I need my map to look like this

( http://earthquake.usgs.gov/eqcenter/pager/us/2007itah/us/5/shakemap.600.jpg)

, where the image data (the yellow and orange bits) that extends out beyond the land boundary is masked by a blue ocean.

Thanks,

Mike

Mike: You can set the background color for the axes instance to blue, and then draw the filled continents over it.

You’ll also have to make your image a masked array, and set the missing values to be transparent.

The basemap examples plotmap_masked.py shows how to do both of these things.

-Jeff

Jeffrey S. Whitaker Phone : (303)497-6313

Meteorologist FAX : (303)497-6449

NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@…120…259…

325 Broadway Office : Skaggs Research Cntr 1D-124

Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg


Michael Hearne

mhearne@…924…

(303) 273-8620

USGS National Earthquake Information Center

1711 Illinois St. Golden CO 80401

Senior Software Engineer

Synergetics, Inc.


Michael Hearne wrote:

Jeff - I looked at that example file, and I think there's a big difference - your etopo base data set is global, and you can plot over the data in the oceans by setting the mask on all pixels less than zero.

My dataset (a map of earthquake shaking) is not global, and actually has NO missing data. I think I need a way to "clip" the data by the land mask - that is, find all of the pixels that are NOT on land, and then mask them off.

Mike:

If it's not global, is it just defined for land points? If so, it can't be a 2-D grid, so you won't be able to plot it with imshow anyway. Can you explain the structure of the data?

Is there an easy way to do this with matplotlib/basemap tools?

Not really. You'll have to define a sea mask for your grid and use that the create a masked array. There is a land-sea mask dataset included in basemap, but it may not match the resolution of your grid.

Regarding my other issue - I used my script to test x/y offset values: [0.05,0.1,0.5,1.0,10] and couldn't see any difference. I'd be more than happy to provide test output, or debugging information...

Just to be clear - these offsets are supposed to move the meridian and/or parallel labels around with respect to the map edge? My actual goal is to get the labels inside the edge of the map (I tried negative numbers to accomplish this, to no effect.)

You need to define an offset as a fraction of the map width - the numbers you are giving are too small to notice any difference. As I said before, try something like -0.01*(m.max-m.min).

On a positive note, I _can_ make solid lines!

Good!

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

I'm afraid that not.
Your problem is that you need to find for each cell of your gridded data
whether you are on land or in an ocean, and mask the corresponding cell in
the latter case. Such a manipulation of geometry is not available in
matplotlib (yet, and will it ever be?).
You could of course start writing your tests. However, that's exactly what OGR
is designed to do, and the couple of functions I sent earlier can be quite
helpful. Once again, OGR is part of the GDAL suite.

···

On Friday 02 November 2007 15:10:15 Michael Hearne wrote:

Is there an easy way to do this with matplotlib/basemap tools?

Jeff -

My data set is actually dynamically generated by a program called ShakeMap. It’s a 2D grid, with an extent usually about 600 kilometers on a side, centered wherever the earthquake happened to be. The ShakeMap program does not know or care that some of the data may be under water, but for display purposed, I do! The grid is also in a geographic projection (latitude/longitude coordinates assumed to be cartesian).

So in this test instance (on a data set near Taiwan), my map width is about 5.91 degrees longitude, and my height is about 5.5 degrees latitude.

If I set xoffset=-0.01*5.91, I get -0.05. This is not noticeably different than the default.

Is the problem that my dataset is not projected?

–Mike

···

On Nov 2, 2007, at 1:33 PM, Jeff Whitaker wrote:

Michael Hearne wrote:

Jeff - I looked at that example file, and I think there’s a big difference - your etopo base data set is global, and you can plot over the data in the oceans by setting the mask on all pixels less than zero.

My dataset (a map of earthquake shaking) is not global, and actually has NO missing data. I think I need a way to “clip” the data by the land mask - that is, find all of the pixels that are NOT on land, and then mask them off.

Mike:

If it’s not global, is it just defined for land points? If so, it can’t be a 2-D grid, so you won’t be able to plot it with imshow anyway. Can you explain the structure of the data?

Is there an easy way to do this with matplotlib/basemap tools?

Not really. You’ll have to define a sea mask for your grid and use that the create a masked array. There is a land-sea mask dataset included in basemap, but it may not match the resolution of your grid.

Regarding my other issue - I used my script to test x/y offset values: [0.05,0.1,0.5,1.0,10] and couldn’t see any difference. I’d be more than happy to provide test output, or debugging information…

Just to be clear - these offsets are supposed to move the meridian and/or parallel labels around with respect to the map edge? My actual goal is to get the labels inside the edge of the map (I tried negative numbers to accomplish this, to no effect.)

You need to define an offset as a fraction of the map width - the numbers you are giving are too small to notice any difference. As I said before, try something like -0.01*(m.max-m.min).

On a positive note, I can make solid lines!

Good!

-Jeff

Jeffrey S. Whitaker Phone : (303)497-6313

Meteorologist FAX : (303)497-6449

NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@…259…

325 Broadway Office : Skaggs Research Cntr 1D-124

Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg


Michael Hearne

mhearne@…120…924…

(303) 273-8620

USGS National Earthquake Information Center

1711 Illinois St. Golden CO 80401

Senior Software Engineer

Synergetics, Inc.


Michael Hearne wrote:

Jeff -

My data set is actually dynamically generated by a program called ShakeMap. It's a 2D grid, with an extent usually about 600 kilometers on a side, centered wherever the earthquake happened to be. The ShakeMap program does not know or care that some of the data may be under water, but for display purposed, I do! The grid is also in a geographic projection (latitude/longitude coordinates assumed to be cartesian).

Michael: You'll either have to create your own land-sea mask, or use Pierre's method with gdal.

So in this test instance (on a data set near Taiwan), my map width is about 5.91 degrees longitude, and my height is about 5.5 degrees latitude.

If I set xoffset=-0.01*5.91, I get -0.05. This is not noticeably different than the default.

Is the problem that my dataset is not projected?

What is m.xmax-m.xmin? (m is the basemap instance, xmin and xmax are instance variables) That's what you have to use - the distance in map projection coordinates, not lat/lon coordinates (although these will be the same if you are using projection='cyl').

-Jeff

···

--Mike
On Nov 2, 2007, at 1:33 PM, Jeff Whitaker wrote:

Michael Hearne wrote:

Jeff - I looked at that example file, and I think there's a big difference - your etopo base data set is global, and you can plot over the data in the oceans by setting the mask on all pixels less than zero.

My dataset (a map of earthquake shaking) is not global, and actually has NO missing data. I think I need a way to "clip" the data by the land mask - that is, find all of the pixels that are NOT on land, and then mask them off.

Mike:

If it's not global, is it just defined for land points? If so, it can't be a 2-D grid, so you won't be able to plot it with imshow anyway. Can you explain the structure of the data?

Is there an easy way to do this with matplotlib/basemap tools?

Not really. You'll have to define a sea mask for your grid and use that the create a masked array. There is a land-sea mask dataset included in basemap, but it may not match the resolution of your grid.

Regarding my other issue - I used my script to test x/y offset values: [0.05,0.1,0.5,1.0,10] and couldn't see any difference. I'd be more than happy to provide test output, or debugging information...

Just to be clear - these offsets are supposed to move the meridian and/or parallel labels around with respect to the map edge? My actual goal is to get the labels inside the edge of the map (I tried negative numbers to accomplish this, to no effect.)

You need to define an offset as a fraction of the map width - the numbers you are giving are too small to notice any difference. As I said before, try something like -0.01*(m.max-m.min).

On a positive note, I _can_ make solid lines!

Good!

-Jeff


Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@…259… <mailto:Jeffrey.S.Whitaker@…259…>
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

------------------------------------------------------
Michael Hearne
mhearne@…924… <mailto:mhearne@…924…>
(303) 273-8620
USGS National Earthquake Information Center
1711 Illinois St. Golden CO 80401
Senior Software Engineer
Synergetics, Inc.
------------------------------------------------------

--
Jeffrey S. Whitaker Phone : (303)497-6313
NOAA/OAR/CDC R/PSD1 FAX : (303)497-6449
325 Broadway Boulder, CO, USA 80305-3328