Enhance Ensemble-Stat to compute neighborhood probability forecasts.

[dwfncar]

This is functionality requested by NOAA/GSD for the DTC VX2 ensemble task.

Reference this Google Doc for details:
https://docs.google.com/document/d/1eepS_dluFguXVT3lzIxfycH64RTi5GVDnXh92onguVk/edit?ts=5bec846c

Also see Craig Schwartz paper:
https://journals.ametsoc.org/doi/full/10.1175/2009WAF2222267.1

Enhance Ensemble-Stat to compute neighborhood probability forecasts and write them the NetCDF output file. Let the user define the neighborhood shape and size to be used. Given a user-defined threshold, support the definition of probabilities in the following ways:
 - frequency of event over the whole cylinder (i.e. for shape=circle) of data.
 - for each grid point, compute a summary value (min, max, mean, median, percentile) and apply the threshold to that value.

*** Question: Can we continue defining the neighborhood size in grid units rather than physical units, like km? ***

In particular, the spatial neighborhood maximum filter is described by Jeff Duda below...

A spatial neighborhood maximum filter: essentially the probability of an event occurring within a radius of a given gridpoint. Example: reflectivity exceeding 40 dBZ. If at any grid point within a neighborhood of an arbitrarily selected radius/size (40 km is commonly used, but it is helpful to do this using a range of neighborhoods, say...10, 20, 30, 40, 50, …, 100 km etc.) reflectivity exceeds 40 dBZ, then the neighborhood probability at the central point of the neighborhood gets set to 1.0 for a given ensemble member. Over a range of ensemble members, you would count the number of ensemble members for which the event occurs. If in 6 of 9 ensemble members there was a grid point within the neighborhood of a given central point where reflectivity exceeded 40 dBZ, then we would set the neighborhood maximum ensemble probability at that central grid point to 0.6667. The output of this technique would be a 2D grid of probabilities (having the same shape as the 2D input grids of ensemble member forecast reflectivities).

[MET-1089] created by johnhg

[JohnHalleyGotway]

This is a feature needed for the 2019 year of the JTTI project. John met with NCAR/MMM on the JTTI project on 20191017. These enhancements are to make neighborhood probabilities available via Ensemble-Stat in MET, per Schwartz et al (2017).

Add new options to the Ensemble-Stat config file for ensemble_flag.nep (neighborhood ensemble probability) and ensemble_flag.nmep (neighborhood max ensemble probability). The idea is pretty straight-forward... just compute the event frequency over the cylinder of data (i.e. multiple ensemble members and spatial area).

Also need to post-process the probabilities by applying a Gaussian smoothing filter.

[JohnHalleyGotway]

Here's the paper about NEPs and NMEPs. Section 2 is all you should need:
https://journals.ametsoc.org/doi/full/10.1175/MWR-D-16-0400.1

[JohnHalleyGotway]

Fixed and merged into develop.

[JohnHalleyGotway]

Jeff Hamilton found that there was a difference in the NMEP output between met-9.0_beta2 and met-9.0_beta3. After investigating, I found that the difference was caused by a bug in beta2 that was fixed in beta3:

I ran ensemble_stat via the debugger and found the reason for the difference in NMEP results between beta2 and beta3. The good news is that beta2 was wrong and beta3 is correct.
Beta3 fixes a bug that was in beta2... and here's the bug:

In ensemble_stat.cc, we apply the smoother the NMEP data:
2595 smooth_dp = smooth_field(prob_dp,
2596 string_to_interpmthd(conf_info.nmep_smooth.method[k].c_str()),
2597 conf_info.nmep_smooth.width[k],
2598 conf_info.nmep_smooth.shape,
2599 conf_info.nmep_smooth.vld_thresh,
2600 conf_info.nmep_smooth.gaussian_dx,
2601 conf_info.nmep_smooth.gaussian_radius);

Here's the definition of that smooth_field() function in data_plane_util.cc:
95 void smooth_field(const DataPlane &dp, DataPlane &smooth_dp,
96 InterpMthd mthd, int width,
97 const GridTemplateFactory::GridTemplates shape,
98 double t, const double gaussian_radius,
99 const double gaussian_dx) {

The bug is the order of the arguments. We pass as input "dx" followed by "radius"... but the function arguments should be "radius" followed by "dx".

In beta3, these options are grouped together in a single struct and the bug is resolved. And that explains the difference.