Calculating potential temperature using the pressure coordinate of a temperature cube

[Luis García Carreras]

Hi,

I am running into issues trying to convert a cube of air_temperature on pressure levels to potential temperature by using the pressure coordinate (the documentation for derived variables uses a separate pressure cube).

I am starting with a cube 'temperature':
<iris 'Cube' of air_temperature / (K) (p: 37; latitude: 5; longitude: 5)>

In order to get the units right, I tried doing:

pressure = temperature.coord('pressure')
pref = iris.coords.AuxCoord(p0, long_name='reference_pressure', units='hPa')
pref.convert_units(pressure.units)
theta = temperature * (pref/pressure)**0.28571429

This gives the error

/apps/python-libs-1.3/lib/iris/coords.pyc in __div__(self, other)
    670
    671     def __div__(self, other):
--> 672         return self.__binary_operator__(other, Coord._MODE_DIV)
    673
    674     def __truediv__(self, other):

/apps/python-libs-1.3/lib/iris/coords.pyc in __binary_operator__(self, other, mode_constant)
    627         if isinstance(other, Coord):
    628             raise iris.exceptions.NotYetImplementedError(
--> 629                 'coord %s coord' % Coord._MODE_SYMBOL[mode_constant])
    630
    631         elif isinstance(other, (int, float, np.number)):

NotYetImplementedError: coord / coord

If I just use the pressure data (i.e. pref.points / temperature.coord('pressure').points) then I get a broadcasting error. I can get around this by manually rebroadcasting this, but I know that multiplying the temperature by its pressure coordinate does automatically broadcast correctly (i.e. temperature * pressure works fine, even though temperature * pressure.points does not).

In summary, is there a way of using a DimCoord like a cube, in order to easily use cube coordinates to calculate other quantities?

[Andrew Dawson]

You should be able to do this calculation. You didn't say what dimensionality p0 has, I'm guessing it is a scalar, and the following assumes that case:

pressure = temperature.coord('pressure')

p0 = 1000.
p0_units = iris.unit.Unit('hPa')
pref = p0_units.convert(p0, pressure.units)
tmp = p0 / pressure
tmp.units = '1'
tmp.points = tmp.points ** 0.28571429
theta = temperature * tmp

* untested so there may be mistakes in here

This is obviously not as simple and nice as what you wanted to do, but I think it works around the issues present in the current version of iris. Firstly instead of using an AuxCoord for pref it just uses a scalar, first making sure it is expressed in the same units as the pressure coordinate. Secondly, you need to manually set the units of (pref / pressure) to '1' (dimensionless). Thirdly I have had to explicity reset the points values to raise them to the power of (R/cp) because iris does not (yet) support raising a coordinate to a power.

I think a few improvements to iris are in order to make this process simpler in future! Let us know how you get on.

Andrew

[Luis García Carreras]

Thanks for the suggestion, it does work (and you're right, p0 is just a scalar). I didn't think that you could convert units in a scalar, so making pref a coordinate does seem unnecessary.

The solution I had come up with before seeing your answer involved making tmp a simple numpy array, and broadcasting to match temperature:

pref = iris.coords.AuxCoord(1000, units='hPa')
pref.convert_units(pressure.units)
tmp = iris.util.broadcast_to_shape(pref / pressure.points, temperature.shape, temperature.coord_dims(pressure))
theta = temperature * (tmp ** 0.28571429)

I imagine both are pretty much equivalent. I thought it could be done more straightforwardly, but actually it hardly matters, just need to write the function once!

Are there any plans to have atmospheric conversions built in somehow (like potential temperature, but also different humidities, etc.)? I can see why it wouldn't be high in the priority list (a nice to have, rather than essential), but it would be very handy if it did exist.

Thanks again!

[Phil Elson]

> Are there any plans to have atmospheric conversions built in somehow (like potential temperature, but also different humidities, etc.)? I can see why it wouldn't be high in the priority list (a nice to have, rather than essential), but it would be very handy if it did exist.

The quick answer is no, but I'd love to have a separate package which did do these kinds of things.

I'd be happy to help set up a package, if others are interested in contributing though.

[Luis García Carreras]

Hi Phil,

sounds great. I'd be happy to help. I've written some code to do temperature/potential temperature conversions as well as between different measures of humidity (the latter don't rely on any iris magic, but at least the output has informative variable names and correct units). I'll try and tidy it up a bit and can email it to you and see what you think - I have no experience with publicly releasing code, so it might not be the prettiest. I'm away the next two weeks, so if I can't get it done this week it will have to be once I'm back.

Cheers!
Luis