Add heights to returned parameters for LCL, LFC, EL and CCL functions

[rdtewx]

What should we add?

The LCL, LFC, EL and CCL functions return pressures and temperatures but not heights, which is useful information too. It would be great to return heights too. A short-term fix would be a documentation update to provide instruction on how to do that leveraging existing functions. Thanks! (BTW, I'm a Python and MetPy noob!)

Reference

No response

[dcamron]

We can't directly calculate height for these quantities without assuming density characteristics of the background atmospheric profile. So we don't make those assumptions for you. We do provide functions that can make those assumptions (or not) for you to explore this yourself, depending on your use-case and available data:

In [1]: from datetime import datetime, UTC
   ...:
   ...: import numpy as np
   ...:
   ...: from metpy.calc import lcl, pressure_to_height_std, thickness_hydrostatic
   ...: from metpy.interpolate import log_interpolate_1d
   ...: from metpy.units import units
   ...: from siphon.simplewebservice.wyoming import WyomingUpperAir
   ...:
   ...: dt = datetime(2011, 4, 27, 12)
   ...: stid = 'BMX'
   ...:
   ...: df = WyomingUpperAir.request_data(dt, stid)
   ...:
   ...: pressure = df['pressure'].values * units.hPa
   ...: height = df['height'].values * units.m
   ...: temperature = df['temperature'].values * units.degC
   ...: dewpoint = df['dewpoint'].values * units.degC
   ...:
   ...: lcl_pressure, lcl_temperature = lcl(pressure[0], temperature[0], dewpoint[0])
   ...: lcl_pressure, lcl_temperature
Out[1]:
(<Quantity(973.328979, 'hectopascal')>,
 <Quantity(15.7658454, 'degree_Celsius')>)

Given this LCL pressure and temperature, we have a few options. We can assume a NOAA Standard Atmosphere and use metpy.calc.pressure_to_height_std:

In [2]: pressure_to_height_std(lcl_pressure).to('m')
Out[2]: <Quantity(337.559412, 'meter')>

We can also estimate the total depth of this layer given more information about the profile by calculating the metpy.calc.thickness_hydrostatic:

In [3]: lcl_layer_pressures = np.append(pressure[pressure >= lcl_pressure], lcl_pressure)
   ...: lcl_layer_temperatures = np.append(temperature[pressure >= lcl_pressure], lcl_temperature)
   ...:
   ...: delta_z = thickness_hydrostatic(lcl_layer_pressures, lcl_layer_temperatures)
   ...: delta_z
Out[3]: <Quantity(126.737276, 'meter')>

Which seems very different. However, our pressure and temperature profiles don't extend to the surface. If we add this thickness to the known height at the beginning of our pressure profile...

In [4]: delta_z + height[0]
Out[4]: <Quantity(304.737276, 'meter')>

it looks pretty similar. Finally, because we do have a full height profile available to us in this case (not always true), we can just estimate the height of our LCL given its nearest pressure / height values via metpy.interpolate.log_interpolate_1d:

In [5]: log_interpolate_1d(lcl_pressure, pressure, height)
Out[5]: <Quantity([305.44692104], 'meter')>

You can see these functions listed and demonstrated in the Reference and Examples sections of the docs, notably sections on sounding calculations. We can always add more. Also, check out the MetPy Cookbook on ProjectPythia.org. I hope this helps.

I'll close this issue shortly unless you have a more specific documentation request or plan to open a documentation PR sometime yourself. Thanks!

[rdtewx]

Understood. I had been using metpy.calc.thickness_hydrostatic which is working well enough but is cumbersome to prep the values the function requires.

On the subject of thickness functions, could that same function be used to estimate/extrapolate the 1000 hPa height when the sounding surface level is above that level? I've been trying but can't get it to work.

[josin77]

This points to another possibility. There are cases when the input data to these functions do not have any pressure information at different heights (like in the case of Temperature and Relative humidity profiles from a Microwave Radiometer).

Right now, to find LCL from such an input data, one option is to estimate pressure levels from the data using metpy.calc.height_to_pressure_std() and then to use the LCL function. But the height_to_pressure function makes use of the hypsometric equation which uses the ideal gas and hydrostatic assumptions.

If we instead could have an LCL function that works based on heights (and not pressure levels), it could be helpful.
For this, we need a parcel_profile function which takes (heights, initial dew point temperature and initial temperature) as inputs.
The LCL could then be returned as the first height where (Tparcel - Tenvironment = 0), and likewise for the other functions mentioned in the original issue.

[dopplershift]

@josin77 The problem with your suggested approach is that to calculate the parcel profile directly as a function of height would require a version of moist_lapse that uses height instead of pressure. This would need an entirely different approach from the current implementation, and not one I'm sure is out there in the literature.

[josin77]

@dopplershift Can we proceed with this approach for the moist lapse rate? :

• Computation of moist lapse rate using
  $Γ_s = Γ_d\dfrac{[1 + (a ⋅ r_s / T )]}{ [1 + (b ⋅ r_s / T ^2)] }$ where $a = 8711 K, b = 1.35\times 10^7 K^2$, and $Γ_d = 9.8 Kkm^{–1}$ (1)
  across different heights. Here, $T$ and $r_s$ as functions of $z$ can be used.

For this, we need $r_s$. To get $r_s$, we can use $r_s = \epsilon \frac{e_s}{P}$ for which $e_s$ in turn can be obtained using metpy.calc.saturation_vapor_pressure() which just takes temperature at different heights as input.

(1) This is Eq.4.38b from Practical Meteorology: An Algebra-based Survey of Atmospheric Science, Version 1.02b, Roland Stull.