add corfidi MCS motion

[wx4stg]

This adds support for Mesoscale Convective System storm motion as described in Corfidi 2003: https://journals.ametsoc.org/view/journals/wefo/18/6/1520-0434_2003_018_0997_cpampf_2_0_co_2.xml

Note that this implementation differs slightly from that of sharppy/SHARPpy's slightly, in that they use a non-pressure weighted wind for the cloud layer winds (this shouldn't have much of an effect on output in most cases) and they take the low-level jet as the mean of 0->1.5km AGL winds. My implementation looks for the maximum below 1.5 km AGL as I believe this is more in line with the original authors' intentions, "In the absence of a distinct low-level speed maximum in the vertical direction, the strongest wind in the lowest 5000 ft (1.5 km) is generally used, in accordance with Bonner (1968)". I'm extremely open to discussion on whether this is actually "better" or not though.

This is also my first contribution, so feel free to edit my documentation/code style/etc.

Thanks!

• Tests added
• Fully documented

[CLAassistant]

All committers have signed the CLA.

[wx4stg]

seems I didn't completely cover this in my unit test, stand by...

[wx4stg]

I'm on the fence as to whether I should be using pressure_weighted_mean or weighted_continuous average for the "cloud-layer winds". The original paper describes "a vector that represents cell advection by the mean cloud-layer wind (with “cloud layer” taken to be the 850–300-hPa layer)1" (Corfidi 2003) but then uses (wind_at_850+wind_at_700+wind_at_500+wind_at_300)/4 to estimate this.

Thinking more on this, this seems to weight the winds at higher pressures more, so maybe I should switch? It doesn't make a huge difference...

[dopplershift]

@wx4stg Just wanted to let you know that I have been looking at this when I have the cycles and been trying to figure out how I want to proceed on the "cloud-layer winds" problem.

[wx4stg]

yeah--no problem, I meant to attend last week's dev call but forgot (oops), but I plan to start attending those once the semester gets going.

[kgoebber]

Okay, so I've dug into this a bit and have attempted to reproduce the values produced in the paper and have not been able to do so (although very close in most circumstances), even with implementing their form of averaging for the cloud layer wind. My best guess is that the LLJ wind is the cause of the differences I'm seeing. I was not able to get the same value they propose for the 16 August 1997 case using either the 850-hPa LLJ method or the maximum wind in the lowest 1500 m.

My thoughts at this point are that the pressure weighted average is likely better. SharpPy averages the winds over the lower 1500 m for the LLJ and also implements non-pressure weighted averaging for both the LLJ and the cloud layer winds.

@mwilson14 Thoughts since you have thought about a lot of sounding indices in the past.

[kgoebber]

Okay, here is the latest of what I have found when working to compare to the table provided in Confide (2003) with the current version of the implemented vector in the PR. In general, it doesn't match the values in the paper, but is close for some soundings, reasonable for many more, and closer to the observed values in other instances. I tried a couple of variations with only using the four levels and using a mean wind for the LLJ and all of those combinations made the errors greater. The only difference that was fairly small and zero in a number of cases was using the 850-hPa wind for the LLJ instead of the max wind in the lowest 1500 m. My best guess as to the differences is that unique LLJ vectors were used for different cases.

I think we should implement a method within the function to be able to specify a particular LLJ u and v-component to override the default behavior. I think it is an open conversation whether to use the maximum in the lowest 1500 m or just the 850-hPa wind for the LLJ. If we do the latter, then we could reduce the need for the height as an extra input variable - though this would not be fully in line with the what the paper describes as a best practice.

Note: I used Wyoming Sounding Archive and had to change the station name for

• 87081912 OMA to OVN
• 89070212 OKC to OUN
• 89071712 OKC to OUN
• 90052800 SIL to LCH
• 91040912 LIT to LZK

