From 8a83ef076698361312883437912bc4a3ae6bddf6 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 15 Apr 2022 12:57:45 -0500 Subject: [PATCH 01/30] first cut at using virtual temperature in parcel profile --- src/metpy/calc/thermo.py | 32 +++++++++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 3 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a85c4d6df4d..a49e16c5de5 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -715,10 +715,20 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which=' units.Quantity(np.nan, temperature.units)) +class parcelPathAssumptions(object): + """ + Holds assumptions made about the parcel path during calculations. + """ + def __init__(self): + self.use_virtual_temperature = True + self.moist_adiabat = 'pseudoadiabatic' + + @exporter.export @preprocess_and_wrap(wrap_like='pressure') @check_units('[pressure]', '[temperature]', '[temperature]') -def parcel_profile(pressure, temperature, dewpoint): +def parcel_profile(pressure, temperature, dewpoint, + assumptions=parcelPathAssumptions()): r"""Calculate the profile a parcel takes through the atmosphere. The parcel starts at `temperature`, and `dewpoint`, lifted up @@ -760,7 +770,8 @@ def parcel_profile(pressure, temperature, dewpoint): Renamed ``dewpt`` parameter to ``dewpoint`` """ - _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint) + _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint, + assumptions=assumptions) return concatenate((t_l, t_u)) @@ -911,7 +922,8 @@ def _check_pressure(pressure): return np.all(pressure[:-1] >= pressure[1:]) -def _parcel_profile_helper(pressure, temperature, dewpoint): +def _parcel_profile_helper(pressure, temperature, dewpoint, + assumptions=parcelPathAssumptions()): """Help calculate parcel profiles. Returns the temperature and pressure, above, below, and including the LCL. The @@ -935,6 +947,13 @@ def _parcel_profile_helper(pressure, temperature, dewpoint): press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl)) temp_lower = dry_lapse(press_lower, temperature) + # Do the virtual temperature correction for the parcel below the LCL + if assumptions.use_virtual_temperature: + # Calculate the relative humidity, mixing ratio, and virtual temperature + rh_lowest = relative_humidity_from_dewpoint(temperature, dewpoint) + rv_lower = mixing_ratio_from_relative_humidity(press_lower[0], temperature, rh_lowest) + temp_lower = virtual_temperature(temp_lower, rv_lower).to(temp_lower.units) + # If the pressure profile doesn't make it to the lcl, we can stop here if _greater_or_close(np.nanmin(pressure), press_lcl): return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units), @@ -955,6 +974,13 @@ def _parcel_profile_helper(pressure, temperature, dewpoint): temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units) temp_upper = temp_upper[::-1][indices] + # Do the virtual temperature correction for the parcel above the LCL + if assumptions.use_virtual_temperature: + + # Calculate the mixing ratio and virtual temperature + rv_upper = mixing_ratio(saturation_vapor_pressure(temp_lcl), press_upper) + temp_upper = virtual_temperature(temp_upper, rv_upper).to(temp_lower.units) + # Return profile pieces return (press_lower[:-1], press_lcl, press_upper[1:], temp_lower[:-1], temp_lcl, temp_upper[1:]) From 7ab2e47b6687430fa4083b0518879c3203eb6253 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Fri, 15 Apr 2022 13:18:20 -0500 Subject: [PATCH 02/30] fix linting with new parcel path class --- src/metpy/calc/thermo.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index a49e16c5de5..e4cef18ff28 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -715,10 +715,9 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which=' units.Quantity(np.nan, temperature.units)) -class parcelPathAssumptions(object): - """ - Holds assumptions made about the parcel path during calculations. - """ +class ParcelPathAssumptions(): + """Holds assumptions made about the parcel path during calculations.""" + def __init__(self): self.use_virtual_temperature = True self.moist_adiabat = 'pseudoadiabatic' @@ -728,7 +727,7 @@ def __init__(self): @preprocess_and_wrap(wrap_like='pressure') @check_units('[pressure]', '[temperature]', '[temperature]') def parcel_profile(pressure, temperature, dewpoint, - assumptions=parcelPathAssumptions()): + assumptions=ParcelPathAssumptions()): r"""Calculate the profile a parcel takes through the atmosphere. The parcel starts at `temperature`, and `dewpoint`, lifted up @@ -923,7 +922,7 @@ def _check_pressure(pressure): def _parcel_profile_helper(pressure, temperature, dewpoint, - assumptions=parcelPathAssumptions()): + assumptions=ParcelPathAssumptions()): """Help calculate parcel profiles. Returns the temperature and pressure, above, below, and including the LCL. The From 1536c9b7fb885458341a539a8adf69fc0b3f1826 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 28 Mar 2023 14:52:01 -0500 Subject: [PATCH 03/30] fix bugs in initial implementation --- src/metpy/calc/thermo.py | 18 +++++------------- tests/calc/test_thermo.py | 16 ++++++++-------- 2 files changed, 13 insertions(+), 21 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index e4cef18ff28..0e4da11ecd6 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -935,6 +935,11 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, Pressure increases between at least two points in your sounding. Using scipy.signal.medfilt may fix this.""" raise InvalidSoundingError(msg) + + if assumptions.use_virtual_temperature: + specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint) + mixing_ratio = mixing_ratio_from_specific_humidity(specific_humidity) + temperature = virtual_temperature(temperature, mixing_ratio) # Find the LCL press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpoint) @@ -946,13 +951,6 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl)) temp_lower = dry_lapse(press_lower, temperature) - # Do the virtual temperature correction for the parcel below the LCL - if assumptions.use_virtual_temperature: - # Calculate the relative humidity, mixing ratio, and virtual temperature - rh_lowest = relative_humidity_from_dewpoint(temperature, dewpoint) - rv_lower = mixing_ratio_from_relative_humidity(press_lower[0], temperature, rh_lowest) - temp_lower = virtual_temperature(temp_lower, rv_lower).to(temp_lower.units) - # If the pressure profile doesn't make it to the lcl, we can stop here if _greater_or_close(np.nanmin(pressure), press_lcl): return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units), @@ -973,12 +971,6 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units) temp_upper = temp_upper[::-1][indices] - # Do the virtual temperature correction for the parcel above the LCL - if assumptions.use_virtual_temperature: - - # Calculate the mixing ratio and virtual temperature - rv_upper = mixing_ratio(saturation_vapor_pressure(temp_lcl), press_upper) - temp_upper = virtual_temperature(temp_upper, rv_upper).to(temp_lower.units) # Return profile pieces return (press_lower[:-1], press_lcl, press_upper[1:], diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index cc8e5827bfc..b509993d29c 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -1261,8 +1261,8 @@ def test_surface_based_cape_cin(array_class): temperature = array_class([22.2, 14.6, 12., 9.4, 7., -38.], units.celsius) dewpoint = array_class([19., -11.2, -10.8, -10.4, -10., -53.2], units.celsius) cape, cin = surface_based_cape_cin(p, temperature, dewpoint) - assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 327.1348028 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -40.0654774 * units('joule / kilogram'), 2) def test_surface_based_cape_cin_with_xarray(): @@ -1285,8 +1285,8 @@ def test_surface_based_cape_cin_with_xarray(): data['temperature'], data['dewpoint'] ) - assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 327.1348028 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -40.0654774 * units('joule / kilogram'), 2) def test_profile_with_nans(): @@ -1326,8 +1326,8 @@ def test_most_unstable_cape_cin_surface(): temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 75.0535446 * units('joule / kilogram'), 2) - assert_almost_equal(mucin, -136.685967 * units('joule / kilogram'), 2) + assert_almost_equal(mucape, 327.1348028 * units('joule / kilogram'), 2) + assert_almost_equal(mucin, -40.0654774 * units('joule / kilogram'), 2) def test_most_unstable_cape_cin(): @@ -1336,8 +1336,8 @@ def test_most_unstable_cape_cin(): temperature = np.array([18.2, 22.2, 17.4, 10., 0., 15]) * units.celsius dewpoint = np.array([19., 19., 14.3, 0., -10., 0.]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 157.11404 * units('joule / kilogram'), 4) - assert_almost_equal(mucin, -31.8406578 * units('joule / kilogram'), 4) + assert_almost_equal(mucape, 174.66467 * units('joule / kilogram'), 4) + assert_almost_equal(mucin, 0 * units('joule / kilogram'), 4) def test_mixed_parcel(): From 7b7f23e6a69956412d3efeea1eb288e64537c76e Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 28 Mar 2023 14:57:34 -0500 Subject: [PATCH 04/30] fix flake8 in thermo --- src/metpy/calc/thermo.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 0e4da11ecd6..b00d2f20274 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -935,7 +935,7 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, Pressure increases between at least two points in your sounding. Using scipy.signal.medfilt may fix this.""" raise InvalidSoundingError(msg) - + if assumptions.use_virtual_temperature: specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint) mixing_ratio = mixing_ratio_from_specific_humidity(specific_humidity) @@ -971,7 +971,6 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, 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:], temp_lower[:-1], temp_lcl, temp_upper[1:]) From 036be157ab580005bb0097446b380908751db6cc Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 28 Mar 2023 19:18:22 -0500 Subject: [PATCH 05/30] Add suggestions from Dan --- src/metpy/calc/thermo.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index b00d2f20274..07ff44ee879 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -936,11 +936,6 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, Using scipy.signal.medfilt may fix this.""" raise InvalidSoundingError(msg) - if assumptions.use_virtual_temperature: - specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint) - mixing_ratio = mixing_ratio_from_specific_humidity(specific_humidity) - temperature = virtual_temperature(temperature, mixing_ratio) - # Find the LCL press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpoint) press_lcl = press_lcl.to(pressure.units) @@ -951,6 +946,14 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl)) temp_lower = dry_lapse(press_lower, temperature) + # Do the virtual temperature correction for the parcel below the LCL + if assumptions.use_virtual_temperature: + # Calculate the relative humidity, mixing ratio, and virtual temperature + rh_lowest = relative_humidity_from_dewpoint(temperature, dewpoint) + rv_lower = mixing_ratio_from_relative_humidity(press_lower[0], temperature, rh_lowest) + temp_lower = virtual_temperature(temp_lower, rv_lower).to(temp_lower.units) + temp_lcl = virtual_temperature(temp_lcl, rv_lower) + # If the pressure profile doesn't make it to the lcl, we can stop here if _greater_or_close(np.nanmin(pressure), press_lcl): return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units), @@ -971,6 +974,13 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units) temp_upper = temp_upper[::-1][indices] + # Do the virtual temperature correction for the parcel above the LCL + if assumptions.use_virtual_temperature: + + # Calculate the mixing ratio and virtual temperature + rv_upper = mixing_ratio(saturation_vapor_pressure(temp_upper), press_upper) + temp_upper = virtual_temperature(temp_upper, rv_upper).to(temp_lower.units) + # Return profile pieces return (press_lower[:-1], press_lcl, press_upper[1:], temp_lower[:-1], temp_lcl, temp_upper[1:]) From 84fccfe1d0c1ecc1394773450149322b4c4e383d Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 28 Mar 2023 19:52:23 -0500 Subject: [PATCH 06/30] compute virtual temperature at the end --- src/metpy/calc/thermo.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 07ff44ee879..c0cfe30ee7c 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -946,14 +946,6 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl)) temp_lower = dry_lapse(press_lower, temperature) - # Do the virtual temperature correction for the parcel below the LCL - if assumptions.use_virtual_temperature: - # Calculate the relative humidity, mixing ratio, and virtual temperature - rh_lowest = relative_humidity_from_dewpoint(temperature, dewpoint) - rv_lower = mixing_ratio_from_relative_humidity(press_lower[0], temperature, rh_lowest) - temp_lower = virtual_temperature(temp_lower, rv_lower).to(temp_lower.units) - temp_lcl = virtual_temperature(temp_lcl, rv_lower) - # If the pressure profile doesn't make it to the lcl, we can stop here if _greater_or_close(np.nanmin(pressure), press_lcl): return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units), @@ -977,7 +969,13 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, # Do the virtual temperature correction for the parcel above the LCL if assumptions.use_virtual_temperature: - # Calculate the mixing ratio and virtual temperature + # Handle below the LCL + rh_lowest = relative_humidity_from_dewpoint(temperature, dewpoint) + rv_lower = mixing_ratio_from_relative_humidity(press_lower[0], temperature, rh_lowest) + temp_lower = virtual_temperature(temp_lower, rv_lower).to(temp_lower.units) + temp_lcl = virtual_temperature(temp_lcl, rv_lower) + + # Handle above the LCL rv_upper = mixing_ratio(saturation_vapor_pressure(temp_upper), press_upper) temp_upper = virtual_temperature(temp_upper, rv_upper).to(temp_lower.units) From 62d49856c94bc02d500b40d12936ced422deb3b1 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 12:44:34 -0500 Subject: [PATCH 07/30] revert to original thermo/tests --- tests/calc/test_thermo.py | 264 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 251 insertions(+), 13 deletions(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index b509993d29c..3850b5bc7a7 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -10,7 +10,7 @@ import xarray as xr from metpy.calc import (brunt_vaisala_frequency, brunt_vaisala_frequency_squared, - brunt_vaisala_period, cape_cin, cross_totals, density, dewpoint, + brunt_vaisala_period, cape_cin, ccl, cross_totals, density, dewpoint, dewpoint_from_relative_humidity, dewpoint_from_specific_humidity, dry_lapse, dry_static_energy, el, equivalent_potential_temperature, exner_function, gradient_richardson_number, InvalidSoundingError, @@ -27,7 +27,7 @@ saturation_equivalent_potential_temperature, saturation_mixing_ratio, saturation_vapor_pressure, showalter_index, specific_humidity_from_dewpoint, specific_humidity_from_mixing_ratio, - static_stability, surface_based_cape_cin, + static_stability, surface_based_cape_cin, sweat_index, temperature_from_potential_temperature, thickness_hydrostatic, thickness_hydrostatic_from_relative_humidity, total_totals_index, vapor_pressure, vertical_totals, vertical_velocity, @@ -35,7 +35,7 @@ virtual_temperature, wet_bulb_temperature) from metpy.calc.thermo import _find_append_zero_crossings from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_nan -from metpy.units import masked_array, units +from metpy.units import is_quantity, masked_array, units def test_relative_humidity_from_dewpoint(): @@ -399,6 +399,129 @@ def test_lcl_nans(): np.nan, 18.82281982535794]) * units.degC) +def test_ccl_basic(): + """First test of CCL calculation. Data: ILX, June 17 2022 00Z.""" + pressure = np.array([993.0, 984.0, 957.0, 948.0, 925.0, 917.0, 886.0, 868.0, 850.0, + 841.0, 813.0, 806.0, 798.0, 738.0, 732.0, 723.0, 716.0, 711.0, + 700.0, 623.0, 621.0, 582.0, 541.0, 500.0, 468.0]) * units.mbar + temperature = np.array([34.6, 33.7, 31.1, 30.1, 27.8, 27.1, 24.3, 22.6, 21.4, + 20.8, 19.6, 19.4, 18.7, 13.0, 13.0, 13.4, 13.5, 13.6, + 13.0, 5.2, 5.0, 1.5, -2.4, -6.7, -10.7]) * units.degC + dewpoint = np.array([19.6, 19.4, 18.7, 18.4, 17.8, 17.5, 16.3, 15.6, 12.4, 10.8, + -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, + -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) * units.degC + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, 763.006048 * units.mbar, 5) + assert_almost_equal(ccl_t, 15.429946 * units.degC, 5) + assert_almost_equal(t_c, 37.991498 * units.degC, 5) + + +def test_ccl_nans(): + """Tests CCL handles nans.""" + pressure = np.array([993.0, 984.0, 957.0, np.nan, 925.0, 917.0, np.nan, 868.0, 850.0, + 841.0, 813.0, 806.0, 798.0, 738.0, 732.0, 723.0, 716.0, 711.0, + 700.0, 623.0, 621.0, 582.0, 541.0, 500.0, 468.0]) * units.mbar + temperature = np.array([34.6, np.nan, 31.1, np.nan, 27.8, 27.1, 24.3, 22.6, 21.4, + 20.8, 19.6, 19.4, 18.7, 13.0, 13.0, 13.4, 13.5, 13.6, + 13.0, 5.2, 5.0, 1.5, -2.4, -6.7, -10.7]) * units.degC + dewpoint = np.array([19.6, 19.4, 18.7, np.nan, 17.8, 17.5, 16.3, 15.6, 12.4, 10.8, + -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, + -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) * units.degC + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, 763.006048 * units.mbar, 5) + assert_almost_equal(ccl_t, 15.429946 * units.degC, 5) + assert_almost_equal(t_c, 37.991498 * units.degC, 5) + + +def test_ccl_unit(): + """Tests CCL pressure and temperature is returned in the correct unit.""" + pressure = (np.array([993.0, 984.0, 957.0, 948.0, 925.0, 917.0, 886.0, 868.0, 850.0, + 841.0, 813.0, 806.0, 798.0, 738.0, 732.0, 723.0, 716.0, 711.0, + 700.0, 623.0, 621.0, 582.0, 541.0, 500.0, 468.0]) * 100) * units.Pa + temperature = (np.array([34.6, 33.7, 31.1, 30.1, 27.8, 27.1, 24.3, 22.6, 21.4, + 20.8, 19.6, 19.4, 18.7, 13.0, 13.0, 13.4, 13.5, 13.6, + 13.0, 5.2, 5.0, 1.5, -2.4, -6.7, -10.7]) + 273.15) * units.kelvin + dewpoint = (np.array([19.6, 19.4, 18.7, 18.4, 17.8, 17.5, 16.3, 15.6, 12.4, 10.8, + -0.4, -3.6, -3.8, -5.0, -6.0, -15.6, -13.2, -11.4, -11.0, + -5.8, -6.2, -14.8, -24.3, -34.7, -38.1]) + 273.15) * units.kelvin + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, (763.006048 * 100) * units.Pa, 3) + assert_almost_equal(ccl_t, (15.429946 + 273.15) * units.kelvin, 3) + assert_almost_equal(t_c, (37.991498 + 273.15) * units.kelvin, 3) + + assert ccl_p.units == pressure.units + assert ccl_t.units == temperature.units + assert t_c.units == temperature.units + + +def test_multiple_ccl(): + """Tests the case where there are multiple CCLs. Data: BUF, May 18 2022 12Z.""" + pressure = np.array([992.0, 990.0, 983.0, 967.0, 950.0, 944.0, 928.0, 925.0, 922.0, + 883.0, 877.7, 858.0, 853.0, 850.0, 835.0, 830.0, 827.0, 826.0, + 813.6, 808.0, 799.0, 784.0, 783.3, 769.0, 760.0, 758.0, 754.0, + 753.0, 738.0, 725.7, 711.0, 704.0, 700.0, 685.0, 672.0, 646.6, + 598.6, 596.0, 587.0, 582.0, 567.0, 560.0, 555.0, 553.3, 537.0, + 526.0, 521.0, 519.0, 515.0, 500.0]) * units.mbar + temperature = np.array([6.8, 6.2, 7.8, 7.6, 7.2, 7.6, 6.6, 6.4, 6.2, 3.2, 2.8, 1.2, + 1.0, 0.8, -0.3, -0.1, 0.4, 0.6, 0.9, 1.0, 0.6, -0.3, -0.3, + -0.7, -1.5, -1.3, 0.2, 0.2, -1.1, -2.1, -3.3, -2.3, -1.7, 0.2, + -0.9, -3.0, -7.3, -7.5, -8.1, -8.3, -9.5, -10.1, -10.7, + -10.8, -12.1, -12.5, -12.7, -12.9, -13.5, -15.5]) * units.degC + dewpoint = np.array([5.1, 5.0, 4.2, 2.7, 2.2, 0.6, -2.4, -2.6, -2.8, -3.8, -3.6, + -3.1, -5.0, -4.2, -1.8, -4.3, -7.6, -6.4, -8.2, -9.0, -10.4, + -9.3, -9.6, -14.7, -11.5, -12.3, -25.8, -25.8, -19.1, -19.6, + -20.3, -42.3, -39.7, -46.8, -46.8, -46.7, -46.5, -46.5, + -52.1, -36.3, -47.5, -30.1, -29.7, -30.4, -37.1, -49.5, + -36.7, -28.9, -28.5, -22.5]) * units.degC + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint) + assert_almost_equal(ccl_p, 680.191653 * units.mbar, 5) + assert_almost_equal(ccl_t, -0.204408 * units.degC, 5) + assert_almost_equal(t_c, 30.8678258 * units.degC, 5) + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, which='bottom') + assert_almost_equal(ccl_p, 886.835325 * units.mbar, 5) + assert_almost_equal(ccl_t, 3.500840 * units.degC, 5) + assert_almost_equal(t_c, 12.5020423 * units.degC, 5) + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, which='all') + assert_array_almost_equal(ccl_p, np.array([886.835325, 680.191653]) * units.mbar, 5) + assert_array_almost_equal(ccl_t, np.array([3.500840, -0.204408]) * units.degC, 5) + assert_array_almost_equal(t_c, np.array([12.5020423, 30.8678258]) * units.degC, 5) + + +def test_ccl_with_ml(): + """Test CCL calculation with a specified mixed-layer depth.""" + pressure = np.array([992.0, 990.0, 983.0, 967.0, 950.0, 944.0, 928.0, 925.0, 922.0, + 883.0, 877.7, 858.0, 853.0, 850.0, 835.0, 830.0, 827.0, 826.0, + 813.6, 808.0, 799.0, 784.0, 783.3, 769.0, 760.0, 758.0, 754.0, + 753.0, 738.0, 725.7, 711.0, 704.0, 700.0, 685.0, 672.0, 646.6, + 598.6, 596.0, 587.0, 582.0, 567.0, 560.0, 555.0, 553.3, 537.0, + 526.0, 521.0, 519.0, 515.0, 500.0]) * units.mbar + temperature = np.array([6.8, 6.2, 7.8, 7.6, 7.2, 7.6, 6.6, 6.4, 6.2, 3.2, 2.8, 1.2, + 1.0, 0.8, -0.3, -0.1, 0.4, 0.6, 0.9, 1.0, 0.6, -0.3, -0.3, + -0.7, -1.5, -1.3, 0.2, 0.2, -1.1, -2.1, -3.3, -2.3, -1.7, 0.2, + -0.9, -3.0, -7.3, -7.5, -8.1, -8.3, -9.5, -10.1, -10.7, + -10.8, -12.1, -12.5, -12.7, -12.9, -13.5, -15.5]) * units.degC + dewpoint = np.array([5.1, 5.0, 4.2, 2.7, 2.2, 0.6, -2.4, -2.6, -2.8, -3.8, -3.6, + -3.1, -5.0, -4.2, -1.8, -4.3, -7.6, -6.4, -8.2, -9.0, -10.4, + -9.3, -9.6, -14.7, -11.5, -12.3, -25.8, -25.8, -19.1, -19.6, + -20.3, -42.3, -39.7, -46.8, -46.8, -46.7, -46.5, -46.5, + -52.1, -36.3, -47.5, -30.1, -29.7, -30.4, -37.1, -49.5, + -36.7, -28.9, -28.5, -22.5]) * units.degC + + ccl_p, ccl_t, t_c = ccl(pressure, temperature, dewpoint, + mixed_layer_depth=500 * units.m, which='all') + + assert_array_almost_equal(ccl_p, np.array( + [850.600930, 784.325312, 737.767377, 648.076147]) * units.mbar, 5) + assert_array_almost_equal(ccl_t, np.array( + [0.840118, -0.280299, -1.118757, -2.875716]) * units.degC, 5) + assert_array_almost_equal(t_c, np.array( + [13.146845, 18.661621, 22.896152, 32.081388]) * units.degC, 5) + + def test_lfc_basic(): """Test LFC calculation.""" levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar @@ -632,7 +755,7 @@ def test_equivalent_potential_temperature_masked(): np.ma.array([311.18586, 313.51781, 315.93971], mask=[False, True, False]), units.kelvin ) - assert isinstance(ept, units.Quantity) + assert is_quantity(ept) assert isinstance(ept.m, np.ma.MaskedArray) assert_array_almost_equal(ept, expected, 3) @@ -656,7 +779,7 @@ def test_saturation_equivalent_potential_temperature_masked(): np.ma.array([335.02750, 338.95813, 343.08740]), units.kelvin ) - assert isinstance(s_ept, units.Quantity) + assert is_quantity(s_ept) assert isinstance(s_ept.m, np.ma.MaskedArray) assert_array_almost_equal(s_ept, expected, 3) @@ -1261,8 +1384,8 @@ def test_surface_based_cape_cin(array_class): temperature = array_class([22.2, 14.6, 12., 9.4, 7., -38.], units.celsius) dewpoint = array_class([19., -11.2, -10.8, -10.4, -10., -53.2], units.celsius) cape, cin = surface_based_cape_cin(p, temperature, dewpoint) - assert_almost_equal(cape, 327.1348028 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -40.0654774 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2) def test_surface_based_cape_cin_with_xarray(): @@ -1285,8 +1408,8 @@ def test_surface_based_cape_cin_with_xarray(): data['temperature'], data['dewpoint'] ) - assert_almost_equal(cape, 327.1348028 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -40.0654774 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2) def test_profile_with_nans(): @@ -1326,8 +1449,8 @@ def test_most_unstable_cape_cin_surface(): temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 327.1348028 * units('joule / kilogram'), 2) - assert_almost_equal(mucin, -40.0654774 * units('joule / kilogram'), 2) + assert_almost_equal(mucape, 75.0535446 * units('joule / kilogram'), 2) + assert_almost_equal(mucin, -136.685967 * units('joule / kilogram'), 2) def test_most_unstable_cape_cin(): @@ -1336,8 +1459,8 @@ def test_most_unstable_cape_cin(): temperature = np.array([18.2, 22.2, 17.4, 10., 0., 15]) * units.celsius dewpoint = np.array([19., 19., 14.3, 0., -10., 0.]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 174.66467 * units('joule / kilogram'), 4) - assert_almost_equal(mucin, 0 * units('joule / kilogram'), 4) + assert_almost_equal(mucape, 157.11404 * units('joule / kilogram'), 4) + assert_almost_equal(mucin, -31.8406578 * units('joule / kilogram'), 4) def test_mixed_parcel(): @@ -1894,6 +2017,40 @@ def test_lcl_grid_surface_lcls(): assert_array_almost_equal(lcl_temperature, temp_truth, 4) +@pytest.fixture() +def index_xarray_data(): + """Create data for testing that index calculations work with xarray data.""" + pressure = xr.DataArray([850., 700., 500.], dims=('isobaric',), attrs={'units': 'hPa'}) + temp = xr.DataArray([[[[296., 295., 294.], [293., 292., 291.]], + [[286., 285., 284.], [283., 282., 281.]], + [[276., 275., 274.], [273., 272., 271.]]]] * units.K, + dims=('time', 'isobaric', 'y', 'x')) + + profile = xr.DataArray([[[[289., 288., 287.], [286., 285., 284.]], + [[279., 278., 277.], [276., 275., 274.]], + [[269., 268., 267.], [266., 265., 264.]]]] * units.K, + dims=('time', 'isobaric', 'y', 'x')) + + dewp = xr.DataArray([[[[294., 293., 292.], [291., 290., 289.]], + [[284., 283., 282.], [281., 280., 279.]], + [[274., 273., 272.], [271., 270., 269.]]]] * units.K, + dims=('time', 'isobaric', 'y', 'x')) + + dirw = xr.DataArray([[[[180., 180., 180.], [180., 180., 180.]], + [[225., 225., 225.], [225., 225., 225.]], + [[270., 270., 270.], [270., 270., 270.]]]] * units.degree, + dims=('time', 'isobaric', 'y', 'x')) + + speed = xr.DataArray([[[[20., 20., 20.], [20., 20., 20.]], + [[25., 25., 25.], [25., 25., 25.]], + [[50., 50., 50.], [50., 50., 50.]]]] * units.knots, + dims=('time', 'isobaric', 'y', 'x')) + + return xr.Dataset({'temperature': temp, 'profile': profile, 'dewpoint': dewp, + 'wind_direction': dirw, 'wind_speed': speed}, + coords={'isobaric': pressure, 'time': ['2020-01-01T00:00Z']}) + + def test_lifted_index(): """Test the Lifted Index calculation.""" pressure = np.array([1014., 1000., 997., 981.2, 947.4, 925., 914.9, 911., @@ -1944,6 +2101,13 @@ def test_lifted_index_500hpa_missing(): assert_almost_equal(li, -7.9176350 * units.delta_degree_Celsius, 1) +def test_lifted_index_xarray(index_xarray_data): + """Test lifted index with a grid of xarray data.""" + result = lifted_index(index_xarray_data.isobaric, index_xarray_data.temperature, + index_xarray_data.profile) + assert_array_almost_equal(result, np.full((1, 1, 2, 3), 7) * units.delta_degC) + + def test_k_index(): """Test the K Index calculation.""" pressure = np.array([1014., 1000., 997., 981.2, 947.4, 925., 914.9, 911., @@ -1968,6 +2132,14 @@ def test_k_index(): assert_almost_equal(ki, 33.5 * units.degC, 2) +def test_k_index_xarray(index_xarray_data): + """Test the K index calculation with a grid of xarray data.""" + result = k_index(index_xarray_data.isobaric, index_xarray_data.temperature, + index_xarray_data.dewpoint) + assert_array_almost_equal(result, + np.array([[[312., 311., 310.], [309., 308., 307.]]]) * units.K) + + def test_gradient_richardson_number(): """Test gradient Richardson number calculation.""" theta = units('K') * np.asarray([254.5, 258.3, 262.2]) @@ -2050,6 +2222,13 @@ def test_total_totals_index(): assert_almost_equal(tt, 45.10 * units.delta_degC, 2) +def test_total_totals_index_xarray(index_xarray_data): + """Test the total totals index calculation with a grid of xarray data.""" + result = total_totals_index(index_xarray_data.isobaric, index_xarray_data.temperature, + index_xarray_data.dewpoint) + assert_array_almost_equal(result, np.full((1, 2, 3), 38.) * units.K) + + def test_vertical_totals(): """Test the Vertical Totals calculation.""" pressure = np.array([1008., 1000., 947., 925., 921., 896., 891., 889., 866., @@ -2069,6 +2248,12 @@ def test_vertical_totals(): assert_almost_equal(vt, 23.70 * units.delta_degC, 2) +def test_vertical_totals_index_xarray(index_xarray_data): + """Test the vertical totals index calculation with a grid of xarray data.""" + result = vertical_totals(index_xarray_data.isobaric, index_xarray_data.temperature) + assert_array_almost_equal(result, np.full((1, 2, 3), 20.) * units.K) + + def test_cross_totals(): """Test the Cross Totals calculation.""" pressure = np.array([1008., 1000., 947., 925., 921., 896., 891., 889., 866., @@ -2094,6 +2279,13 @@ def test_cross_totals(): assert_almost_equal(ct, 21.40 * units.delta_degC, 2) +def test_cross_totals_index_xarray(index_xarray_data): + """Test the cross totals index calculation with a grid of xarray data.""" + result = cross_totals(index_xarray_data.isobaric, index_xarray_data.temperature, + index_xarray_data.dewpoint) + assert_array_almost_equal(result, np.full((1, 2, 3), 18.) * units.K) + + 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., @@ -2157,3 +2349,49 @@ def test_parcel_profile_with_lcl_as_dataset_duplicates(): profile = parcel_profile_with_lcl_as_dataset(pressure, temperature, dewpoint) xr.testing.assert_allclose(profile, truth, atol=1e-5) + + +def test_sweat_index(): + """Test the SWEAT Index calculation.""" + pressure = np.array([1008., 1000., 947., 925., 921., 896., 891., 889., 866., + 858., 850., 835., 820., 803., 733., 730., 700., 645., + 579., 500., 494., 466., 455., 441., 433., 410., 409., + 402., 400., 390., 388., 384., 381., 349., 330., 320., + 306., 300., 278., 273., 250., 243., 208., 200., 196., + 190., 179., 159., 151., 150., 139.]) * units.hPa + temperature = np.array([27.4, 26.4, 22.9, 21.4, 21.2, 20.7, 20.6, 21.2, 19.4, + 19.1, 18.8, 17.8, 17.4, 16.3, 11.4, 11.2, 10.2, 6.1, + 0.6, -4.9, -5.5, -8.5, -9.9, -11.7, -12.3, -13.7, -13.8, + -14.9, -14.9, -16.1, -16.1, -16.9, -17.3, -21.7, -24.5, -26.1, + -28.3, -29.5, -33.1, -34.2, -39.3, -41., -50.2, -52.5, -53.5, + -55.2, -58.6, -65.2, -68.1, -68.5, -72.5]) * units.degC + dewpoint = np.array([24.9, 24.6, 22., 20.9, 20.7, 14.8, 13.6, 12.2, 16.8, + 16.6, 16.5, 15.9, 13.6, 13.2, 11.3, 11.2, 8.6, 4.5, + -0.8, -8.1, -9.5, -12.7, -12.7, -12.8, -13.1, -24.7, -24.4, + -21.9, -24.9, -36.1, -31.1, -26.9, -27.4, -33., -36.5, -47.1, + -31.4, -33.5, -40.1, -40.8, -44.1, -45.6, -54., -56.1, -56.9, + -58.6, -61.9, -68.4, -71.2, -71.6, -77.2]) * units.degC + speed = np.array([0., 3., 10., 12., 12., 14., 14., 14., 12., + 12., 12., 12., 11., 11., 12., 12., 10., 10., + 8., 5., 4., 1., 0., 3., 5., 10., 10., + 11., 11., 13., 14., 14., 15., 23., 23., 24., + 24., 24., 26., 27., 28., 30., 25., 24., 26., + 28., 33., 29., 32., 26., 26.]) * units.knot + direction = np.array([0., 170., 200., 205., 204., 200., 197., 195., 180., + 175., 175., 178., 181., 185., 160., 160., 165., 165., + 203., 255., 268., 333., 0., 25., 40., 83., 85., + 89., 90., 100., 103., 107., 110., 90., 88., 87., + 86., 85., 85., 85., 60., 55., 60., 50., 46., + 40., 45., 35., 50., 50., 50.]) * units.degree + + sweat = sweat_index(pressure, temperature, dewpoint, speed, direction) + assert_almost_equal(sweat, 227., 2) + + +def test_sweat_index_xarray(index_xarray_data): + """Test the SWEAT index calculation with a grid of xarray data.""" + result = sweat_index(index_xarray_data.isobaric, index_xarray_data.temperature, + index_xarray_data.dewpoint, index_xarray_data.wind_speed, + index_xarray_data.wind_direction) + assert_array_almost_equal(result, np.array([[[[490.2, 478.2, 466.2], + [454.2, 442.2, 430.2]]]])) \ No newline at end of file From 911a9410129e86a404437691749323a62bddeac7 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 12:45:20 -0500 Subject: [PATCH 08/30] add reverted thermo module --- src/metpy/calc/thermo.py | 1105 ++++++++++++++++++++++++++++++++++---- 1 file changed, 1006 insertions(+), 99 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index c0cfe30ee7c..8b3dd99e74b 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -44,6 +44,13 @@ def relative_humidity_from_dewpoint(temperature, dewpoint): `pint.Quantity` Relative humidity + Examples + -------- + >>> from metpy.calc import relative_humidity_from_dewpoint + >>> from metpy.units import units + >>> relative_humidity_from_dewpoint(25 * units.degC, 12 * units.degC).to('percent') + + See Also -------- saturation_vapor_pressure @@ -220,6 +227,14 @@ def dry_lapse(pressure, temperature, reference_pressure=None, vertical_dim=0): `pint.Quantity` The parcel's resulting temperature at levels given by ``pressure`` + Examples + -------- + >>> from metpy.calc import dry_lapse + >>> from metpy.units import units + >>> plevs = [1000, 925, 850, 700] * units.hPa + >>> dry_lapse(plevs, 15 * units.degC).to('degC') + + See Also -------- moist_lapse : Calculate parcel temperature assuming liquid saturation processes @@ -277,6 +292,15 @@ def moist_lapse(pressure, temperature, reference_pressure=None): `pint.Quantity` The resulting parcel temperature at levels given by `pressure` + Examples + -------- + >>> from metpy.calc import moist_lapse + >>> from metpy.units import units + >>> plevs = [925, 850, 700, 500, 300, 200] * units.hPa + >>> moist_lapse(plevs, 5 * units.degC).to('degC') + + See Also -------- dry_lapse : Calculate parcel temperature assuming dry adiabatic processes @@ -401,6 +425,13 @@ def lcl(pressure, temperature, dewpoint, max_iters=50, eps=1e-5): eps : float, optional The desired relative error in the calculated value, defaults to 1e-5. + Examples + -------- + >>> from metpy.calc import lcl + >>> from metpy.units import units + >>> lcl(943 * units.hPa, 33 * units.degC, 28 * units.degC) + (, ) + See Also -------- parcel_profile @@ -442,6 +473,109 @@ def _lcl_iter(p, p0, w, t): return lcl_p, globals()['dewpoint']._nounit(vapor_pressure._nounit(lcl_p, w)) +@exporter.export +@preprocess_and_wrap() +@check_units('[pressure]', '[temperature]', '[temperature]') +def ccl(pressure, temperature, dewpoint, height=None, mixed_layer_depth=None, which='top'): + r"""Calculate the convective condensation level (CCL) and convective temperature. + + This function is implemented directly based on the definition of the CCL, + as in [USAF1990]_, and finding where the ambient temperature profile intersects + the line of constant mixing ratio starting at the surface, using the surface dewpoint + or the average dewpoint of a shallow layer near the surface. + + Parameters + ---------- + pressure : `pint.Quantity` + Atmospheric pressure profile + + temperature : `pint.Quantity` + Temperature at the levels given by `pressure` + + dewpoint : `pint.Quantity` + Dewpoint at the levels given by `pressure` + + height : `pint.Quantity`, optional + Atmospheric heights at the levels given by `pressure`. + Only needed when specifying a mixed layer depth as a height. + + mixed_layer_depth : `pint.Quantity`, optional + The thickness of the mixed layer as a pressure or height above the bottom + of the layer (default None). + + which: str, optional + Pick which CCL value to return; must be one of 'top', 'bottom', or 'all'. + 'top' returns the lowest-pressure CCL (default), + 'bottom' returns the highest-pressure CCL, + 'all' returns every CCL in a `Pint.Quantity` array. + + Returns + ------- + `pint.Quantity` + CCL Pressure + + `pint.Quantity` + CCL Temperature + + `pint.Quantity` + Convective Temperature + + See Also + -------- + lcl, lfc, el + + Notes + ----- + Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). + Since this function returns scalar values when given a profile, this will return Pint + Quantities even when given xarray DataArray profiles. + + Examples + -------- + >>> import metpy.calc as mpcalc + >>> from metpy.units import units + >>> pressure = [993, 957, 925, 886, 850, 813, 798, 732, 716, 700] * units.mbar + >>> temperature = [34.6, 31.1, 27.8, 24.3, 21.4, 19.6, 18.7, 13, 13.5, 13] * units.degC + >>> dewpoint = [19.6, 18.7, 17.8, 16.3, 12.4, -0.4, -3.8, -6, -13.2, -11] * units.degC + >>> ccl_p, ccl_t, t_c = mpcalc.ccl(pressure, temperature, dewpoint) + >>> ccl_p, t_c + (, ) + """ + pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint) + + # If the mixed layer is not defined, take the starting dewpoint to be the + # first element of the dewpoint array and calculate the corresponding mixing ratio. + if mixed_layer_depth is None: + p_start, dewpoint_start = pressure[0], dewpoint[0] + vapor_pressure_start = saturation_vapor_pressure(dewpoint_start) + r_start = mixing_ratio(vapor_pressure_start, p_start) + + # Else, calculate the mixing ratio of the mixed layer. + else: + vapor_pressure_profile = saturation_vapor_pressure(dewpoint) + r_profile = mixing_ratio(vapor_pressure_profile, pressure) + r_start = mixed_layer(pressure, r_profile, height=height, + depth=mixed_layer_depth)[0] + + # rt_profile is the temperature-pressure profile with a fixed mixing ratio + rt_profile = globals()['dewpoint'](vapor_pressure(pressure, r_start)) + + x, y = find_intersections(pressure, rt_profile, temperature, + direction='increasing', log_x=True) + + # In the case of multiple CCLs, select which to return + if which == 'top': + x, y = x[-1], y[-1] + elif which == 'bottom': + x, y = x[0], y[0] + elif which not in ['top', 'bottom', 'all']: + raise ValueError(f'Invalid option for "which": {which}. Valid options are ' + '"top", "bottom", and "all".') + + x, y = x.to(pressure.units), y.to(temperature.units) + return x, y, dry_lapse(pressure[0], y, x).to(temperature.units) + + @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[temperature]', '[temperature]', '[temperature]') @@ -489,6 +623,30 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi `pint.Quantity` LFC temperature, or array of same if which='all' + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, lfc + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # calculate LFC + >>> lfc(p, T, Td) + (, ) + See Also -------- parcel_profile @@ -674,6 +832,32 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which=' `pint.Quantity` EL temperature, or array of same if which='all' + Examples + -------- + >>> from metpy.calc import el, dewpoint_from_relative_humidity, parcel_profile + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compute parcel profile temperature + >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') + >>> # calculate EL + >>> el(p, T, Td, prof) + (, ) + See Also -------- parcel_profile @@ -715,19 +899,10 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which=' units.Quantity(np.nan, temperature.units)) -class ParcelPathAssumptions(): - """Holds assumptions made about the parcel path during calculations.""" - - def __init__(self): - self.use_virtual_temperature = True - self.moist_adiabat = 'pseudoadiabatic' - - @exporter.export @preprocess_and_wrap(wrap_like='pressure') @check_units('[pressure]', '[temperature]', '[temperature]') -def parcel_profile(pressure, temperature, dewpoint, - assumptions=ParcelPathAssumptions()): +def parcel_profile(pressure, temperature, dewpoint): r"""Calculate the profile a parcel takes through the atmosphere. The parcel starts at `temperature`, and `dewpoint`, lifted up @@ -751,6 +926,35 @@ def parcel_profile(pressure, temperature, dewpoint, `pint.Quantity` The parcel's temperatures at the specified pressure levels + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, parcel_profile + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # computer parcel temperature + >>> parcel_profile(p, T[0], Td[0]).to('degC') + + See Also -------- lcl, moist_lapse, dry_lapse, parcel_profile_with_lcl, parcel_profile_with_lcl_as_dataset @@ -769,8 +973,7 @@ def parcel_profile(pressure, temperature, dewpoint, Renamed ``dewpt`` parameter to ``dewpoint`` """ - _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint, - assumptions=assumptions) + _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint) return concatenate((t_l, t_u)) @@ -814,6 +1017,36 @@ def parcel_profile_with_lcl(pressure, temperature, dewpoint): The parcel profile temperatures at all of the levels in the returned pressures array, including the LCL + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, parcel_profile_with_lcl + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # computer parcel temperature + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> p_wLCL, T_wLCL, Td_wLCL, prof_wLCL = parcel_profile_with_lcl(p, T, Td) + >>> print(f'Shape of original pressure array: {p.shape}') + Shape of original pressure array: (30,) + >>> print(f'Shape of pressure array from function: {p_wLCL.shape}') + Shape of pressure array from function: (31,) + >>> print(p == p_wLCL) + False + See Also -------- lcl, moist_lapse, dry_lapse, parcel_profile, parcel_profile_with_lcl_as_dataset @@ -921,8 +1154,7 @@ def _check_pressure(pressure): return np.all(pressure[:-1] >= pressure[1:]) -def _parcel_profile_helper(pressure, temperature, dewpoint, - assumptions=ParcelPathAssumptions()): +def _parcel_profile_helper(pressure, temperature, dewpoint): """Help calculate parcel profiles. Returns the temperature and pressure, above, below, and including the LCL. The @@ -966,19 +1198,6 @@ def _parcel_profile_helper(pressure, temperature, dewpoint, temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units) temp_upper = temp_upper[::-1][indices] - # Do the virtual temperature correction for the parcel above the LCL - if assumptions.use_virtual_temperature: - - # Handle below the LCL - rh_lowest = relative_humidity_from_dewpoint(temperature, dewpoint) - rv_lower = mixing_ratio_from_relative_humidity(press_lower[0], temperature, rh_lowest) - temp_lower = virtual_temperature(temp_lower, rv_lower).to(temp_lower.units) - temp_lcl = virtual_temperature(temp_lcl, rv_lower) - - # Handle above the LCL - rv_upper = mixing_ratio(saturation_vapor_pressure(temp_upper), press_upper) - temp_upper = virtual_temperature(temp_upper, rv_upper).to(temp_lower.units) - # Return profile pieces return (press_lower[:-1], press_lcl, press_upper[1:], temp_lower[:-1], temp_lcl, temp_upper[1:]) @@ -1016,6 +1235,13 @@ def vapor_pressure(pressure, mixing_ratio): `pint.Quantity` Ambient water vapor (partial) pressure in the same units as ``pressure`` + Examples + -------- + >>> from metpy.calc import vapor_pressure + >>> from metpy.units import units + >>> vapor_pressure(988 * units.hPa, 18 * units('g/kg')).to('hPa') + + See Also -------- saturation_vapor_pressure, dewpoint @@ -1050,6 +1276,13 @@ def saturation_vapor_pressure(temperature): `pint.Quantity` Saturation water vapor (partial) pressure + Examples + -------- + >>> from metpy.calc import saturation_vapor_pressure + >>> from metpy.units import units + >>> saturation_vapor_pressure(25 * units.degC).to('hPa') + + See Also -------- vapor_pressure, dewpoint @@ -1089,6 +1322,12 @@ def dewpoint_from_relative_humidity(temperature, relative_humidity): `pint.Quantity` Dewpoint temperature + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity + >>> from metpy.units import units + >>> dewpoint_from_relative_humidity(10 * units.degC, 50 * units.percent) + .. versionchanged:: 1.0 Renamed ``rh`` parameter to ``relative_humidity`` @@ -1119,6 +1358,13 @@ def dewpoint(vapor_pressure): `pint.Quantity` Dewpoint temperature + Examples + -------- + >>> from metpy.calc import dewpoint + >>> from metpy.units import units + >>> dewpoint(22 * units.hPa) + + See Also -------- dewpoint_from_relative_humidity, saturation_vapor_pressure, vapor_pressure @@ -1175,6 +1421,13 @@ def mixing_ratio(partial_press, total_press, molecular_weight_ratio=mpconsts.nou `pint.Quantity` The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g) + Examples + -------- + >>> from metpy.calc import mixing_ratio + >>> from metpy.units import units + >>> mixing_ratio(25 * units.hPa, 1000 * units.hPa).to('g/kg') + + See Also -------- saturation_mixing_ratio, vapor_pressure @@ -1217,6 +1470,13 @@ def saturation_mixing_ratio(total_press, temperature): `pint.Quantity` Saturation mixing ratio, dimensionless + Examples + -------- + >>> from metpy.calc import saturation_mixing_ratio + >>> from metpy.units import units + >>> saturation_mixing_ratio(983 * units.hPa, 25 * units.degC).to('g/kg') + + Notes ----- This function is a straightforward implementation of the equation given in many places, @@ -1273,6 +1533,13 @@ def equivalent_potential_temperature(pressure, temperature, dewpoint): `pint.Quantity` Equivalent potential temperature of the parcel + Examples + -------- + >>> from metpy.calc import equivalent_potential_temperature + >>> from metpy.units import units + >>> equivalent_potential_temperature(850*units.hPa, 20*units.degC, 18*units.degC) + + Notes ----- [Bolton1980]_ formula for Theta-e is used, since according to @@ -1340,6 +1607,13 @@ def saturation_equivalent_potential_temperature(pressure, temperature): `pint.Quantity` Saturation equivalent potential temperature of the parcel + Examples + -------- + >>> from metpy.calc import saturation_equivalent_potential_temperature + >>> from metpy.units import units + >>> saturation_equivalent_potential_temperature(500 * units.hPa, -20 * units.degC) + + Notes ----- [Bolton1980]_ formula for Theta-e is used (for saturated case), since according to @@ -1395,6 +1669,13 @@ def virtual_temperature( `pint.Quantity` Corresponding virtual temperature of the parcel + Examples + -------- + >>> from metpy.calc import virtual_temperature + >>> from metpy.units import units + >>> virtual_temperature(283 * units.K, 12 * units('g/kg')) + + Notes ----- .. math:: T_v = T \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})} @@ -1441,6 +1722,13 @@ def virtual_potential_temperature(pressure, temperature, mixing_ratio, `pint.Quantity` Corresponding virtual potential temperature of the parcel + Examples + -------- + >>> from metpy.calc import virtual_potential_temperature + >>> from metpy.units import units + >>> virtual_potential_temperature(500 * units.hPa, -15 * units.degC, 1 * units('g/kg')) + + Notes ----- .. math:: \Theta_v = \Theta \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})} @@ -1471,7 +1759,7 @@ def density(pressure, temperature, mixing_ratio, molecular_weight_ratio=mpconsts Total atmospheric pressure temperature: `pint.Quantity` - Air temperature + Air temperature (or the virtual temperature if the mixing ratio is set to 0) mixing_ratio : `pint.Quantity` Mass mixing ratio (dimensionless) @@ -1486,9 +1774,16 @@ def density(pressure, temperature, mixing_ratio, molecular_weight_ratio=mpconsts `pint.Quantity` Corresponding density of the parcel + Examples + -------- + >>> from metpy.calc import density + >>> from metpy.units import units + >>> density(1000 * units.hPa, 10 * units.degC, 24 * units('g/kg')) + + Notes ----- - .. math:: \rho = \frac{p}{R_dT_v} + .. math:: \rho = \frac{\epsilon p\,(1+w)}{R_dT\,(w+\epsilon)} .. versionchanged:: 1.0 Renamed ``mixing`` parameter to ``mixing_ratio`` @@ -1527,6 +1822,14 @@ def relative_humidity_wet_psychrometric(pressure, dry_bulb_temperature, wet_bulb `pint.Quantity` Relative humidity + Examples + -------- + >>> from metpy.calc import relative_humidity_wet_psychrometric + >>> from metpy.units import units + >>> relative_humidity_wet_psychrometric(1000 * units.hPa, 19 * units.degC, + ... 10 * units.degC).to('percent') + + See Also -------- psychrometric_vapor_pressure_wet, saturation_vapor_pressure @@ -1581,6 +1884,18 @@ def psychrometric_vapor_pressure_wet(pressure, dry_bulb_temperature, wet_bulb_te `pint.Quantity` Vapor pressure + Examples + -------- + >>> from metpy.calc import psychrometric_vapor_pressure_wet, saturation_vapor_pressure + >>> from metpy.units import units + >>> vp = psychrometric_vapor_pressure_wet(958 * units.hPa, 25 * units.degC, + ... 12 * units.degC) + >>> print(f'Vapor Pressure: {vp:.2f}') + Vapor Pressure: 628.15 pascal + >>> rh = (vp / saturation_vapor_pressure(25 * units.degC)).to('percent') + >>> print(f'RH: {rh:.2f}') + RH: 19.83 percent + See Also -------- saturation_vapor_pressure @@ -1628,7 +1943,7 @@ def mixing_ratio_from_relative_humidity(pressure, temperature, relative_humidity temperature: `pint.Quantity` Air temperature - relative_humidity: array_like + relative_humidity: array-like The relative humidity expressed as a unitless ratio in the range [0, 1]. Can also pass a percentage if proper units are attached. @@ -1637,6 +1952,16 @@ def mixing_ratio_from_relative_humidity(pressure, temperature, relative_humidity `pint.Quantity` Mixing ratio (dimensionless) + Examples + -------- + >>> from metpy.calc import mixing_ratio_from_relative_humidity + >>> from metpy.units import units + >>> p = 1000. * units.hPa + >>> T = 28.1 * units.degC + >>> rh = .65 + >>> mixing_ratio_from_relative_humidity(p, T, rh).to('g/kg') + + See Also -------- relative_humidity_from_mixing_ratio, saturation_mixing_ratio @@ -1684,6 +2009,14 @@ def relative_humidity_from_mixing_ratio(pressure, temperature, mixing_ratio): `pint.Quantity` Relative humidity + Examples + -------- + >>> from metpy.calc import relative_humidity_from_mixing_ratio + >>> from metpy.units import units + >>> relative_humidity_from_mixing_ratio(1013.25 * units.hPa, + ... 30 * units.degC, 18/1000).to('percent') + + See Also -------- mixing_ratio_from_relative_humidity, saturation_mixing_ratio @@ -1721,6 +2054,15 @@ def mixing_ratio_from_specific_humidity(specific_humidity): `pint.Quantity` Mixing ratio + Examples + -------- + >>> from metpy.calc import mixing_ratio_from_specific_humidity + >>> from metpy.units import units + >>> sh = [4.77, 12.14, 6.16, 15.29, 12.25] * units('g/kg') + >>> mixing_ratio_from_specific_humidity(sh).to('g/kg') + + See Also -------- mixing_ratio, specific_humidity_from_mixing_ratio @@ -1754,6 +2096,13 @@ def specific_humidity_from_mixing_ratio(mixing_ratio): `pint.Quantity` Specific humidity + Examples + -------- + >>> from metpy.calc import specific_humidity_from_mixing_ratio + >>> from metpy.units import units + >>> specific_humidity_from_mixing_ratio(19 * units('g/kg')) + + See Also -------- mixing_ratio, mixing_ratio_from_specific_humidity @@ -1796,6 +2145,14 @@ def relative_humidity_from_specific_humidity(pressure, temperature, specific_hum `pint.Quantity` Relative humidity + Examples + -------- + >>> from metpy.calc import relative_humidity_from_specific_humidity + >>> from metpy.units import units + >>> relative_humidity_from_specific_humidity(1013.25 * units.hPa, + ... 30 * units.degC, 18/1000).to('percent') + + See Also -------- relative_humidity_from_mixing_ratio @@ -1861,6 +2218,32 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' `pint.Quantity` Convective Inhibition (CIN) + Examples + -------- + >>> from metpy.calc import cape_cin, dewpoint_from_relative_humidity, parcel_profile + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compture parcel temperature + >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') + >>> # calculate surface based CAPE/CIN + >>> cape_cin(p, T, Td, prof) + (, ) + See Also -------- lfc, el @@ -2021,6 +2404,31 @@ def most_unstable_parcel(pressure, temperature, dewpoint, height=None, bottom=No `pint.Quantity` Pressure, temperature, and dewpoint of most unstable parcel in the profile + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, most_unstable_parcel + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # find most unstable parcel of depth 50 hPa + >>> most_unstable_parcel(p, T, Td, depth=50*units.hPa) + (, , + , 0) + integer Index of the most unstable parcel in the given profile @@ -2058,14 +2466,26 @@ def isentropic_interpolation(levels, pressure, temperature, *args, vertical_dim= Parameters ---------- - levels : array + levels : array-like One-dimensional array of desired potential temperature surfaces - pressure : array + pressure : array-like One-dimensional array of pressure levels - temperature : array + temperature : array-like Array of temperature + + args : array-like, optional + Any additional variables will be interpolated to each isentropic level. + + Returns + ------- + list + List with pressure at each isentropic level, followed by each additional + argument interpolated to isentropic coordinates. + + Other Parameters + ---------------- vertical_dim : int, optional The axis corresponding to the vertical in the temperature array, defaults to 0. @@ -2080,17 +2500,8 @@ def isentropic_interpolation(levels, pressure, temperature, *args, vertical_dim= The desired absolute error in the calculated value, defaults to 1e-6. bottom_up_search : bool, optional - Controls whether to search for levels bottom-up, or top-down. Defaults to - True, which is bottom-up search. - - args : array, optional - Any additional variables will be interpolated to each isentropic level - - Returns - ------- - list - List with pressure at each isentropic level, followed by each additional - argument interpolated to isentropic coordinates. + Controls whether to search for levels bottom-up (starting at lower indices), + or top-down (starting at higher indices). Defaults to True, which is bottom-up search. See Also -------- @@ -2122,14 +2533,11 @@ def _isen_iter(iter_log_p, isentlevs_nd, ka, a, b, pok): fp = exner * (ka * t - a) return iter_log_p - (f / fp) - # Get dimensions in temperature - ndim = temperature.ndim - # Convert units pressure = pressure.to('hPa') temperature = temperature.to('kelvin') - slices = [np.newaxis] * ndim + slices = [np.newaxis] * temperature.ndim slices[vertical_dim] = slice(None) slices = tuple(slices) pressure = units.Quantity(np.broadcast_to(pressure[slices].magnitude, temperature.shape), @@ -2139,7 +2547,7 @@ def _isen_iter(iter_log_p, isentlevs_nd, ka, a, b, pok): sort_pressure = np.argsort(pressure.m, axis=vertical_dim) sort_pressure = np.swapaxes(np.swapaxes(sort_pressure, 0, vertical_dim)[::-1], 0, vertical_dim) - sorter = broadcast_indices(pressure, sort_pressure, ndim, vertical_dim) + sorter = broadcast_indices(sort_pressure, temperature.shape, vertical_dim) levs = pressure[sorter] tmpk = temperature[sorter] @@ -2235,8 +2643,8 @@ def isentropic_interpolation_as_dataset( eps : float, optional The desired absolute error in the calculated value, defaults to 1e-6. bottom_up_search : bool, optional - Controls whether to search for levels bottom-up, or top-down. Defaults to - True, which is bottom-up search. + Controls whether to search for levels bottom-up (starting at lower indices), + or top-down (starting at higher indices). Defaults to True, which is bottom-up search. Returns ------- @@ -2404,6 +2812,30 @@ def most_unstable_cape_cin(pressure, temperature, dewpoint, **kwargs): `pint.Quantity` Most unstable Convective Inhibition (CIN) + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, most_unstable_cape_cin + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # calculate most unstbale CAPE/CIN + >>> most_unstable_cape_cin(p, T, Td) + (, ) + See Also -------- cape_cin, most_unstable_parcel, parcel_profile @@ -2457,6 +2889,29 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs): `pint.Quantity` Mixed-layer Convective INhibition (CIN) + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, mixed_layer_cape_cin + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 25.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .75, .56, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> mixed_layer_cape_cin(p, T, Td, depth=50 * units.hPa) + (, ) + See Also -------- cape_cin, mixed_parcel, parcel_profile @@ -2532,6 +2987,31 @@ def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None, `pint.Quantity` Dewpoint of the mixed parcel + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, mixed_parcel + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # find the mixed parcel of depth 50 hPa + >>> mixed_parcel(p, T, Td, depth=50 * units.hPa) + (, , + ) + Notes ----- Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). @@ -2609,6 +3089,30 @@ def mixed_layer(pressure, *args, height=None, bottom=None, depth=None, interpola `pint.Quantity` The mixed value of the data variable + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, mixed_layer + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # find mixed layer T and Td of depth 50 hPa + >>> mixed_layer(p, T, Td, depth=50 * units.hPa) + [, ] + Notes ----- Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). @@ -2656,6 +3160,13 @@ def dry_static_energy(height, temperature): `pint.Quantity` Dry static energy + Examples + -------- + >>> from metpy.calc import dry_static_energy + >>> from metpy.units import units + >>> dry_static_energy(1000 * units.meters, 8 * units.degC) + + See Also -------- montgomery_streamfunction @@ -2702,6 +3213,13 @@ def moist_static_energy(height, temperature, specific_humidity): `pint.Quantity` Moist static energy + Examples + -------- + >>> from metpy.calc import moist_static_energy + >>> from metpy.units import units + >>> moist_static_energy(1000 * units.meters, 8 * units.degC, 8 * units('g/kg')) + + Notes ----- .. math:: \text{moist static energy} = c_{pd} T + gz + L_v q @@ -2774,6 +3292,44 @@ def thickness_hydrostatic(pressure, temperature, mixing_ratio=None, `pint.Quantity` The thickness of the layer in meters + Examples + -------- + >>> import metpy.calc as mpcalc + >>> from metpy.units import units + >>> temperature = [278, 275, 270] * units.kelvin + >>> pressure = [950, 925, 900] * units.millibar + >>> mpcalc.thickness_hydrostatic(pressure, temperature) + + + >>> bottom, depth = 950 * units.millibar, 25 * units.millibar + >>> mpcalc.thickness_hydrostatic(pressure, temperature, bottom=bottom, depth=depth) + + + To include the mixing ratio in the calculation: + + >>> r = [0.005, 0.006, 0.002] * units.dimensionless + >>> mpcalc.thickness_hydrostatic(pressure, temperature, mixing_ratio=r, + ... bottom=bottom, depth=depth) + + + Compute the 1000-500 hPa Thickness + + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # specify a layer + >>> layer = (p <= 1000 * units.hPa) & (p >= 500 * units.hPa) + >>> # compute the hydrostatic thickness + >>> mpcalc.thickness_hydrostatic(p[layer], T[layer]) + + See Also -------- thickness_hydrostatic_from_relative_humidity, pressure_to_height_std, virtual_temperature @@ -2878,6 +3434,31 @@ def thickness_hydrostatic_from_relative_humidity(pressure, temperature, relative `pint.Quantity` The thickness of the layer in meters + Examples + -------- + >>> from metpy.calc import thickness_hydrostatic_from_relative_humidity + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> ip1000_500 = (p <= 1000 * units.hPa) & (p >= 500 * units.hPa) + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # compute hydrostatic thickness from RH + >>> thickness_hydrostatic_from_relative_humidity(p[ip1000_500], + ... T[ip1000_500], + ... rh[ip1000_500]) + + See Also -------- thickness_hydrostatic, pressure_to_height_std, virtual_temperature, @@ -3081,6 +3662,13 @@ def wet_bulb_temperature(pressure, temperature, dewpoint): `pint.Quantity` Wet-bulb temperature + Examples + -------- + >>> from metpy.calc import wet_bulb_temperature + >>> from metpy.units import units + >>> wet_bulb_temperature(993 * units.hPa, 32 * units.degC, 15 * units.degC) + + See Also -------- lcl, moist_lapse @@ -3143,6 +3731,30 @@ def static_stability(pressure, temperature, vertical_dim=0): `pint.Quantity` The profile of static stability + Examples + -------- + >>> from metpy.calc import static_stability + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # Static Stability Parameter + >>> static_stability(p, T).to('m^2 s^-2 Pa^-2') + .. versionchanged:: 1.0 Renamed ``axis`` parameter ``vertical_dim`` @@ -3182,6 +3794,12 @@ def dewpoint_from_specific_humidity(pressure, temperature, specific_humidity): `pint.Quantity` Dew point temperature + Examples + -------- + >>> from metpy.calc import dewpoint_from_specific_humidity + >>> from metpy.units import units + >>> dewpoint_from_specific_humidity(1000 * units.hPa, 10 * units.degC, 5 * units('g/kg')) + .. versionchanged:: 1.0 Changed signature from ``(specific_humidity, temperature, pressure)`` @@ -3233,6 +3851,13 @@ def vertical_velocity_pressure(w, pressure, temperature, mixing_ratio=0): `pint.Quantity` Vertical velocity in terms of pressure (in Pascals / second) + Examples + -------- + >>> from metpy.calc import vertical_velocity_pressure + >>> from metpy.units import units + >>> vertical_velocity_pressure(0.5 * units('cm/s'), 700 * units.hPa, 5 * units.degC) + + See Also -------- density, vertical_velocity @@ -3285,6 +3910,13 @@ def vertical_velocity(omega, pressure, temperature, mixing_ratio=0): `pint.Quantity` Vertical velocity in terms of height (in meters / second) + Examples + -------- + >>> from metpy.calc import vertical_velocity + >>> from metpy.units import units + >>> vertical_velocity(-15 * units('Pa/s'), 700 * units.hPa, 5 * units.degC) + + See Also -------- density, vertical_velocity_pressure @@ -3313,6 +3945,12 @@ def specific_humidity_from_dewpoint(pressure, dewpoint): `pint.Quantity` Specific humidity + Examples + -------- + >>> from metpy.calc import specific_humidity_from_dewpoint + >>> from metpy.units import units + >>> specific_humidity_from_dewpoint(988 * units.hPa, 15 * units.degC).to('g/kg') + .. versionchanged:: 1.0 Changed signature from ``(dewpoint, pressure)`` @@ -3327,18 +3965,20 @@ def specific_humidity_from_dewpoint(pressure, dewpoint): @exporter.export -@preprocess_and_wrap() +@add_vertical_dim_from_xarray +@preprocess_and_wrap(broadcast=('pressure', 'temperature', 'parcel_profile')) @check_units('[pressure]', '[temperature]', '[temperature]') -def lifted_index(pressure, temperature, parcel_profile): +def lifted_index(pressure, temperature, parcel_profile, vertical_dim=0): """Calculate Lifted Index from the pressure temperature and parcel profile. Lifted index formula derived from [Galway1956]_ and referenced by [DoswellSchultz2006]_: - LI = T500 - Tp500 + + .. math:: LI = T500 - Tp500 where: - T500 is the measured temperature at 500 hPa - Tp500 is the temperature of the lifted parcel at 500 hPa + * :math:`T500` is the measured temperature at 500 hPa + * :math:`Tp500` is the temperature of the lifted parcel at 500 hPa Calculation of the lifted index is defined as the temperature difference between the observed 500 hPa temperature and the temperature of a parcel lifted from the @@ -3356,36 +3996,68 @@ def lifted_index(pressure, temperature, parcel_profile): parcel_profile : `pint.Quantity` Temperature profile of the parcel + vertical_dim : int, optional + The axis corresponding to vertical, defaults to 0. Automatically determined from + xarray DataArray arguments. + Returns ------- `pint.Quantity` Lifted Index + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, lifted_index, parcel_profile + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compute the parcel temperatures from surface parcel + >>> prof = parcel_profile(p, T[0], Td[0]) + >>> # calculate the LI + >>> lifted_index(p, T, prof) + + """ # find the measured temperature and parcel profile temperature at 500 hPa. t500, tp500 = interpolate_1d(units.Quantity(500, 'hPa'), - pressure, temperature, parcel_profile) + pressure, temperature, parcel_profile, axis=vertical_dim) # calculate the lifted index. - return t500 - tp500.to(units.degC) + return t500 - tp500 @exporter.export -@preprocess_and_wrap() +@add_vertical_dim_from_xarray +@preprocess_and_wrap(broadcast=('pressure', 'temperature', 'dewpoint')) @check_units('[pressure]', '[temperature]', '[temperature]') -def k_index(pressure, temperature, dewpoint): +def k_index(pressure, temperature, dewpoint, vertical_dim=0): """Calculate K Index from the pressure temperature and dewpoint. K Index formula derived from [George1960]_: - K = (T850 - T500) + Td850 - (T700 - Td700) + + .. math:: K = (T850 - T500) + Td850 - (T700 - Td700) where: - T850 is the temperature at 850 hPa - T700 is the temperature at 700 hPa - T500 is the temperature at 500 hPa - Td850 is the dewpoint at 850 hPa - Td700 is the dewpoint at 700 hPa + * :math:`T850` is the temperature at 850 hPa + * :math:`T700` is the temperature at 700 hPa + * :math:`T500` is the temperature at 500 hPa + * :math:`Td850` is the dewpoint at 850 hPa + * :math:`Td700` is the dewpoint at 700 hPa Calculation of the K Index is defined as the temperature difference between the static instability between 850 hPa and 500 hPa, add with the moisture @@ -3402,16 +4074,43 @@ def k_index(pressure, temperature, dewpoint): dewpoint : `pint.Quantity` Dewpoint temperature corresponding to pressure + vertical_dim : int, optional + The axis corresponding to vertical, defaults to 0. Automatically determined from + xarray DataArray arguments. + Returns ------- `pint.Quantity` K Index + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, k_index + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> k_index(p, T, Td) + + """ - # Find temperature and dewpoint at 850, 700 and 500 hPa - (t850, t700, t500), (td850, td700, _) = interpolate_1d(units.Quantity([850, 700, 500], - 'hPa'), pressure, temperature, - dewpoint) + # Find temperature and dewpoint at 850, 700, and 500 hPa + (t850, t700, t500), (td850, td700, _) = interpolate_1d( + units.Quantity([850, 700, 500], 'hPa'), pressure, temperature, dewpoint, + axis=vertical_dim) # Calculate k index. return ((t850 - t500) + td850 - (t700 - td700)).to(units.degC) @@ -3475,7 +4174,7 @@ def gradient_richardson_number(height, potential_temperature, u, v, vertical_dim def scale_height(temperature_bottom, temperature_top): r"""Calculate the scale height of a layer. - .. math:: H = \frac{R_d \overline{T}}{g} + .. math:: H = \frac{R_d \overline{T}}{g} This function assumes dry air, but can be used with the virtual temperature to account for moisture. @@ -3492,6 +4191,14 @@ def scale_height(temperature_bottom, temperature_top): ------- `pint.Quantity` Scale height of layer + + Examples + -------- + >>> from metpy.calc import scale_height + >>> from metpy.units import units + >>> scale_height(20 * units.degC, -50 * units.degC) + + """ t_bar = 0.5 * (temperature_bottom + temperature_top) return (mpconsts.nounit.Rd * t_bar) / mpconsts.nounit.g @@ -3504,28 +4211,54 @@ def showalter_index(pressure, temperature, dewpoint): """Calculate Showalter Index. Showalter Index derived from [Galway1956]_: - SI = T500 - Tp500 + + .. math:: SI = T500 - Tp500 where: - T500 is the measured temperature at 500 hPa - Tp500 is the temperature of the parcel at 500 hPa when lifted from 850 hPa + + * :math:`T500` is the measured temperature at 500 hPa + * :math:`Tp500` is the temperature of the parcel at 500 hPa when lifted from 850 hPa Parameters ---------- - pressure : `pint.Quantity` - Atmospheric pressure, in order from highest to lowest pressure + pressure : `pint.Quantity` + Atmospheric pressure, in order from highest to lowest pressure - temperature : `pint.Quantity` - Ambient temperature corresponding to ``pressure`` + temperature : `pint.Quantity` + Ambient temperature corresponding to ``pressure`` - dewpoint : `pint.Quantity` - Ambient dew point temperatures corresponding to ``pressure`` + dewpoint : `pint.Quantity` + Ambient dew point temperatures corresponding to ``pressure`` Returns ------- `pint.Quantity` Showalter index + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, showalter_index + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compute the showalter index + >>> showalter_index(p, T, Td) + + """ # find the measured temperature and dew point temperature at 850 hPa. t850, td850 = interpolate_1d(units.Quantity(850, 'hPa'), pressure, temperature, dewpoint) @@ -3541,19 +4274,21 @@ def showalter_index(pressure, temperature, dewpoint): @exporter.export -@preprocess_and_wrap() +@add_vertical_dim_from_xarray +@preprocess_and_wrap(broadcast=('pressure', 'temperature', 'dewpoint')) @check_units('[pressure]', '[temperature]', '[temperature]') -def total_totals_index(pressure, temperature, dewpoint): +def total_totals_index(pressure, temperature, dewpoint, vertical_dim=0): """Calculate Total Totals Index from the pressure temperature and dewpoint. Total Totals Index formula derived from [Miller1972]_: - TT = (T850 + Td850) - (2 * T500) + + .. math:: TT = (T850 + Td850) - (2 * T500) where: - T850 is the temperature at 850 hPa - T500 is the temperature at 500 hPa - Td850 is the dewpoint at 850 hPa + * :math:`T850` is the temperature at 850 hPa + * :math:`T500` is the temperature at 500 hPa + * :math:`Td850` is the dewpoint at 850 hPa Calculation of the Total Totals Index is defined as the temperature at 850 hPa plus the dewpoint at 850 hPa, minus twice the temperature at 500 hPa. This index consists of @@ -3570,33 +4305,64 @@ def total_totals_index(pressure, temperature, dewpoint): dewpoint : `pint.Quantity` Dewpoint temperature corresponding to pressure + vertical_dim : int, optional + The axis corresponding to vertical, defaults to 0. Automatically determined from + xarray DataArray arguments. + Returns ------- `pint.Quantity` Total Totals Index + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, total_totals_index + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compute the TT index + >>> total_totals_index(p, T, Td) + + """ # Find temperature and dewpoint at 850 and 500 hPa. (t850, t500), (td850, _) = interpolate_1d(units.Quantity([850, 500], 'hPa'), - pressure, temperature, dewpoint) + pressure, temperature, dewpoint, + axis=vertical_dim) # Calculate total totals index. return (t850 - t500) + (td850 - t500) @exporter.export -@preprocess_and_wrap() +@add_vertical_dim_from_xarray +@preprocess_and_wrap(broadcast=('pressure', 'temperature')) @check_units('[pressure]', '[temperature]') -def vertical_totals(pressure, temperature): +def vertical_totals(pressure, temperature, vertical_dim=0): """Calculate Vertical Totals from the pressure and temperature. Vertical Totals formula derived from [Miller1972]_: - VT = T850 - T500 + + .. math:: VT = T850 - T500 where: - T850 is the temperature at 850 hPa - T500 is the temperature at 500 hPa + * :math:`T850` is the temperature at 850 hPa + * :math:`T500` is the temperature at 500 hPa Calculation of the Vertical Totals is defined as the temperature difference between 850 hPa and 500 hPa. This is a part of the Total Totals Index. @@ -3609,33 +4375,57 @@ def vertical_totals(pressure, temperature): temperature : `pint.Quantity` Temperature corresponding to pressure + vertical_dim : int, optional + The axis corresponding to vertical, defaults to 0. Automatically determined from + xarray DataArray arguments. + Returns ------- `pint.Quantity` Vertical Totals + Examples + -------- + >>> from metpy.calc import vertical_totals + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # compute vertical totals index + >>> vertical_totals(p, T) + + """ # Find temperature at 850 and 500 hPa. (t850, t500) = interpolate_1d(units.Quantity([850, 500], 'hPa'), - pressure, temperature) + pressure, temperature, axis=vertical_dim) # Calculate vertical totals. return t850 - t500 @exporter.export -@preprocess_and_wrap() +@add_vertical_dim_from_xarray +@preprocess_and_wrap(broadcast=('pressure', 'temperature', 'dewpoint')) @check_units('[pressure]', '[temperature]', '[temperature]') -def cross_totals(pressure, temperature, dewpoint): +def cross_totals(pressure, temperature, dewpoint, vertical_dim=0): """Calculate Cross Totals from the pressure temperature and dewpoint. Cross Totals formula derived from [Miller1972]_: - CT = Td850 - T500 + + .. math:: CT = Td850 - T500 where: - Td850 is the dewpoint at 850 hPa - T500 is the temperature at 500 hPa + * :math:`Td850` is the dewpoint at 850 hPa + * :math:`T500` is the temperature at 500 hPa Calculation of the Cross Totals is defined as the difference between dewpoint at 850 hPa and temperature at 500 hPa. This is a part of the Total Totals Index. @@ -3651,15 +4441,132 @@ def cross_totals(pressure, temperature, dewpoint): dewpoint : `pint.Quantity` Dewpoint temperature corresponding to pressure + vertical_dim : int, optional + The axis corresponding to vertical, defaults to 0. Automatically determined from + xarray DataArray arguments. + Returns ------- `pint.Quantity` Cross Totals + Examples + -------- + >>> from metpy.calc import dewpoint_from_relative_humidity, cross_totals + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate dewpoint + >>> Td = dewpoint_from_relative_humidity(T, rh) + >>> # compute the cross totals index + >>> cross_totals(p, T, Td) + + """ # Find temperature and dewpoint at 850 and 500 hPa (_, t500), (td850, _) = interpolate_1d(units.Quantity([850, 500], 'hPa'), - pressure, temperature, dewpoint) + pressure, temperature, dewpoint, axis=vertical_dim) # Calculate vertical totals. return td850 - t500 + + +@exporter.export +@add_vertical_dim_from_xarray +@preprocess_and_wrap(broadcast=('pressure', 'temperature', 'dewpoint', 'speed', 'direction')) +@check_units('[pressure]', '[temperature]', '[temperature]', '[speed]') +def sweat_index(pressure, temperature, dewpoint, speed, direction, vertical_dim=0): + """Calculate SWEAT Index. + + SWEAT Index derived from [Miller1972]_: + + .. math:: SWEAT = 12Td_{850} + 20(TT - 49) + 2f_{850} + f_{500} + 125(S + 0.2) + + where: + + * :math:`Td_{850}` is the dewpoint at 850 hPa; the first term is set to zero + if :math:`Td_{850}` is negative. + * :math:`TT` is the total totals index; the second term is set to zero + if :math:`TT` is less than 49 + * :math:`f_{850}` is the wind speed at 850 hPa + * :math:`f_{500}` is the wind speed at 500 hPa + * :math:`S` is the shear term: :math:`sin{(dd_{850} - dd_{500})}`, where + :math:`dd_{850}` and :math:`dd_{500}` are the wind directions at 850 hPa and 500 hPa, + respectively. It is set to zero if any of the following conditions are not met: + + 1. :math:`dd_{850}` is between 130 - 250 degrees + 2. :math:`dd_{500}` is between 210 - 310 degrees + 3. :math:`dd_{500} - dd_{850} > 0` + 4. both the wind speeds are greater than or equal to 15 kts + + Calculation of the SWEAT Index consists of a low-level moisture, instability, + and the vertical wind shear (both speed and direction). This index aim to + determine the likeliness of severe weather and tornadoes. + + Parameters + ---------- + pressure : `pint.Quantity` + Pressure level(s), in order from highest to lowest pressure + + temperature : `pint.Quantity` + Temperature corresponding to pressure + + dewpoint : `pint.Quantity` + Dewpoint temperature corresponding to pressure + + speed : `pint.Quantity` + Wind speed corresponding to pressure + + direction : `pint.Quantity` + Wind direction corresponding to pressure + + vertical_dim : int, optional + The axis corresponding to vertical, defaults to 0. Automatically determined from + xarray DataArray arguments. + + Returns + ------- + `pint.Quantity` + SWEAT Index + + """ + # Find dewpoint at 850 hPa. + td850 = interpolate_1d(units.Quantity(850, 'hPa'), pressure, dewpoint, axis=vertical_dim) + + # Find total totals index. + tt = total_totals_index(pressure, temperature, dewpoint, vertical_dim=vertical_dim) + + # Find wind speed and direction at 850 and 500 hPa + (f850, f500), (dd850, dd500) = interpolate_1d(units.Quantity([850, 500], 'hPa'), + pressure, speed, direction, + axis=vertical_dim) + + # First term is set to zero if Td850 is negative + first_term = 12 * np.clip(td850.m_as('degC'), 0, None) + + # Second term is set to zero if TT is less than 49 + second_term = 20 * np.clip(tt.m_as('degC') - 49, 0, None) + + # Shear term is set to zero if any of four conditions are not met + required = ((units.Quantity(130, 'deg') <= dd850) & (dd850 <= units.Quantity(250, 'deg')) + & (units.Quantity(210, 'deg') <= dd500) & (dd500 <= units.Quantity(310, 'deg')) + & (dd500 - dd850 > 0) + & (f850 >= units.Quantity(15, 'knots')) + & (f500 >= units.Quantity(15, 'knots'))) + shear_term = np.atleast_1d(125 * (np.sin(dd500 - dd850) + 0.2)) + shear_term[~required] = 0 + + # Calculate sweat index. + return first_term + second_term + (2 * f850.m) + f500.m + shear_term \ No newline at end of file From b6aa39ab2f6faf07a5bcde552f0c7b22de67b2de Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 13:03:24 -0500 Subject: [PATCH 09/30] add missing sweat index calc --- tests/calc/test_thermo.py | 46 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 38f054f8be3..6d2d3d167f1 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -2349,3 +2349,49 @@ def test_parcel_profile_with_lcl_as_dataset_duplicates(): profile = parcel_profile_with_lcl_as_dataset(pressure, temperature, dewpoint) xr.testing.assert_allclose(profile, truth, atol=1e-5) + + +def test_sweat_index(): + """Test the SWEAT Index calculation.""" + pressure = np.array([1008., 1000., 947., 925., 921., 896., 891., 889., 866., + 858., 850., 835., 820., 803., 733., 730., 700., 645., + 579., 500., 494., 466., 455., 441., 433., 410., 409., + 402., 400., 390., 388., 384., 381., 349., 330., 320., + 306., 300., 278., 273., 250., 243., 208., 200., 196., + 190., 179., 159., 151., 150., 139.]) * units.hPa + temperature = np.array([27.4, 26.4, 22.9, 21.4, 21.2, 20.7, 20.6, 21.2, 19.4, + 19.1, 18.8, 17.8, 17.4, 16.3, 11.4, 11.2, 10.2, 6.1, + 0.6, -4.9, -5.5, -8.5, -9.9, -11.7, -12.3, -13.7, -13.8, + -14.9, -14.9, -16.1, -16.1, -16.9, -17.3, -21.7, -24.5, -26.1, + -28.3, -29.5, -33.1, -34.2, -39.3, -41., -50.2, -52.5, -53.5, + -55.2, -58.6, -65.2, -68.1, -68.5, -72.5]) * units.degC + dewpoint = np.array([24.9, 24.6, 22., 20.9, 20.7, 14.8, 13.6, 12.2, 16.8, + 16.6, 16.5, 15.9, 13.6, 13.2, 11.3, 11.2, 8.6, 4.5, + -0.8, -8.1, -9.5, -12.7, -12.7, -12.8, -13.1, -24.7, -24.4, + -21.9, -24.9, -36.1, -31.1, -26.9, -27.4, -33., -36.5, -47.1, + -31.4, -33.5, -40.1, -40.8, -44.1, -45.6, -54., -56.1, -56.9, + -58.6, -61.9, -68.4, -71.2, -71.6, -77.2]) * units.degC + speed = np.array([0., 3., 10., 12., 12., 14., 14., 14., 12., + 12., 12., 12., 11., 11., 12., 12., 10., 10., + 8., 5., 4., 1., 0., 3., 5., 10., 10., + 11., 11., 13., 14., 14., 15., 23., 23., 24., + 24., 24., 26., 27., 28., 30., 25., 24., 26., + 28., 33., 29., 32., 26., 26.]) * units.knot + direction = np.array([0., 170., 200., 205., 204., 200., 197., 195., 180., + 175., 175., 178., 181., 185., 160., 160., 165., 165., + 203., 255., 268., 333., 0., 25., 40., 83., 85., + 89., 90., 100., 103., 107., 110., 90., 88., 87., + 86., 85., 85., 85., 60., 55., 60., 50., 46., + 40., 45., 35., 50., 50., 50.]) * units.degree + + sweat = sweat_index(pressure, temperature, dewpoint, speed, direction) + assert_almost_equal(sweat, 227., 2) + + +def test_sweat_index_xarray(index_xarray_data): + """Test the SWEAT index calculation with a grid of xarray data.""" + result = sweat_index(index_xarray_data.isobaric, index_xarray_data.temperature, + index_xarray_data.dewpoint, index_xarray_data.wind_speed, + index_xarray_data.wind_direction) + assert_array_almost_equal(result, np.array([[[[490.2, 478.2, 466.2], + [454.2, 442.2, 430.2]]]])) From 0d9f1fd3ebf5068ce5a2fc52376c6b099af5b634 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 14:04:46 -0500 Subject: [PATCH 10/30] add new virtual temperature from dewpoint function --- src/metpy/calc/thermo.py | 61 +++++++++++++++++++++++++++++++++++++++ tests/calc/test_thermo.py | 12 +++++++- 2 files changed, 72 insertions(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 7d241d926df..68ba0ed85f7 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1688,6 +1688,67 @@ def virtual_temperature( / (molecular_weight_ratio * (1 + mixing_ratio))) +@exporter.export +@preprocess_and_wrap(wrap_like='temperature', broadcast=('pressure', + 'temperature', + 'dewpoint')) +@check_units('[pressure]', '[temperature]', '[temperature]') +def virtual_temperature_from_dewpoint( + pressure, temperature, dewpoint, molecular_weight_ratio=mpconsts.nounit.epsilon +): + r"""Calculate virtual temperature. + + This calculation must be given an air parcel's temperature and mixing ratio. + The implementation uses the formula outlined in [Hobbs2006]_ pg.80. + + Parameters + ---------- + pressure: `pint.Quantity` + Total atmospheric pressure + + temperature: `pint.Quantity` + Air temperature + + dewpoint : `pint.Quantity` + Dewpoint temperature + + molecular_weight_ratio : `pint.Quantity` or float, optional + The ratio of the molecular weight of the constituent gas to that assumed + for air. Defaults to the ratio for water vapor to dry air. + (:math:`\epsilon\approx0.622`) + + Returns + ------- + `pint.Quantity` + Corresponding virtual temperature of the parcel + + Examples + -------- + >>> from metpy.calc import virtual_temperature_from_dewpoint + >>> from metpy.units import units + >>> virtual_temperature_from_dewpoint(1000 * units.hPa, + 30 * units.degC, + 25 * units.degC) + + + Notes + ----- + .. math:: T_v = T \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})} + + .. versionchanged:: 1.0 + Renamed ``mixing`` parameter to ``mixing_ratio`` + + """ + # Convert the dewpoint to specific humidity + specific_humidity = specific_humidity_from_dewpoint(pressure, dewpoint) + + # Convert the specific humidity to mixing ratio + mixing_ratio = mixing_ratio_from_specific_humidity(specific_humidity) + + # Calculate virtual temperature with given parameters + return virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio) + + @exporter.export @preprocess_and_wrap( wrap_like='temperature', diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 6d2d3d167f1..4e821247c7a 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -32,7 +32,8 @@ thickness_hydrostatic_from_relative_humidity, total_totals_index, vapor_pressure, vertical_totals, vertical_velocity, vertical_velocity_pressure, virtual_potential_temperature, - virtual_temperature, wet_bulb_temperature) + virtual_temperature, virtual_temperature_from_dewpoint, + wet_bulb_temperature) from metpy.calc.thermo import _find_append_zero_crossings from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_nan from metpy.units import is_quantity, masked_array, units @@ -792,6 +793,15 @@ def test_virtual_temperature(): assert_almost_equal(tv, 288.2796 * units.kelvin, 3) +def test_virtual_temperature_from_dewpoint(): + """Test virtual temperature calculation.""" + p = 1000 * units.hPa + t = 30 * units.degC + td = 25 * units.degC + tv = virtual_temperature_from_dewpoint(p, t, td) + assert_almost_equal(tv, 33.6740 * units.degC, 3) + + def test_virtual_potential_temperature(): """Test virtual potential temperature calculation.""" p = 999. * units.mbar From d393916471d6d70e6271a744b40b20059a9d18ae Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 14:52:38 -0500 Subject: [PATCH 11/30] add new function to overrides --- docs/_templates/overrides/metpy.calc.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index 8efa5adee2a..c5920a81505 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -59,6 +59,7 @@ Moist Thermodynamics vertical_velocity_pressure virtual_potential_temperature virtual_temperature + virtual_temperature_from_dewpoint wet_bulb_temperature From 22aeec62a54accd7d5b38818a3f470ce3a6d37dd Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 15:28:11 -0500 Subject: [PATCH 12/30] use virtual temperature for CAPE --- src/metpy/calc/thermo.py | 9 ++-- tests/calc/test_thermo.py | 88 +++++++++++++++++++-------------------- 2 files changed, 50 insertions(+), 47 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 68ba0ed85f7..6c94a8d1233 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1726,9 +1726,7 @@ def virtual_temperature_from_dewpoint( -------- >>> from metpy.calc import virtual_temperature_from_dewpoint >>> from metpy.units import units - >>> virtual_temperature_from_dewpoint(1000 * units.hPa, - 30 * units.degC, - 25 * units.degC) + >>> virtual_temperature_from_dewpoint(1000 * units.hPa, 30 * units.degC, 25 * units.degC) Notes @@ -2339,6 +2337,11 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' """ pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, dewpoint, parcel_profile) + + # Convert the temperature/parcel profile to virtual temperature + temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint) + parcel_profile = virtual_temperature_from_dewpoint(pressure, parcel_profile, dewpoint) + # Calculate LFC limit of integration lfc_pressure, _ = lfc(pressure, temperature, dewpoint, parcel_temperature_profile=parcel_profile, which=which_lfc) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 4e821247c7a..01a8810ac7d 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -1066,8 +1066,8 @@ def test_cape_cin(): dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 75.05354 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -89.890078 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 75.229846 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -90.004823 * units('joule / kilogram'), 2) def test_cape_cin_no_el(): @@ -1077,8 +1077,8 @@ def test_cape_cin_no_el(): dewpoint = np.array([19., -11.2, -10.8, -10.4]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC') cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 0.08610409 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -89.8900784 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 0.08381723 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -90.0048229 * units('joule / kilogram'), 2) def test_cape_cin_no_lfc(): @@ -1394,8 +1394,8 @@ def test_surface_based_cape_cin(array_class): temperature = array_class([22.2, 14.6, 12., 9.4, 7., -38.], units.celsius) dewpoint = array_class([19., -11.2, -10.8, -10.4, -10., -53.2], units.celsius) cape, cin = surface_based_cape_cin(p, temperature, dewpoint) - assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 75.2298456 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -137.133362 * units('joule / kilogram'), 2) def test_surface_based_cape_cin_with_xarray(): @@ -1418,8 +1418,8 @@ def test_surface_based_cape_cin_with_xarray(): data['temperature'], data['dewpoint'] ) - assert_almost_equal(cape, 75.0535446 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -136.685967 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 75.2298456 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -137.133362 * units('joule / kilogram'), 2) def test_profile_with_nans(): @@ -1459,8 +1459,8 @@ def test_most_unstable_cape_cin_surface(): temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 75.0535446 * units('joule / kilogram'), 2) - assert_almost_equal(mucin, -136.685967 * units('joule / kilogram'), 2) + assert_almost_equal(mucape, 75.2298456 * units('joule / kilogram'), 2) + assert_almost_equal(mucin, -137.133362 * units('joule / kilogram'), 2) def test_most_unstable_cape_cin(): @@ -1469,8 +1469,8 @@ def test_most_unstable_cape_cin(): temperature = np.array([18.2, 22.2, 17.4, 10., 0., 15]) * units.celsius dewpoint = np.array([19., 19., 14.3, 0., -10., 0.]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 157.11404 * units('joule / kilogram'), 4) - assert_almost_equal(mucin, -31.8406578 * units('joule / kilogram'), 4) + assert_almost_equal(mucape, 157.354834 * units('joule / kilogram'), 4) + assert_almost_equal(mucin, -32.0881909 * units('joule / kilogram'), 4) def test_mixed_parcel(): @@ -1490,8 +1490,8 @@ def test_mixed_layer_cape_cin(multiple_intersections): """Test the calculation of mixed layer cape/cin.""" pressure, temperature, dewpoint = multiple_intersections mlcape, mlcin = mixed_layer_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mlcape, 987.7323 * units('joule / kilogram'), 2) - assert_almost_equal(mlcin, -20.6727628 * units('joule / kilogram'), 2) + assert_almost_equal(mlcape, 987.7085876 * units('joule / kilogram'), 2) + assert_almost_equal(mlcin, -20.80235875 * units('joule / kilogram'), 2) def test_mixed_layer(): @@ -1806,10 +1806,10 @@ def test_multiple_lfcs_simple(multiple_intersections): 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.3534595 * units.mbar, 3) - assert_almost_equal(lfc_temp_top, 4.8848999 * units.degC, 3) - assert_almost_equal(lfc_pressure_bottom, 884.14790 * units.mbar, 3) - assert_almost_equal(lfc_temp_bottom, 13.95707016 * units.degC, 3) + assert_almost_equal(lfc_pressure_top, 705.3534497 * units.mbar, 3) + assert_almost_equal(lfc_temp_top, 4.884899 * units.degC, 3) + assert_almost_equal(lfc_pressure_bottom, 884.1478935 * units.mbar, 3) + assert_almost_equal(lfc_temp_bottom, 13.95706975 * units.degC, 3) assert_almost_equal(len(lfc_pressure_all), 2, 0) @@ -1817,8 +1817,8 @@ def test_multiple_lfs_wide(multiple_intersections): """Test 'wide' LFC for sounding with multiple LFCs.""" levels, temperatures, dewpoints = multiple_intersections lfc_pressure_wide, lfc_temp_wide = lfc(levels, temperatures, dewpoints, which='wide') - assert_almost_equal(lfc_pressure_wide, 705.3534595 * units.hPa, 3) - assert_almost_equal(lfc_temp_wide, 4.8848999 * units.degC, 3) + assert_almost_equal(lfc_pressure_wide, 705.3534497 * units.hPa, 3) + assert_almost_equal(lfc_temp_wide, 4.88489902 * units.degC, 3) def test_invalid_which(multiple_intersections): @@ -1842,10 +1842,10 @@ def test_multiple_els_simple(multiple_intersections): 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.151466 * units.mbar, 3) - assert_almost_equal(el_temp_top, -56.81015490 * units.degC, 3) - assert_almost_equal(el_pressure_bottom, 849.7998947 * units.mbar, 3) - assert_almost_equal(el_temp_bottom, 12.4233265 * units.degC, 3) + assert_almost_equal(el_pressure_top, 228.151524 * units.mbar, 3) + assert_almost_equal(el_temp_top, 228.151524 * units.degC, 3) + assert_almost_equal(el_pressure_bottom, 849.7998957 * units.mbar, 3) + assert_almost_equal(el_temp_bottom, 12.42268288 * units.degC, 3) assert_almost_equal(len(el_pressure_all), 2, 0) @@ -1853,24 +1853,24 @@ def test_multiple_el_wide(multiple_intersections): """Test 'wide' EL for sounding with multiple ELs.""" levels, temperatures, dewpoints = multiple_intersections el_pressure_wide, el_temp_wide = el(levels, temperatures, dewpoints, which='wide') - assert_almost_equal(el_pressure_wide, 228.151466 * units.hPa, 3) - assert_almost_equal(el_temp_wide, -56.81015490 * units.degC, 3) + assert_almost_equal(el_pressure_wide, 228.151524 * units.hPa, 3) + assert_almost_equal(el_temp_wide, -56.81015357 * units.degC, 3) def test_muliple_el_most_cape(multiple_intersections): """Test 'most_cape' EL for sounding with multiple ELs.""" levels, temperatures, dewpoints = multiple_intersections el_pressure_wide, el_temp_wide = el(levels, temperatures, dewpoints, which='most_cape') - assert_almost_equal(el_pressure_wide, 228.151466 * units.hPa, 3) - assert_almost_equal(el_temp_wide, -56.81015490 * units.degC, 3) + assert_almost_equal(el_pressure_wide, 228.151524 * units.hPa, 3) + assert_almost_equal(el_temp_wide, -56.81015356 * units.degC, 3) def test_muliple_lfc_most_cape(multiple_intersections): """Test 'most_cape' LFC for sounding with multiple LFCs.""" levels, temperatures, dewpoints = multiple_intersections lfc_pressure_wide, lfc_temp_wide = lfc(levels, temperatures, dewpoints, which='most_cape') - assert_almost_equal(lfc_pressure_wide, 705.3534595 * units.hPa, 3) - assert_almost_equal(lfc_temp_wide, 4.8848999 * units.degC, 3) + assert_almost_equal(lfc_pressure_wide, 705.35344968 * units.hPa, 3) + assert_almost_equal(lfc_temp_wide, 4.88489902 * units.degC, 3) def test_el_lfc_most_cape_bottom(): @@ -1884,9 +1884,9 @@ def test_el_lfc_most_cape_bottom(): lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints, which='most_cape') el_pressure, el_temp = el(levels, temperatures, dewpoints, which='most_cape') assert_almost_equal(lfc_pressure, 900.73235 * units.hPa, 3) - assert_almost_equal(lfc_temp, 14.672512 * units.degC, 3) - assert_almost_equal(el_pressure, 849.7998947 * units.hPa, 3) - assert_almost_equal(el_temp, 12.4233265 * units.degC, 3) + assert_almost_equal(lfc_temp, 14.6717153 * units.degC, 3) + assert_almost_equal(el_pressure, 849.79989566 * units.hPa, 3) + assert_almost_equal(el_temp, 12.42268288 * units.degC, 3) def test_cape_cin_top_el_lfc(multiple_intersections): @@ -1894,8 +1894,8 @@ def test_cape_cin_top_el_lfc(multiple_intersections): levels, temperatures, dewpoints = multiple_intersections parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='top') - assert_almost_equal(cape, 1258.94592 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -97.752333 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 1259.157926 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -97.9763618 * units('joule / kilogram'), 3) def test_cape_cin_bottom_el_lfc(multiple_intersections): @@ -1903,8 +1903,8 @@ def test_cape_cin_bottom_el_lfc(multiple_intersections): levels, temperatures, dewpoints = multiple_intersections parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_el='bottom') - assert_almost_equal(cape, 2.18798 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -8.16221 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 2.20076511 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -8.22115922 * units('joule / kilogram'), 3) def test_cape_cin_wide_el_lfc(multiple_intersections): @@ -1913,8 +1913,8 @@ def test_cape_cin_wide_el_lfc(multiple_intersections): parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='wide', which_el='wide') - assert_almost_equal(cape, 1258.9459 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -97.752333 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 1259.157926 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -97.9763618 * units('joule / kilogram'), 3) def test_cape_cin_custom_profile(): @@ -1924,7 +1924,7 @@ def test_cape_cin_custom_profile(): dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) + 5 * units.delta_degC cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 1440.463208696 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 1442.0101452 * units('joule / kilogram'), 2) assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 2) @@ -2010,7 +2010,7 @@ def test_cape_cin_value_error(): -35.9, -26.7, -37.7, -43.1, -33.9, -40.9, -46.1, -34.9, -33.9, -33.7, -33.3, -42.5, -50.3, -49.7, -49.5, -58.3, -61.3]) * units.degC cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint) - expected_cape, expected_cin = 2007.040698 * units('joules/kg'), 0.0 * units('joules/kg') + expected_cape, expected_cin = 1973.828926 * units('joules/kg'), 0.0 * units('joules/kg') assert_almost_equal(cape, expected_cape, 3) assert_almost_equal(cin, expected_cin, 3) @@ -2083,7 +2083,7 @@ def test_lifted_index(): -57.5]) * units.degC parcel_prof = parcel_profile(pressure, temperature[0], dewpoint[0]) li = lifted_index(pressure, temperature, parcel_prof) - assert_almost_equal(li, -7.9176350 * units.delta_degree_Celsius, 2) + assert_almost_equal(li, -7.9115691 * units.delta_degree_Celsius, 2) def test_lifted_index_500hpa_missing(): @@ -2108,7 +2108,7 @@ def test_lifted_index_500hpa_missing(): -57.5]) * units.degC parcel_prof = parcel_profile(pressure, temperature[0], dewpoint[0]) li = lifted_index(pressure, temperature, parcel_prof) - assert_almost_equal(li, -7.9176350 * units.delta_degree_Celsius, 1) + assert_almost_equal(li, -7.9806301 * units.delta_degree_Celsius, 1) def test_lifted_index_xarray(index_xarray_data): @@ -2204,7 +2204,7 @@ def test_showalter_index(): -48.9, -50.2, -51.5, -53.3, -55.5, -55.9]), 'degC') result = showalter_index(pressure, temps, dewp) - assert_almost_equal(result, units.Quantity(7.6024, 'delta_degC'), 4) + assert_almost_equal(result, units.Quantity(7.6028, 'delta_degC'), 4) def test_total_totals_index(): From 1af6af1ba91c536f922c44ac6f9f687fdd061e19 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 15:33:26 -0500 Subject: [PATCH 13/30] fix last few failing tests --- tests/calc/test_thermo.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 01a8810ac7d..2160f85bbcc 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -666,8 +666,8 @@ def test_sensitive_sounding(): assert_almost_equal(lfc_temp, 20.498 * units.degC, 2) pos, neg = surface_based_cape_cin(p, t, td) - assert_almost_equal(pos, 0.1115 * units('J/kg'), 3) - assert_almost_equal(neg, -6.0866 * units('J/kg'), 3) + assert_almost_equal(pos, 0.0 * units('J/kg'), 3) + assert_almost_equal(neg, 0.0 * units('J/kg'), 3) def test_lfc_sfc_precision(): @@ -1843,7 +1843,7 @@ def test_multiple_els_simple(multiple_intersections): 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.151524 * units.mbar, 3) - assert_almost_equal(el_temp_top, 228.151524 * units.degC, 3) + assert_almost_equal(el_temp_top, -56.810153566 * units.degC, 3) assert_almost_equal(el_pressure_bottom, 849.7998957 * units.mbar, 3) assert_almost_equal(el_temp_bottom, 12.42268288 * units.degC, 3) assert_almost_equal(len(el_pressure_all), 2, 0) @@ -2204,7 +2204,7 @@ def test_showalter_index(): -48.9, -50.2, -51.5, -53.3, -55.5, -55.9]), 'degC') result = showalter_index(pressure, temps, dewp) - assert_almost_equal(result, units.Quantity(7.6028, 'delta_degC'), 4) + assert_almost_equal(result, units.Quantity(7.6024, 'delta_degC'), 4) def test_total_totals_index(): From 994b00bcd0a970c682c149811e902b3f08263702 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 15:36:49 -0500 Subject: [PATCH 14/30] fix linting issues --- src/metpy/calc/thermo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 6c94a8d1233..f05e83eca55 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2337,11 +2337,11 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' """ pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, dewpoint, parcel_profile) - + # Convert the temperature/parcel profile to virtual temperature temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint) parcel_profile = virtual_temperature_from_dewpoint(pressure, parcel_profile, dewpoint) - + # Calculate LFC limit of integration lfc_pressure, _ = lfc(pressure, temperature, dewpoint, parcel_temperature_profile=parcel_profile, which=which_lfc) From fb7b8e6db27abcccbe4bd4dbaedb9561e53d01bd Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 16:10:25 -0500 Subject: [PATCH 15/30] update docstring values --- src/metpy/calc/thermo.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index f05e83eca55..d208efe9f8b 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1727,7 +1727,7 @@ def virtual_temperature_from_dewpoint( >>> from metpy.calc import virtual_temperature_from_dewpoint >>> from metpy.units import units >>> virtual_temperature_from_dewpoint(1000 * units.hPa, 30 * units.degC, 25 * units.degC) - + Notes ----- @@ -2301,7 +2301,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') >>> # calculate surface based CAPE/CIN >>> cape_cin(p, T, Td, prof) - (, ) + (, ) See Also -------- @@ -2898,7 +2898,7 @@ def most_unstable_cape_cin(pressure, temperature, dewpoint, **kwargs): >>> Td = dewpoint_from_relative_humidity(T, rh) >>> # calculate most unstbale CAPE/CIN >>> most_unstable_cape_cin(p, T, Td) - (, ) + (, ) See Also -------- @@ -2974,7 +2974,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs): >>> # calculate dewpoint >>> Td = dewpoint_from_relative_humidity(T, rh) >>> mixed_layer_cape_cin(p, T, Td, depth=50 * units.hPa) - (, ) + (, ) See Also -------- From be33d31777d9ba92dd9fc9e5affc8d2ccfe4c150 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 3 Apr 2023 17:24:18 -0500 Subject: [PATCH 16/30] fix units for virtual temperature --- src/metpy/calc/thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index d208efe9f8b..b64a30f30a8 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1727,7 +1727,7 @@ def virtual_temperature_from_dewpoint( >>> from metpy.calc import virtual_temperature_from_dewpoint >>> from metpy.units import units >>> virtual_temperature_from_dewpoint(1000 * units.hPa, 30 * units.degC, 25 * units.degC) - + Notes ----- From 7356935042bee9665921d330e4084d9a81d65024 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Wed, 12 Apr 2023 11:12:35 -0500 Subject: [PATCH 17/30] Fix use of parcel moisture profile --- src/metpy/calc/thermo.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index b64a30f30a8..e23c022c8d1 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2301,7 +2301,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') >>> # calculate surface based CAPE/CIN >>> cape_cin(p, T, Td, prof) - (, ) + (, ) See Also -------- @@ -2338,9 +2338,29 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, dewpoint, parcel_profile) + # Calculate the lcl to use for determining the moisture of the parcel profile + pressure_lcl, _ = lcl(pressure[0], temperature[0], dewpoint[0]) + below_lcl = pressure < pressure_lcl + above_lcl = pressure >= pressure_lcl + + # Calculate the dewpoint of the lower part of the parcel profile + sfc_specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint[0]) + mixing_below_lcl = mixing_ratio_from_specific_humidity(sfc_specific_humidity)\ + * np.ones(len(pressure[below_lcl])) + + # Check for the whether the lcl is in the pressure profile, if not, end there + if _greater_or_close(np.nanmin(pressure), pressure_lcl): + parcel_mixing_ratio = mixing_below_lcl + else: + # Use the assumption that the parcel profile is saturated above the LCL and merge + mixing_above_lcl = mixing_ratio_from_relative_humidity(pressure[above_lcl], + temperature[above_lcl], + 1.) + parcel_mixing_ratio = concatenate([mixing_below_lcl, mixing_above_lcl]) + # Convert the temperature/parcel profile to virtual temperature temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint) - parcel_profile = virtual_temperature_from_dewpoint(pressure, parcel_profile, dewpoint) + parcel_profile = virtual_temperature(parcel_profile, parcel_mixing_ratio) # Calculate LFC limit of integration lfc_pressure, _ = lfc(pressure, temperature, dewpoint, From d8c78acb2f0e1fbd46cfaf69d88476b031cd551e Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Wed, 12 Apr 2023 22:05:01 -0500 Subject: [PATCH 18/30] Add fixed booleans --- src/metpy/calc/thermo.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index e23c022c8d1..7d5da0e2f08 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2301,7 +2301,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') >>> # calculate surface based CAPE/CIN >>> cape_cin(p, T, Td, prof) - (, ) + (, ) See Also -------- @@ -2338,10 +2338,9 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' pressure, temperature, dewpoint, parcel_profile = _remove_nans(pressure, temperature, dewpoint, parcel_profile) - # Calculate the lcl to use for determining the moisture of the parcel profile pressure_lcl, _ = lcl(pressure[0], temperature[0], dewpoint[0]) - below_lcl = pressure < pressure_lcl - above_lcl = pressure >= pressure_lcl + below_lcl = pressure > pressure_lcl + above_lcl = pressure <= pressure_lcl # Calculate the dewpoint of the lower part of the parcel profile sfc_specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint[0]) From 404147faa974e0837432abbf5e1207f111ad0c36 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 13 Apr 2023 10:32:16 -0500 Subject: [PATCH 19/30] fix failing tests --- tests/calc/test_thermo.py | 56 +++++++++++++++++++-------------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 2160f85bbcc..60e28887e61 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -1066,8 +1066,8 @@ def test_cape_cin(): dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 75.229846 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -90.004823 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 215.056976 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -9.94798721 * units('joule / kilogram'), 2) def test_cape_cin_no_el(): @@ -1077,8 +1077,8 @@ def test_cape_cin_no_el(): dewpoint = np.array([19., -11.2, -10.8, -10.4]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC') cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 0.08381723 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -90.0048229 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 12.74623773 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -9.947987213 * units('joule / kilogram'), 2) def test_cape_cin_no_lfc(): @@ -1394,8 +1394,8 @@ def test_surface_based_cape_cin(array_class): temperature = array_class([22.2, 14.6, 12., 9.4, 7., -38.], units.celsius) dewpoint = array_class([19., -11.2, -10.8, -10.4, -10., -53.2], units.celsius) cape, cin = surface_based_cape_cin(p, temperature, dewpoint) - assert_almost_equal(cape, 75.2298456 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -137.133362 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 215.05697634 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -33.0633599455 * units('joule / kilogram'), 2) def test_surface_based_cape_cin_with_xarray(): @@ -1418,8 +1418,8 @@ def test_surface_based_cape_cin_with_xarray(): data['temperature'], data['dewpoint'] ) - assert_almost_equal(cape, 75.2298456 * units('joule / kilogram'), 2) - assert_almost_equal(cin, -137.133362 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 215.056976346 * units('joule / kilogram'), 2) + assert_almost_equal(cin, -33.0633599455 * units('joule / kilogram'), 2) def test_profile_with_nans(): @@ -1459,8 +1459,8 @@ def test_most_unstable_cape_cin_surface(): temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 75.2298456 * units('joule / kilogram'), 2) - assert_almost_equal(mucin, -137.133362 * units('joule / kilogram'), 2) + assert_almost_equal(mucape, 215.056976346 * units('joule / kilogram'), 2) + assert_almost_equal(mucin, -33.0633599455 * units('joule / kilogram'), 2) def test_most_unstable_cape_cin(): @@ -1469,8 +1469,8 @@ def test_most_unstable_cape_cin(): temperature = np.array([18.2, 22.2, 17.4, 10., 0., 15]) * units.celsius dewpoint = np.array([19., 19., 14.3, 0., -10., 0.]) * units.celsius mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mucape, 157.354834 * units('joule / kilogram'), 4) - assert_almost_equal(mucin, -32.0881909 * units('joule / kilogram'), 4) + assert_almost_equal(mucape, 173.749389796 * units('joule / kilogram'), 4) + assert_almost_equal(mucin, -20.968278741 * units('joule / kilogram'), 4) def test_mixed_parcel(): @@ -1482,16 +1482,16 @@ def test_mixed_parcel(): dewpoint, depth=250 * units.hPa) assert_almost_equal(parcel_pressure, 959. * units.hPa, 6) - assert_almost_equal(parcel_temperature, 28.7401463 * units.degC, 6) - assert_almost_equal(parcel_dewpoint, 7.1534658 * units.degC, 6) + assert_almost_equal(parcel_temperature, 28.7401463 * units.degC, 2) + assert_almost_equal(parcel_dewpoint, 7.15346585151 * units.degC, 2) def test_mixed_layer_cape_cin(multiple_intersections): """Test the calculation of mixed layer cape/cin.""" pressure, temperature, dewpoint = multiple_intersections mlcape, mlcin = mixed_layer_cape_cin(pressure, temperature, dewpoint) - assert_almost_equal(mlcape, 987.7085876 * units('joule / kilogram'), 2) - assert_almost_equal(mlcin, -20.80235875 * units('joule / kilogram'), 2) + assert_almost_equal(mlcape, 1132.706800436 * units('joule / kilogram'), 2) + assert_almost_equal(mlcin, -13.4809966289 * units('joule / kilogram'), 2) def test_mixed_layer(): @@ -1883,10 +1883,10 @@ def test_el_lfc_most_cape_bottom(): -6.9, -9.5, -12., -14.6, -15.8]) * units.degC lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints, which='most_cape') el_pressure, el_temp = el(levels, temperatures, dewpoints, which='most_cape') - assert_almost_equal(lfc_pressure, 900.73235 * units.hPa, 3) - assert_almost_equal(lfc_temp, 14.6717153 * units.degC, 3) - assert_almost_equal(el_pressure, 849.79989566 * units.hPa, 3) - assert_almost_equal(el_temp, 12.42268288 * units.degC, 3) + assert_almost_equal(lfc_pressure, 714.358842 * units.hPa, 3) + assert_almost_equal(lfc_temp, 5.415421178 * units.degC, 3) + assert_almost_equal(el_pressure, 658.61917328 * units.hPa, 3) + assert_almost_equal(el_temp, 1.951095056 * units.degC, 3) def test_cape_cin_top_el_lfc(multiple_intersections): @@ -1894,8 +1894,8 @@ def test_cape_cin_top_el_lfc(multiple_intersections): levels, temperatures, dewpoints = multiple_intersections parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='top') - assert_almost_equal(cape, 1259.157926 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -97.9763618 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 1345.18884959 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -35.179268355 * units('joule / kilogram'), 3) def test_cape_cin_bottom_el_lfc(multiple_intersections): @@ -1903,8 +1903,8 @@ def test_cape_cin_bottom_el_lfc(multiple_intersections): levels, temperatures, dewpoints = multiple_intersections parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_el='bottom') - assert_almost_equal(cape, 2.20076511 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -8.22115922 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 4.57262449 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -5.9471237534 * units('joule / kilogram'), 3) def test_cape_cin_wide_el_lfc(multiple_intersections): @@ -1913,8 +1913,8 @@ def test_cape_cin_wide_el_lfc(multiple_intersections): parcel_prof = parcel_profile(levels, temperatures[0], dewpoints[0]).to('degC') cape, cin = cape_cin(levels, temperatures, dewpoints, parcel_prof, which_lfc='wide', which_el='wide') - assert_almost_equal(cape, 1259.157926 * units('joule / kilogram'), 3) - assert_almost_equal(cin, -97.9763618 * units('joule / kilogram'), 3) + assert_almost_equal(cape, 1345.1888496 * units('joule / kilogram'), 3) + assert_almost_equal(cin, -35.179268355 * units('joule / kilogram'), 3) def test_cape_cin_custom_profile(): @@ -1924,7 +1924,7 @@ def test_cape_cin_custom_profile(): dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) + 5 * units.delta_degC cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) - assert_almost_equal(cape, 1442.0101452 * units('joule / kilogram'), 2) + assert_almost_equal(cape, 1650.61208729 * units('joule / kilogram'), 2) assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 2) @@ -2010,7 +2010,7 @@ def test_cape_cin_value_error(): -35.9, -26.7, -37.7, -43.1, -33.9, -40.9, -46.1, -34.9, -33.9, -33.7, -33.3, -42.5, -50.3, -49.7, -49.5, -58.3, -61.3]) * units.degC cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint) - expected_cape, expected_cin = 1973.828926 * units('joules/kg'), 0.0 * units('joules/kg') + expected_cape, expected_cin = 2098.688061 * units('joules/kg'), 0.0 * units('joules/kg') assert_almost_equal(cape, expected_cape, 3) assert_almost_equal(cin, expected_cin, 3) From c8d13ab04e20a0a8ce2a2de35f039aed4da4f150 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 13 Apr 2023 10:39:10 -0500 Subject: [PATCH 20/30] fix line continuation --- src/metpy/calc/thermo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 7d5da0e2f08..d522ed4d3d5 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2344,8 +2344,8 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' # Calculate the dewpoint of the lower part of the parcel profile sfc_specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint[0]) - mixing_below_lcl = mixing_ratio_from_specific_humidity(sfc_specific_humidity)\ - * np.ones(len(pressure[below_lcl])) + mixing_below_lcl = mixing_ratio_from_specific_humidity(sfc_specific_humidity) * np.ones( + len(pressure[below_lcl])) # Check for the whether the lcl is in the pressure profile, if not, end there if _greater_or_close(np.nanmin(pressure), pressure_lcl): From dc8ec57804988cf7b9c30f90206ae1fee45532d3 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 13 Apr 2023 14:01:18 -0500 Subject: [PATCH 21/30] fix the docstring values --- src/metpy/calc/thermo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index d522ed4d3d5..f697bec2428 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2301,7 +2301,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' >>> prof = parcel_profile(p, T[0], Td[0]).to('degC') >>> # calculate surface based CAPE/CIN >>> cape_cin(p, T, Td, prof) - (, ) + (, ) See Also -------- @@ -2917,7 +2917,7 @@ def most_unstable_cape_cin(pressure, temperature, dewpoint, **kwargs): >>> Td = dewpoint_from_relative_humidity(T, rh) >>> # calculate most unstbale CAPE/CIN >>> most_unstable_cape_cin(p, T, Td) - (, ) + (, ) See Also -------- @@ -2993,7 +2993,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs): >>> # calculate dewpoint >>> Td = dewpoint_from_relative_humidity(T, rh) >>> mixed_layer_cape_cin(p, T, Td, depth=50 * units.hPa) - (, ) + (, ) See Also -------- From 370e091a8636f80b3fc4fa9687b44c8e7b2bc106 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Thu, 13 Apr 2023 17:23:38 -0500 Subject: [PATCH 22/30] update docstrings to include virtual temperature --- src/metpy/calc/thermo.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index f697bec2428..d5d5c016921 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2311,9 +2311,11 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' ----- Formula adopted from [Hobbs1977]_. - .. math:: \text{CAPE} = -R_d \int_{LFC}^{EL} (T_{parcel} - T_{env}) d\text{ln}(p) + .. math:: \text{CAPE} = -R_d \int_{LFC}^{EL} + (T_{{v}_{parcel}} - T_{{v}_{env}}) d\text{ln}(p) - .. math:: \text{CIN} = -R_d \int_{SFC}^{LFC} (T_{parcel} - T_{env}) d\text{ln}(p) + .. math:: \text{CIN} = -R_d \int_{SFC}^{LFC} + (T_{{v}_{parcel}} - T_{{v}_{env}}) d\text{ln}(p) * :math:`CAPE` is convective available potential energy @@ -2323,8 +2325,8 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' * :math:`SFC` is the level of the surface or beginning of parcel path * :math:`R_d` is the gas constant * :math:`g` is gravitational acceleration - * :math:`T_{parcel}` is the parcel temperature - * :math:`T_{env}` is environment temperature + * :math:`T_{{v}_{parcel}}` is the parcel virtual temperature + * :math:`T_{{v}_{env}}` is environment virtual temperature * :math:`p` is atmospheric pressure Only functions on 1D profiles (not higher-dimension vertical cross sections or grids). From 9015cb23f548ed49394043b749862c5ae96835ae Mon Sep 17 00:00:00 2001 From: Max Grover Date: Tue, 18 Apr 2023 16:40:17 -0500 Subject: [PATCH 23/30] Update src/metpy/calc/thermo.py Co-authored-by: Ryan May --- src/metpy/calc/thermo.py | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index d5d5c016921..0d0b5b11305 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2344,20 +2344,10 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' below_lcl = pressure > pressure_lcl above_lcl = pressure <= pressure_lcl - # Calculate the dewpoint of the lower part of the parcel profile - sfc_specific_humidity = specific_humidity_from_dewpoint(pressure[0], dewpoint[0]) - mixing_below_lcl = mixing_ratio_from_specific_humidity(sfc_specific_humidity) * np.ones( - len(pressure[below_lcl])) - - # Check for the whether the lcl is in the pressure profile, if not, end there - if _greater_or_close(np.nanmin(pressure), pressure_lcl): - parcel_mixing_ratio = mixing_below_lcl - else: - # Use the assumption that the parcel profile is saturated above the LCL and merge - mixing_above_lcl = mixing_ratio_from_relative_humidity(pressure[above_lcl], - temperature[above_lcl], - 1.) - parcel_mixing_ratio = concatenate([mixing_below_lcl, mixing_above_lcl]) + # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated + # based on the temperature above the LCL + parcel_mixing_ratio = np.where(below_lcl, saturation_mixing_ratio(press, dewp), + saturation_mixing_ratio(press, temp)) # Convert the temperature/parcel profile to virtual temperature temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint) From 1a011932c1c40a2228239ca1dfd5dd91bf52b17a Mon Sep 17 00:00:00 2001 From: Max Grover Date: Wed, 19 Apr 2023 07:08:42 -0500 Subject: [PATCH 24/30] Remove above LCL --- src/metpy/calc/thermo.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 0d0b5b11305..8d34643bd04 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2342,7 +2342,6 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' pressure_lcl, _ = lcl(pressure[0], temperature[0], dewpoint[0]) below_lcl = pressure > pressure_lcl - above_lcl = pressure <= pressure_lcl # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated # based on the temperature above the LCL From 87f718a76b18588f2dc4d0d6f360752dbf2ddca8 Mon Sep 17 00:00:00 2001 From: Max Grover Date: Wed, 19 Apr 2023 07:14:54 -0500 Subject: [PATCH 25/30] Fix variable names in numpy.where --- src/metpy/calc/thermo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 8d34643bd04..61fbb4f563f 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2345,8 +2345,8 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated # based on the temperature above the LCL - parcel_mixing_ratio = np.where(below_lcl, saturation_mixing_ratio(press, dewp), - saturation_mixing_ratio(press, temp)) + parcel_mixing_ratio = np.where(below_lcl, saturation_mixing_ratio(presssure, dewpoint), + saturation_mixing_ratio(presssure, temperature)) # Convert the temperature/parcel profile to virtual temperature temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint) From 39e26f8fb4fcf847a441aff0dc83b8267fe00d8b Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Wed, 19 Apr 2023 07:28:14 -0500 Subject: [PATCH 26/30] fix misspelled pressure --- src/metpy/calc/thermo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 61fbb4f563f..f19dc542dd8 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2345,8 +2345,8 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom' # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated # based on the temperature above the LCL - parcel_mixing_ratio = np.where(below_lcl, saturation_mixing_ratio(presssure, dewpoint), - saturation_mixing_ratio(presssure, temperature)) + parcel_mixing_ratio = np.where(below_lcl, saturation_mixing_ratio(pressure, dewpoint), + saturation_mixing_ratio(pressure, temperature)) # Convert the temperature/parcel profile to virtual temperature temperature = virtual_temperature_from_dewpoint(pressure, temperature, dewpoint) From 408600703e7a282b53f72816a6f8549c4d4b7725 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Wed, 19 Apr 2023 07:43:25 -0500 Subject: [PATCH 27/30] fix mlcin --- src/metpy/calc/thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index f19dc542dd8..f056eaaa8be 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -2984,7 +2984,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs): >>> # calculate dewpoint >>> Td = dewpoint_from_relative_humidity(T, rh) >>> mixed_layer_cape_cin(p, T, Td, depth=50 * units.hPa) - (, ) + (, ) See Also -------- From 01d94bdbced4ee38127ea88b42c4b79b64996404 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 24 Apr 2023 09:47:51 -0500 Subject: [PATCH 28/30] update sensitive sounding test --- tests/calc/test_thermo.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 60e28887e61..572cfc2b6e2 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -655,19 +655,19 @@ def test_sensitive_sounding(): 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, + t = units.Quantity([25.1, 24.5, 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.422 * units.mbar, 2) - assert_almost_equal(lfc_temp, 20.498 * units.degC, 2) + assert_almost_equal(lfc_pressure, 952.8445 * units.mbar, 2) + assert_almost_equal(lfc_temp, 20.94469 * units.degC, 2) pos, neg = surface_based_cape_cin(p, t, td) - assert_almost_equal(pos, 0.0 * units('J/kg'), 3) - assert_almost_equal(neg, 0.0 * units('J/kg'), 3) + assert_almost_equal(pos, 0.106791 * units('J/kg'), 3) + assert_almost_equal(neg, -282.620677 * units('J/kg'), 3) def test_lfc_sfc_precision(): From 612eeba0390ed1231a35c06f5fdee71d0e4c4085 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 24 Apr 2023 10:06:50 -0500 Subject: [PATCH 29/30] add suggestion from ryan --- src/metpy/calc/thermo.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index f056eaaa8be..ffd6723cfbc 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1737,11 +1737,8 @@ def virtual_temperature_from_dewpoint( Renamed ``mixing`` parameter to ``mixing_ratio`` """ - # Convert the dewpoint to specific humidity - specific_humidity = specific_humidity_from_dewpoint(pressure, dewpoint) - - # Convert the specific humidity to mixing ratio - mixing_ratio = mixing_ratio_from_specific_humidity(specific_humidity) + # Convert dewpoint to mixing ratio + mixing_ratio = saturation_mixing_ratio(pressure, temperature) # Calculate virtual temperature with given parameters return virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio) From cf9e3878644c6ed997af5b8853ec62a4778b85bf Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Mon, 24 Apr 2023 13:10:36 -0500 Subject: [PATCH 30/30] use dewpoint instead of temperature --- src/metpy/calc/thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index ffd6723cfbc..bd91134eca1 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -1738,7 +1738,7 @@ def virtual_temperature_from_dewpoint( """ # Convert dewpoint to mixing ratio - mixing_ratio = saturation_mixing_ratio(pressure, temperature) + mixing_ratio = saturation_mixing_ratio(pressure, dewpoint) # Calculate virtual temperature with given parameters return virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio)