From e79c0d976d43eefa1c96d914e4bcf267b46130b5 Mon Sep 17 00:00:00 2001 From: Drew Camron Date: Tue, 5 Apr 2022 01:25:59 -0500 Subject: [PATCH 1/3] Handle duplicates in moist profile (fixes #2309) Repeat values breaks solve_ivp in moist_lapse. Without breaking existing handling of values close to reference_pressure, remove duplicate pressures for moist_lapse input, then re-add results at duplicate levels profile to match pre-1.1 behavior. --- src/metpy/calc/thermo.py | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index e9e7dae1cb8..79f635c304d 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -748,11 +748,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`` @@ -809,7 +815,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 @@ -934,9 +944,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:], From 57e6c3c1821e7debfad0089a8879f8b103ccf9bc Mon Sep 17 00:00:00 2001 From: Drew Camron Date: Tue, 5 Apr 2022 01:33:01 -0500 Subject: [PATCH 2/3] Test repeat pressure handling in moist profile Externally tested with MetPy 1.0.1. --- tests/calc/test_thermo.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index df53327214a..a0deabedf64 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -2092,3 +2092,23 @@ 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) From 83e541f1c5d631fa661ac411ee91013fcb777d9a Mon Sep 17 00:00:00 2001 From: Drew Camron Date: Wed, 6 Apr 2022 13:59:58 -0500 Subject: [PATCH 3/3] Add test for duplicates in profile `_as_dataset` Make sure duplicate handling doesn't interrupt dataset construction and proper alignment of results. --- tests/calc/test_thermo.py | 45 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index a0deabedf64..cc8e5827bfc 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -2112,3 +2112,48 @@ def test_parcel_profile_drop_duplicates(): 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)