feat: implement altimeter and MSLP inverse functions with round-trip …

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
 ---------

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]')

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