basemap: Mask the ocean [SEC=UNCLASSIFIED]

Listers,

I'm using basemap to plot randomly sampled values (x,y,z) through hexbin. This produces a very nice result. Some sample code is:

···

----------
import numpy as np
from numpy.random import seed
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.mlab import griddata

ll_lat = -38.39477 # extent of area of interest
ll_lon = 144.54767
ur_lat = -37.51642
ur_lon = 145.67144

num_points = 100 # sample points

# create random sampling over the area of interest
seed(0)
data = np.ones((3, num_points))
data[0,:] *= ll_lon + np.random.random((num_points))*(ur_lon-ll_lon)
data[1,:] *= ll_lat + np.random.random((num_points))*(ur_lat-ll_lat)
data[2,:] *= np.random.random((num_points))*10000

# plot the data
fig = plt.figure()
ax = fig.add_subplot(111)
m = Basemap(projection='cyl', llcrnrlat=ll_lat, urcrnrlat=ur_lat,
            llcrnrlon=ll_lon, urcrnrlon=ur_lon, resolution='f',
            suppress_ticks=False, area_thresh=0.5)
plt.hexbin(data[0,:], data[1,:], data[2,:], zorder=3)
m.fillcontinents(color=(0.8,0.8,0.8,0), zorder=1)
m.drawcoastlines(linewidth=0.25, color='k', zorder=2)
plt.show()
----------

This contrived example shows a sparse set of hexagons on both land and ocean. I would like the hexagons over the ocean to be hidden. I can make the ones on land disappear by changing the 'zorder' parameter of .hexbin() to 0. However I have found no way of doing the inverse and hiding hexagons over the ocean.

Using drawlsmask() is too crude at a 5-minute resolution.

Any ideas?

Thanks,
Ross

Ross,

one way is to mask (or remove) ocean points using the _geoslib module provided with basemap.
When you create a Basemap instance, you can retrieve all its polygons land (continents and islands) with “mymap.coastpolygons”.

Thay are stored as numpy arrays, and you can convert them to _geoslib.Polygon objects :

poly = _geoslib.Polygon(N.asarray(coastalpoly).T)

Then you loop over all Polygons and all (x,y) points and test :

good_point = _geoslib.Point((x,y)).within(poly)

Thanks to this method, you can choose you optimal resolution.
You can even compute the intersection of you hexagons with coastal polygons using .intersection() and .area (instead of simply checking if the center is inside) and then reject points depending the fraction of the cell covered by land (or ocean).

···

On Mon, Nov 2, 2009 at 8:07 AM, <Ross.Wilson@…2849…> wrote:

Listers,

I’m using basemap to plot randomly sampled values (x,y,z) through hexbin. This produces a very nice result. Some sample code is:


import numpy as np

from numpy.random import seed

import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

from matplotlib.mlab import griddata

ll_lat = -38.39477 # extent of area of interest

ll_lon = 144.54767

ur_lat = -37.51642

ur_lon = 145.67144

num_points = 100 # sample points

create random sampling over the area of interest

seed(0)

data = np.ones((3, num_points))

data[0,:] = ll_lon + np.random.random((num_points))(ur_lon-ll_lon)

data[1,:] = ll_lat + np.random.random((num_points))(ur_lat-ll_lat)

data[2,:] *= np.random.random((num_points))*10000

plot the data

fig = plt.figure()

ax = fig.add_subplot(111)

m = Basemap(projection=‘cyl’, llcrnrlat=ll_lat, urcrnrlat=ur_lat,

        llcrnrlon=ll_lon, urcrnrlon=ur_lon, resolution='f',

        suppress_ticks=False, area_thresh=0.5)

plt.hexbin(data[0,:], data[1,:], data[2,:], zorder=3)

m.fillcontinents(color=(0.8,0.8,0.8,0), zorder=1)

m.drawcoastlines(linewidth=0.25, color=‘k’, zorder=2)

plt.show()


This contrived example shows a sparse set of hexagons on both land and ocean. I would like the hexagons over the ocean to be hidden. I can make the ones on land disappear by changing the ‘zorder’ parameter of .hexbin() to 0. However I have found no way of doing the inverse and hiding hexagons over the ocean.

Using drawlsmask() is too crude at a 5-minute resolution.

Any ideas?

Thanks,

Ross


Come build with us! The BlackBerry® Developer Conference in SF, CA

is the only developer event you need to attend this year. Jumpstart your

developing skills, take BlackBerry mobile applications to market and stay

