Calculation fixes

docs/_templates/overrides/metpy.calc.rst

@@ -79,6 +79,7 @@ calc
       most_unstable_cape_cin
       most_unstable_parcel
       parcel_profile
+      parcel_profile_with_lcl
       significant_tornado
       storm_relative_helicity
       supercell_composite
@@ -126,7 +127,7 @@ calc
       brunt_vaisala_period
       friction_velocity
       tke
-   
+
 
    Mathematical Functions
    ----------------------
@@ -173,7 +174,7 @@ calc
       reduce_point_density
       resample_nn_1d
       smooth_gaussian
-      
+
 
    Deprecated
    ----------
@@ -182,7 +183,7 @@ calc
 
    .. autosummary::
       :toctree: ./
-      
+
       get_wind_components
       get_wind_dir
       get_wind_speed

metpy/calc/kinematics.py

@@ -571,8 +571,14 @@ def storm_relative_helicity(u, v, heights, depth, bottom=0 * units.m,
     int_layers = (storm_relative_u[1:] * storm_relative_v[:-1] -
                   storm_relative_u[:-1] * storm_relative_v[1:])
 
+    # Need to manually check for masked value because sum() on masked array with non-default
+    # mask will return a masked value rather than 0. See numpy/numpy#11736
     positive_srh = int_layers[int_layers.magnitude > 0.].sum()
+    if np.ma.is_masked(positive_srh):
+        positive_srh = 0.0 * units('meter**2 / second**2')
     negative_srh = int_layers[int_layers.magnitude < 0.].sum()
+    if np.ma.is_masked(negative_srh):
+        negative_srh = 0.0 * units('meter**2 / second**2')
 
     return (positive_srh.to('meter ** 2 / second ** 2'),
             negative_srh.to('meter ** 2 / second ** 2'),

metpy/calc/tests/test_kinematics.py

@@ -443,6 +443,23 @@ def test_storm_relative_helicity_agl():
     assert_almost_equal(total_srh, 300. * units('meter ** 2 / second ** 2 '), 6)
 
 
+def test_storm_relative_helicity_masked():
+    """Test that srh does not return masked values."""
+    h = np.ma.array([20.72, 234.85, 456.69, 683.21])
+    u = np.ma.array([-2.32, -3.23, 0.736, 9.07])
+    v = np.ma.array([8.31, 13.57, 25.56, 30.55])
+    u = np.ma.array(np.zeros((4,)))
+    v = np.zeros_like(u)
+    pos, neg, com = storm_relative_helicity(units.knot * u, units.knot * v, units.meter * h,
+                                            depth=500 * units.meter,
+                                            storm_u=15.77463015050421 * units('m/s'),
+                                            storm_v=21.179437759755647 * units('m/s'))
+
+    assert not np.ma.is_masked(pos)
+    assert not np.ma.is_masked(neg)
+    assert not np.ma.is_masked(com)
+
+
 def test_absolute_vorticity_asym():
     """Test absolute vorticity calculation with a complicated field."""
     u = np.array([[2, 4, 8], [0, 2, 2], [4, 6, 8]]) * units('m/s')

metpy/calc/tests/test_thermo.py

@@ -16,7 +16,7 @@
                         mixed_parcel, mixing_ratio, mixing_ratio_from_relative_humidity,
                         mixing_ratio_from_specific_humidity, moist_lapse,
                         moist_static_energy, most_unstable_cape_cin, most_unstable_parcel,
-                        parcel_profile, potential_temperature,
+                        parcel_profile, parcel_profile_with_lcl, potential_temperature,
                         psychrometric_vapor_pressure_wet,
                         relative_humidity_from_dewpoint,
                         relative_humidity_from_mixing_ratio,
@@ -140,6 +140,25 @@ def test_parcel_profile():
     assert_array_almost_equal(prof, true_prof, 2)
 
 
+def test_parcel_profile_lcl():
+    """Test parcel profile with lcl calculation."""
+    p = np.array([1004., 1000., 943., 928., 925., 850., 839., 749., 700., 699.]) * units.hPa
+    t = np.array([24.2, 24., 20.2, 21.6, 21.4, 20.4, 20.2, 14.4, 13.2, 13.]) * units.degC
+    td = np.array([21.9, 22.1, 19.2, 20.5, 20.4, 18.4, 17.4, 8.4, -2.8, -3.0]) * units.degC
+
+    true_prof = np.array([297.35, 297.01, 294.5, 293.48, 292.92, 292.81, 289.79, 289.32,
+                          285.15, 282.59, 282.53]) * units.kelvin
+    true_p = np.insert(p.m, 2, 970.699) * units.mbar
+    true_t = np.insert(t.m, 2, 22.047) * units.degC
+    true_td = np.insert(td.m, 2, 20.609) * units.degC
+
+    pressure, temp, dewp, prof = parcel_profile_with_lcl(p, t, td)
+    assert_almost_equal(pressure, true_p, 3)
+    assert_almost_equal(temp, true_t, 3)
+    assert_almost_equal(dewp, true_td, 3)
+    assert_array_almost_equal(prof, true_prof, 2)
+
+
 def test_parcel_profile_saturated():
     """Test parcel_profile works when LCL in levels (issue #232)."""
     levels = np.array([1000., 700., 500.]) * units.mbar
@@ -287,14 +306,36 @@ def test_lfc_equals_lcl():
     levels = np.array([912., 905.3, 874.4, 850., 815.1, 786.6, 759.1,
                        748., 732.2, 700., 654.8]) * units.mbar
     temperatures = np.array([29.4, 28.7, 25.2, 22.4, 19.4, 16.8,
-                             14.3, 13.2, 12.6, 11.4, 7.1]) * units.celsius
+                             14.0, 13.2, 12.6, 11.4, 7.1]) * units.celsius
     dewpoints = np.array([18.4, 18.1, 16.6, 15.4, 13.2, 11.4, 9.6,
                           8.8, 0., -18.6, -22.9]) * units.celsius
     lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints)
     assert_almost_equal(lfc_pressure, 777.0333 * units.mbar, 2)
     assert_almost_equal(lfc_temp, 15.8714 * units.celsius, 2)
 
 
+def test_sensitive_sounding():
+    """Test quantities for a sensitive sounding (#902)."""
+    # This sounding has a very small positive area in the low level. It's only captured
+    # properly if the parcel profile includes the LCL, otherwise it breaks LFC and CAPE
+    p = units.Quantity([1004., 1000., 943., 928., 925., 850., 839., 749., 700., 699.,
+                        603., 500., 404., 400., 363., 306., 300., 250., 213., 200.,
+                        176., 150.], 'hectopascal')
+    t = units.Quantity([24.2, 24., 20.2, 21.6, 21.4, 20.4, 20.2, 14.4, 13.2, 13., 6.8, -3.3,
+                        -13.1, -13.7, -17.9, -25.5, -26.9, -37.9, -46.7, -48.7, -52.1, -58.9],
+                       'degC')
+    td = units.Quantity([21.9, 22.1, 19.2, 20.5, 20.4, 18.4, 17.4, 8.4, -2.8, -3.0, -15.2,
+                         -20.3, -29.1, -27.7, -24.9, -39.5, -41.9, -51.9, -60.7, -62.7, -65.1,
+                         -71.9], 'degC')
+    lfc_pressure, lfc_temp = lfc(p, t, td)
+    assert_almost_equal(lfc_pressure, 947.476 * units.mbar, 2)
+    assert_almost_equal(lfc_temp, 20.498 * units.degC, 2)
+
+    pos, neg = surface_based_cape_cin(p, t, td)
+    assert_almost_equal(pos, 0.112 * units('J/kg'), 3)
+    assert_almost_equal(neg, -6.075 * units('J/kg'), 3)
+
+
 def test_lfc_sfc_precision():
     """Test LFC when there are precision issues with the parcel path."""
     levels = np.array([839., 819.4, 816., 807., 790.7, 763., 736.2,
@@ -796,7 +837,7 @@ def test_surface_based_cape_cin():
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     cape, cin = surface_based_cape_cin(p, temperature, dewpoint)
     assert_almost_equal(cape, 58.0368212 * units('joule / kilogram'), 6)
-    assert_almost_equal(cin, -89.8073512 * units('joule / kilogram'), 6)
+    assert_almost_equal(cin, -136.597240 * units('joule / kilogram'), 6)
 
 
 def test_most_unstable_cape_cin_surface():

metpy/calc/thermo.py

@@ -350,8 +350,11 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None):
     """
     # Default to surface parcel if no profile or starting pressure level is given
     if parcel_temperature_profile is None:
-        parcel_temperature_profile = parcel_profile(pressure, temperature[0],
-                                                    dewpt[0]).to('degC')
+        new_stuff = parcel_profile_with_lcl(pressure, temperature, dewpt)
+        pressure, temperature, _, parcel_temperature_profile = new_stuff
+        temperature = temperature.to('degC')
+        parcel_temperature_profile = parcel_temperature_profile.to('degC')
+
     # The parcel profile and data have the same first data point, so we ignore
     # that point to get the real first intersection for the LFC calculation.
     x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:],
@@ -413,8 +416,11 @@ def el(pressure, temperature, dewpt, parcel_temperature_profile=None):
     """
     # Default to surface parcel if no profile or starting pressure level is given
     if parcel_temperature_profile is None:
-        parcel_temperature_profile = parcel_profile(pressure, temperature[0],
-                                                    dewpt[0]).to('degC')
+        new_stuff = parcel_profile_with_lcl(pressure, temperature, dewpt)
+        pressure, temperature, _, parcel_temperature_profile = new_stuff
+        temperature = temperature.to('degC')
+        parcel_temperature_profile = parcel_temperature_profile.to('degC')
+
     # If the top of the sounding parcel is warmer than the environment, there is no EL
     if parcel_temperature_profile[-1] > temperature[-1]:
         return np.nan * pressure.units, np.nan * temperature.units
@@ -456,27 +462,99 @@ def parcel_profile(pressure, temperature, dewpt):
     --------
     lcl, moist_lapse, dry_lapse
 
+    """
+    _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpt)
+    return concatenate((t_l, t_u))
+
+
+@exporter.export
+@preprocess_xarray
+@check_units('[pressure]', '[temperature]', '[temperature]')
+def parcel_profile_with_lcl(pressure, temperature, dewpt):
+    r"""Calculate the profile a parcel takes through the atmosphere.
+
+    The parcel starts at `temperature`, and `dewpt`, lifted up
+    dry adiabatically to the LCL, and then moist adiabatically from there.
+    `pressure` specifies the pressure levels for the profile. This function returns
+    a profile that includes the LCL.
+
+    Parameters
+    ----------
+    pressure : `pint.Quantity`
+        The atmospheric pressure level(s) of interest. The first entry should be the starting
+        point pressure.
+    temperature : `pint.Quantity`
+        The atmospheric temperature at the levels in `pressure`. The first entry should be the
+        starting point temperature.
+    dewpt : `pint.Quantity`
+        The atmospheric dew point at the levels in `pressure`. The first entry should be the
+        starting dew point.
+
+    Returns
+    -------
+    pressure : `pint.Quantity`
+        The parcel profile pressures, which includes the specified levels and the LCL
+    ambient_temperature : `pint.Quantity`
+        The atmospheric temperature values, including the value interpolated to the LCL level
+    ambient_dew_point : `pint.Quantity`
+        The atmospheric dew point values, including the value interpolated to the LCL level
+    profile_temperature : `pint.Quantity`
+        The parcel profile temperatures at all of the levels in the returned pressures array,
+        including the LCL.
+
+    See Also
+    --------
+    lcl, moist_lapse, dry_lapse, parcel_profile
+
+    """
+    p_l, p_lcl, p_u, t_l, t_lcl, t_u = _parcel_profile_helper(pressure, temperature[0],
+                                                              dewpt[0])
+    new_press = concatenate((p_l, p_lcl, p_u))
+    prof_temp = concatenate((t_l, t_lcl, t_u))
+    new_temp = _insert_lcl_level(pressure, temperature, p_lcl)
+    new_dewp = _insert_lcl_level(pressure, dewpt, p_lcl)
+    return new_press, new_temp, new_dewp, prof_temp
+
+
+def _parcel_profile_helper(pressure, temperature, dewpt):
+    """Help calculate parcel profiles.
+
+    Returns the temperature and pressure, above, below, and including the LCL. The
+    other calculation functions decide what to do with the pieces.
+
     """
     # Find the LCL
-    lcl_pressure, _ = lcl(pressure[0], temperature, dewpt)
-    lcl_pressure = lcl_pressure.to(pressure.units)
+    press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpt)
+    press_lcl = press_lcl.to(pressure.units)
 
     # Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the
     # LCL is included in the levels. It's slightly redundant in that case, but simplifies
     # the logic for removing it later.
-    press_lower = concatenate((pressure[pressure >= lcl_pressure], lcl_pressure))
-    t1 = dry_lapse(press_lower, temperature)
+    press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl))
+    temp_lower = dry_lapse(press_lower, temperature)
 
     # If the pressure profile doesn't make it to the lcl, we can stop here
-    if _greater_or_close(np.nanmin(pressure), lcl_pressure.m):
-        return t1[:-1]
+    if _greater_or_close(np.nanmin(pressure), press_lcl.m):
+        return (press_lower[:-1], press_lcl, np.array([]) * press_lower.units,
+                temp_lower[:-1], temp_lcl, np.array([]) * temp_lower.units)
 
     # Find moist pseudo-adiabatic profile starting at the LCL
-    press_upper = concatenate((lcl_pressure, pressure[pressure < lcl_pressure]))
-    t2 = moist_lapse(press_upper, t1[-1]).to(t1.units)
+    press_upper = concatenate((press_lcl, pressure[pressure < press_lcl]))
+    temp_upper = moist_lapse(press_upper, temp_lower[-1]).to(temp_lower.units)
+
+    # Return profile pieces
+    return (press_lower[:-1], press_lcl, press_upper[1:],
+            temp_lower[:-1], temp_lcl, temp_upper[1:])
+
+
+def _insert_lcl_level(pressure, temperature, lcl_pressure):
+    """Insert the LCL pressure into the profile."""
+    interp_temp = interpolate_1d(lcl_pressure, pressure, temperature)
 
-    # Return LCL *without* the LCL point
-    return concatenate((t1[:-1], t2[1:]))
+    # Pressure needs to be increasing for searchsorted, so flip it and then convert
+    # the index back to the original array
+    loc = pressure.size - pressure[::-1].searchsorted(lcl_pressure)
+    return np.insert(temperature.m, loc, interp_temp.m) * temperature.units
 
 
 @exporter.export
@@ -1571,8 +1649,8 @@ def surface_based_cape_cin(pressure, temperature, dewpoint):
     cape_cin, parcel_profile
 
     """
-    profile = parcel_profile(pressure, temperature[0], dewpoint[0])
-    return cape_cin(pressure, temperature, dewpoint, profile)
+    p, t, td, profile = parcel_profile_with_lcl(pressure, temperature, dewpoint)
+    return cape_cin(p, t, td, profile)
 
 
 @exporter.export