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