Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

does MetPy's heat_index match NWS heat_index calculation by Rothfusz 1990? #1096

Closed
akrherz opened this issue Jul 16, 2019 · 13 comments · Fixed by #1098
Closed

does MetPy's heat_index match NWS heat_index calculation by Rothfusz 1990? #1096

akrherz opened this issue Jul 16, 2019 · 13 comments · Fixed by #1098
Labels
Area: Calc Pertains to calculations Type: Enhancement Enhancement to existing functionality
Milestone

Comments

@akrherz
Copy link
Contributor

akrherz commented Jul 16, 2019

The present MetPy heat_index calculation has a doc of

The implementation uses the formula outlined in [Rothfusz1990]. This equation is a multi-variable least-squares regression of the values obtained in [Steadman1979].

But the Weather Prediction Center website on the topic shows a number of corrections and then contains the confusing verbiage

The Rothfusz regression is not valid for extreme temperature and relative humidity conditions beyond the range of data considered by Steadman.

So if for nothing else than to document this situation, does MetPy need an additional heat index formulation? and/or the adjustments included for certain temperature and humidity ranges?

@akrherz
Copy link
Contributor Author

akrherz commented Jul 17, 2019

AWIPS II python code includes adjustments. The plot thickens!

@akrherz
Copy link
Contributor Author

akrherz commented Jul 17, 2019

GEMPAK's prheat.f

So the question I am attempting to resolve now is if the 40% RH threshold is right for the Rothfusz equation used by MetPy or that was only meant to be applied to Steadman.

@akrherz
Copy link
Contributor Author

akrherz commented Jul 17, 2019

So an attempt to summarize my two cents.

It is important to understand that there is no exact heat index equation, but my opinion is that it would be good to closely match what the National Weather Service does.

e6dcd58 introduced the heat_index equation to MetPy back on 19 Dec 2008 and I am sure @dopplershift has a great recollection of his coding thoughts that day :) This code implemented the Rothfusz 1990 multiple regression of the tables found in Steadman 1979. The MetPy code comments denoted that this regression is undefined for temperatures below 80F or relative humidities below 40%.

So the initial question, is where did those two constraints come from? The terse Rothfusz 1990 paper does not seem to mention them. My assumption is that the GEMPAK code was consulted for this and it contains the following code comment:

The Rothfusz regression is optimal for TMPF > ~80 and RELH > ~40%

But the present day GEMPAK code, at least, only checks the TMPF requirement. It also has a curious TMPF < 40F requirement, which may be some glorious typo / long standing bug? The AWIPS II code also doesn't have a 40% requirement.

So I am wondering what to do here? :) I think the path forward, that mostly matches what NWS does.

  • Drop the 40% relative humidity condition for masking.
  • Leave the 80F threshold for the mask_undefined flag.
  • Implement the corrections to Rothfusz found in GEMPAK and AWIPS II.

I am willing to PR this and would love to be told how I have messed up this assessment! :)

@dopplershift
Copy link
Member

I think your assessment is sound. The only context I can add for the original commit is that the driving force behind adding the heat_index calculation was to mimic the meteograms produced by the Oklahoma Mesonet. The only hypothesis I have at this point is that back then the Mesonet site listed that limitation back then--currently they don't list the RH limitation. Actually, I just stumbled upon:

image

which caps out at 40%. It's entirely possible that's the basis I was using (it would seem to be consistent with the undefined <80F as well). However, my feelings now are that we have no scholarly evidence, or any citation whatsoever, for an exact RH>40% limitation, so I believe we should drop it.

In searching for documentation on RH, I found the following:

I think we should do what's documented as the NWS approach, since that's objectively the most consistent with the Steadman table. I actually don't see any justification in the papers of cutting off below 80, only changing formulas. So I think the path forward is:

  • Implement the documented formula/algorithm from the papers/site
  • Cite all of those resources

I don't see any mention in those resources about the 80F limit, but I do see practical use of it, like GEMPAK and the Oklahoma Mesonet. So should we keep it as an option with mask_undefined?

@akrherz
Copy link
Contributor Author

akrherz commented Jul 17, 2019

Thanks @dopplershift for the thorough comment. I wonder if the 80F 40% threshold is a back-handed way to always keep the heat_index >= air temperature. Do you see it as a problem that this proposed change would allow the heat_index to return values less than the temperature?

I am fine with keeping the current default logic of masking heat_index < 80F by default.

@dopplershift
Copy link
Member

I don't have a problem with heat_index returning values less than the temperature--the original table has such values.

@zbruick
Copy link
Contributor

zbruick commented Jul 17, 2019

And I don't have an off-hand example of this either, but I've seen NWS have heat indices below temperature forecasts in the grids too, so I think this is pretty common.

@akrherz
Copy link
Contributor Author

akrherz commented Jul 18, 2019

As I continue to review this confusing fun, it is unclear to me that what I PR'd last night is fully replicating what the NWS does. I also suspect there is no canonical NWS implementation of heat index :) The Anderson 2013 paper seems to imply that the canonical implementation is the javascript code on the WPC website, lol.

The AWIPS II code has a hard coded setting of the heat_index to the temperature when the temperature is less than 80F. GEMPAK does not have this.

En-lieu of something more authoritative, I think we are good to go with this change.

@akrherz
Copy link
Contributor Author

akrherz commented Jul 19, 2019

Yo MetPy, feast your eyes on this, bah ha ha

new_heatindex

Oh noes, that plot does not exactly match values with plot posted above! For example, 80F and 100% RH, Metpy patched has 89 and NWS chart has 87.... and what does official WPC javascript say: 89! We win. The official chart is not right, lol.

I need to get a life.

Gist with code

@akrherz
Copy link
Contributor Author

akrherz commented Jul 19, 2019

Might as well show what patched MetPy now computes for Heat Index by Dew Point.
hi

@dopplershift
Copy link
Member

That's amazing. I would love a PR that adds that as an example. Would prefer it without seaborn, but might even be worth it with seaborn.

@akrherz
Copy link
Contributor Author

akrherz commented Jul 20, 2019

Thanks for the kind words @dopplershift , there is no need to introduce a seaborn requirement just for this plot. I'll just convert the heatmap usage into pure matplotlib code and make a PR.

@akrherz
Copy link
Contributor Author

akrherz commented Jul 22, 2019

@dopplershift I updated the gist to remove the seaborn usage. I'm still interested in adding this to the examples, but I can't get sphinx-gallery installed on my conda-forge env today, sigh. Will try later.

ahuang11 pushed a commit to ahuang11/MetPy that referenced this issue Aug 30, 2019
 - adds additional corrections used by the NWS operational heat index
 - removes default lower bound of 40% Relative Humidity for undefined HI
 - updates tests as this change produces different test results
 - closes Unidata#1096
@dopplershift dopplershift added Area: Calc Pertains to calculations Type: Enhancement Enhancement to existing functionality labels Sep 27, 2019
@dopplershift dopplershift added this to the 0.11 milestone Sep 27, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Area: Calc Pertains to calculations Type: Enhancement Enhancement to existing functionality
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants