Merge pull request #595 from jrleeman/precip_water_cleanup

Cleanup precipitable water

docs/references.rst

@@ -94,9 +94,6 @@ References
            composite and significant tornado parameters. Preprints, 22nd Conf. on Severe Local
            Storms, Hyannis, MA, Amer. Meteor. Soc.
 
-.. [Tsonis2008] Tsonis, A. A., 2008: An introduction to atmospheric thermodynamics.
-           Cambridge Univ. Press, Cambridge.
-
 .. [Ziv1994] Ziv, B., and Alpert, P., 1994: Isobaric to Isentropic Interpolation Errors
            and Implication to Potential Vorticity Analysis. *J. Appl. Meteor.*, **33**,
            694-703, doi:`10.1175/1520-0450(1994)033%3C0694:ITIIEA%3E2.0.CO;2

metpy/calc/indices.py

@@ -15,39 +15,48 @@
 
 @exporter.export
 @check_units('[temperature]', '[pressure]', '[pressure]')
-def precipitable_water(dewpt, pressure, top=400 * units('hPa')):
+def precipitable_water(dewpt, pressure, bottom=None, top=None):
     r"""Calculate precipitable water through the depth of a sounding.
 
-    Default layer depth is sfc-400 hPa. Formula used is:
+    Formula used is:
 
-    .. math:: \frac{1}{pg} \int\limits_0^d x \,dp
+    .. math::  -\frac{1}{\rho_l g} \int\limits_{p_\text{bottom}}^{p_\text{top}} r dp
 
-    from [Tsonis2008]_, p. 170.
+    from [Salby1996]_, p. 28.
 
     Parameters
     ----------
     dewpt : `pint.Quantity`
         Atmospheric dewpoint profile
     pressure : `pint.Quantity`
         Atmospheric pressure profile
+    bottom: `pint.Quantity`, optional
+        Bottom of the layer, specified in pressure. Defaults to None (highest pressure).
     top: `pint.Quantity`, optional
-        The top of the layer, specified in pressure. Defaults to 400 hPa.
+        The top of the layer, specified in pressure. Defaults to None (lowest pressure).
 
     Returns
     -------
     `pint.Quantity`
         The precipitable water in the layer
 
     """
-    sort_inds = np.argsort(pressure[::-1])
+    # Sort pressure and dewpoint to be in decreasing pressure order (increasing height)
+    sort_inds = np.argsort(pressure)[::-1]
     pressure = pressure[sort_inds]
     dewpt = dewpt[sort_inds]
 
-    pres_layer, dewpt_layer = get_layer(pressure, dewpt, depth=pressure[0] - top)
+    if top is None:
+        top = np.nanmin(pressure) * pressure.units
+
+    if bottom is None:
+        bottom = np.nanmax(pressure) * pressure.units
+
+    pres_layer, dewpt_layer = get_layer(pressure, dewpt, bottom=bottom, depth=bottom - top)
 
     w = mixing_ratio(saturation_vapor_pressure(dewpt_layer), pres_layer)
-    # Since pressure is in decreasing order, pw will be the negative of what we want.
-    # Thus the *-1
+
+    # Since pressure is in decreasing order, pw will be the opposite sign of that expected.
     pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units) /
                 (g * rho_l))
     return pw.to('millimeters')

metpy/calc/tests/test_indices.py

@@ -21,7 +21,20 @@ def test_precipitable_water():
     """Test precipitable water with observed sounding."""
     with UseSampleData():
         data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC', source='wyoming')
-    pw = precipitable_water(data.variables['dewpoint'][:], data.variables['pressure'][:])
+    pw = precipitable_water(data.variables['dewpoint'][:], data.variables['pressure'][:],
+                            top=400 * units.hPa)
+    truth = (0.8899441949243486 * units('inches')).to('millimeters')
+    assert_array_equal(pw, truth)
+
+
+def test_precipitable_water_no_bounds():
+    """Test precipitable water with observed sounding and no bounds given."""
+    with UseSampleData():
+        data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC', source='wyoming')
+    dewpoint = data.variables['dewpoint'][:]
+    pressure = data.variables['pressure'][:]
+    inds = pressure >= 400 * units.hPa
+    pw = precipitable_water(dewpoint[inds], pressure[inds])
     truth = (0.8899441949243486 * units('inches')).to('millimeters')
     assert_array_equal(pw, truth)
 

metpy/calc/tools.py

@@ -323,6 +323,10 @@ def _get_bound_pressure_height(pressure, bound, heights=None, interpolate=True):
         The bound pressure and height.
 
     """
+    sort_inds = np.argsort(pressure)[::-1]
+    pressure = pressure[sort_inds]
+    if heights is not None:
+        heights = heights[sort_inds]
     # Bound is given in pressure
     if bound.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
         # If the bound is in the pressure data, we know the pressure bound exactly
@@ -538,7 +542,7 @@ def get_layer(pressure, *args, **kwargs):
 
     # If the bottom is not specified, make it the surface pressure
     if bottom is None:
-        bottom = pressure[0]
+        bottom = np.nanmax(pressure) * pressure.units
 
     bottom_pressure, bottom_height = _get_bound_pressure_height(pressure, bottom,
                                                                 heights=heights,