ahead of the curve. Join us from November 9 - 12, 2009. Register now!

http://p.sf.net/sfu/devconference


Matplotlib-users mailing list

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users


Stephane Raynaud

Stephane Raynaud wrote:

Ross,

one way is to mask (or remove) ocean points using the _geoslib module provided with basemap.
When you create a Basemap instance, you can retrieve all its polygons land (continents and islands) with "mymap.coastpolygons".
Thay are stored as numpy arrays, and you can convert them to _geoslib.Polygon objects :

poly = _geoslib.Polygon(N.asarray(coastalpoly).T)

Then you loop over all Polygons and all (x,y) points and test :

good_point = _geoslib.Point((x,y)).within(poly)

Thanks to this method, you can choose you optimal resolution.
You can even compute the intersection of you hexagons with coastal polygons using .intersection() and .area (instead of simply checking if the center is inside) and then reject points depending the fraction of the cell covered by land (or ocean).

Following Stephane's excellent suggestion, here's a prototype Basemap method that checks to see if a point is on land or over water. Ross - if you find it useful I'll include it in the next release. Note that it will be slow for lots of points or large map regions.

-Jeff

hexbin.py (1.56 KB)

···

On Mon, Nov 2, 2009 at 8:07 AM, <Ross.Wilson@…2808… > <mailto:Ross.Wilson@…2808…>> wrote:

    Listers,

    I'm using basemap to plot randomly sampled values (x,y,z) through
    hexbin. This produces a very nice result. Some sample code is:
    ----------
    import numpy as np
    from numpy.random import seed
    import matplotlib.pyplot as plt
    from mpl_toolkits.basemap import Basemap
    from matplotlib.mlab import griddata

    ll_lat = -38.39477 # extent of area of interest
    ll_lon = 144.54767
    ur_lat = -37.51642
    ur_lon = 145.67144

    num_points = 100 # sample points

    # create random sampling over the area of interest
    seed(0)
    data = np.ones((3, num_points))
    data[0,:] *= ll_lon + np.random.random((num_points))*(ur_lon-ll_lon)
    data[1,:] *= ll_lat + np.random.random((num_points))*(ur_lat-ll_lat)
    data[2,:] *= np.random.random((num_points))*10000

    # plot the data
    fig = plt.figure()
    ax = fig.add_subplot(111)
    m = Basemap(projection='cyl', llcrnrlat=ll_lat, urcrnrlat=ur_lat,
               llcrnrlon=ll_lon, urcrnrlon=ur_lon, resolution='f',
               suppress_ticks=False, area_thresh=0.5)
    plt.hexbin(data[0,:], data[1,:], data[2,:], zorder=3)
    m.fillcontinents(color=(0.8,0.8,0.8,0), zorder=1)
    m.drawcoastlines(linewidth=0.25, color='k', zorder=2)
    plt.show()
    ----------

    This contrived example shows a sparse set of hexagons on both land
    and ocean. I would like the hexagons over the ocean to be hidden.
     I can make the ones on land disappear by changing the 'zorder'
    parameter of .hexbin() to 0. However I have found no way of doing
    the inverse and hiding hexagons over the ocean.

    Using drawlsmask() is too crude at a 5-minute resolution.

    Any ideas?

    Thanks,
    Ross

    ------------------------------------------------------------------------------
    Come build with us! The BlackBerry(R) Developer Conference in SF, CA
    is the only developer event you need to attend this year.
    Jumpstart your
    developing skills, take BlackBerry mobile applications to market
    and stay
    ahead of the curve. Join us from November 9 - 12, 2009. Register now!
    http://p.sf.net/sfu/devconference
    _______________________________________________
    Matplotlib-users mailing list
    Matplotlib-users@lists.sourceforge.net
    <mailto:Matplotlib-users@lists.sourceforge.net>
    https://lists.sourceforge.net/lists/listinfo/matplotlib-users

--
Stephane Raynaud
------------------------------------------------------------------------

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
------------------------------------------------------------------------

_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
  
--
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-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Stephane Raynaud wrote:

Ross,

one way is to mask (or remove) ocean points using the _geoslib module
provided with basemap.
When you create a Basemap instance, you can retrieve all its polygons
land (continents and islands) with "mymap.coastpolygons".
Thay are stored as numpy arrays, and you can convert them to
_geoslib.Polygon objects :

poly = _geoslib.Polygon(N.asarray(coastalpoly).T)

Then you loop over all Polygons and all (x,y) points and test :

good_point = _geoslib.Point((x,y)).within(poly)

