Add Mean Sea Level Pressure Reduction

docs/api/references.rst

@@ -142,6 +142,12 @@ References
 
 .. [NWS10201] 2017: `National Weather Service Instruction 10-201 <../_static/NWS_10-201.pdf>`_.
 
+
+.. [Pauley1998] Pauley, P., 1998: An example of uncertainty in sea level pressure reduction.
+           *Wea. Forecasting.*, 13, 833–850,
+           doi: `10.1175/1520-0434(1998)013%3C0833:AEOUIS%3E2.0.CO;2
+           <https://doi.org/10.1175/1520-0434(1998)013%3C0833:AEOUIS%3E2.0.CO;2>`_.
+
 .. [Philips1957] Philips, N. A., 1957: A coordinate system having some special
            advantages for numerical forecasting. *J. Meteor.*, **14**, 184-185,
            doi:`10.1175/1520-0469(1957)014%3C0184:ACSHSS%3E2.0.CO;2

src/metpy/calc/basic.py

@@ -1202,3 +1202,56 @@ def _check_radians(value, max_radians=2 * np.pi):
         warnings.warn('Input over {} radians. '
                       'Ensure proper units are given.'.format(np.nanmax(max_radians)))
     return value
+
+
+@exporter.export
+@preprocess_and_wrap(wrap_like='pressure')
+@check_units('[pressure]', '[temperature]', '[length]', '[pressure]')
+def mean_sea_level_pressure(pressure, temperature, geopotential_height, vapor_pressure):
+    """Reduce surface pressure to mean sea level.
+
+    Parameters
+    ----------
+    pressure : `pint.Quantity`
+        Station pressure.
+
+    temperature : `pint.Quantity`
+        Station temperature.
+
+    geopotential_height : `pint.Quantity`
+        Station geopotential height.
+
+    vapor_pressure : `pint.Quantity`
+        Station water vapor pressure.
+
+    Returns
+    -------
+    `pint.Quantity`
+        Pressure reduced to mean sea level.
+
+    Notes
+    -----
+    This function is an implementation of the method used in [NWS1963]_, Ch. 7. The humidity
+    correction term has been fitted to be continuous and the local station correction has been
+    ignored. See also [Pauley1998]_.
+    """
+    # define function for empirical humidity correction
+    def humidity_correction(height):
+        """Calculate humidity correction term.
+
+        Fitted from table of empirically derived values found in NWS Manual of Barometry.
+        """
+        const = np.array([5.04969979e-6, 1.30236277, 1.96926239e-1])
+        return (const[0] * (height.to('meter').m**const[1]) + const[2]) * units('degF / hPa')
+
+    # define constants
+    k = 29.28980 * units('meter / kelvin')
+    lapse_rate = 6.5 * units('kelvin / kilometer')
+
+    # calculate mean temperature for extrapolation
+    mean_temp = temperature.to('kelvin') + \
+        ((lapse_rate * geopotential_height) / 2.).to('kelvin') + \
+        (vapor_pressure * humidity_correction(geopotential_height)).to('kelvin')
+
+    # calculate mslp from hypsometric equation
+    return pressure * np.exp(geopotential_height / (k * mean_temp))

tests/calc/test_basic.py

@@ -12,9 +12,9 @@
                         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)
+                        mean_sea_level_pressure, 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)
 from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_array_equal)
 from metpy.units import units
 
@@ -796,3 +796,23 @@ def test_altimeter_to_sea_level_pressure_hpa():
     res = altimeter_to_sea_level_pressure(altim, elev, temp)
     truth = 1016.246 * units.hectopascal
     assert_almost_equal(res, truth, 3)
+
+
+def test_mean_sea_level_pressure():
+    """Test mean sea level pressure."""
+    p = 1013.25 * units('hPa')
+    T = 20. * units.degC
+    es = 3 * units('hPa')
+    z = 0. * units.meter
+    p0 = mean_sea_level_pressure(p, T, z, es)
+    assert_almost_equal(p, p0, 7)
+
+
+def test_mean_sea_level_pressure_1000m():
+    """Test mean sea level pressure with height above sea level."""
+    p = 920. * units('hPa')
+    T = 20. * units.degC
+    es = 1 * units('hPa')
+    z = 3270. * units.feet
+    p0 = mean_sea_level_pressure(p, T, z, es)
+    assert_almost_equal(1011.9606088890467 * units('hPa'), p0, 7)