Date	Station	Forecast wdir/spd	Obs wdir/spd	MetPy wdir/spd	Diff w/ MetPy wdir	Diff w/ MetPy spd	FCST-Obs wdir	FCST-Obs spd
83052000	VCT	235/60	260/55	229/58	-31	3	-25	5
83052000	DRT	245/45	265/40	248/57	-17	17	-20	5
83070112	PIA	270/43	290/45	242/34	-48	-11	-20	-2
83071912	BIS	290/57	300/50	264/62	-36	12	-10	7
83072000	GRB	300/54	310/50	296/59	-14	9	-10	4
83072200	WAL	315/53	315/45	308/72	-7	27	0	8
87061912	AMA	295/30	315/30	270/26	-45	-4	-20	0
87070412	UMN	280/30	290/30	282/30	-8	0	-10	0
87070512	UMN	270/32	270/35	306/26	36	-9	0	-3
87071112	STC	225/48	230/44	224/27	-6	-17	-5	4
87072012	STC	225/41	290/35	240/49	-50	14	-65	6
87081912	OVN	290/55	320/30	284/63	-36	33	-30	25
87082000	UMN	310/40	330/30	307/42	-23	12	-20	10
87111600	LCH	260/07	260/10	231/54	-29	44	0	-3
88040600	SLO	255/49	270/40	252/49	-18	9	-15	9
88051000	HTS	280/72	270/30	278/69	8	39	10	42
88071700	ACY	325/24	315/22	340/36	25	14	10	2
88081600	WAL	310/25	300/30	326/26	26	-4	10	-5
88081800	ACY	330/57	330/55	334/60	4	5	0	2
89022112	CHS	230/51	230/37	242/56	12	19	0	14
89050500	SEP	290/53	310/65	288/51	-22	-14	-20	-12
89050600	CHS	240/45	230/45	258/50	28	5	10	0
89061700	WAL	215/46	220/45	225/45	5	-0	-5	1
89070212	OUN	340/50	340/50	337/54	-3	4	0	0
89071712	OUN	300/32	315/35	300/26	-15	-9	-15	-3
89080600	HTS	270/30	290/25	299/18	9	-7	-20	5
89080700	AMA	330/25	350/25	326/29	-24	4	-20	0
89112100	IAD	285/60	285/55	289/86	4	31	0	5
90021012	AHN	240/57	270/50	226/82	-44	32	-30	7
90052800	LCH	250/30	270/30	277/67	7	37	-20	0
90062000	TOP	290/40	270/45	281/45	11	0	20	-5
90082612	GRB	260/27	300/30	241/31	-59	1	-40	-3
91040912	LZK	240/47	240/40	249/42	9	2	0	7
91050500	GGG	240/38	250/30	234/41	-16	11	-10	8
91050600	AYS	270/35	270/40	272/33	2	-7	0	-5
91050700	ACY	240/31	270/40	230/59	-40	19	-30	-9
91052900	FNT	290/37	300/35	296/39	-4	4	-10	2
91070800	FNT	260/50	270/45	258/49	-12	4	-10	5
92070300	FNT	300/48	290/40	301/43	11	3	10	8
92070300	UMN	250/47	270/45	241/43	-29	-2	-20	2
93060412	PAH	270/50	280/52	256/41	-24	-11	-10	-2
93060500	HAT	285/67	300/60	274/67	-26	7	-15	7
93080100	PAH	330/38	320/40	322/25	2	-15	10	-2
93080100	TOP	315/35	320/35	305/44	-15	9	-5	0
94041800	PIA	300/55	280/55	312/69	32	14	20	0
96050512	LZK	270/40	270/70	252/34	-18	-36	0	-30
96052112	CHH	270/50	270/60	263/48	-7	-12	0	-10
97081612	DTX	265/40	275/45	255/51	-20	6	-10	-5

[mwilson14]

This implementation makes sense to me! I agree that the pressure-weighted wind is likely better and that having an option in the function to specify a LLJ vector is a good idea. Also, sorry for missing that first comment I was tagged in, that was the weekend I moved to Boulder so things were a bit chaotic.

[kgoebber]

@wx4stg, okay so I think we have a final recommendation on implementation for this PR!