Thanks to this method, you can choose you optimal resolution.
You can even compute the intersection of you hexagons with coastal
polygons using .intersection() and .area (instead of simply checking
if the center is inside) and then reject points depending the fraction
of the cell covered by land (or ocean).

Following Stephane's excellent suggestion, here's a prototype Basemap
method that checks to see if a point is on land or over water. Ross -
if you find it useful I'll include it in the next release. Note that it
will be slow for lots of points or large map regions.

-Jeff

···

-----Original Message-----
From: Jeff Whitaker [mailto:jswhit@…146…]
Sent: Tuesday, 3 November 2009 02:53
To: Stephane Raynaud
Cc: Wilson Ross; matplotlib-users@lists.sourceforge.net
Subject: Re: [Matplotlib-users] basemap: Mask the ocean [SEC=UNCLASSIFIED]

--------------------------

Yes, Stephane's approach is nice and Jeff has nicely encapsulated the approach. I'll put that into my bag of tricks!

However it doesn't quite do what I want. My data does not have any points in the ocean; the hex binning creates hexagonal patches that extend out into the ocean. As a physicist I say "that's a representation artifact" and leave it at that, but my end-users want that 'bleed' into the ocean removed. My argument that they are losing data falls on deaf ears.

Here's an even more contrived example that improves on my poor previous attempt to explain the problem:
--------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

ll_lat = -38.10 # extent of area of interest
ll_lon = 145.06
ur_lat = -38.00
ur_lon = 145.16

# create data points all on land
fudge = ((ll_lon+0.05, ur_lat-0.00, 1000),
         (ll_lon+0.05, ur_lat-0.01, 2000),
         (ll_lon+0.05, ur_lat-0.02, 3000),
         (ll_lon+0.05, ur_lat-0.03, 4000),
         (ll_lon+0.04, ur_lat-0.025, 1000),
         (ll_lon+0.043, ur_lat-0.036, 10000),
         (ll_lon+0.047, ur_lat-0.041, 20000),
         (ll_lon+0.07, ur_lat-0.07, 4000),
         (ll_lon+0.08, ur_lat-0.08, 3000),
         (ll_lon+0.09, ur_lat-0.09, 2000),
         (ll_lon+0.10, ur_lat-0.10, 1000))

data = np.ones((3, len(fudge)))
for (i, (lon, lat, val)) in enumerate(fudge):
    data[0,i] = lon
    data[1,i] = lat
    data[2,i] = val

# plot the data
fig = plt.figure()
m = Basemap(projection='cyl', llcrnrlat=ll_lat, urcrnrlat=ur_lat,
                   llcrnrlon=ll_lon, urcrnrlon=ur_lon, resolution='f')
plt.hexbin(data[0,:], data[1,:], data[2,:], gridsize=5)
m.drawcoastlines(linewidth=0.5, color='k')
plt.show()
--------------------------------------------------------

None of the data points are on land, as can be seen if you change the 'gridsize' parameter to 100 or so (land is to the NE). However, when gridsize=5, the red hexagon in the middle extends out into the ocean. My users want that hexagon and others like it clipped at the coastline.

Crazy people, users.

Ross

Maybe using path clipping you can plot only over land like this :

···

import numpy as np
from mpl_toolkits.basemap import Basemap

import pylab as plt
import matplotlib.patches as patches
import matplotlib.transforms as transforms

ll_lat = -38.10 # extent of area of interest
ll_lon = 145.06
ur_lat = -38.00
ur_lon = 145.16

create data points all on land

fudge = ((ll_lon+0.05, ur_lat-0.00, 1000),
(ll_lon+0.05, ur_lat-0.01, 2000),
(ll_lon+0.05, ur_lat-0.02, 3000),
(ll_lon+0.05, ur_lat-0.03, 4000),

    (ll_lon+0.04,  ur_lat-0.025,  1000),
    (ll_lon+0.043, ur_lat-0.036,  10000),
    (ll_lon+0.047, ur_lat-0.041,  20000),
    (ll_lon+0.07,  ur_lat-0.07,   4000),
    (ll_lon+0.08,  ur_lat-0.08,   3000),

    (ll_lon+0.09,  ur_lat-0.09,   2000),
    (ll_lon+0.10,  ur_lat-0.10,   1000))

data = np.ones((3, len(fudge)))
for (i, (lon, lat, val)) in enumerate(fudge):
data[0,i] = lon
data[1,i] = lat

data[2,i] = val

plot the data

