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.wil...@ga.gov.au> wrote:
>
>
> -----Original Message-----
> From: Jeff Whitaker [mailto:jsw...@fastmail.fm]
> 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
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
