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
On Mon, Nov 2, 2009 at 8:07 AM, <ross.wil...@ga.gov.au
<mailto:ross.wil...@ga.gov.au>> 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.whita...@noaa.gov
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
import numpy as np
from numpy import ma
from numpy.random import seed
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, _geoslib
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
class Basemap2(Basemap):
def is_land(self,xpt,ypt):
landpt = False
n = 0
for x,y in self.coastpolygons:
type = self.coastpolygontypes[n]
if type in [1,3]: # land or island in lake
b = np.asarray([x,y]).T
poly = _geoslib.Polygon(b)
landpt = _geoslib.Point((xpt,ypt)).within(poly)
if landpt: break
n = n + 1
return landpt
# plot the data
fig = plt.figure()
m = Basemap2(projection='cyl', llcrnrlat=ll_lat, urcrnrlat=ur_lat,
llcrnrlon=ll_lon, urcrnrlon=ur_lon, resolution='f',
suppress_ticks=False, area_thresh=0.5)
for npt in range(data.shape[1]):
if not m.is_land(data[0,npt],data[1,npt]):
data[:,npt] = 1.e30
data = ma.masked_values(data,1.e30)
plt.hexbin(data[0,:], data[1,:], data[2,:],zorder=2)
m.fillcontinents(color=(0.8,0.8,0.8,0))
m.drawcoastlines(linewidth=0.25, color='k')
plt.show()
------------------------------------------------------------------------------
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
