From 10bbda76e491f5d31973b091369cb882173f74d7 Mon Sep 17 00:00:00 2001 From: zbruick Date: Wed, 18 Sep 2019 16:26:24 -0600 Subject: [PATCH 1/2] Return input temperature units for LCL/LFC/EL. Disallow LFC when it and EL both occur below LCL. --- metpy/calc/tests/test_thermo.py | 50 ++++++++++++++++++++++++++++++++- metpy/calc/thermo.py | 20 +++++++------ 2 files changed, 61 insertions(+), 9 deletions(-) diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py index a16c4121b02..59cbb4e9bfc 100644 --- a/metpy/calc/tests/test_thermo.py +++ b/metpy/calc/tests/test_thermo.py @@ -1,4 +1,4 @@ -# Copyright (c) 2008,2015,2016,2017,2018 MetPy Developers. +# Copyright (c) 2008,2015,2016,2017,2018,2019 MetPy Developers. # Distributed under the terms of the BSD 3-Clause License. # SPDX-License-Identifier: BSD-3-Clause """Test the `thermo` module.""" @@ -260,6 +260,16 @@ def test_lcl(): assert_almost_equal(lcl_temperature, 17.676 * units.degC, 2) +def test_lcl_kelvin(): + """Test LCL temperature is returned as Kelvin, if temperature is Kelvin.""" + temperature = 273.09723 * units.kelvin + lcl_pressure, lcl_temperature = lcl(1017.16 * units.mbar, temperature, + 264.5351 * units.kelvin) + assert_almost_equal(lcl_pressure, 889.416 * units.mbar, 2) + assert_almost_equal(lcl_temperature, 262.827 * units.kelvin, 2) + assert lcl_temperature.units == temperature.units + + def test_lcl_convergence(): """Test LCL calculation convergence failure.""" with pytest.raises(RuntimeError): @@ -276,6 +286,19 @@ def test_lfc_basic(): assert_almost_equal(lfc_temp, 9.705 * units.celsius, 2) +def test_lfc_kelvin(): + """Test that LFC temperature returns Kelvin if Kelvin is provided.""" + pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar + temperature = (np.array([22.2, 14.6, 12., 9.4, 7., -49.] + ) + 273.15) * units.kelvin + dewpoint = (np.array([19., -11.2, -10.8, -10.4, -10., -53.2] + ) + 273.15) * units.kelvin + lfc_pressure, lfc_temp = lfc(pressure, temperature, dewpoint) + assert_almost_equal(lfc_pressure, 727.415 * units.mbar, 2) + assert_almost_equal(lfc_temp, 9.705 * units.degC, 2) + assert lfc_temp.units == temperature.units + + def test_lfc_ml(): """Test Mixed-Layer LFC calculation.""" levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar @@ -499,6 +522,17 @@ def test_el(): assert_almost_equal(el_temperature, -11.7027 * units.degC, 3) +def test_el_kelvin(): + """Test that EL temperature returns Kelvin if Kelvin is provided.""" + levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar + temperatures = (np.array([22.2, 14.6, 12., 9.4, 7., -38.]) + 273.15) * units.kelvin + dewpoints = (np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) + 273.15) * units.kelvin + el_pressure, el_temp = el(levels, temperatures, dewpoints) + assert_almost_equal(el_pressure, 470.4075 * units.mbar, 3) + assert_almost_equal(el_temp, -11.7027 * units.degC, 3) + assert el_temp.units == temperatures.units + + def test_el_ml(): """Test equilibrium layer calculation for a mixed parcel.""" levels = np.array([959., 779.2, 751.3, 724.3, 700., 400., 269.]) * units.mbar @@ -535,6 +569,20 @@ def test_no_el_multi_crossing(): assert_nan(el_temperature, temperatures.units) +def test_lfc_and_el_below_lcl(): + """Test that LFC and EL are returned as NaN if both are below LCL.""" + dewpoint = [264.5351, 261.13443, 259.0122, 252.30063, 248.58017, 242.66582] * units.kelvin + temperature = [273.09723, 268.40173, 263.56207, 260.257, 256.63538, + 252.91345] * units.kelvin + pressure = [1017.16, 950, 900, 850, 800, 750] * units.hPa + el_pressure, el_temperature = el(pressure, temperature, dewpoint) + lfc_pressure, lfc_temperature = lfc(pressure, temperature, dewpoint) + assert_nan(lfc_pressure, pressure.units) + assert_nan(lfc_temperature, temperature.units) + assert_nan(el_pressure, pressure.units) + assert_nan(el_temperature, temperature.units) + + def test_el_lfc_equals_lcl(): """Test equilibrium layer calculation when the lfc equals the lcl.""" levels = np.array([912., 905.3, 874.4, 850., 815.1, 786.6, 759.1, 748., diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py index 00628858faa..4c110add9ba 100644 --- a/metpy/calc/thermo.py +++ b/metpy/calc/thermo.py @@ -1,4 +1,4 @@ -# Copyright (c) 2008,2015,2016,2017,2018 MetPy Developers. +# Copyright (c) 2008,2015,2016,2017,2018,2019 MetPy Developers. # Distributed under the terms of the BSD 3-Clause License. # SPDX-License-Identifier: BSD-3-Clause """Contains a collection of thermodynamic calculations.""" @@ -359,7 +359,7 @@ def _lcl_iter(p, p0, w, t): fp = so.fixed_point(_lcl_iter, pressure.m, args=(pressure.m, w, temperature), xtol=eps, maxiter=max_iters) lcl_p = fp * pressure.units - return lcl_p, dewpoint(vapor_pressure(lcl_p, w)) + return lcl_p, dewpoint(vapor_pressure(lcl_p, w)).to(temperature.units) @exporter.export @@ -406,8 +406,7 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta if parcel_temperature_profile is None: 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') + parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) if dewpt_start is None: dewpt_start = dewpt[0] @@ -446,8 +445,14 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta idx = x < this_lcl[0] # LFC height < LCL height, so set LFC = LCL if not any(idx): - x, y = this_lcl - return x, y + el_pres, _ = find_intersections(pressure[1:], parcel_temperature_profile[1:], + temperature[1:], direction='decreasing', + log_x=True) + if el_pres > this_lcl[0]: + return np.nan * pressure.units, np.nan * temperature.units + else: + x, y = this_lcl + return x, y # Otherwise, find all LFCs that exist above the LCL # What is returned depends on which flag as described in the docstring else: @@ -508,8 +513,7 @@ def el(pressure, temperature, dewpt, parcel_temperature_profile=None, which='top if parcel_temperature_profile is None: 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') + parcel_temperature_profile = parcel_temperature_profile.to(temperature.units) # If the top of the sounding parcel is warmer than the environment, there is no EL if parcel_temperature_profile[-1] > temperature[-1]: From 60ca79087c85c2e3380ab9fed660c346d5c418fc Mon Sep 17 00:00:00 2001 From: zbruick Date: Thu, 19 Sep 2019 09:29:06 -0600 Subject: [PATCH 2/2] Reduce return statements in `lfc` --- metpy/calc/thermo.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py index 4c110add9ba..feb6666265c 100644 --- a/metpy/calc/thermo.py +++ b/metpy/calc/thermo.py @@ -435,10 +435,10 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta mask = pressure < this_lcl[0] if np.all(_less_or_close(parcel_temperature_profile[mask], temperature[mask])): # LFC doesn't exist - return np.nan * pressure.units, np.nan * temperature.units + x, y = np.nan * pressure.units, np.nan * temperature.units else: # LFC = LCL x, y = this_lcl - return x, y + return x, y # LFC exists. Make sure it is no lower than the LCL else: @@ -448,11 +448,11 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta el_pres, _ = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:], direction='decreasing', log_x=True) - if el_pres > this_lcl[0]: - return np.nan * pressure.units, np.nan * temperature.units + if np.min(el_pres) > this_lcl[0]: + x, y = np.nan * pressure.units, np.nan * temperature.units else: x, y = this_lcl - return x, y + return x, y # Otherwise, find all LFCs that exist above the LCL # What is returned depends on which flag as described in the docstring else: