atmqty Examples: Transforming AtmQty Objects

Description

This web page provides examples of using transform_AQobj to interpolate the variables in the domain of one AtmQty object to the domain described by another AtmQty object. In other words, transform_AQobj transforms one AtmQty object into another, by means of interpolation. The data used in these examples comes from the ECMWF ERA-40 reanalysis. To avoid presenting too many examples, results from only a subset of all the possible uses for transform_AQobj is shown. Additional output from examples illustrating other uses are provided; comparison to all the examples operates as a "unit test" of the function.

A link to the full source code needed and the data file used in these examples is at the bottom of this page. For more information on the transform_AQobj command, type help(transform_AQobj) or see the pydoc generated documentation.

Import Statements and Data

All examples below use the following import statements and data:

>>> import atmqty
>>> import cdms
>>> import MV
>>> import Numeric as N

The data is from a netCDF data file of fields from the ECMWF ERA-40 reanalysis for 15 Jan 2001 (0Z) on pressure levels. We read in the data using netCDF utilities from cdms and convert the arrays to Numeric arrays. Vectors lon and lat are longitude and latitude, respectively, in degrees, and lev is pressure levels, in hPa. We transpose the arrays so that the first index cycles through latitude, the second through longitude, and the third through levels. We also reverse the latitude order so that latitude increases with increasing index (in the dataset latitudes start from 90°N and progress south):

>>> f = cdms.open('era40_2001-01-15.nc')

>>> T = MV.filled(f('t', squeeze=1))
>>> u = MV.filled(f('u', squeeze=1))
>>> v = MV.filled(f('v', squeeze=1))
>>> q = MV.filled(f('q', squeeze=1))

>>> T_obj = f['t']
>>> lon = T_obj.getLongitude()[:]
>>> lat = T_obj.getLatitude()[::-1]
>>> lev = T_obj.getLevel()[:].astype(N.Float)

>>> T = N.transpose(T, axes=(1,2,0))[::-1,:,:]
>>> u = N.transpose(u, axes=(1,2,0))[::-1,:,:]
>>> v = N.transpose(v, axes=(1,2,0))[::-1,:,:]
>>> q = N.transpose(q, axes=(1,2,0))[::-1,:,:]

Temperature is in K, winds are in m/s, and specific humidity in kg/kg. The levels are at 250, 300, 400, 500, 600, 700, and 850 hPa. Longitude and latitude values are separated every 2.5 degrees.

We create an AtmQty object (x) for the above ECMWF domain and fill it with the above atmospheric quantities. Using those quantities, we compute potential temperature and altitude, storing it in the object:

>>> x = atmqty.AtmQty( lat, lon, lev \
...                  , lev_type="pressure", missing=1e+29 )
>>> x.add_qty({'T':T, 'u':u, 'v':v, 'q':q})
>>> x.theta()
>>> x.alt()

The zonal wind at 500 hPa is obtained by the AtmQty method get_qty:

>>> u500 = x.get_qty('u')[:,:,3]

and looks like:

The 500 hPa temperature looks like:

Examples

Horizontal and Pressure Coordinate Interpolation

Here we interpolate the quantities in the AtmQty object x onto a new horizontal domain (lat2, lon2) as well as a new set of pressure levels (lev2). The horizontal domain has longitude and latitude spacings of 4 degrees. Of the 4 new pressure levels, one of the levels (900 hPa) is outside the range of levels in the original dataset in object x:

>>> lev2 = N.array([900., 800., 500., 350.])
>>> lat2 = N.arange(20)*4.0 - 30.0
>>> lon2 = N.arange(50)*4.0 + 20.0

We create an AtmQty object (x2) to specify this new domain:

>>> x2 = atmqty.AtmQty( lat2, lon2, lev2 \
...                   , lev_type="pressure", missing=1e+29 )

but do not add any atmospheric quantities to it. We need to create x2 because transform_AQobj reads in the description of the target domain from an AtmQty object. We interpolate all values in x into the domain specified by x2 by:

>>> x2 = atmqty.transform_AQobj.transform_AQobj(x, x2)

The 500 hPa zonal wind in x2 then is:

This looks very similar to the 500 hPa zonal wind in x, except a little coarser. Note that all atmospheric quantities at 900 hPa (except pressure, of course) have the value of 1e+29, since transform_AQobj does not interpolate outside of the original domain.

Interpolation to Isentropic Levels

We specify 4 isentropic levels, at 300, 315, 330, and 350 K. We define a new AtmQty object with the same latitude-longitude grid as x:

>>> lev3 = N.array([300., 315., 330., 350.])
>>> x3 = atmqty.AtmQty( lat, lon, lev3 \
...                   , lev_type="isentropic", missing=1e+14 )

Using transform_AQobj we transform x into x3. After that, we also calculate pressure on the x3 isentropic surfaces:

>>> x3 = atmqty.transform_AQobj.transform_AQobj(x, x3)
>>> x3.press()

The zonal wind on the 315 K isentropic level is:

While the pressure on that (315 K) surface is:

As an example of the vertical interpolation, we examine a vertical profile at latitude 45 S, longitude 87.5 E:

The abscissa is temperature in K, the asterisk are the values in pressure coordinates from x, and the diamonds are the x3 values interpolated from x, on isentropic levels. Note that at this x3 location, the domain in x only successfully interpolates non-missing values to 3 of the isentropic levels.

Change the Missing Value Value

In the x3 example above, the missing value value is also changed from the value used in x (1e+29) to 1e+14, by setting the missing keyword in the x3 instantiation call to 1e+14. The transform_AQobj changes all the instances of missing values in x to this new value, and uses the new value if an element becomes undefined during the interpolation so that all missing values in x3 are designated by 1e+14.

Full Code and Data File For Examples

• test_trans.py: Note that the code fragments in the examples on this page given above are slightly different from the full code in this file, since the code in test_trans.py generates the plots shown.
• era40_2001-01-15.nc (1.5 Mb): Please see here for the restrictions ECMWF has as to how this data can be used.