BUG: Fix edge cases in profiles calculations (Fixes #902)

By not including the LCL in the profiles used in calculations, some odd
corner cases result where, e.g., the intersection points with the
ambient profile change. It also has an impact in cape/cin calculations.
Refactor profile handling to explicitly work with a profile that
includes the LCL and expose this as an API function.

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/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
@@ -295,6 +314,28 @@ def test_lfc_equals_lcl():
     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,

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