1. Add in an input option for a u and v LLJ. This will allow a more nimble method for end users who might know or want to specify a certain LLJ instead of using the default option that is implemented in the function.
2. Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.
3. Put a second example in the function that demonstrates finding the maximum LLJ in the lowest 1500 m and putting that LLJ value in the optional variables.

Let @dcamron or I know if you would like any help or want to punt any of these updates to us if you don't have the cycles at the midpoint of the current semester.

[wx4stg]

@kgoebber thanks for the review! I have a few questions..

Add in an input option for a u and v LLJ. This will allow a more nimble method for end users who might know or want to specify a certain LLJ instead of using the default option that is implemented in the function.

Ok, just to clarify, I'll add keyword args, something like llj_u and llj_v, which are None by default. If the function is called without specifying, then proceed to find LLJ as maximum below 850 hPa?

Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.

I agree with the idea of using pressure coordinate height for this, but this introduces a bit of an edge case: what if the sounding is from something like Denver where the surface pressure occasionally dips below 850 hPa? Raise ValueError? High altitude MCSes are... certainly not very typical, but just curious about how to handle this.

Let dcamron or I know if you would like any help or want to punt any of these updates to us if you don't have the cycles at the midpoint of the current semester.

Hilariously, these updates could not have come at a better time, we just finished midterms so I have a little breathing room this week :)

[kgoebber]

@kgoebber thanks for the review! I have a few questions..

Add in an input option for a u and v LLJ. This will allow a more nimble method for end users who might know or want to specify a certain LLJ instead of using the default option that is implemented in the function.

Ok, just to clarify, I'll add keyword args, something like llj_u and llj_v, which are None by default. If the function is called without specifying, then proceed to find LLJ as maximum below 850 hPa?

Yep! So in the end it should be something like...

def corfidi_storm_motion(pressure, u, v, *, u_llj, v_llj):

@dopplershift, @dcamron thoughts on naming of LLJ u and v components?

Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.

I agree with the idea of using pressure coordinate height for this, but this introduces a bit of an edge case: what if the sounding is from something like Denver where the surface pressure occasionally dips below 850 hPa? Raise ValueError? High altitude MCSes are... certainly not very typical, but just curious about how to handle this.

I think raise a value error would be good and prompt to specify an LLJ u and v component. @dopplershift @dcamron thoughts?

Let dcamron or I know if you would like any help or want to punt any of these updates to us if you don't have the cycles at the midpoint of the current semester.

Hilariously, these updates could not have come at a better time, we just finished midterms so I have a little breathing room this week :)

Excellent! Hope midterms went well!

[dopplershift]

def corfidi_storm_motion(pressure, u, v, *, u_llj, v_llj):

@dopplershift, @dcamron thoughts on naming of LLJ u and v components?

While we usually try to avoid terse/abbreviated names so that the function is more self-documenting, u_low_level_jet seems a bit verbose and ugly for something we are requiring as a keyword argument, so I think I'm ok with eg. u_llj--we just of course have to spell out "low-level jet" when documenting what these are.

Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.

I think raise a value error would be good and prompt to specify an LLJ u and v component. @dopplershift @dcamron thoughts?

Agreed. Something like ValueError("Unable to automatically determine low-level jet winds, please pass them in manually with u_llj and v_llj.)

[wx4stg]

@kgoebber I think this covers the requested changes and provides tests, if you still have the code from your earlier comparison, you can try it again.

strangely I can't get the bots to run..

[kgoebber]

There is a conflict with the main branch. There are a couple of ways to resolve the conflict. In general it is best of you update your main MetPy branch

git checkout main

git pull upstream main

git checkout <branch>

git rebase -i main

The last command is the rebase in an interactive mode. It will notify you if there are conflicts that are not resolved automatically that might need your intervention. Once you have rebased, then you can push the branch updates as per usual.

[wx4stg]

I see. I wasn't sure if y'all would want to handle that or not. thanks.

[dopplershift]

Thanks for working with us on getting this in @wx4stg !