Merge pull request #621 from jthielen/thickness-hydrostatic

Add a Thickness Calculation

metpy/calc/tests/test_thermo.py

@@ -19,7 +19,7 @@
                         saturation_mixing_ratio,
                         saturation_vapor_pressure,
                         specific_humidity_from_mixing_ratio,
-                        surface_based_cape_cin, vapor_pressure,
+                        surface_based_cape_cin, thickness_hydrostatic, vapor_pressure,
                         virtual_potential_temperature, virtual_temperature)
 from metpy.calc.thermo import _find_append_zero_crossings
 from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_nan
@@ -744,3 +744,39 @@ def test_moist_static_energy():
     """Tests the moist static energy calculation."""
     mse = moist_static_energy(1000 * units.m, 25 * units.degC, 0.012 * units.dimensionless)
     assert_almost_equal(mse, 339.4594 * units('kJ/kg'), 6)
+
+
+def test_thickness_hydrostatic():
+    """Tests the thickness calculation for a moist layer."""
+    pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa
+    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.degC
+    mixing = np.array([0.01458, 0.00209, 0.00224, 0.00240, 0.00256, 0.00010])
+    thickness = thickness_hydrostatic(pressure, temperature, mixing=mixing)
+    assert_almost_equal(thickness, 9892.07 * units.m, 2)
+
+
+def test_thickness_hydrostatic_subset():
+    """Tests the thickness calculation with a subset of the moist layer."""
+    pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa
+    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.degC
+    mixing = np.array([0.01458, 0.00209, 0.00224, 0.00240, 0.00256, 0.00010])
+    thickness = thickness_hydrostatic(pressure, temperature, mixing=mixing,
+                                      bottom=850 * units.hPa, depth=150 * units.hPa)
+    assert_almost_equal(thickness, 1630.81 * units.m, 2)
+
+
+def test_thickness_hydrostatic_isothermal():
+    """Tests the thickness calculation for a dry isothermal layer at 0 degC."""
+    pressure = np.arange(1000, 500 - 1e-10, -10) * units.hPa
+    temperature = np.zeros_like(pressure) * units.degC
+    thickness = thickness_hydrostatic(pressure, temperature)
+    assert_almost_equal(thickness, 5542.12 * units.m, 2)
+
+
+def test_thickness_hydrostatic_isothermal_subset():
+    """Tests the thickness calculation for a dry isothermal layer subset at 0 degC."""
+    pressure = np.arange(1000, 500 - 1e-10, -10) * units.hPa
+    temperature = np.zeros_like(pressure) * units.degC
+    thickness = thickness_hydrostatic(pressure, temperature, bottom=850 * units.hPa,
+                                      depth=350 * units.hPa)
+    assert_almost_equal(thickness, 4242.68 * units.m, 2)

metpy/calc/thermo.py

@@ -1508,3 +1508,75 @@ def moist_static_energy(heights, temperature, specific_humidity):
 
     """
     return (dry_static_energy(heights, temperature) + Lv * specific_humidity).to('kJ/kg')
+
+
+@exporter.export
+@check_units('[pressure]', '[temperature]')
+def thickness_hydrostatic(pressure, temperature, **kwargs):
+    r"""Calculate the thickness of a layer via the hypsometric equation.
+
+    This thickness calculation uses the pressure and temperature profiles (and optionally
+    mixing ratio) via the hypsometric equation with virtual temperature adjustment
+
+    .. math:: Z_2 - Z_1 = -\frac{R_d}{g} \int_{p_1}^{p_2} T_v d\ln p,
+
+    which is based off of Equation 3.24 in [Hobbs2006]_.
+
+    This assumes a hydrostatic atmosphere.
+
+    Layer bottom and depth specified in pressure.
+
+    Parameters
+    ----------
+    pressure : `pint.Quantity`
+        Atmospheric pressure profile
+    temperature : `pint.Quantity`
+        Atmospheric temperature profile
+    mixing : `pint.Quantity`, optional
+        Profile of dimensionless mass mixing ratio. If none is given, virtual temperature
+        is simply set to be the given temperature.
+    molecular_weight_ratio : `pint.Quantity` or float, optional
+        The ratio of the molecular weight of the constituent gas to that assumed
+        for air. Defaults to the ratio for water vapor to dry air.
+        (:math:`\epsilon\approx0.622`).
+    bottom : `pint.Quantity`, optional
+        The bottom of the layer in pressure. Defaults to the first observation.
+    depth : `pint.Quantity`, optional
+        The depth of the layer in hPa. Defaults to the full profile if bottom is not given,
+        and 100 hPa if bottom is given.
+
+    Returns
+    -------
+    `pint.Quantity`
+        The thickness of the layer in meters.
+
+    See Also
+    --------
+    pressure_to_height_std, virtual_temperature
+
+    """
+    mixing = kwargs.pop('mixing', None)
+    molecular_weight_ratio = kwargs.pop('molecular_weight_ratio', epsilon)
+    bottom = kwargs.pop('bottom', None)
+    depth = kwargs.pop('depth', None)
+
+    # Get the data for the layer, conditional upon bottom/depth being specified and mixing
+    # ratio being given
+    if bottom is None and depth is None:
+        if mixing is None:
+            layer_p, layer_virttemp = pressure, temperature
+        else:
+            layer_p = pressure
+            layer_virttemp = virtual_temperature(temperature, mixing, molecular_weight_ratio)
+    else:
+        if mixing is None:
+            layer_p, layer_virttemp = get_layer(pressure, temperature, bottom=bottom,
+                                                depth=depth)
+        else:
+            layer_p, layer_temp, layer_w = get_layer(pressure, temperature, mixing,
+                                                     bottom=bottom, depth=depth)
+            layer_virttemp = virtual_temperature(layer_temp, layer_w, molecular_weight_ratio)
+
+    # Take the integral (with unit handling) and return the result in meters
+    return (- Rd / g * np.trapz(layer_virttemp.to('K'), x=np.log(layer_p / units.hPa)) *
+            units.K).to('m')