fig = plt.figure()
m = Basemap(projection=‘cyl’, llcrnrlat=ll_lat, urcrnrlat=ur_lat,
llcrnrlon=ll_lon, urcrnrlon=ur_lon, resolution=‘f’)

+CHANGES

coll = plt.hexbin(data[0,:], data[1,:], data[2,:], gridsize=5)
ax = plt.gca()
for n, poly in enumerate(m.coastpolygons):
type = m.coastpolygontypes[n]
if type in [1,3]:
p = patches.Polygon(np.asarray(poly).T)

    p.set_color('none')
    ax.add_patch(p)
    m.set_axes_limits(ax=ax)
    coll.set_clip_path(p)

-CHANGES

plt.show()

On Wed, Nov 4, 2009 at 1:17 AM, <Ross.Wilson@…2849…> wrote:

-----Original Message-----

From: Jeff Whitaker [mailto:jswhit@…1836…46…]

Sent: Tuesday, 3 November 2009 02:53

To: Stephane Raynaud

Cc: Wilson Ross; matplotlib-users@lists.sourceforge.net

Subject: Re: [Matplotlib-users] basemap: Mask the ocean [SEC=UNCLASSIFIED]

Stephane Raynaud wrote:

Ross,

one way is to mask (or remove) ocean points using the _geoslib module

provided with basemap.

When you create a Basemap instance, you can retrieve all its polygons

land (continents and islands) with “mymap.coastpolygons”.

Thay are stored as numpy arrays, and you can convert them to

_geoslib.Polygon objects :

poly = _geoslib.Polygon(N.asarray(coastalpoly).T)

Then you loop over all Polygons and all (x,y) points and test :

good_point = _geoslib.Point((x,y)).within(poly)

Thanks to this method, you can choose you optimal resolution.

You can even compute the intersection of you hexagons with coastal

polygons using .intersection() and .area (instead of simply checking

if the center is inside) and then reject points depending the fraction

of the cell covered by land (or ocean).

Following Stephane’s excellent suggestion, here’s a prototype Basemap

method that checks to see if a point is on land or over water. Ross -

if you find it useful I’ll include it in the next release. Note that it

will be slow for lots of points or large map regions.

-Jeff


Yes, Stephane’s approach is nice and Jeff has nicely encapsulated the approach. I’ll put that into my bag of tricks!

However it doesn’t quite do what I want. My data does not have any points in the ocean; the hex binning creates hexagonal patches that extend out into the ocean. As a physicist I say “that’s a representation artifact” and leave it at that, but my end-users want that ‘bleed’ into the ocean removed. My argument that they are losing data falls on deaf ears.

Here’s an even more contrived example that improves on my poor previous attempt to explain the problem:


import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

ll_lat = -38.10 # extent of area of interest

ll_lon = 145.06

ur_lat = -38.00

ur_lon = 145.16

create data points all on land

fudge = ((ll_lon+0.05, ur_lat-0.00, 1000),

     (ll_lon+0.05,  ur_lat-0.01,   2000),

     (ll_lon+0.05,  ur_lat-0.02,   3000),

     (ll_lon+0.05,  ur_lat-0.03,   4000),

     (ll_lon+0.04,  ur_lat-0.025,  1000),

     (ll_lon+0.043, ur_lat-0.036,  10000),

     (ll_lon+0.047, ur_lat-0.041,  20000),

     (ll_lon+0.07,  ur_lat-0.07,   4000),

     (ll_lon+0.08,  ur_lat-0.08,   3000),

     (ll_lon+0.09,  ur_lat-0.09,   2000),

     (ll_lon+0.10,  ur_lat-0.10,   1000))

data = np.ones((3, len(fudge)))

for (i, (lon, lat, val)) in enumerate(fudge):

data[0,i] = lon

data[1,i] = lat

data[2,i] = val

plot the data

fig = plt.figure()

m = Basemap(projection=‘cyl’, llcrnrlat=ll_lat, urcrnrlat=ur_lat,

llcrnrlon=ll_lon, urcrnrlon=ur_lon, resolution=‘f’)

plt.hexbin(data[0,:], data[1,:], data[2,:], gridsize=5)

m.drawcoastlines(linewidth=0.5, color=‘k’)

plt.show()


None of the data points are on land, as can be seen if you change the ‘gridsize’ parameter to 100 or so (land is to the NE). However, when gridsize=5, the red hexagon in the middle extends out into the ocean. My users want that hexagon and others like it clipped at the coastline.

Crazy people, users.

Ross


Stephane Raynaud