Menu
▾
▴
SF.net SVN: matplotlib:[6825] trunk/toolkits/basemap
|
From: <js...@us...> - 2009-01-25 14:29:37
|
Revision: 6825
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=6825&view=rev
Author: jswhit
Date: 2009-01-25 14:29:31 +0000 (Sun, 25 Jan 2009)
Log Message:
-----------
added maskoceans function and example.
Modified Paths:
--------------
trunk/toolkits/basemap/Changelog
trunk/toolkits/basemap/MANIFEST.in
trunk/toolkits/basemap/examples/README
trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
Added Paths:
-----------
trunk/toolkits/basemap/examples/maskoceans.py
Modified: trunk/toolkits/basemap/Changelog
===================================================================
--- trunk/toolkits/basemap/Changelog 2009-01-25 12:56:39 UTC (rev 6824)
+++ trunk/toolkits/basemap/Changelog 2009-01-25 14:29:31 UTC (rev 6825)
@@ -1,4 +1,5 @@
version 0.99.4 (not yet released)
+ * added 'maskoceans' function.
* update pupynere to version 1.0.8 (supports writing large files).
* added more informative error message in readshapefile when
one of the shapefile components can't be found.
Modified: trunk/toolkits/basemap/MANIFEST.in
===================================================================
--- trunk/toolkits/basemap/MANIFEST.in 2009-01-25 12:56:39 UTC (rev 6824)
+++ trunk/toolkits/basemap/MANIFEST.in 2009-01-25 14:29:31 UTC (rev 6825)
@@ -70,6 +70,7 @@
include examples/plotprecip.py
include examples/nws_precip_conus_20061222.nc
include examples/NetCDFFile_tst.py
+include examples/maskoceans.py
include examples/README
include lib/mpl_toolkits/__init__.py
include lib/mpl_toolkits/basemap/__init__.py
Modified: trunk/toolkits/basemap/examples/README
===================================================================
--- trunk/toolkits/basemap/examples/README 2009-01-25 12:56:39 UTC (rev 6824)
+++ trunk/toolkits/basemap/examples/README 2009-01-25 14:29:31 UTC (rev 6825)
@@ -123,3 +123,5 @@
save_background.py shows how to save a map background and reuse it in another
figure (without having to redraw coastlines).
+
+maskoceans.py shows how to mask 'wet' areas on a plot.
Added: trunk/toolkits/basemap/examples/maskoceans.py
===================================================================
--- trunk/toolkits/basemap/examples/maskoceans.py (rev 0)
+++ trunk/toolkits/basemap/examples/maskoceans.py 2009-01-25 14:29:31 UTC (rev 6825)
@@ -0,0 +1,46 @@
+from mpl_toolkits.basemap import Basemap, shiftgrid, maskoceans, interp
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.mlab as mlab
+
+topodatin = mlab.load('etopo20data.gz')
+lonsin = mlab.load('etopo20lons.gz')
+latsin = mlab.load('etopo20lats.gz')
+
+# shift data so lons go from -180 to 180 instead of 20 to 380.
+topoin,lons1 = shiftgrid(180.,topodatin,lonsin,start=False)
+lats1 = latsin
+
+fig=plt.figure()
+# setup basemap
+m=Basemap(resolution='l',projection='lcc',lon_0=-100,lat_0=40,width=8.e6,height=6.e6)
+lons, lats = np.meshgrid(lons1,lats1)
+x, y = m(lons, lats)
+# interpolate land/sea mask to topo grid, mask ocean values.
+topo = maskoceans(lons, lats, topoin, inlands=False)
+# make contour plot (ocean values will be masked)
+#CS=m.contourf(x,y,topo,np.arange(-300,3001,50),cmap=plt.cm.jet,extend='both')
+im=m.pcolormesh(x,y,topo,cmap=plt.cm.jet,vmin=-300,vmax=3000)
+# draw coastlines.
+m.drawcoastlines()
+plt.title('ETOPO data with marine areas masked (original grid)')
+
+fig=plt.figure()
+# interpolate topo data to higher resolution grid (to better match
+# the land/sea mask).
+nlats = 3*topoin.shape[0]
+nlons = 3*topoin.shape[1]
+lons = np.linspace(-180,180,nlons)
+lats = np.linspace(-90,90,nlats)
+lons, lats = np.meshgrid(lons, lats)
+topo = interp(topoin,lons1,lats1,lons,lats,order=1)
+x, y = m(lons, lats)
+# interpolate land/sea mask to topo grid, mask ocean values.
+topo = maskoceans(lons, lats, topo, inlands=False)
+# make contour plot (ocean values will be masked)
+#CS=m.contourf(x,y,topo,np.arange(-300,3001,50),cmap=plt.cm.jet,extend='both')
+im=m.pcolormesh(x,y,topo,cmap=plt.cm.jet,vmin=-300,vmax=3000)
+# draw coastlines.
+m.drawcoastlines()
+plt.title('ETOPO data with marine areas masked (data on finer grid)')
+plt.show()
Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py 2009-01-25 12:56:39 UTC (rev 6824)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py 2009-01-25 14:29:31 UTC (rev 6825)
@@ -8,6 +8,8 @@
:func:`interp`: bilinear interpolation between rectilinear grids.
+:func:`maskoceans`: mask 'wet' points of an input array.
+
:func:`shiftgrid`: shifts global lat/lon grids east or west.
:func:`addcyclic`: Add cyclic (wraparound) point in longitude.
@@ -3113,20 +3115,8 @@
# read in.
if self.lsmask is None:
# read in land/sea mask.
- lsmaskf = open(os.path.join(basemap_datadir,'5minmask.bin'),'rb')
- nlons = 4320; nlats = nlons/2
- delta = 360./float(nlons)
- lsmask = np.reshape(np.fromstring(lsmaskf.read(),np.uint8),(nlats,nlons))
- lsmask_lons = np.arange(-180,180.,delta)
- lsmask_lats = np.arange(-90.,90+0.5*delta,delta)
- # add cyclic point in longitude
- lsmask, lsmask_lons = addcyclic(lsmask, lsmask_lons)
- nlons = nlons + 1; nlats = nlats + 1
- # add North Pole point (assumed water)
- tmparr = np.zeros((nlats,nlons),lsmask.dtype)
- tmparr[0:nlats-1,0:nlons] = lsmask
- lsmask = tmparr
- lsmaskf.close()
+ lsmask_lons, lsmask_lats, lsmask = _readlsmask()
+
# instance variable lsmask is set on first invocation,
# it contains the land-sea mask interpolated to the native
# projection grid. Further calls to drawlsmask will not
@@ -3950,3 +3940,53 @@
"""
cdftime = netcdftime.utime(units,calendar=calendar)
return cdftime.date2num(dates)
+
+def maskoceans(lonsin,latsin,datain,inlands=False):
+ """
+ mask data (``datain``), defined on a grid with latitudes ``latsin``
+ longitudes ``lonsin`` so that points over water will not be plotted.
+
+ .. tabularcolumns:: |l|L|
+
+ ============== ====================================================
+ Arguments Description
+ ============== ====================================================
+ lonsin, latsin rank-2 arrays containing longitudes and latitudes of
+ grid.
+ datain rank-2 input array on grid defined by ``lonsin`` and
+ ``latsin``.
+ inlands if False, mask only ocean points. If True, mask
+ ocean points and points over inland water bodies.
+ Default False.
+ ============== ====================================================
+
+ returns a masked array the same shape as datain with "wet" points masked.
+ """
+ # read in land/sea mask.
+ lsmask_lons, lsmask_lats, lsmask = _readlsmask()
+ # nearest-neighbor interpolation to output grid.
+ lsmasko = interp(lsmask,lsmask_lons,lsmask_lats,lonsin,latsin,masked=True,order=0)
+ # mask input data.
+ if inlands: # mask inland water bodies.
+ mask = np.logical_or(lsmasko==0,lsmasko==2)
+ else: # mask just marine areas.
+ mask = lsmasko == 0
+ return ma.masked_array(datain,mask=mask)
+
+def _readlsmask():
+ # read in land/sea mask.
+ lsmaskf = open(os.path.join(basemap_datadir,'5minmask.bin'),'rb')
+ nlons = 4320; nlats = nlons/2
+ delta = 360./float(nlons)
+ lsmask = np.reshape(np.fromstring(lsmaskf.read(),np.uint8),(nlats,nlons))
+ lsmask_lons = np.arange(-180,180.,delta)
+ lsmask_lats = np.arange(-90.,90+0.5*delta,delta)
+ # add cyclic point in longitude
+ lsmask, lsmask_lons = addcyclic(lsmask, lsmask_lons)
+ nlons = nlons + 1; nlats = nlats + 1
+ # add North Pole point (assumed water)
+ tmparr = np.zeros((nlats,nlons),lsmask.dtype)
+ tmparr[0:nlats-1,0:nlons] = lsmask
+ lsmask = tmparr
+ lsmaskf.close()
+ return lsmask_lons, lsmask_lats, lsmask
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
