Add altimeter and sea level pressure calculations from station pressure

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

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)
+    <Quantity(1014.36653, 'hectopascal')>
+
+    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)
+    <Quantity(1013.28978, 'hectopascal')>
+
+    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.
 

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