Skip to content

Commit

Permalink
look for maximum below 850hPa, allow u_llj and v_llj override
Browse files Browse the repository at this point in the history
  • Loading branch information
@wx4stg
wx4stg committed Oct 30, 2023
1 parent 3be4f1d commit d0d86ae
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 36 deletions.
60 changes: 42 additions & 18 deletions src/metpy/calc/indices.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -339,15 +339,15 @@ def bunkers_storm_motion(pressure, u, v, height):

@exporter.export
@preprocess_and_wrap()
@check_units('[pressure]', '[speed]', '[speed]', '[length]')
def corfidi_storm_motion(pressure, u, v, height):
@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 1.5 km AGL)
(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
Expand All @@ -364,8 +364,11 @@ def corfidi_storm_motion(pressure, u, v, height):
v : `pint.Quantity`
Full profile of the V-component of the wind
height : `pint.Quantity`
Full profile of height
u_llj : `pint.Quantity`, optional
U-component of low-level jet
v_llj : `pint.Quantity`, optional
V-component of low-level jet
Returns
-------
Expand All @@ -380,13 +383,22 @@ def corfidi_storm_motion(pressure, u, v, height):
>>> from metpy.calc import corfidi_storm_motion, wind_components
>>> from metpy.units import units
>>> p = [1000, 925, 850, 700, 500, 400] * units.hPa
>>> h = [250, 700, 1500, 3100, 5720, 7120] * units.meters
>>> 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, h)
>>> corfidi_storm_motion(p, u, v)
(<Quantity([16.4274315 7.75758388], 'knot')>,
<Quantity([36.32782655 35.21132283], 'knot')>)
>>> # 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)
(<Quantity([4.90039505 1.47297683], 'knot')>,
<Quantity([24.8007901 28.92671577], 'knot')>)
Notes
-----
Expand All @@ -395,17 +407,29 @@ def corfidi_storm_motion(pressure, u, v, height):
Quantities even when given xarray DataArray profiles.
"""
# remove nans from input
pressure, u, v, height = _remove_nans(pressure, u, v, height)
# convert height to AGL
height = height - height[0]
# 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')

# 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(height <= units.Quantity(1500, 'meter'))
llj_index = np.argmax(wind_magnitude[:lowlevel_index])
llj_inverse = units.Quantity.from_list((-u[llj_index], -v[llj_index]))
# 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
Expand All @@ -417,7 +441,7 @@ def corfidi_storm_motion(pressure, u, v, height):
depth = bottom - pressure[-1]
else:
depth = units.Quantity(550, 'hectopascal')
cloud_layer_winds = mean_pressure_weighted(pressure, u, v, height=height,
cloud_layer_winds = mean_pressure_weighted(pressure, u, v,
bottom=bottom,
depth=depth)

Expand Down
48 changes: 30 additions & 18 deletions tests/calc/test_indices.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,33 +180,45 @@ 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'],
data['height']))
truth_full = [4.38681947, -26.16100158,
22.11242453, -15.67399095] * units('kt')
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)
motion_no_bottom = concatenate(corfidi_storm_motion(data['pressure'][6:],
data['u_wind'][6:], data['v_wind'][6:],
data['height'][6:]))
truth_no_bottom = [5.70216267, -29.03497766,
24.74311094, -21.42194311] * units('kt')
assert_almost_equal(motion_no_bottom.flatten(), truth_no_bottom, 8)

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'][6:], data['u_wind'][6:], data['v_wind'][6:])

motion_no_top = concatenate(corfidi_storm_motion(data['pressure'][:37],
data['u_wind'][:37], data['v_wind'][:37],
data['height'][:37]))
truth_no_top = [4.1892675, -25.20826346,
21.71732059, -13.76851471] * units('kt')
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)

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,
data['height']))
truth_with_nans = [6.6604286, -26.30560547,
23.79507672, -16.0832665] * units('kt')
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)

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_bulk_shear():
Expand Down

0 comments on commit d0d86ae

Please sign in to comment.