atmqty Examples: Transforming AtmQty Objects (Additional)

Description and Import Statements

This web page provides additional examples of using transform_AQobj to interpolate the variables in the domain of one AtmQty object to the domain described by another AtmQty object. The examples here use synthetic data, the same data that are used in some of the docstring examples in transform_AQobj.

The basic examples page was kept limited to avoid there being too much clutter. That page, plus this one, when compared to function output, operate 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:

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

Synthetic Data

First we define a global domain and create the AtmQty object (x) for that domain. The vertical coordinate is pressure levels in hPa and latitude and longitude are in degrees:

>>> nx = 73
>>> ny = 37
>>> lon = N.arange(nx) * 5.0 - 180.0
>>> lat = N.arange(ny) * 5.0 - 90.0
>>> lev = N.array([850., 700., 500., 300., 250.])

>>> x = atmqty.AtmQty( lat, lon, lev \
...                  , lev_type="pressure", missing=1e+29 )

The synthetic data for temperature, zonal wind, and meridional wind is created with the following code. Temperature is in K and winds are in m/s. Note that we make some of the values missing:

>>> lat_2d = N.reshape( N.repeat(lat, nx), (ny, nx) )
>>> lon_2d = N.reshape( N.repeat(N.reshape(lon,(1,nx)), ny) \
...                   , (ny, nx) )

>>> u_2d = N.sin(lon_2d*N.pi/180.) * \
...        N.sin((lat_2d-90.)*N.pi/30.)
>>> v_2d = N.sin((lon_2d-30.)*N.pi/180.) * \
...        N.sin((lat_2d-90.)*N.pi/30.)
>>> N.putmask( u_2d, N.where(N.absolute(u_2d) < 1e-10, 1, 0), 0. )
>>> N.putmask( v_2d, N.where(N.absolute(v_2d) < 1e-10, 1, 0), 0. )

>>> T_2d = N.outerproduct( N.sin(lat*N.pi/360.) \
...                     , N.cos(lon*N.pi/360.) ) * 10. + 257.3
>>> T = N.zeros(x.qty_shape, typecode=N.Float)
>>> u = N.zeros(x.qty_shape, typecode=N.Float)
>>> v = N.zeros(x.qty_shape, typecode=N.Float)

>>> T[:,:,0] = T_2d
>>> u[:,:,0] = u_2d
>>> v[:,:,0] = v_2d

>>> for ilev in range(1,len(lev)):
...    T[:,:,ilev] = T[:,:,ilev-1] - N.sin(ilev*4.0/N.pi)*5.0
...    u[:,:,ilev] = u[:,:,ilev-1] + N.sin(ilev*4.0/N.pi)*4.0 + 3.0
...    v[:,:,ilev] = v[:,:,ilev-1] + N.sin(ilev*4.0/N.pi)*2.0

>>> T[23:30,17,3] = 1e+29
>>> u[10:18,23:32,1] = 1e+29
>>> v[17:24,62:70,2] = 1e+29

We add T, u, and v to the AtmQty object x and calculate potential temperature and altitude:

>>> x.add_qty({'T':T, 'u':u, 'v':v})
>>> x.theta()
>>> x.alt(T_surf=T_2d)

The temperature, zonal wind, and meridional wind at 500 hPa look like the following:

Additional Examples Using Synthetic Data

Horizontal Coordinate Interpolation To Smaller Sub-Domain

Here we interpolate the quantities in the AtmQty object x onto a new horizontal domain (lat2, lon2) while keeping the same pressure levels as x. The new AtmQty object is x2:

>>> lat2 = N.arange(40)*2.0 - 20.0
>>> lon2 = N.arange(50)*4.0 - 90.0
>>> x2 = atmqty.AtmQty( lat2, lon2, lev \
...                   , lev_type="pressure", missing=1e+29 )
>>> 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.

Horizontal Coordinate Interpolation To Larger Domain

Here we interpolate the quantities in the AtmQty object x2 onto the domain of x. The new AtmQty object is x2a:

>>> x2a = atmqty.AtmQty( lat, lon, lev \
...                    , lev_type="pressure", missing=1e+29 )
>>> x2a = atmqty.transform_AQobj.transform_AQobj(x2, x2a)

The 500 hPa zonal wind in x2a then is:

This illustrates how because the domain in x2 is a sub-set of the global domain, when you interpolate it to the global domain (when you create x2a), all x2a points outside the x2 domain are set to missing.

Horizontal and Pressure Coordinate Interpolation

Here we interpolate the quantities in the AtmQty object x onto the horizontal domain of x2 and a new set of pressure levels. The new AtmQty object is x3:

>>> lev4 = N.array([900., 800., 600., 500.])
>>> x3 = atmqty.AtmQty( lat2, lon2, lev4 \
...                   , lev_type="pressure", missing=1e+17 )
>>> x3 = atmqty.transform_AQobj.transform_AQobj(x, x3)

The 500 hPa zonal wind in x3 should look exactly the same as the equivalent field in x2:

In this example we also changed the missing value value from 1e+29 (in x) to 1e+17 (in x3).

Pressure to Isentropic Levels Interpolation

We interpolate x from pressure levels to the following isentropic levels. The new AtmQty object is x4:

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

We compare zonal wind in x (asterisk) and x4 (diamond) for a vertical section at latitude 20 N and longitude 75 E:

Full Code and Data File For Examples

• test_trans_add.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_add.py generates the plots shown.