0

I want to calculate the layer thickness from ECMWF data.

What I have are the variables T (within one layer) and relhum (within one layer) and pressure (both at layer interfaces, and at layer midpoints).

I have trouble understanding the pressure argument to the metpy.calc.thickness_hydrostatic_from_relative_humidity function:

When I look at the formula given in the docstring,

$$ Z_2 - Z_1 = -\frac{R_d}{g} \int_{p_1}^{p_2} T_v d\ln p $$

it seems to me as if the function would return the LHS, i.e., $Z_2 - Z_1$.

However, it also seems to me as if the pressure, temperature, relative_humidity arguments all have to have the same dimensions.

I find this confusing: In order to get the thickness of the layer between $Z_1$ and $Z_2$, I would expect to have to input both pressures $p_1$ and $p_2$ as well. However, the temperatures are usually defined at full levels, i.e., I have one less temperature than I have pressures.

For example, in order to calculate the thickness of the bottom layer I would expect to give pressure at surface, pressure at top of layer, and temperature (and humidity) within the layer. But when I try to use the function like this, I get a 

ValueError: operands could not be broadcast together with shapes (361,1440,79) (361,1440,78)

Please help me understand how to use this function properly.

andreas-h
  • 10,679
  • 18
  • 60
  • 78

1 Answers1

1

However, it also seems to me as if the pressure, temperature, relative_humidity arguments all have to have the same dimensions.

[...]

For example, in order to calculate the thickness of the bottom layer I would expect to give pressure at surface, pressure at top of layer, and temperature (and humidity) within the layer.

Indeed, contrary to your expectation, the pressure, temperature, and relative_humidity variables all need to be defined on the same levels, since, as given MetPy's documentation, thickness_hydrostatic_from_relative_humidity operates on atmospheric profiles. Note that in MetPy's terminology, a "layer" is over any set of vertical levels, not just between two levels.

The need for data defined on common levels can also be seen in the calculation itself, since integration is performed via the trapezoidal rule. For example, in the minimal case of a layer with two levels of data, the formula used by MetPy

given integration via the trapezoidal rule, reduces to

In practice, if you indeed have "offset" or "staggered" levels of data, where your temperature and relative_humidity have a different vertical coordinate than your pressure, you will need to interpolation your data to be on common levels to use thickness_hydrostatic_from_relative_humidity (or any of MetPy's other profile calculations). A simple destaggering approach that may be applicable in your case, with a pressure array of shape (361,1440,79) and temperature array of shape (361,1440,78), could look like

pressure_destaggered = (pressure[..., :-1] + pressure[..., 1:]) / 2

There are two caveats with this approach, however:

  • thickness_hydrostatic_from_relative_humidity was only designed to work on (1D) profiles of atmospheric data, not 3D grids, so exercise caution with your result (specifying the optional bottom and depth arguments are likely to fail, and the default integration axis may or may not be correct)
  • this destaggering will not work if you only have two pressure levels/one temperature level, since MetPy needs at least two levels to define a layer. In this case, you will need to approximate the temperature/relative_humidity at the two pressure levels in some fashion (such as assuming constant temperature/relative_humidity between the two levels).
Jon Thielen
  • 479
  • 2
  • 7