From 3d99a3444d7140a97883c8cd002775ea4710d576 Mon Sep 17 00:00:00 2001 From: Marco Fronzi Date: Mon, 10 Aug 2026 16:17:23 +1000 Subject: [PATCH] Add altimeter_from_station_pressure and sea_level_pressure_from_station_pressure Adds the analytic inverse of altimeter_to_station_pressure and exposes the hypsometric reduction from station pressure to sea-level pressure. Fixes #3369 --- docs/_templates/overrides/metpy.calc.rst | 2 + src/metpy/calc/basic.py | 112 +++++++++++++++++++++++ tests/calc/test_basic.py | 68 ++++++++++++-- 3 files changed, 176 insertions(+), 6 deletions(-) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index 926a82b58f3..02a1d62d6f5 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -190,10 +190,12 @@ Standard Atmosphere .. autosummary:: :toctree: ./ + altimeter_from_station_pressure altimeter_to_sea_level_pressure altimeter_to_station_pressure height_to_pressure_std pressure_to_height_std + sea_level_pressure_from_station_pressure Smoothing --------- diff --git a/src/metpy/calc/basic.py b/src/metpy/calc/basic.py index 7f516b03a9b..74dfe0e1006 100644 --- a/src/metpy/calc/basic.py +++ b/src/metpy/calc/basic.py @@ -1296,6 +1296,118 @@ def altimeter_to_sea_level_pressure(altimeter_value, height, temperature): return psfc * np.exp(height / h) +@exporter.export +@preprocess_and_wrap(wrap_like='pressure') +@check_units('[pressure]', '[length]') +def altimeter_from_station_pressure(pressure, height): + r"""Convert station pressure to an altimeter setting. + + This is the inverse of `altimeter_to_station_pressure`. It is useful when a station + reports its actual (station) pressure but an altimeter setting is required, e.g. for + comparison against METAR observations. The definitions of altimeter setting and station + pressure are taken from [Smithsonian1951]_. A standard atmosphere [NOAA1976]_ is assumed. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure at the station elevation + + height : `pint.Quantity` + Elevation of the station measuring pressure + + Returns + ------- + `pint.Quantity` + The altimeter setting + + Examples + -------- + >>> from metpy.calc import altimeter_from_station_pressure + >>> from metpy.units import units + >>> altimeter_from_station_pressure(900. * units.hPa, 1000. * units.m) + + + See Also + -------- + altimeter_to_station_pressure, sea_level_pressure_from_station_pressure + + Notes + ----- + Inverting Equation 1 and Equation 3 of the Smithsonian Handbook (1951) p. 269, as + implemented in `altimeter_to_station_pressure`, gives + + .. math:: A_{mb} = \left[\left(p_{mb} - 0.3\right)^n + + \frac{p_{0}^n a H_{b}}{T_{0}}\right]^\frac{1}{n} + + where :math:`p_{0}` = 1013.25 mb is standard sea-level pressure, :math:`T_{0}` = 288 K is + standard sea-level temperature, :math:`a` is the standard atmosphere lapse rate + :math:`6.5^{\circ}C. km.^{-1}`, :math:`H_{b}` is the station elevation, and + :math:`n = \frac{a R_{d}}{g} = 0.190284`. + + """ + # N-Value + n = (mpconsts.Rd * gamma / mpconsts.g).to_base_units() + + return ((pressure - units.Quantity(0.3, 'hPa')) ** n + + (p0.to(pressure.units) ** n * gamma * height) / t0) ** (1 / n) + + +@exporter.export +@preprocess_and_wrap(wrap_like='pressure') +@check_units('[pressure]', '[length]', '[temperature]') +def sea_level_pressure_from_station_pressure(pressure, height, temperature): + r"""Convert station pressure to sea-level pressure. + + Sea-level pressure is a pressure value obtained by the theoretical reduction of barometric + pressure to sea level. It is assumed that the atmosphere extends to sea level below the + station and that the properties of the atmosphere are related to conditions observed at + the station [Smithsonian1951]_. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure at the station elevation + + height : `pint.Quantity` + Elevation of the station measuring pressure + + temperature : `pint.Quantity` + Temperature at the station + + Returns + ------- + `pint.Quantity` + The sea-level pressure + + Examples + -------- + >>> from metpy.calc import sea_level_pressure_from_station_pressure + >>> from metpy.units import units + >>> sea_level_pressure_from_station_pressure(900. * units.hPa, 1000. * units.m, + ... 15. * units.degC) + + + See Also + -------- + altimeter_to_sea_level_pressure, altimeter_from_station_pressure + + Notes + ----- + This function is implemented using the following equations from Wallace and Hobbs (1977), + as described in more detail in `altimeter_to_sea_level_pressure`: + + .. math:: p_{sealevel} = p_{station} exp\left(\frac{\Delta z}{H}\right) + + where :math:`\Delta z` is the station elevation and :math:`H = \frac{R_{d}T}{g}` is the + scale height. + + """ + # Calculate the scale height + h = mpconsts.Rd * temperature / mpconsts.g + + return pressure * np.exp(height / h) + + def _check_radians(value, max_radians=2 * np.pi): """Input validation of values that could be in degrees instead of radians. diff --git a/tests/calc/test_basic.py b/tests/calc/test_basic.py index dc8770ef4d7..e557f0762f7 100644 --- a/tests/calc/test_basic.py +++ b/tests/calc/test_basic.py @@ -9,12 +9,14 @@ import xarray as xr from metpy.calc import (add_height_to_pressure, add_pressure_to_height, - altimeter_to_sea_level_pressure, altimeter_to_station_pressure, - apparent_temperature, coriolis_parameter, geopotential_to_height, - heat_index, height_to_geopotential, height_to_pressure_std, - pressure_to_height_std, sigma_to_pressure, smooth_circular, - smooth_gaussian, smooth_n_point, smooth_rectangular, smooth_window, - wind_components, wind_direction, wind_speed, windchill, zoom_xarray) + altimeter_from_station_pressure, altimeter_to_sea_level_pressure, + altimeter_to_station_pressure, apparent_temperature, + coriolis_parameter, geopotential_to_height, heat_index, + height_to_geopotential, height_to_pressure_std, pressure_to_height_std, + sea_level_pressure_from_station_pressure, sigma_to_pressure, + smooth_circular, smooth_gaussian, smooth_n_point, smooth_rectangular, + smooth_window, wind_components, wind_direction, wind_speed, windchill, + zoom_xarray) from metpy.cbook import get_test_data from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_array_equal from metpy.units import units @@ -814,6 +816,60 @@ def test_altimeter_to_sea_level_pressure_hpa(array_type): assert_array_almost_equal(res, truth, 3) +def test_altimeter_from_station_pressure_hpa(): + """Test converting station pressure to altimeter setting with hectopascals.""" + station = 900. * units.hPa + elev = 1000. * units.m + res = altimeter_from_station_pressure(station, elev) + truth = 1014.3665 * units.hectopascal + assert_almost_equal(res, truth, 3) + + +def test_altimeter_from_station_pressure_inhg(): + """Test converting station pressure to altimeter setting with inches of mercury.""" + station = 26.57 * units.inHg + elev = 1000. * units.m + res = altimeter_from_station_pressure(station, elev).to('inHg') + truth = 29.9465 * units.inHg + assert_almost_equal(res, truth, 3) + + +def test_altimeter_from_station_pressure_sea_level(): + """Test that at sea level the altimeter setting is the station pressure less 0.3 hPa.""" + res = altimeter_from_station_pressure(1000. * units.hPa, 0. * units.m) + truth = 999.7 * units.hectopascal + assert_almost_equal(res, truth, 4) + + +def test_altimeter_station_pressure_round_trip(array_type): + """Test that the altimeter and station pressure conversions invert each other.""" + mask = [False, True, False, True] + altim = array_type([965., 982., 1013.25, 1040.], 'hectopascal', mask=mask) + elev = 1000. * units.m + station = altimeter_to_station_pressure(altim, elev) + res = altimeter_from_station_pressure(station, elev) + assert_array_almost_equal(res, altim, 6) + + +def test_sea_level_pressure_from_station_pressure(): + """Test converting station pressure to sea level pressure.""" + res = sea_level_pressure_from_station_pressure(900. * units.hPa, 1000. * units.m, + 15. * units.degC) + truth = 1013.2898 * units.hectopascal + assert_almost_equal(res, truth, 3) + + +def test_sea_level_pressure_from_station_pressure_consistency(): + """Test consistency with the existing altimeter-based sea level pressure routine.""" + altim = 29.92 * units.inHg + elev = 1000. * units.m + temp = 15. * units.degC + station = altimeter_to_station_pressure(altim, elev) + res = sea_level_pressure_from_station_pressure(station, elev, temp).to('hPa') + truth = altimeter_to_sea_level_pressure(altim, elev, temp).to('hPa') + assert_almost_equal(res, truth, 8) + + def test_zoom_xarray(): """Test zoom_xarray on 2D DataArray.""" data = xr.open_dataset(get_test_data('GFS_test.nc', False))