|
| 1 | +--- |
| 2 | +layout: default |
| 3 | +title: Masking_Data_Land-Sea_masks Tutorial |
| 4 | +--- |
| 5 | + |
| 6 | +# Masking_Data_Land-Sea_masks Tutorial |
| 7 | +[download iPython Notebook](Masking_Data_Land-Sea_masks.ipynb) |
| 8 | + |
| 9 | +# In this example we will show how to generate masks, including lans/sea masks |
| 10 | + |
| 11 | +Notebook can be accessed [here](notebook.ipynb) |
| 12 | + |
| 13 | +## Preparing the notebook |
| 14 | + |
| 15 | + |
| 16 | +```python |
| 17 | +import requests |
| 18 | +r = requests.get("https://uvcdat.llnl.gov/cdat/sample_data/clt.nc",stream=True) |
| 19 | +with open("clt.nc","wb") as f: |
| 20 | + for chunk in r.iter_content(chunk_size=1024): |
| 21 | + if chunk: # filter local_filename keep-alive new chunks |
| 22 | + f.write(chunk) |
| 23 | + |
| 24 | +# and load data |
| 25 | +import cdms2 |
| 26 | +f = cdms2.open("clt.nc") |
| 27 | +clt = f("clt", time=slice(0,1), squeeze=1) # Get first month |
| 28 | +u = f("u", level=slice(0,1), squeeze=1) |
| 29 | +v = f("v", level=slice(0,1), squeeze=1) |
| 30 | +clt = clt.regrid(u.getGrid(), regridTool="regrid2") # Put data on same grid |
| 31 | + |
| 32 | +# computes wind speed |
| 33 | +import MV2 |
| 34 | +speed = MV2.sqrt(u**2+v**2) |
| 35 | +print "Max speed:", speed.max() |
| 36 | +print "Mean speed:",speed.mean() |
| 37 | +print "Min speed:",speed.min() |
| 38 | + |
| 39 | +# Prepare graphics |
| 40 | +import vcs |
| 41 | +x=vcs.init() |
| 42 | +``` |
| 43 | + |
| 44 | + Max speed: 68.9132 |
| 45 | + Mean speed: 16.2591233086 |
| 46 | + Min speed: 0.0611087 |
| 47 | + |
| 48 | + |
| 49 | +## Value based masks |
| 50 | + |
| 51 | + |
| 52 | +```python |
| 53 | +# Let's mask out area where wind speed is greater than twice mean |
| 54 | +mask = MV2.greater(speed,speed.mean()*2.) |
| 55 | + |
| 56 | +# Let's apply this mask |
| 57 | +clt_masked = MV2.masked_where(mask,clt) |
| 58 | +x.plot(clt_masked) |
| 59 | +``` |
| 60 | + |
| 61 | + /Users/doutriaux1/anaconda2/envs/2.12-nox/lib/python2.7/site-packages/vcs/VTKPlots.py:998: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future. |
| 62 | + Check the NumPy 1.11 release notes for more information. |
| 63 | + data[:] = numpy.ma.masked_invalid(data, numpy.nan) |
| 64 | + /Users/doutriaux1/anaconda2/envs/2.12-nox/lib/python2.7/site-packages/numpy/ma/core.py:6385: MaskedArrayFutureWarning: In the future the default for ma.maximum.reduce will be axis=0, not the current None, to match np.maximum.reduce. Explicitly pass 0 or None to silence this warning. |
| 65 | + return self.reduce(a) |
| 66 | + /Users/doutriaux1/anaconda2/envs/2.12-nox/lib/python2.7/site-packages/numpy/ma/core.py:6385: MaskedArrayFutureWarning: In the future the default for ma.minimum.reduce will be axis=0, not the current None, to match np.minimum.reduce. Explicitly pass 0 or None to silence this warning. |
| 67 | + return self.reduce(a) |
| 68 | + |
| 69 | + |
| 70 | + |
| 71 | + |
| 72 | + |
| 73 | + |
| 74 | + |
| 75 | + |
| 76 | + |
| 77 | +## Land-sea Masks |
| 78 | + |
| 79 | +### Generating a landsea mask on any grid |
| 80 | + |
| 81 | +Conveniently CDAT can generate masks for you (for regular grids only). |
| 82 | + |
| 83 | +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. |
| 84 | + |
| 85 | +more on the technique used can be read [here](https://www-pcmdi.llnl.gov/publications/pdf/58.pdf) |
| 86 | + |
| 87 | + |
| 88 | +```python |
| 89 | +import cdutil |
| 90 | +mask = cdutil.generateLandSeaMask(clt) |
| 91 | +mask2 = MV2.where(mask._mask,1.,mask) # Not needed for cdutil versions >= 2.12.2017.9.25 |
| 92 | +mask2.setAxisList(mask.getAxisList()) # Not needed for cdutil versions >= 2.12.2017.9.25 |
| 93 | +clt_masked = MV2.masked_where(mask2,clt) |
| 94 | +x.clear() |
| 95 | +x.plot(clt_masked) |
| 96 | +``` |
| 97 | + |
| 98 | + |
| 99 | + |
| 100 | + |
| 101 | + |
| 102 | + |
| 103 | + |
| 104 | + |
| 105 | +### Surface type by region masks |
| 106 | + |
| 107 | +CDAT also provide capabilities to mask regions. Original regions and their numbers come from [EzGet](http://github.com/uv-cdat/ezget) |
| 108 | + |
| 109 | +The function requires both a land/sea mask and a file reporting "regions", default "region" mask is as follow: |
| 110 | +<img src="colorgeog.png"> |
| 111 | + |
| 112 | +Regions tables is: |
| 113 | +<img src="table.png"> |
| 114 | + |
| 115 | + |
| 116 | + |
| 117 | +```python |
| 118 | +regions, guess = cdutil.generateSurfaceTypeByRegionMask(mask2*100., verbose=False) |
| 119 | +``` |
| 120 | + |
| 121 | + Done : | | 0.00Done : ## | 4.76Done : #### | 9.Done : ###### | 14.29Done : ######## | 19.Done : ########## | 23.Done : ############ | 28.Done : ############## | 33.33Done : ################ | 38.10Done : ################## | 42.86Done : #################### | 47.62Done : ##################### | 52.Done : ####################### | 57.Done : ######################### | 61.Done : ########################### | 66.Done : ############################# | 71.Done : ############################### | 76.Done : ################################# | 80.95Done : ################################### | 85.71Done : ##################################### | 90.Done : ####################################### | 95.24Done : ########################################| 100.00% |
| 122 | + |
| 123 | + |
| 124 | + |
| 125 | +```python |
| 126 | +# let's take a look |
| 127 | +x.clear() |
| 128 | +x.plot(regions) |
| 129 | +``` |
| 130 | + |
| 131 | + |
| 132 | + |
| 133 | + |
| 134 | + |
| 135 | + |
| 136 | + |
| 137 | + |
| 138 | + |
| 139 | +```python |
| 140 | +# Now let's extract the indian ocean which according to table are area 205 and 206 |
| 141 | +ind1 = MV2.equal(regions,205) |
| 142 | +ind2 = MV2.equal(regions,206) |
| 143 | +indian_ocean = MV2.logical_or(ind1,ind2) |
| 144 | + |
| 145 | +clt_indian_ocean = MV2.masked_where(MV2.logical_not(indian_ocean),clt) |
| 146 | +x.clear() |
| 147 | +x.plot(clt_indian_ocean(longitude=(15,150),latitude=(-90,35)),ratio="autot") |
| 148 | +``` |
| 149 | + |
| 150 | + |
| 151 | + |
| 152 | + |
| 153 | + |
| 154 | + |
| 155 | + |
0 commit comments