diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index 926a82b58f3..0ec2da72b4f 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -194,6 +194,8 @@ Standard Atmosphere altimeter_to_station_pressure height_to_pressure_std pressure_to_height_std + station_pressure_to_altimeter + station_to_sea_level_pressure Smoothing --------- diff --git a/src/metpy/calc/basic.py b/src/metpy/calc/basic.py index 7f516b03a9b..9675bb410f2 100644 --- a/src/metpy/calc/basic.py +++ b/src/metpy/calc/basic.py @@ -1223,6 +1223,87 @@ def altimeter_to_station_pressure(altimeter_value, height): + units.Quantity(0.3, 'hPa')) +@exporter.export +@preprocess_and_wrap(wrap_like='station_pressure') +@check_units('[pressure]', '[length]') +def station_pressure_to_altimeter(station_pressure, height): + r"""Convert station pressure to the altimeter setting. + + This function calculates the altimeter setting, which is the pressure value + to which an aircraft altimeter scale is set so that it indicates the + altitude above mean sea-level. This is the mathematical inverse of + altimeter_to_station_pressure. + + Parameters + ---------- + station_pressure : `pint.Quantity` + The atmospheric pressure at the designated station elevation + + height: `pint.Quantity` + Elevation of the station measuring pressure + + Returns + ------- + `pint.Quantity` + The altimeter setting value (in. Hg or hPa) + + Notes + ----- + This function is implemented by inverting the equations from the + Smithsonian Handbook (1951) p. 269 used in altimeter_to_station_pressure. + + The 0.3 hPa offset is subtracted from the station pressure before + performing the exponentiation to ensure mathematical consistency with + the original forward formula. + """ + # Calculate the Poisson constant n (approx. 0.190284) + # n = (Rd * gamma) / g + n = (mpconsts.Rd * gamma / mpconsts.g).to_base_units() + + # 1. Subtract the 0.3 hPa offset from station pressure first + # 2. Raise to power of n + # 3. Add the height correction term + # 4. Raise the entire result to 1/n + return (((station_pressure - units.Quantity(0.3, 'hPa')) ** n + + ((p0.to(station_pressure.units) ** n * gamma * height) / t0)) ** (1 / n)) + + +@exporter.export +@preprocess_and_wrap(wrap_like='station_pressure') +@check_units('[pressure]', '[length]', '[temperature]') +def station_to_sea_level_pressure(station_pressure, height, temperature): + r"""Convert station pressure to sea level pressure (MSLP). + + This function calculates the Mean Sea Level Pressure (MSLP) using the + station pressure, the station elevation, and the current temperature. + This follows the requirements for atmospheric pressure reduction. + + Parameters + ---------- + station_pressure : `pint.Quantity` + The atmospheric pressure measured at the station + + height : `pint.Quantity` + The elevation of the station + + temperature : `pint.Quantity` + The air temperature at the station + + Returns + ------- + `pint.Quantity` + The calculated sea level pressure + + Notes + ----- + The formula used is a variation of the hypsometric equation: + .. math:: P_{slp} = P_s \exp \left( \frac{g \cdot H}{R_d \cdot \bar{T}} \right) + """ + # Rd is the gas constant for dry air, g is acceleration due to gravity + # These are typically pulled from metpy.constants + return station_pressure * np.exp((mpconsts.g * height) / (mpconsts.Rd * temperature)) + + @exporter.export @preprocess_and_wrap(wrap_like='altimeter_value') @check_units('[pressure]', '[length]', '[temperature]') diff --git a/tests/calc/test_basic.py b/tests/calc/test_basic.py index dc8770ef4d7..a1300cf538a 100644 --- a/tests/calc/test_basic.py +++ b/tests/calc/test_basic.py @@ -14,6 +14,7 @@ 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, + station_pressure_to_altimeter, station_to_sea_level_pressure, 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 @@ -790,6 +791,59 @@ def test_altimeter_to_station_pressure_hpa(array_type): assert_array_almost_equal(res, truth, 3) +def test_station_pressure_to_altimeter_inhg(): + """Test the station pressure to altimeter function with inches of mercury.""" + station_pressure = 950.96498 * units.hectopascal + elev = 500 * units.m + res = station_pressure_to_altimeter(station_pressure, elev) + truth = 29.8 * units.inHg + assert_almost_equal(res, truth, 3) + + +def test_station_pressure_to_altimeter_hpa(array_type): + """Test the station pressure to altimeter function with hectopascals.""" + mask = [False, True, False, True] + station_pressure = array_type( + [784.262996, 838.651657, 896.037821, 954.639265], 'hectopascal', mask=mask + ) + elev = array_type([2000., 1500., 1000., 500.], 'meter', mask=mask) + res = station_pressure_to_altimeter(station_pressure, elev) + truth = array_type([1000., 1005., 1010., 1013.], 'hectopascal', mask=mask) + assert_array_almost_equal(res, truth, 3) + + +def test_station_to_sea_level_pressure_basic(): + """Test MSLP calculation with standard values.""" + p_station = 950.96498 * units.hPa + elev = 500 * units.m + temp = 30 * units.degC + # Based on the inverse of the existing test truth + expected_mslp = 1006.089 * units.hPa + + res = station_to_sea_level_pressure(p_station, elev, temp) + assert_almost_equal(res, expected_mslp, 2) + + +def test_altimeter_round_trip(array_type): + """Test that moving forward and backward returns the original value.""" + # Use a range of values to ensure the 0.3 hPa offset is handled correctly everywhere + altim_start = array_type([29.92, 30.00, 31.00], 'inHg') + elev = array_type([0, 1000, 2000], 'meter') + + # Forward then Backward + intermediate_p = altimeter_to_station_pressure(altim_start, elev) + final_altim = station_pressure_to_altimeter(intermediate_p, elev) + + assert_array_almost_equal(altim_start, final_altim, 4) + + +def test_station_pressure_to_altimeter_invalid_units(): + """Verify that the decorator catches incorrect unit types.""" + with pytest.raises(ValueError): + # Passing meters where a pressure unit is expected + station_pressure_to_altimeter(units.Quantity(100, 'm'), units.Quantity(100, 'm')) + + def test_altimiter_to_sea_level_pressure_inhg(): """Test the altimeter to sea level pressure function with inches of mercury.""" altim = 29.8 * units.inHg