`atmqty`

Examples: Transforming `AtmQty`

Objects

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.

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:

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.

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.

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.

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