In this example we will show how to generate masks, including lans/sea masks

Notebook can be accessed here

Preparing the notebook

In [1]:
import requests
r = requests.get("https://cdat.llnl.gov/cdat/sample_data/clt.nc",stream=True)
with open("clt.nc","wb") as f:
    for chunk in r.iter_content(chunk_size=1024):
        if chunk:  # filter local_filename keep-alive new chunks
            f.write(chunk)

# and load data
import cdms2
f = cdms2.open("clt.nc")
clt = f("clt", time=slice(0,1), squeeze=1) # Get first month
u = f("u", level=slice(0,1), squeeze=1)
v = f("v", level=slice(0,1), squeeze=1)
clt = clt.regrid(u.getGrid(), regridTool="regrid2") # Put data on same grid

# computes wind speed
import MV2
speed = MV2.sqrt(u**2+v**2)
print "Max speed:", speed.max()
print "Mean speed:",speed.mean()
print "Min speed:",speed.min()

# Prepare graphics
import vcs
x=vcs.init()
Max speed: 68.91321
Mean speed: 16.25912330863402
Min speed: 0.061108682

Value based masks

In [2]:
# Let's mask out area where wind speed is greater than twice mean
mask = MV2.greater(speed,speed.mean()*2.)

# Let's apply this mask
clt_masked = MV2.masked_where(mask,clt)
x.plot(clt_masked)
/Users/doutriaux1/anaconda2/envs/cdat/lib/python2.7/site-packages/vtk/util/numpy_support.py:135: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.
  assert not numpy.issubdtype(z.dtype, complex), \
Out[2]:

Land-sea Masks

Generating a landsea mask on any grid

Conveniently CDAT can generate masks for you (for regular grids only).

The observed data set used here as the basis for creating realistic model land/sea masks was obtained from the U.S. Navy on a 1/6 degree longitude-latitude grid.

more on the technique used can be read here

In [3]:
import cdutil
mask = cdutil.generateLandSeaMask(clt)
mask2 = MV2.where(mask._mask,1.,mask) # Not needed for cdutil versions >= 2.12.2017.9.25
mask2.setAxisList(mask.getAxisList()) # Not needed for cdutil versions >= 2.12.2017.9.25
clt_masked = MV2.masked_where(mask2,clt)
x.clear()
x.plot(clt_masked)
Out[3]:

Surface type by region masks

CDAT also provide capabilities to mask regions. Original regions and their numbers come from EzGet

The function requires both a land/sea mask and a file reporting "regions", default "region" mask is as follow:

Regions tables is:

In [4]:
regions, guess = cdutil.generateSurfaceTypeByRegionMask(mask2*100., verbose=False)
Done : |                                       | 0.00Done : ##                                      | 4.76Done : ####                                    | 9.52Done : ######                                  | 14.29Done : ########                                | 19.Done : ##########                              | 23.81Done : ############                            | 28.57Done : ##############                          | 33.33Done : ################                        | 38.10Done : ##################                      | 42.Done : ####################                    | 47.62Done : #####################                   | 52.38Done : #######################                 | 57.14Done : #########################               | 61.90Done : ###########################             | 66.67Done : #############################           | 71.43Done : ###############################         | 76.19Done : #################################       | 80.95Done : ###################################     | 85.71Done : #####################################   | 90.48Done : ####################################### | 95.24Done : ########################################| 100.00%
In [5]:
# let's take a look
x.clear()
x.plot(regions)
Out[5]:
In [6]:
# Now let's extract the indian ocean which according to table are area 205 and 206
ind1 = MV2.equal(regions,205)
ind2 = MV2.equal(regions,206)
indian_ocean = MV2.logical_or(ind1,ind2)

clt_indian_ocean = MV2.masked_where(MV2.logical_not(indian_ocean),clt)
x.clear()
x.plot(clt_indian_ocean(longitude=(15,150),latitude=(-90,35)),ratio="autot")
Out[6]: