Merge pull request #2412 from dcamron/parcel-profile-duplicate

Drop repeat pressures in moist region of `parcel_profile`

src/metpy/calc/thermo.py

@@ -744,11 +744,17 @@ def parcel_profile(pressure, temperature, dewpoint):
 
     See Also
     --------
-    lcl, moist_lapse, dry_lapse
+    lcl, moist_lapse, dry_lapse, parcel_profile_with_lcl, parcel_profile_with_lcl_as_dataset
 
     Notes
     -----
     Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
+    Duplicate pressure levels return duplicate parcel temperatures. Consider preprocessing
+    low-precision, high frequency profiles with tools like `scipy.medfilt`,
+    `pandas.drop_duplicates`, or `numpy.unique`.
+
+    Will only return Pint Quantities, even when given xarray DataArray profiles. To
+    obtain a xarray Dataset instead, use `parcel_profile_with_lcl_as_dataset` instead.
 
     .. versionchanged:: 1.0
        Renamed ``dewpt`` parameter to ``dewpoint``
@@ -805,7 +811,11 @@ def parcel_profile_with_lcl(pressure, temperature, dewpoint):
     Notes
     -----
     Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
-    Also, will only return Pint Quantities, even when given xarray DataArray profiles. To
+    Duplicate pressure levels return duplicate parcel temperatures. Consider preprocessing
+    low-precision, high frequency profiles with tools like `scipy.medfilt`,
+    `pandas.drop_duplicates`, or `numpy.unique`.
+
+    Will only return Pint Quantities, even when given xarray DataArray profiles. To
     obtain a xarray Dataset instead, use `parcel_profile_with_lcl_as_dataset` instead.
 
     .. versionchanged:: 1.0
@@ -930,9 +940,20 @@ def _parcel_profile_helper(pressure, temperature, dewpoint):
         return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units),
                 temp_lower[:-1], temp_lcl, units.Quantity(np.array([]), temp_lower.units))
 
-    # Find moist pseudo-adiabatic profile starting at the LCL
+    # Establish profile above LCL
     press_upper = concatenate((press_lcl, pressure[pressure < press_lcl]))
-    temp_upper = moist_lapse(press_upper, temp_lower[-1]).to(temp_lower.units)
+
+    # Remove duplicate pressure values from remaining profile. Needed for solve_ivp in
+    # moist_lapse. unique will return remaining values sorted ascending.
+    unique, indices, counts = np.unique(press_upper.m, return_inverse=True, return_counts=True)
+    unique = units.Quantity(unique, press_upper.units)
+    if np.any(counts > 1):
+        warnings.warn(f'Duplicate pressure(s) {unique[counts > 1]:~P} provided. '
+                      'Output profile includes duplicate temperatures as a result.')
+
+    # Find moist pseudo-adiabatic profile starting at the LCL, reversing above sorting
+    temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units)
+    temp_upper = temp_upper[::-1][indices]
 
     # Return profile pieces
     return (press_lower[:-1], press_lcl, press_upper[1:],

tests/calc/test_thermo.py

@@ -2092,3 +2092,68 @@ def test_cross_totals():
 
     ct = cross_totals(pressure, temperature, dewpoint)
     assert_almost_equal(ct, 21.40 * units.delta_degC, 2)
+
+
+def test_parcel_profile_drop_duplicates():
+    """Test handling repeat pressures in moist region of profile."""
+    pressure = np.array([962., 951., 937.9, 925., 908., 905.7, 894., 875.,
+                         41.3, 40.8, 37., 36.8, 32., 30., 27.7, 27.7, 26.4]) * units.hPa
+
+    temperature = units.Quantity(19.6, 'degC')
+
+    dewpoint = units.Quantity(18.6, 'degC')
+
+    truth = np.array([292.75, 291.78965331, 291.12778784, 290.61996294,
+                      289.93681828, 289.84313902, 289.36183185, 288.5626898,
+                      135.46280886, 134.99220142, 131.27369084, 131.07055878,
+                      125.93977169, 123.63877507, 120.85291224, 120.85291224,
+                      119.20448296]) * units.kelvin
+
+    with pytest.warns(UserWarning, match='Duplicate pressure'):
+        profile = parcel_profile(pressure, temperature, dewpoint)
+        assert_almost_equal(profile, truth, 5)
+
+
+def test_parcel_profile_with_lcl_as_dataset_duplicates():
+    """Test that parcel profile dataset creation works with duplicate pressures in profile."""
+    pressure = np.array(
+        [951., 951., 937.9, 925., 908., 30., 27.7, 27.7, 26.4, 25.1]
+    ) * units.hPa
+
+    temperature = np.array(
+        [20., 20., 19.5, 19., 18.6, -58.5, -58.1, -58.1, -57.2, -56.2]
+    ) * units.degC
+
+    dewpoint = np.array(
+        [19.4, 19.4, 19., 18.6, 18.3, -73.5, -75.1, -75.1, -77., -78.8]
+    ) * units.degC
+
+    truth = xr.Dataset(
+        {
+            'ambient_temperature': (
+                ('isobaric',),
+                np.insert(temperature.m, 2, 19.679237747615478) * units.degC
+            ),
+            'ambient_dew_point': (
+                ('isobaric',),
+                np.insert(dewpoint.m, 2, 19.143390198092384) * units.degC
+            ),
+            'parcel_temperature': (
+                ('isobaric',),
+                [
+                    293.15, 293.15, 292.40749167, 292.22841462, 291.73069653, 291.06139433,
+                    125.22698955, 122.40534065, 122.40534065, 120.73573642, 119.0063293
+                ] * units.kelvin
+            )
+        },
+        coords={
+            'isobaric': (
+                'isobaric',
+                np.insert(pressure.m, 2, 942.6)
+            )
+        }
+    )
+
+    profile = parcel_profile_with_lcl_as_dataset(pressure, temperature, dewpoint)
+
+    xr.testing.assert_allclose(profile, truth, atol=1e-5)