Apply NMEP to observation field

[RogerHar]

In the paper that introduced NMEP, Schwartz and Sobash (2017), section 2b(3) says "When objectively verifying NMEPs against observations, for consistency, application of neighborhoods to observations must follow NMEP procedures to determine observed event occurrence". So you should verify the NMEP probababilities generated from an ensemble against a binary field indicating whether the event was observed anywhere within the neighborhood surrounding each gridpoint (assuming gridded observations/analyses).

What's the best way to do that in MET 11? NMEP is explained in the Gen-Ens-Prod chapter (section 9.2.1), though I'm a bit confused to see that the NMEP config settings are in section 5.3.1 EnsembleStatConfig of the Configuration File Overview. Would I be right in thinking they only apply to the ensemble, not the observations? I guess I could treat each observation field as an ensemble of size 1 and run them through Gen-Ens-Prod too separately. Or is there a better way?

[JohnHalleyGotway]

Hi Roger, sorry for the delay in responding on this one. I did discuss with this Christina Kalb (@CPKalb) who's used NMEP in some METplus ensemble verification. And I also sent this discussion link to Craig Schwartz (@weather4evr). Michelle Harrold (@michelleharrold) and Jeff Beck (@JeffBeck-NOAA) have also thought through similar questions. I'm hoping they can address the scientific considerations when comparing model-derived NMEP fields to observations.

But I can speak to the configuration options.

First, thanks for pointing out that NMEP appears in both Gen-Ens-Prod section 9.2.1 and Ensemble-Stat section 5.3.1. While we did indeed remove the ens dictionary from the Ensemble-Stat config file in MET version 11.0.0, we (perhaps not unsurprisingly) failed to clarify the change in that section of the config file. I'll write up a documentation issue to clarify those docs for the next release.

The NMEP output is indeed only available from the Gen-Ens-Prod tool. So it can only be derived from the ensemble data, not the observations. I always have to think very carefully about the NEP and NMEP values. Their definition requires careful consideration when deciding how to compare them to observations. For one ensemble evaluation project with @michelleharrold and @JeffBeck-NOAA, I believe they chose to skip NMEP because of the complications it introduced.

NEP is the fraction of the grid points at which the specified event occurs, computed across the spatial neighborhood and ensemble members. So the denominator is the number of neighborhood grid points times the number of members.

NMEP is the fraction of the ensemble members for which the specified event occurs somewhere in the spatial neighborhood. So the denominator is simply the number of ensemble members.

The User's Guide notes that these ensemble products can be evaluated as probabilities in Grid-Stat or Point-Stat.

But the big question is what, if any, additional processing should be done to the observations to make for a "fair comparison"?
I'll note there are multiple ways you could evaluate the NEP and NMEP output.

1. You could apply MET's probabilistic verification methods in Point-Stat and Grid-Stat and compute reliability, resolution, and so on. If you interpret NEP and NMEP as simply being ensemble-derived probabilities, perhaps no additional processing of the observations is required?

2. You could run Grid-Stat to compute fractions skill score. But I believe that's where things get confusing. We decided it was reasonable for NEP but not for NMEP. For NEP, I believe we set the nbrhd.field = OBS; so that the fractional coverage computations were only applied to the observation data, and we used the raw NEP input values as-is. Of course, we made sure to configure Grid-Stat to use the same threshold and neighborhood size as was used to create the NEP field in Gen-Ens-Prod in the first place.

@RogerHar I do see that we discussed something similar in this issue comment. At the time, it was decided that computing a true ensemble-based fractions skill score was not widely done and not yet needed in MET.

Hopefully @weather4evr and @michelleharrold can weigh in on these details.

[JohnHalleyGotway]

@RogerHar I wanted to add another comment about verifying NEP and NMEP probabilities using Grid-Stat. The NEP value that's computed by Gen-Ens-Prod is the probability that the event (i.e. precip > 0) will occur somewhere within the neighborhood. As such, verifying it against the raw observation value at each grid point, as I suggested above, may not be best.

Another option to consider is "smoothing" the gridded observation field by replacing the value at each grid point with the maximum of the values in the surrounding neighborhood.

For example, the following interp setting defines the value at each grid point as the maximum in a circle with diameter 5 centered on that point.

interp = {
   field      = OBS;
   vld_thresh = 1.0;
   type = [ 
      { method = MAX; width = 5; shape = CIRCLE; }
   ];
}

I am not saying that this is the "correct" thing to do... just pointing it out as a potentially useful configuration option.

[JeffBeck-NOAA]

@RogerHar I wanted to add another comment about verifying NEP and NMEP probabilities using Grid-Stat. The NEP value that's computed by Gen-Ens-Prod is the probability that the event (i.e. precip > 0) will occur somewhere within the neighborhood. As such, verifying it against the raw observation value at each grid point, as I suggested above, may not be best.

Hi @JohnHalleyGotway, did you mean to say "The NMEP value that's computed by Gen-Ens-Prod..."?

[RogerHar]

Ah, thanks John, while it's probably not the only 'correct' thing to do, that does sound like it's exactly what I want to do and attempted to describe above. I didn't know I could process the obs field in that way - in fact I'd never noticed I could use 'interp' within tools other than Point-Stat. So I can run Gen-Ens-Prod on my ensemble fields to get an 'NMEP field', then run this and my observation field through Grid-Stat with an 'interp' setting similar to the one you suggest here. Excellent, I'll mark this as the answer.

[JeffBeck-NOAA]

The NMEP output is indeed only available from the Gen-Ens-Prod tool. So it can only be derived from the ensemble data, not the observations. I always have to think very carefully about the NEP and NMEP values. Their definition requires careful consideration when deciding how to compare them to observations. For one ensemble evaluation project with @michelleharrold and @JeffBeck-NOAA, I believe they chose to skip NMEP because of the complications it introduced.

Yes, we decided to forgo verification of NMEP for now, due to uncertainty on how best to process the observations. However, @JohnHalleyGotway's suggestion to "replace the observation value at each grid point with the maximum of the values in the surrounding neighborhood" sounds like a valid "apples-to-apples" comparison between obs and NMEP. @michelleharrold, @weather4evr, thoughts?

[weather4evr]

I believe what @RogerHar marked as correct is indeed the right answer, so long as the model field is already an NMEP (presumably produced by Gen-Ens-Prod). Then, both the observations and forecast probabilities are "the likelihood of an event occurring within X km of a point". Care would need to be taken to ensure the neighborhood size and exceedance thresholds are kept constant between Gen-Ens-Prod and grid_stat.

[weather4evr]

Also, I've been working on a way to do all this with python embedding, such that the Gen-Ens-Prod step can be removed altogether...more on this later.

[michelleharrold]

I don't have much to add (but this has been a great thread to summarize several internal DTC conversations we have had throughout the last year!) -- Our testing and evaluation decided to only use the NMEP output from gen-ens-prod (no vx component) because we were wavering on how to best implement the obs processing. I think @JohnHalleyGotway provided a great path forward for obs processing in grid-stat to match output in gen-ens-prod, and I am excited to seethe Python embedding solution @weather4evr is working on.