Mask a variable on pressure levels where data are below surface

This tutorial shows how to mask data

In some case earth system model output are interpolated to som standard pressure levels, this can lead to data being generated bellow the ground

This tutorial shows how to mask these erroneous data. Data courtesy of Jerry Potter

The CDAT software was developed by LLNL. This tutorial was written by Charles Doutriaux. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Download the Jupyter Notebook

Prepare Notebook

Back to Top

In [1]:
from __future__ import print_function
import cdms2
import numpy
import vcs
import requests
import os

for filename in ['ps.nc','ta.nc']:
    if not os.path.exists(filename):
        r = requests.get("https://cdat.llnl.gov/cdat/sample_data/notebooks/{}".format(filename), stream=True)
        with open(filename,"wb") as f:
            for chunk in r.iter_content(chunk_size=1024):
                if chunk:  # filter local_filename keep-alive new chunks
                    f.write(chunk)

Load the data

Back to Top

In [2]:
PS_file = cdms2.open("ps.nc")
TA_file = cdms2.open("ta.nc")

ps = PS_file("ps")
ta = TA_file("ta")

print("Shaped: ps {} and ta {}".format(ps.shape,ta.shape))
Shaped: ps (1, 320, 640) and ta (1, 37, 145, 288)

Regrid surface pressure data to target grid

Back to Top

In [3]:
ps = ps.regrid(ta.getGrid())
print("ps new shape:",ps.shape)
/software/anaconda53/envs/cdms2/lib/python3.7/site-packages/cdms2/avariable.py:1347: Warning: 
avariable.regrid: We chose regridTool = esmf for you among the following choices:
   Tools ->    'regrid2' (old behavior)
               'esmf' (conserve, patch, linear) or
               'libcf' (linear)
  warnings.warn(message, Warning)
/software/anaconda53/envs/cdms2/lib/python3.7/site-packages/cdms2/avariable.py:1354: Warning: 
avariable.regrid: We chose regridMethod = linear for you among the following choices:
    'conserve' or 'linear' or 'patch'
  warnings.warn(message, Warning)
ps new shape: (1, 145, 288)

Mask Data Below Surface

Back to Top

In [4]:
# Loop through levels and mask where pressure is less than ps

levels = ta.getLevel()
print("Levels are:", levels[:])
for i,level in enumerate(levels):
    low = numpy.less(ps,level)
    ta[:,i] = numpy.ma.masked_where(low,ta[:,i])
Levels are: [100000.  97500.  95000.  92500.  90000.  87500.  85000.  82500.  80000.
  77500.  75000.  70000.  65000.  60000.  55000.  50000.  45000.  40000.
  35000.  30000.  25000.  22500.  20000.  17500.  15000.  12500.  10000.
   7000.   5000.   3000.   2000.   1000.    700.    500.    300.    200.
    100.]

Plot Data

Back to Top

In [5]:
x = vcs.init(bg=True,geometry=(600,400))
x.plot(ta[:,0])
Out[5]: