Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix pressure_to_height and reverse inaccurate after 11km #2084

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix pressure_to_height and reverse inaccurate after 11km #2084

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions docs/api/references.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,10 @@ References
Evidence and climate significance. *Geophys. Res. Lett.*, **38**, L01706,
`doi:10.1029/2010GL045777 <https://doi.org/10.1029/2010GL045777>`_.

.. [Kraus2004] H. Kraus, 2004: Die Atmosphaere der Erde,
Springer, 470pp., Sections II.1.4. and II.6.1.2.,
doi:`10.1007/3-540-35017-9 <https://doi.org/10.1007/3-540-35017-9>`_.

.. [Lackmann2011] Lackmann, G., 2011: *Midlatitude Synoptic Meteorology*. Amer. Meteor. Soc.,
345 pp.

Expand Down
123 changes: 106 additions & 17 deletions src/metpy/calc/basic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,16 @@
p0 = units.Quantity(1013.25, 'hPa')
gamma = units.Quantity(6.5, 'K/km')

# Values according to the 1976 U.S. Standard atmosphere [NOAA1976]_.
# List of tuples (height, temperature, pressure, temperature gradient)
_STANDARD_ATMOSPHERE = list(
zip(units.Quantity([0, 11, 20, 32, 47, 51, 71], 'km'),
units.Quantity([288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65], 'K'),
units.Quantity([101325, 22632.1, 5474.89, 868.019, 110.906, 66.9389, 3.95642], 'Pa'),
units.Quantity([0.0065, 0, -0.001, -0.0028, 0, 0.0028, float('NAN')], 'K/m')))

_HEIGHT, _TEMPERATURE, _PRESSURE, _TEMPERATURE_GRADIENT = 0, 1, 2, 3


@exporter.export
@preprocess_and_wrap(wrap_like='u')
Expand Down Expand Up @@ -414,27 +424,67 @@
@preprocess_and_wrap(wrap_like='pressure')
@check_units('[pressure]')
def pressure_to_height_std(pressure):
r"""Convert pressure data to height using the U.S. standard atmosphere [NOAA1976]_.
r"""Convert pressure to height (km).

The implementation uses the formula outlined in [Hobbs1977]_ pg.60-61.
Conversion of pressure to height (km) with the
hydrostatic equation, according to the profile of the
1976 U.S. Standard atmosphere [NOAA1976]_.
Reference [Kraus2004]_.

Parameters
----------
pressure : `pint.Quantity`
pressure : `pint.Quantity` or `xarray.DataArray`
Atmospheric pressure

Returns
-------
`pint.Quantity`
Corresponding height value(s)
`pint.Quantity` or `xarray.DataArray`
Corresponding height value(s) (kilometers)

Notes
-----
.. math:: Z = \frac{T_0}{\Gamma}[1-\frac{p}{p_0}^\frac{R\Gamma}{g}]

.. math:: Z = \begin{cases}
Z_0 + \frac{T_0 - T_0 \cdot \exp\left(\frac{\Gamma \cdot R}
{g\cdot\log(\frac{p}{p0})}\right)}{\Gamma}&\Gamma \neq 0
\\Z_0 - \frac{R \cdot T_0}{g \cdot \log(\frac{p}{p_0})} &\text{else}
\end{cases}
"""
return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(
mpconsts.Rd * gamma / mpconsts.g))
is_array = hasattr(pressure.magnitude, '__len__')
if not is_array:
pressure = units.Quantity([pressure.magnitude], pressure.units)

# Initialize the return array.
is_dask = not hasattr(pressure.magnitude, 'fill')
if not is_dask:
z = units.Quantity(np.full_like(pressure, np.nan), 'km')
else:
# full_like for Dask seems to be broken, use an alternative means of filling with nan
# TypeError: ones_like() got an unexpected keyword argument 'subok'
z = units.Quantity(np.empty_like(pressure), 'km')
z += np.nan

for i, ((z0, t0, p0, gamma), (_z1, _t1, p1, _)) in enumerate(zip(
_STANDARD_ATMOSPHERE[:-1], _STANDARD_ATMOSPHERE[1:])):
p1 = _STANDARD_ATMOSPHERE[i + 1][_PRESSURE]
indices = (pressure > p1) & (pressure <= p0)

if i == 0:
indices |= (pressure >= p0)
if is_dask:
indices = indices.compute()

if gamma != 0:
z[indices] = (z0 + 1. / gamma * (
t0 - t0 * np.exp(gamma * mpconsts.Rd / mpconsts.g * np.log(
pressure[indices] / p0)))).to(units.km)
else:
z[indices] = (z0 - (mpconsts.Rd * t0) / mpconsts.g * np.log(pressure[indices] / p0
)).to(units.km)

if np.isnan(z).any():
raise ValueError('Height to pressure conversion not implemented for z > 71km')

Check warning on line 485 in src/metpy/calc/basic.py

View check run for this annotation

Codecov / codecov/patch

src/metpy/calc/basic.py#L485 

Added line #L485 was not covered by tests

return z if is_array else z[0]


@exporter.export
Expand Down Expand Up @@ -564,26 +614,65 @@
@preprocess_and_wrap(wrap_like='height')
@check_units('[length]')
def height_to_pressure_std(height):
r"""Convert height data to pressures using the U.S. standard atmosphere [NOAA1976]_.
r"""Convert height to pressure (hPa).

The implementation inverts the formula outlined in [Hobbs1977]_ pg.60-61.
Conversion of height to pressure (hPa) with the
hydrostatic equation, according to the profile of the
1976 U.S. Standard atmosphere [NOAA1976]_.
Reference [Kraus2004]_.

Parameters
----------
height : `pint.Quantity`
height : `pint.Quantity` or `xarray.DataArray`
Atmospheric height

Returns
-------
`pint.Quantity`
Corresponding pressure value(s)
`pint.Quantity` or `xarray.DataArray`
Corresponding pressure value(s) (hPa)

Notes
-----
.. math:: p = p_0 e^{\frac{g}{R \Gamma} \text{ln}(1-\frac{Z \Gamma}{T_0})}

.. math:: p = \begin{cases}
p_0 \cdot \left[\frac{T_0 - \Gamma \cdot (Z - Z_0)}{T_0}\right]^
{\frac{g}{\Gamma \cdot R}} &\Gamma \neq 0
\\p_0 \cdot \exp\left(\frac{-g \cdot (Z - Z_0)}{R \cdot T_0}\right) &\text{else}
\end{cases}
"""
return p0 * (1 - (gamma / t0) * height) ** (mpconsts.g / (mpconsts.Rd * gamma))
is_array = hasattr(height.magnitude, '__len__')
if not is_array:
height = units.Quantity([height.magnitude], height.units)

# Initialize the return array.
is_dask = not hasattr(height.magnitude, 'fill')
if not is_dask:
p = units.Quantity(np.full_like(height, np.nan), 'hPa')
else:
# full_like for Dask seems to be broken, use an alternative means of filling with nan
# TypeError: ones_like() got an unexpected keyword argument 'subok'
p = units.Quantity(np.empty_like(height), 'hPa')
p += np.nan

for i, ((z0, t0, p0, gamma), (z1, _t1, _p1, _)) in enumerate(zip(
_STANDARD_ATMOSPHERE[:-1], _STANDARD_ATMOSPHERE[1:])):
indices = (height >= z0) & (height < z1)

if i == 0:
indices |= height < z0
if is_dask:
indices = indices.compute()

if gamma != 0:
p[indices] = (p0 * ((t0 - gamma * (height[indices] - z0)) / t0) ** (
mpconsts.g / (gamma * mpconsts.Rd))).to(units.hPa)
else:
p[indices] = (p0 * np.exp(-mpconsts.g * (height[indices] - z0) / (mpconsts.Rd * t0)
)).to(units.hPa)

if np.isnan(p).any():
raise ValueError('Height to pressure conversion not implemented for z > 71km')

Check warning on line 673 in src/metpy/calc/basic.py

View check run for this annotation

Codecov / codecov/patch

src/metpy/calc/basic.py#L673 

Added line #L673 was not covered by tests

return p if is_array else p[0]


@exporter.export
Expand Down
8 changes: 4 additions & 4 deletions tests/calc/test_basic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -328,7 +328,7 @@ def test_pressure_to_heights_basic(array_type):
mask = [False, True, False, True]
pressures = array_type([975.2, 987.5, 956., 943.], 'mbar', mask=mask)
heights = pressure_to_height_std(pressures)
values = array_type([321.5, 216.5, 487.6, 601.7], 'meter', mask=mask)
values = array_type([321.6579, 216.5843, 487.8407, 601.9013], 'meter', mask=mask)
assert_array_almost_equal(heights, values, 1)


Expand All @@ -343,7 +343,7 @@ def test_heights_to_pressure_basic(array_type):

def test_pressure_to_heights_units():
"""Test that passing non-mbar units works."""
assert_almost_equal(pressure_to_height_std(29 * units.inHg), 262.8498 * units.meter, 3)
assert_almost_equal(pressure_to_height_std(29 * units.inHg), 262.9867 * units.meter, 3)


def test_coriolis_force(array_type):
Expand All @@ -362,7 +362,7 @@ def test_add_height_to_pressure(array_type):
pressure_in = array_type([1000., 900., 800.], 'hPa', mask=mask)
height = array_type([877.17421094, 500., 300.], 'meter', mask=mask)
pressure_out = add_height_to_pressure(pressure_in, height)
truth = array_type([900., 846.725, 770.666], 'hPa', mask=mask)
truth = array_type([900.0464, 846.7529, 770.6813], 'hPa', mask=mask)
assert_array_almost_equal(pressure_out, truth, 2)


Expand All @@ -372,7 +372,7 @@ def test_add_pressure_to_height(array_type):
height_in = array_type([110.8286757, 250., 500.], 'meter', mask=mask)
pressure = array_type([100., 200., 300.], 'hPa', mask=mask)
height_out = add_pressure_to_height(height_in, pressure)
truth = array_type([987.971601, 2114.957, 3534.348], 'meter', mask=mask)
truth = array_type([988.4233, 2115.9024, 3535.8360], 'meter', mask=mask)
assert_array_almost_equal(height_out, truth, 3)


Expand Down
64 changes: 32 additions & 32 deletions tests/calc/test_calc_tools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -217,39 +217,39 @@ def get_bounds_data():

@pytest.mark.parametrize('pressure, bound, hgts, interp, expected', [
(get_bounds_data()[0], 900 * units.hPa, None, True,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 900 * units.hPa, None, False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 870 * units.hPa, None, True,
(870 * units.hPa, 1.2665298 * units.kilometer)),
(870 * units.hPa, 1.2671 * units.kilometer)),
(get_bounds_data()[0], 870 * units.hPa, None, False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(get_bounds_data()[0], 0.9880028 * units.kilometer, None, True,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(get_bounds_data()[0], 0.9880028 * units.kilometer, None, False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(get_bounds_data()[0], 1.2665298 * units.kilometer, None, True,
(870 * units.hPa, 1.2665298 * units.kilometer)),
(get_bounds_data()[0], 1.2665298 * units.kilometer, None, False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 0.9885 * units.kilometer, None, True,
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 0.9885 * units.kilometer, None, False,
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 1.2671 * units.kilometer, None, True,
(870 * units.hPa, 1.2671 * units.kilometer)),
(get_bounds_data()[0], 1.2671 * units.kilometer, None, False,
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 900 * units.hPa, get_bounds_data()[1], True,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 900 * units.hPa, get_bounds_data()[1], False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 870 * units.hPa, get_bounds_data()[1], True,
(870 * units.hPa, 1.2643214 * units.kilometer)),
(870 * units.hPa, 1.2649 * units.kilometer)),
(get_bounds_data()[0], 870 * units.hPa, get_bounds_data()[1], False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(get_bounds_data()[0], 0.9880028 * units.kilometer, get_bounds_data()[1], True,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(get_bounds_data()[0], 0.9880028 * units.kilometer, get_bounds_data()[1], False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(get_bounds_data()[0], 1.2665298 * units.kilometer, get_bounds_data()[1], True,
(870.9869087 * units.hPa, 1.2665298 * units.kilometer)),
(get_bounds_data()[0], 1.2665298 * units.kilometer, get_bounds_data()[1], False,
(900 * units.hPa, 0.9880028 * units.kilometer)),
(get_bounds_data()[0], 0.98800289 * units.kilometer, get_bounds_data()[1], True,
(900 * units.hPa, 0.9880028 * units.kilometer))
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 0.9885 * units.kilometer, get_bounds_data()[1], True,
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 0.9885 * units.kilometer, get_bounds_data()[1], False,
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 1.2671 * units.kilometer, get_bounds_data()[1], True,
(870.9869087 * units.hPa, 1.2671 * units.kilometer)),
(get_bounds_data()[0], 1.2671 * units.kilometer, get_bounds_data()[1], False,
(900 * units.hPa, 0.9885 * units.kilometer)),
(get_bounds_data()[0], 0.98859 * units.kilometer, get_bounds_data()[1], True,
(900 * units.hPa, 0.9885 * units.kilometer))
])
def test_get_bound_pressure_height(pressure, bound, hgts, interp, expected):
"""Test getting bounds in layers with various parameter combinations."""
Expand Down Expand Up @@ -317,9 +317,9 @@ def test_get_layer_float32_no_heights():
p_l, u_l, v_l = get_layer(p, u, v, depth=1000 * units.meter)
assert_array_equal(p_l[:-1], p[:-1])
assert_array_almost_equal(u_l[:-1], u[:-1], 7)
assert_almost_equal(u_l[-1], 3.0455916 * units('m/s'), 4)
assert_almost_equal(u_l[-1], 3.03666 * units('m/s'), 4)
assert_array_almost_equal(v_l[:-1], v[:-1], 7)
assert_almost_equal(v_l[-1], 20.2149378 * units('m/s'), 4)
assert_almost_equal(v_l[-1], 20.21832 * units('m/s'), 4)
assert p_l.dtype == p.dtype
assert u_l.dtype == u.dtype
assert v_l.dtype == v.dtype
Expand Down Expand Up @@ -355,8 +355,8 @@ def layer_test_data():
(layer_test_data()[0], layer_test_data()[1], None, None, 150 * units.hPa, False,
(np.array([1000, 900]) * units.hPa, np.array([25.0, 16.666666]) * units.degC)),
(layer_test_data()[0], layer_test_data()[1], None, 2 * units.km, 3 * units.km, True,
(np.array([794.85264282, 700., 600., 540.01696548]) * units.hPa,
np.array([7.93049516, 0., -8.33333333, -13.14758845]) * units.degC))
(np.array([794.9484, 700., 600., 540.1925]) * units.hPa,
np.array([7.938, 0., -8.33333333, -13.133]) * units.degC))
])
def test_get_layer(pressure, variable, heights, bottom, depth, interp, expected):
"""Test get_layer functionality."""
Expand Down Expand Up @@ -385,8 +385,8 @@ def test_get_layer_masked():
p = units.Quantity(np.ma.array([1000, 500, 400]), 'hPa')
u = units.Quantity(np.arange(3), 'm/s')
p_layer, u_layer = get_layer(p, u, depth=units.Quantity(6000, 'm'))
true_p_layer = units.Quantity([1000., 500., 464.4742], 'hPa')
true_u_layer = units.Quantity([0., 1., 1.3303], 'm/s')
true_p_layer = units.Quantity([1000., 500., 464.6738], 'hPa')
true_u_layer = units.Quantity([0., 1., 1.3284], 'm/s')
assert_array_almost_equal(p_layer, true_p_layer, 4)
assert_array_almost_equal(u_layer, true_u_layer, 4)

Expand Down