diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index fc2836cfdd1..f7a78e3ba58 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -71,6 +71,7 @@ Soundings bulk_shear bunkers_storm_motion + corfidi_storm_motion cape_cin ccl critical_angle diff --git a/docs/api/references.rst b/docs/api/references.rst index 5a3796a5cbe..7191be729fa 100644 --- a/docs/api/references.rst +++ b/docs/api/references.rst @@ -39,6 +39,12 @@ References doi:`10.1175/1520-0434(2000)015\<0061:PSMUAN\>2.0.CO;2 2.0.CO;2>`_. +.. [Corfidi2003] Corfidi, S. F., 2003: Cold Pools and MCS Propagation: + Forecasting the Motion of Downwind-Developing MCSs. + *Wea. Forecasting.*, **18**, 997–1017, + doi:`10.1175/1520-0434(2003)018\<0997:CPAMPF\>2.0.CO;2 + 2.0.CO;2>`_. + .. [CODATA2018] Tiesinga, Eite, Peter J. Mohr, David B. Newell, and Barry N. Taylor, 2020: The 2018 CODATA Recommended Values of the Fundamental Physical Constants (Web Version 8.1). Database developed by J. Baker, M. Douma, and S. Kotochigova. diff --git a/src/metpy/calc/indices.py b/src/metpy/calc/indices.py index e7560015ca2..7b3c349f668 100644 --- a/src/metpy/calc/indices.py +++ b/src/metpy/calc/indices.py @@ -4,6 +4,7 @@ """Contains calculation of various derived indices.""" import numpy as np +from .basic import wind_speed from .thermo import mixing_ratio, saturation_vapor_pressure from .tools import _remove_nans, get_layer from .. import constants as mpconsts @@ -339,6 +340,124 @@ def bunkers_storm_motion(pressure, u, v, height): return right_mover, left_mover, wind_mean +@exporter.export +@preprocess_and_wrap() +@check_units('[pressure]', '[speed]', '[speed]', u_llj='[speed]', v_llj='[speed]') +def corfidi_storm_motion(pressure, u, v, *, u_llj=None, v_llj=None): + r"""Calculate upwind- and downwind-developing MCS storm motions using the Corfidi method. + + Method described by ([Corfidi2003]_): + + * Cloud-layer winds, estimated as the mean of the wind from 850 hPa to 300 hPa + * Convergence along the cold pool, defined as the negative of the low-level jet + (estimated as maximum in winds below 850 hPa unless specified in u_llj and v_llj) + * The vector sum of the above is used to describe upwind propagating MCSes + * Downwind propagating MCSes are taken as propagating along the sum of the cloud-layer wind + and negative of the low level jet, therefore the cloud-level winds are doubled + to re-add storm motion + + Parameters + ---------- + pressure : `pint.Quantity` + Pressure from full profile + + u : `pint.Quantity` + Full profile of the U-component of the wind + + v : `pint.Quantity` + Full profile of the V-component of the wind + + u_llj : `pint.Quantity`, optional + U-component of low-level jet + + v_llj : `pint.Quantity`, optional + V-component of low-level jet + + Returns + ------- + upwind_prop: (`pint.Quantity`, `pint.Quantity`) + Scalar U- and V- components of Corfidi upwind propagating MCS motion + + downwind_prop: (`pint.Quantity`, `pint.Quantity`) + Scalar U- and V- components of Corfidi downwind propagating MCS motion + + Examples + -------- + >>> from metpy.calc import corfidi_storm_motion, wind_components + >>> from metpy.units import units + >>> p = [1000, 925, 850, 700, 500, 400] * units.hPa + >>> wdir = [165, 180, 190, 210, 220, 250] * units.degree + >>> speed = [5, 15, 20, 30, 50, 60] * units.knots + >>> u, v = wind_components(speed, wdir) + >>> corfidi_storm_motion(p, u, v) + (, + ) + >>> # Example with manually specified low-level jet as max wind below 1500m + >>> import numpy as np + >>> h = [100, 250, 700, 1500, 3100, 5720] * units.meters + >>> lowest1500_index = np.argmin(h <= units.Quantity(1500, 'meter')) + >>> llj_index = np.argmax(speed[:lowest1500_index]) + >>> llj_u, llj_v = u[llj_index], v[llj_index] + >>> corfidi_storm_motion(p, u, v, u_llj=llj_u, v_llj=llj_v) + (, + ) + + + 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. + + """ + # If user specifies only one component of low-level jet, raise + if (u_llj is None) ^ (v_llj is None): + raise ValueError('Must specify both u_llj and v_llj or neither') + + # remove nans from input + pressure, u, v = _remove_nans(pressure, u, v) + + # If LLJ specified, use that + if u_llj is not None and v_llj is not None: + # find inverse of low-level jet + llj_inverse = units.Quantity.from_list((-1 * u_llj, -1 * v_llj)) + # If pressure values contain values below 850 hPa, find low-level jet + elif np.max(pressure) >= units.Quantity(850, 'hectopascal'): + # convert u/v wind to wind speed and direction + wind_magnitude = wind_speed(u, v) + # find inverse of low-level jet + lowlevel_index = np.argmin(pressure >= units.Quantity(850, 'hectopascal')) + llj_index = np.argmax(wind_magnitude[:lowlevel_index]) + llj_inverse = units.Quantity.from_list((-u[llj_index], -v[llj_index])) + # If LLJ not specified and maximum pressure value is above 850 hPa, raise + else: + raise ValueError('Must specify low-level jet or ' + 'specify pressure values below 850 hPa') + + # cloud layer mean wind + # don't select outside bounds of given data + if pressure[0] < units.Quantity(850, 'hectopascal'): + bottom = pressure[0] + else: + bottom = units.Quantity(850, 'hectopascal') + if pressure[-1] > units.Quantity(300, 'hectopascal'): + depth = bottom - pressure[-1] + else: + depth = units.Quantity(550, 'hectopascal') + cloud_layer_winds = mean_pressure_weighted(pressure, u, v, + bottom=bottom, + depth=depth) + + cloud_layer_winds = units.Quantity.from_list(cloud_layer_winds) + + # calculate corfidi vectors + upwind = cloud_layer_winds + llj_inverse + + downwind = 2 * cloud_layer_winds + llj_inverse + + return upwind, downwind + + @exporter.export @preprocess_and_wrap() @check_units('[pressure]', '[speed]', '[speed]') diff --git a/tests/calc/test_indices.py b/tests/calc/test_indices.py index 77c5f637d1d..1a5ab52ebcd 100644 --- a/tests/calc/test_indices.py +++ b/tests/calc/test_indices.py @@ -9,7 +9,7 @@ import pytest import xarray as xr -from metpy.calc import (bulk_shear, bunkers_storm_motion, critical_angle, +from metpy.calc import (bulk_shear, bunkers_storm_motion, corfidi_storm_motion, critical_angle, mean_pressure_weighted, precipitable_water, significant_tornado, supercell_composite, weighted_continuous_average) from metpy.testing import (assert_almost_equal, assert_array_almost_equal, get_upper_air_data, @@ -176,6 +176,67 @@ def test_bunkers_motion(): assert_almost_equal(motion.flatten(), truth, 8) +def test_corfidi_motion(): + """Test corfidi MCS motion with observed sounding.""" + data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC') + motion_full = concatenate(corfidi_storm_motion(data['pressure'], + data['u_wind'], data['v_wind'])) + truth_full = [20.60174457, -22.38741441, + 38.32734963, -11.90040377] * units('kt') + assert_almost_equal(motion_full.flatten(), truth_full, 8) + + +def test_corfidi_motion_override_llj(): + """Test corfidi MCS motion with overridden LLJ.""" + data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC') + motion_override = concatenate(corfidi_storm_motion(data['pressure'], + data['u_wind'], data['v_wind'], + u_llj=0 * units('kt'), + v_llj=0 * units('kt'))) + truth_override = [17.72560506, 10.48701063, + 35.45121012, 20.97402126] * units('kt') + assert_almost_equal(motion_override.flatten(), truth_override, 8) + + with pytest.raises(ValueError): + corfidi_storm_motion(data['pressure'], data['u_wind'], + data['v_wind'], u_llj=10 * units('kt')) + + with pytest.raises(ValueError): + corfidi_storm_motion(data['pressure'], data['u_wind'], + data['v_wind'], v_llj=10 * units('kt')) + + +def test_corfidi_corfidi_llj_unaivalable(): + """Test corfidi MCS motion where the LLJ is unailable.""" + data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC') + with pytest.raises(ValueError): + corfidi_storm_motion(data['pressure'][6:], data['u_wind'][6:], data['v_wind'][6:]) + + +def test_corfidi_corfidi_cloudlayer_trimmed(): + """Test corfidi MCS motion where sounding does not include the entire cloud layer.""" + data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC') + motion_no_top = concatenate(corfidi_storm_motion(data['pressure'][:37], + data['u_wind'][:37], data['v_wind'][:37])) + truth_no_top = [20.40419260, -21.43467629, + 37.93224569, -9.99492754] * units('kt') + assert_almost_equal(motion_no_top.flatten(), truth_no_top, 8) + + +def test_corfidi_motion_with_nans(): + """Test corfidi MCS motion with observed sounding with nans.""" + data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC') + u_with_nans = data['u_wind'] + u_with_nans[6:10] = np.nan + v_with_nans = data['v_wind'] + v_with_nans[6:10] = np.nan + motion_with_nans = concatenate(corfidi_storm_motion(data['pressure'], + u_with_nans, v_with_nans)) + truth_with_nans = [20.01078763, -22.65208606, + 37.14543575, -12.42974709] * units('kt') + assert_almost_equal(motion_with_nans.flatten(), truth_with_nans, 8) + + def test_bunkers_motion_with_nans(): """Test Bunkers storm motion with observed sounding.""" data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')