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