diff --git a/metpy/calc/tests/test_thermo.py b/metpy/calc/tests/test_thermo.py index 2bb67054b6d..b50c5f163c1 100644 --- a/metpy/calc/tests/test_thermo.py +++ b/metpy/calc/tests/test_thermo.py @@ -1201,6 +1201,69 @@ def test_lfc_not_below_lcl(): assert_almost_equal(lfc_temp, 6.4992871 * units.celsius, 6) +def test_multiple_lfcs(): + """Test sounding with multiple LFCs. + + If which='top', return lowest-pressure LFC. + If which='all', return all LFCs + + """ + levels = np.array([966., 937.2, 925., 904.6, 872.6, 853., 850., 836., 821., 811.6, 782.3, + 754.2, 726.9, 700., 648.9, 624.6, 601.1, 595., 587., 576., 555.7, + 534.2, 524., 500., 473.3, 400., 384.5, 358., 343., 308.3, 300., 276., + 273., 268.5, 250., 244.2, 233., 200.]) * units.mbar + temperatures = np.array([18.2, 16.8, 16.2, 15.1, 13.3, 12.2, 12.4, 14., 14.4, + 13.7, 11.4, 9.1, 6.8, 4.4, -1.4, -4.4, -7.3, -8.1, + -7.9, -7.7, -8.7, -9.8, -10.3, -13.5, -17.1, -28.1, -30.7, + -35.3, -37.1, -43.5, -45.1, -49.9, -50.4, -51.1, -54.1, -55., + -56.7, -57.5]) * units.degC + dewpoints = np.array([16.9, 15.9, 15.5, 14.2, 12.1, 10.8, 8.6, 0., -3.6, -4.4, + -6.9, -9.5, -12., -14.6, -15.8, -16.4, -16.9, -17.1, -27.9, -42.7, + -44.1, -45.6, -46.3, -45.5, -47.1, -52.1, -50.4, -47.3, -57.1, + -57.9, -58.1, -60.9, -61.4, -62.1, -65.1, -65.6, + -66.7, -70.5]) * units.degC + lfc_pressure_top, lfc_temp_top = lfc(levels, temperatures, dewpoints) + lfc_pressure_bottom, lfc_temp_bottom = lfc(levels, temperatures, dewpoints, + which='bottom') + lfc_pressure_all, _ = lfc(levels, temperatures, dewpoints, which='all') + assert_almost_equal(lfc_pressure_top, 705.5170051 * units.mbar, 6) + assert_almost_equal(lfc_temp_top, 4.8922235 * units.degC, 6) + assert_almost_equal(lfc_pressure_bottom, 884.3290151 * units.mbar, 6) + assert_almost_equal(lfc_temp_bottom, 13.9597571 * units.degC, 6) + assert_almost_equal(len(lfc_pressure_all), 2, 0) + + +def test_multiple_els(): + """Test sounding with multiple ELs. + + If which='top', return lowest-pressure EL. + If which='all', return all ELs + + """ + levels = np.array([966., 937.2, 925., 904.6, 872.6, 853., 850., 836., 821., 811.6, 782.3, + 754.2, 726.9, 700., 648.9, 624.6, 601.1, 595., 587., 576., 555.7, + 534.2, 524., 500., 473.3, 400., 384.5, 358., 343., 308.3, 300., 276., + 273., 268.5, 250., 244.2, 233., 200.]) * units.mbar + temperatures = np.array([18.2, 16.8, 16.2, 15.1, 13.3, 12.2, 12.4, 14., 14.4, + 13.7, 11.4, 9.1, 6.8, 4.4, -1.4, -4.4, -7.3, -8.1, + -7.9, -7.7, -8.7, -9.8, -10.3, -13.5, -17.1, -28.1, -30.7, + -35.3, -37.1, -43.5, -45.1, -49.9, -50.4, -51.1, -54.1, -55., + -56.7, -57.5]) * units.degC + dewpoints = np.array([16.9, 15.9, 15.5, 14.2, 12.1, 10.8, 8.6, 0., -3.6, -4.4, + -6.9, -9.5, -12., -14.6, -15.8, -16.4, -16.9, -17.1, -27.9, -42.7, + -44.1, -45.6, -46.3, -45.5, -47.1, -52.1, -50.4, -47.3, -57.1, + -57.9, -58.1, -60.9, -61.4, -62.1, -65.1, -65.6, + -66.7, -70.5]) * units.degC + el_pressure_top, el_temp_top = el(levels, temperatures, dewpoints) + el_pressure_bottom, el_temp_bottom = el(levels, temperatures, dewpoints, which='bottom') + el_pressure_all, _ = el(levels, temperatures, dewpoints, which='all') + assert_almost_equal(el_pressure_top, 228.3671032 * units.mbar, 6) + assert_almost_equal(el_temp_top, -56.8123126 * units.degC, 6) + assert_almost_equal(el_pressure_bottom, 849.7958931 * units.mbar, 6) + assert_almost_equal(el_temp_bottom, 12.4233265 * units.degC, 6) + assert_almost_equal(len(el_pressure_all), 2, 0) + + def test_cape_cin_custom_profile(): """Test the CAPE and CIN calculation with a custom profile passed to LFC and EL.""" p = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar diff --git a/metpy/calc/thermo.py b/metpy/calc/thermo.py index fae1a08d1ce..a36dfe9746e 100644 --- a/metpy/calc/thermo.py +++ b/metpy/calc/thermo.py @@ -362,7 +362,8 @@ def _lcl_iter(p, p0, w, t): @exporter.export @preprocess_xarray @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') -def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_start=None): +def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_start=None, + which='top'): r"""Calculate the level of free convection (LFC). This works by finding the first intersection of the ideal parcel path and @@ -382,11 +383,14 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta dewpt_start: `pint.Quantity`, optional The dewpoint of the parcel for which to calculate the LFC. Defaults to the surface dewpoint. + which: str, optional + Pick which LFC to return. Options are 'top', 'bottom', and 'all'. + Default is the 'top' (lowest pressure) LFC. Returns ------- `pint.Quantity` - The LFC pressure and temperature + The LFC pressure and temperature, or arrays of same if which='all' See Also -------- @@ -439,17 +443,30 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta if not any(idx): x, y = this_lcl return x, y - # Otherwise, make select first candidate LFC above the LCL + # Otherwise, find all LFCs that exist above the LCL + # What is returned depends on which flag as described in the docstring else: - x = x[idx] - y = y[idx] - return x[0], y[0] + x, y = _multiple_el_lfc_options(x, y, idx, which) + return x, y + + +def _multiple_el_lfc_options(intersect_pressures, intersect_temperatures, valid_x, + which): + """Choose which ELs and LFCs to return from a sounding.""" + p_list, t_list = intersect_pressures[valid_x], intersect_temperatures[valid_x] + if which == 'all': + x, y = p_list, t_list + elif which == 'bottom': + x, y = p_list[0], t_list[0] + elif which == 'top': + x, y = p_list[-1], t_list[-1] + return x, y @exporter.export @preprocess_xarray @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') -def el(pressure, temperature, dewpt, parcel_temperature_profile=None): +def el(pressure, temperature, dewpt, parcel_temperature_profile=None, which='top'): r"""Calculate the equilibrium level. This works by finding the last intersection of the ideal parcel path and @@ -467,11 +484,14 @@ def el(pressure, temperature, dewpt, parcel_temperature_profile=None): parcel_temperature_profile: `pint.Quantity`, optional The parcel temperature profile from which to calculate the EL. Defaults to the surface parcel profile. + which: str, optional + Pick which EL to return. Options are 'top', 'bottom', and 'all'. + Default is the 'top' (lowest pressure) EL. Returns ------- `pint.Quantity, pint.Quantity` - The EL pressure and temperature + The EL pressure and temperature, or lists of same if which='all' See Also -------- @@ -491,10 +511,12 @@ def el(pressure, temperature, dewpt, parcel_temperature_profile=None): # Otherwise the last intersection (as long as there is one, and it's not # below the LCL) is the EL - x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:]) + x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:], + direction='decreasing') lcl_p, _ = lcl(pressure[0], temperature[0], dewpt[0]) + idx = x < lcl_p if len(x) > 0 and x[-1] < lcl_p: - return x[-1], y[-1] + return _multiple_el_lfc_options(x, y, idx, which) else: return np.nan * pressure.units, np.nan * temperature.units