What is the general workflow to verify ensemble probability?

[lesnow93]

Hello again METplus forum. I went over the ensemble_stat tutorial material (videos, use cases, etc.), yet it's still not obvious what the general workflow should be for taking ensemble member output and verifying the probabilities. I think the warn-on-forecast use case example is closest to what I'm looking to do, but I wanted to get the perspective of those experienced in using ensemble_stat.

Let's go with the most obvious "working" example of accumulated precipitation from an ensemble modeling system consisting of a series of individual deterministic forecasts. Should we first generate the [neighborhood] ensemble relative frequency/probability before running it through ensemble_stat, or does ensemble_stat take care of that [neighborhood] probability generation internally during verification?

Second, is the general work flow for a config file as follows: pcp_combine of precip for each ensemble member to obtain desired accumulation intervals, run ensemble_stat, then run grid_stat to verify probability output stats created by ensemble_stat? And is grid_stat output then used to examine the various ensemble verification metrics (e.g., ROC/AUC, reliability, spread-skill, etc.)?

[JohnHalleyGotway]

@lesnow93, I see you have some questions about verifying ensemble-derived probabilities in METplus.

The current version of MET is 10.1.0 (part of the METplus version 4.1.0 coordinated release). We are currently working on version MET 11.0.0 (part of the METplus version 5.0.0 coordinated release. It's due out in December 2022, and changes in that version impact the answer to your question.

The biggest changes to be aware of are:

1. MET version 10.1.0 added the gen_ens_prod tool which replaces the ensemble product generation step of ensemble-stat (see Create the new Gen-Ens-Prod tool for ensemble product generation. MET#1904).
2. MET version 11.0.0 (in 12/2022) will REMOVE that ensemble product generation step from ensemble-stat (see Remove ensemble post-processing from the Ensemble-Stat tool. MET#1908). We kept it between 10.1.0 and 11.0.0 to give users time to adjust.
3. MET version 11.0.0 (in 12/2022) will ENHANCE ensemble_stat to compute probabilistic statistics on the fly without needing to call grid_stat or point_stat (see Enhance Ensemble-Stat to compute probabilistic statistics for user-defined or climatology-based thresholds. MET#1259).

So the workflow for verifying the probability of precip really depends on what version you're using.

MET version 10.1.0:

1. Run pcp_combine on each of the ensemble members, as needed, to compute the accumulation interval of interest.
   2a. OPTION: Run ensemble_stat on the pcp_combine output to verify the ensemble against gridded and/or point observations. This gives you ranked histograms and continuous ensemble stats, like CRPS. It also writes derived ensemble probabilities to a NetCDF file.
   2b. OPTION: If you don't care about the ensemble stats in 2a, just run gen_ens_prod instead of ensemble_stat. That writes derived ensemble probabilities to a NetCDF file.
2. Run point_stat and/or grid_stat to compare the probabilities from 2a or 2b to point/gridded observations. They writes probabilistic line types, like PCT, PSTD, PRC, and PJC.

MET version 11.0.0 (or the latest development version for it, MET-11.0.0-beta3):

1. SAME AS ABOVE: Run pcp_combine on each of the ensemble members, as needed, to extract the accumulation interval of interest.

Next, you have a choice...

If you care about ensemble statistics and only want simple probability of exceedance probabilities, run ensemble_stat to compute ranked histograms (RHIST line type), continuous ensemble stats (ECNT line type), AND ALSO compute probabilistic statistics to the PCT, PSTD, PRC, PJC line types. And you're done.

The prob_cat_thresh config option for ensemble_stat defines the precip thresholds for which you want probabilities (e.g. >0 or >25.4, something like that).

The prob_pct_thresh config option for ensemble_stat defines how to populate an Nx2 probabilistic contingency table (e.g. PCT output line).

If you want fancier probabilities, like Neighborhood Maximum Ensemble Probabilities (NMEP), run gen_ens_prod on the pcp_combine output to generate them. Then run point_stat and/or grid_stat to compare those probabilities to point/gridded observations. They write probabilistic line types, like PCT, PSTD, PRC, and PJC.

So you can see the answer has changed with the addition of gen_ens_prod and changes to ensemble_stat.

Hope that helps!

[lesnow93]

Thanks so much John H-G for your detailed outline of the workflow. This is exactly what I need to get me going! I can't wait for the MET 11 / METplus 5.0.0 coordinated release in December (looming deadline), so I'll get started with the MET 10.1.0/METplus 4.1.0 installation.

I did view the tutorial info about peeling off the generation of derived ensemble products into gen_ens_prod. Your reply was a helpful reminder of that transition to the next coordinated release. So I will set up my config to generate the ensemble verif stats via 1/2a/3, i.e. pcp_combine on each member, ensemble_stat, then grid_stat to compared against the gridded obs precip data.

I definitely like the new workflow for METplus v5.0.0. Stripping away point_stat/grid_stat in the final step to verify ensemble probability makes a lot of sense.

-Jonathan

[lesnow93]

Since this thread is still open on ensemble_stat, here is a reply on my attempt to get pcp_combine+ensemble_stat working. The format field for precip accumulation interval is not properly substituted by ensemble_stat and as a result, metplus fails because it's looking for the wrong-named input files. Specifically, 'A{level?fmt=%HH}' (which should be A06, A12, etc. and works fine in pcp_combine and grid_stat) is being set to 'A00' for each instance of the FCST_ENSEMBLE_STAT_INPUT_TEMPLATE and OBS_ENSEMBLE_STAT_GRID_INPUT_TEMPLATE.

FYI, I was successful in running ensemble_stat and creating *_ens.nc and *_orank.nc files by hard-coding in 'A06' into the observation and forecast ensemble_stat template environment variables. However, the other files created by ensemble_stat are all of 0 size, so maybe I have a separate issue going on with 0-sized .stat, ecnt, rps, rhist, ssvar, phist, and relp.txt files?

I attached the config and log files for your reference.

Cheers,
Jonathan

hiwat_EnsembleStat_APCP_v1.conf.txt
metplus_201803291800_EnsStat_v1.log

[JohnHalleyGotway]

@lesnow93 while I know a lot about the MET tools themselves, @georgemccabe knows the most about setting up the METplus configuration files properly. I see the following in the config file you attached:

# The settings below aren't working for me. ensemble_stat looks for "A00", instead of "A06", "A12", "A24".
OBS_ENSEMBLE_STAT_VAR1_NAME = APCP
#OBS_ENSEMBLE_STAT_VAR1_LEVELS = A06, A12, A24
#OBS_ENSEMBLE_STAT_VAR1_LEVELS = A06
OBS_ENSEMBLE_STAT_VAR1_LEVELS = "(*,*)"

@georgemccabe do you have any advice for @lesnow93 about processing multiple accumulation intervals in Ensemble-Stat? I believe that for each of 3 accumulation intervals, he want to:

1. Run pcp_combine to compute that accumulation interval.
   • Run -add to add up hourly observations.
   • Run -subtract for each ensemble member to compute the accumulation.
2. Run Ensemble-Stat on the pcp_combine output from 1. to compute statistics.

So he's looping over multiple ensemble members as well as multiple accumulation intervals.

@georgemccabe any advice for how to best handle this situation?

[georgemccabe]

Hi @lesnow93, it looks like 00 is being substituted for {level?fmt=%HH} because you have FCST_ENSEMBLE_STAT_VAR1_LEVELS = "(*,*)". While this is the correct notation for accessing the NetCDF level from the MET output file generated by PCPCombine, it will not properly substitute the accumulation you need in the filename template substitution. MET writes NetCDF files with fields that contain a level attribute that matches the input accumulation from the GRIB input, so you should be able to set:

FCST_ENSEMBLE_STAT_VAR1_NAME = APCP
FCST_ENSEMBLE_STAT_VAR1_LEVELS = A06

It looks like you had tried that before, but it is now commented out. What error message do you get when you set it that way?

[lesnow93]

Hi George,

Yes, unfortunately, I did try the settings as you suggested above (FCST_ENSEMBLE_STAT_VAR1_LEVELS = A06) and I had the same result. So I next tried FCST_ENSEMBLE_STAT_VAR1_LEVELS = "(,)" based on the WoF use case (I believe), but that didn't work either. The substitution becomes 'A00' each time. The only time it works is if I hard-code the template to be 'A06' as opposed to '{level?fmt=%HH}'

-Jonathan

[georgemccabe]

Hi @lesnow93, I did some testing locally and it looks like there is a bug in the EnsembleStat wrapper where the level is not set properly for forecast input templates. It looks like it works fine in the obs templates, but an error is thrown if level is referenced in the output template. I created issue #1910 to fix this. Thank you for bringing this to our attention!

Will you need a bugfix release to METplus v4.1 that includes this fix? The fix will also be included in the upcoming METplus v5.0.0 release.

[georgemccabe]

Hi @lesnow93, Pull request #1911 fixed this issue on the main_v4.1 branch, so you can obtain this fix by pulling the latest changes from that branch.

[lesnow93]

Well, it looks like I still have errors in using the configurable level forecast field. More on that later after I do a bit more debugging.

In the meantime, I'm having trouble understanding the difference between "prob_cat_thresh" and "prob_pct_thresh". I think that prob_cat_thresh is a list of thresholds of the variable of interest over which we will be evaluating [probabilistic] verification results (typically the same as the thresholds used to create ensemble relative frequencies), whereas prob_pct_thresh is a list of probabilities (or size of probability bin widths) for constructing ROC/AUC stats, for example.

I can't seem to find any METplus use cases that provide a clear example of configuring both of these entries. There is an explanation provided in the MET user's guide on page 246 of the PDF, but I think there may be a typo with the example for prob_pct_thresh. About 2/3 of the way down the page, the example shows creating probability bin widths of 0.25 each using "prob_cat_thresh = [ ==0.25 ];", but I think this should say "prob_pct_thresh = [ ==0.25 ];".

Could you let me know if my interpretation of these 2 parameters is correct and if that is a typo in the MET user's guide?

thanks much as always,
Jonathan

[JohnHalleyGotway]

@lesnow93, yes indeed, that was a typo you found. I've corrected it in the main_v10.1, latest, and develop version of the docs. You can see the updated version of the ensemble_stat configuration file section.

Thanks for letting us know! And please let us know about any other typos you find, or any other documentation that is not clear.

Thanks!

[lesnow93]

Thanks for that confirmation, John H-G (and George). I have a couple additional problems/questions for you:

(1) I'm still struggling with how prob_pct_thresh should be set. In my metplus config, I tried setting "ENSEMBLE_STAT_PROB_PCT_THRESH = ==0.10", which seemed to work in using probability thresholds every 10% for the prc (Receiver Operating Characteristic) output file type. I had 11 total thresholds and accompanying PODY and POFD, ranging from 0 to 1, every 0.1. However, the PODY and POFD results were not what I expected. I thought that by definition, the probability=0 threshold should result in PODY and POFD both = 1, and probability=1 threshold should produce a PODY and POFD=0. But, in the output files, I have some fractional value for the prob=0 threshold. Perhaps you can see what might be wrong in my ensemble_stat config file?

(2) Additionally, I have my ensemble_stat METplus config file set up to process a list of 3 accumulation thresholds (A06, A12, and A24), however, the ensemble_stat output files only include A06 accumulation interval.

I attached the METplus config, log, and sample _prc.txt files.

Thanks again for the help,
-Jonathan

ensemble_stat_HIWAT_ensemble_APCP_2018032918_F18_20180330_120000V_prc.txt
metplus_201803291800_EnsStat_v2.log
hiwat_EnsembleStat_APCP_v2.conf.txt

[JohnHalleyGotway]

Jonathan,

You're right that when creating a ROC curve the first point is always (0,0) and the last point is always (1,1). When we create ROC curves, we automatically add those points and then use the values in PRC line to fill in the points in-between. The PCT table DOES NOT compute those (0,0) and (1,1) values numerically. Let me explain using a simpler case of prob_pct_thresh = ==0.25. This defines a 4x2 probabilistic contingency table. For each observation, the event either is or is not occurring. And for each obs, we interpolate the gridded forecast probabilities to that location. That interpolated probability value is in the range [0, 1], and it's placed in one of the probability bins as described below, or see this section of the documentation:

Category 1: 0.00 <= p <0.25
Category 2: 0.25 <= p <0.50
Category 3: 0.50 <= p <0.75
Category 4: 0.75 <= p <=1.00

So there is no explicit bin for being exactly equal to 0. And there is no explicit bin for being exactly equal to 1.0. Each of the 4 bins above includes all probabilities in the defined range.

I'll defer to @georgemccabe to consider why your METplus output only contains 6-hour accumulations.

[georgemccabe]

Hi @lesnow93,

The EnsembleStat wrapper only substitutes the first level/accumulation in the filename templates, even though all of the accumulations are defined in the field information. You are likely only getting the 6 hr accumulation because it is only reading those files.

For now, you should be able to replace {level?fmt=%HH} with * in your forecast input template to read all of the accumulations. Keep in mind that if you have other similarly named files in the directory for other accumulations that you are not processing, then it will include those. This may not be an issue since the field information defines the accumulations to process, so the other information would not be read. I'm not sure if the extra files will affect the valid file/field thresholds, so that may need to be adjusted. You can try this and review the log output to check if everything is being processed as it should.

A more robust solution would be to make enhancements to the wrapper logic to handle this situation in a better way. There are 2 improvements that I could see useful here:

1. Other wrappers like GridStat have an option GRID_STAT_ONCE_PER_FIELD which loops over the field information and calls the tool once for each field name/level. This allows you to process each accumulation separately and write output files that are named accordingly to include the accumulation output in separate files. This option does not currently exist for EnsembleStat wrapper, but it could be implemented.

2. Instead of just using the first field name/level in the filename template substitution, we could adjust the logic to check if you have included the level filename template tag in an input template and expand it to find all files that match each of the level/accumulations you are processing. This would behave similarly to using a wildcard for the level as I mentioned above, but it would ensure that only the requested accumulations are processed if you have other accumulations in the directory.

Let me know if either or both of those suggested enhancements would be help in your case. Also let me know if you are not able to get the desired results using the wildcard character.

[lesnow93]

Perfect. Thanks for the clarification! Will v11 of MET be able to verify the NEP and NMEP within the ensemble_stat run, or will that still require one to run point_stat/grid_stat separately afterward to verify the neighborhood probabilities?

[JohnHalleyGotway]

You will run Gen-Ens-Prod to generate the NEP and NMEP output and then verify them using Point-Stat and/or Grid-Stat.

FYI, we do create development versions of our releases every 6 weeks or so which we call "beta" releases. The current one is beta4 and we're working away on beta5 right now. The development releases can be found on the DTC website and GitHub releases.

You're welcome to download/compile/test them out and provide feedback. However, note that the format of the output may change from version to version. So once the official version is released, please switch to using that. We do not provide support for development releases, other than fixing any bugs that are found in them.

[JohnHalleyGotway]

Just wanted to follow up on this discussion to see if questions remain. If not, please feel free to select the "Mark As Answer" button for one of the previous responses, and then I'll lock it to prevent future posts. If I don't hear otherwise in the next week or so, I'll just go ahead and mark it as being answered.

Thanks!

[lesnow93]

I'm still batch-processing my ensemble runs through the equivalent of Gen-Ens-Prod and Ensemble-Stat (still using ES for both NEP/NMEP and grid point by point scores). I haven't yet gotten to trying to run grid_stat on the NMEP output. I'll try doing that this week on a single daily ensemble run to test out that process. If you could keep this thread open for a bit longer, that'll be great!

[JohnHalleyGotway]

Hi @lesnow93, any followup on this discussion or can we mark it as being answered and create a new discussion if new questions arise?

[lesnow93]

Just go ahead and mark as answered. Thanks! -JonC From: John Halley Gotway ***@***.***> Sent: Tuesday, January 24, 2023 3:13 PM To: dtcenter/METplus ***@***.***> Cc: Case, Jonathan (MSFC-ST11)[ENSCO INC] ***@***.***>; Mention ***@***.***> Subject: [EXTERNAL] Re: [dtcenter/METplus] What is the general workflow to verify ensemble probability? (Discussion #1890) Hi @lesnow93<https://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Flesnow93&data=05%7C01%7CJonathan.Case-1%40nasa.gov%7Ce0803082a0de47796a0908dafe4fcfc4%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C638101916001176669%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=p5A2YK7XxrDwLpALqJICXqas%2FrtU7rDZMZtu%2B55waUM%3D&reserved=0>, any followup on this discussion or can we mark it as being answered and create a new discussion if new questions arise? - Reply to this email directly, view it on GitHub<https://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdtcenter%2FMETplus%2Fdiscussions%2F1890%23discussioncomment-4770879&data=05%7C01%7CJonathan.Case-1%40nasa.gov%7Ce0803082a0de47796a0908dafe4fcfc4%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C638101916001176669%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=OvQiPFovlfHjV3JJTYqgf5FXlwvXHgHM9DZX4fvvSHQ%3D&reserved=0>, or unsubscribe<https://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABITRCJD4BYJQUFY7MEMS6TWUBAWVANCNFSM6AAAAAARNNHICY&data=05%7C01%7CJonathan.Case-1%40nasa.gov%7Ce0803082a0de47796a0908dafe4fcfc4%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C638101916001176669%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=6Mz2pX3ZcIbxj2qE4XyD3cyo2wowV7VCjxkwAcqtfcI%3D&reserved=0>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>