Feature 1583 es hira (#2056)

met/data/config/ConfigConstants

@@ -88,6 +88,7 @@ AW_MEAN        = 20;
 GAUSSIAN       = 21;
 MAXGAUSS       = 22;
 GEOG_MATCH     = 23;
+HIRA           = 24;
 
 // Interpolation types
 NONE           = 1;

met/docs/Users_Guide/config_options.rst

@@ -610,7 +610,7 @@ using the following entries:
     * MAXGAUSS    to compute the maximum value in the neighborhood
       and apply a Gaussian smoother to the result
 
-    The BEST and GEOG_MATCH interpolation options are not valid for regridding.
+    The BEST, GEOG_MATCH, and HIRA options are not valid for regridding.
 
 * The "width" entry specifies a regridding width, when applicable.
   - width = 4;    To regrid using a 4x4 box or circle with diameter 4.
@@ -1689,7 +1689,10 @@ This dictionary may include the following entries:
     * MAXGAUSS    for the maximum value followed by a Gaussian smoother
 
     * GEOG_MATCH  for the nearest grid point where the land/sea mask
-      and geography criteria are satisfied.
+      and geography criteria are satisfied
+
+    * HIRA        for all neighborhood points to define a spatial
+      ensemble (only in Ensemble-Stat)
 
     The BUDGET, FORCE, GAUSSIAN, and MAXGAUSS methods are not valid for
     interpolating to point locations. For grid-to-grid comparisons, the

met/docs/Users_Guide/ensemble-stat.rst

@@ -17,6 +17,8 @@ Ensemble forecasts derived from a set of deterministic ensemble members
 
 Ensemble forecasts are often created as a set of deterministic forecasts. The ensemble members are rarely used separately. Instead, they can be combined in various ways to produce a forecast. MET can combine the ensemble members into some type of summary forecast according to user specifications. Ensemble means are the most common, and can be paired with the ensemble variance or spread. Maximum, minimum and other summary values are also available, with details in the practical information section.
 
+Typically an ensemble is constructed by selecting a single forecast value from each member for each observation. When the High Resolution Assessment (HiRA) interpolation method is chosen, all of the nearby neighborhood points surrounding each observation from each member are used. Therefore, processing an N-member ensemble using a HiRA neighborhood of size M produces ensemble output with size N*M. This approach fully leverages information from all nearby grid points to evaluate the ensemble quality.
+
 The ensemble relative frequency is the simplest method for turning a set of deterministic forecasts into something resembling a probability forecast. MET will create the ensemble relative frequency as the proportion of ensemble members forecasting some event. For example, if 5 out of 10 ensemble members predict measurable precipitation at a grid location, then the ensemble relative frequency of precipitation will be :math:`5/10=0.5`. If the ensemble relative frequency is calibrated (unlikely) then this could be thought of as a probability of precipitation.
 
 The neighborhood ensemble probability (NEP) and neighborhood maximum ensemble probability (NMEP) methods are described in :ref:`Schwartz and Sobash (2017) <Schwartz-2017>`. They are an extension of the ensemble relative frequencies described above. The NEP value is computed by averaging the relative frequency of the event within the neighborhood over all ensemble members. The NMEP value is computed as the fraction of ensemble members for which the event is occurring somewhere within the surrounding neighborhood. The NMEP output is typically smoothed using a Gaussian kernel filter. The neighborhood sizes and smoothing options can be customized in the configuration file.
@@ -161,6 +163,8 @@ ____________________
 
 The configuration options listed above are common to many MET tools and are described in :numref:`config_options`.
 
+Note that the **HIRA** interpolation method is only supported in Ensemble-Stat.
+
 _____________________
 
 .. code-block:: none

met/scripts/config/EnsembleStatConfig

@@ -249,6 +249,10 @@ interp = {
       {
          method = NEAREST;
          width  = 1;
+      },
+      {
+         method = HIRA;
+         width  = 2;
       }
    ];
 }

met/scripts/config/STATAnalysisConfig

@@ -84,7 +84,7 @@ jobs = [
    "-job aggregate      -line_type GRAD                                          -vx_mask DTC165   -vx_mask DTC166 -by FCST_VAR      -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_GRAD.stat",
    "-job aggregate      -line_type ISC                        -fcst_thresh >0.0  -vx_mask TILE_TOT                 -fcst_var APCP_12 -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_ISC.stat",
    "-job aggregate      -line_type RHIST                      -obtype MC_PCP     -vx_mask HUC4_1605 -vx_mask HUC4_1803 -vx_mask HUC4_1804 -vx_mask HUC4_1805 -vx_mask HUC4_1806 -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_RHIST.stat",
-   "-job aggregate_stat -line_type ORANK -out_line_type RHIST -obtype ADPSFC     -vx_mask HUC4_1605 -vx_mask HUC4_1803 -vx_mask HUC4_1804 -vx_mask HUC4_1805 -vx_mask HUC4_1806 -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_stat_ORANK_RHIST.stat"
+   "-job aggregate_stat -line_type ORANK -out_line_type RHIST -obtype ADPSFC     -vx_mask HUC4_1605,HUC4_1803,HUC4_1804,HUC4_1805,HUC4_1806 -by INTERP_MTHD,INTERP_PNTS -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_stat_ORANK_RHIST.stat"
 ];
 
 ////////////////////////////////////////////////////////////////////////////////

met/src/basic/vx_config/config_util.cc

@@ -141,9 +141,10 @@ RegridInfo::RegridInfo() {
 void RegridInfo::validate() {
 
    // Check for unsupported regridding options
-   if(method == InterpMthd_Best ||
+   if(method == InterpMthd_Best       ||
       method == InterpMthd_Geog_Match ||
-      method == InterpMthd_Gaussian) {
+      method == InterpMthd_Gaussian   ||
+      method == InterpMthd_HiRA) {
       mlog << Error << "\nRegridInfo::validate() -> "
            << "\"" << interpmthd_to_string(method)
            << "\" not valid for regridding, only interpolating.\n\n";
@@ -2272,6 +2273,7 @@ InterpMthd int_to_interpmthd(int i) {
    else if(i == conf_const.lookup_int(interpmthd_gaussian_str))    m = InterpMthd_Gaussian;
    else if(i == conf_const.lookup_int(interpmthd_maxgauss_str))    m = InterpMthd_MaxGauss;
    else if(i == conf_const.lookup_int(interpmthd_geog_match_str))  m = InterpMthd_Geog_Match;
+   else if(i == conf_const.lookup_int(interpmthd_hira_str))        m = InterpMthd_HiRA;
    else {
       mlog << Error << "\nconf_int_to_interpmthd() -> "
            << "Unexpected value of " << i

met/src/basic/vx_util/interp_mthd.cc

@@ -44,6 +44,7 @@ ConcatString interpmthd_to_string(const InterpMthd m) {
       case(InterpMthd_Gaussian):    out = interpmthd_gaussian_str;    break;
       case(InterpMthd_MaxGauss):    out = interpmthd_maxgauss_str;    break;
       case(InterpMthd_Geog_Match):  out = interpmthd_geog_match_str;  break;
+      case(InterpMthd_HiRA):        out = interpmthd_hira_str;        break;
 
       case(InterpMthd_None):
       default:                      out = interpmthd_none_str;        break;
@@ -77,6 +78,7 @@ InterpMthd string_to_interpmthd(const char *mthd_str) {
    else if(strcmp(mthd_str, interpmthd_gaussian_str  )  == 0) m = InterpMthd_Gaussian;
    else if(strcmp(mthd_str, interpmthd_maxgauss_str  )  == 0) m = InterpMthd_MaxGauss;
    else if(strcmp(mthd_str, interpmthd_geog_match_str)  == 0) m = InterpMthd_Geog_Match;
+   else if(strcmp(mthd_str, interpmthd_hira_str)        == 0) m = InterpMthd_HiRA;
    else                                                       m = InterpMthd_None;
 
    return(m);

met/src/basic/vx_util/interp_mthd.h

@@ -42,7 +42,8 @@ enum InterpMthd {
    InterpMthd_Lower_Left,
    InterpMthd_Gaussian,
    InterpMthd_MaxGauss,
-   InterpMthd_Geog_Match
+   InterpMthd_Geog_Match,
+   InterpMthd_HiRA
 };
 
 //
@@ -69,6 +70,7 @@ static const char interpmthd_lower_left_str[]  = "LOWER_LEFT";
 static const char interpmthd_gaussian_str[]    = "GAUSSIAN";
 static const char interpmthd_maxgauss_str[]    = "MAXGAUSS";
 static const char interpmthd_geog_match_str[]  = "GEOG_MATCH";
+static const char interpmthd_hira_str[]        = "HIRA";
 
 ///////////////////////////////////////////////////////////////////////////////
 

met/src/basic/vx_util/interp_util.cc

@@ -1020,6 +1020,7 @@ double compute_sfc_interp(const DataPlane &dp,
          mlog << Error << "\ncompute_sfc_interp() -> "
               << "unsupported interpolation method encountered: "
               << interpmthd_to_string(mthd) << "(" << mthd << ")\n\n";
+         exit(1);
    }
 
    delete gt;

met/src/libcode/vx_data2d/var_info.h

@@ -258,9 +258,9 @@ inline int          VarInfo::accum_attr()     const { return(SetAttrAccum);    }
 //
 
 struct InputInfo {
-    VarInfo * var_info;         // Variable information to read
-    int file_index;             // Index in file_list of file to read
-    StringArray * file_list;    // Array of files (unallocated)
+   VarInfo * var_info;      // Variable information to read
+   int file_index;          // Index in file_list of file to read
+   StringArray * file_list; // Array of files (unallocated)
 };
 
 ////////////////////////////////////////////////////////////////////////
@@ -270,33 +270,37 @@ struct InputInfo {
 //
 class EnsVarInfo {
 
-private:
-    vector<InputInfo> inputs;     // Vector of InputInfo
-    VarInfo * ctrl_info;          // Field info for control member
-public:
-    EnsVarInfo();
-    ~EnsVarInfo();
-    EnsVarInfo(const EnsVarInfo &);
+   private:
+      vector<InputInfo> inputs; // Vector of InputInfo
+      VarInfo * ctrl_info;      // Field info for control member
 
-    void clear();
-    void assign(const EnsVarInfo &);
+   public:
+      EnsVarInfo();
+      ~EnsVarInfo();
+      EnsVarInfo(const EnsVarInfo &);
+
+      void clear();
+      void assign(const EnsVarInfo &);
 
-    void add_input(InputInfo);
-    int inputs_n();
+      void add_input(InputInfo);
+      int inputs_n();
 
-    void set_ctrl(VarInfo *);
-    VarInfo * get_ctrl(int);
+      void set_ctrl(VarInfo *);
+      VarInfo * get_ctrl(int);
 
-    // Get VarInfo from first InputInfo if requested VarInfo is NULL
-    VarInfo * get_var_info(int index=0);
-    ConcatString get_file(int index=0);
-    int get_file_index(int index=0);
+      // Get VarInfo from first InputInfo if requested VarInfo is NULL
+      VarInfo * get_var_info(int index=0);
+      ConcatString get_file(int index=0);
+      int get_file_index(int index=0);
+
+      ConcatString nc_var_str;      // Ensemble variable name strings
+      ThreshArray cat_ta;           // Ensemble categorical thresholds
+      ConcatString raw_magic_str;   // Magic string w/o var substitution
 
-    ConcatString nc_var_str;      // Ensemble variable name strings
-    ThreshArray cat_ta;           // Ensemble categorical thresholds
-    ConcatString raw_magic_str;   // Magic string w/o var substitution
 };
 
+////////////////////////////////////////////////////////////////////////
+
 ConcatString raw_magic_str(Dictionary i_edict, GrdFileType file_type);
 
 ///////////////////////////////////////////////////////////////////////////////

met/src/libcode/vx_statistics/ens_stats.cc

@@ -177,8 +177,8 @@ void ECNTInfo::clear() {
 
    othresh.clear();
    n_ens      = n_pair      = 0;
-   crps_emp   = crpscl_emp  = crpss_emp   = bad_data_double;
-   crps_gaus  = crpscl_gaus = crpss_gaus  = bad_data_double;
+   crps_emp   = crpscl_emp  = crpss_emp  = bad_data_double;
+   crps_gaus  = crpscl_gaus = crpss_gaus = bad_data_double;
    ign        = bad_data_double;
    me         = rmse       = spread      = bad_data_double;
    me_oerr    = rmse_oerr  = spread_oerr = bad_data_double;
@@ -228,6 +228,9 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
    // Store the number of ensemble members
    n_ens = pd.n_ens;
 
+   // HiRA data has no ensemble mean
+   bool skip_mean = pd.mn_na.has(bad_data_double);
+
    // Compute empirical CRPS scores
    crps_emp   = pd.crps_emp_na.wmean(pd.wgt_na);
    crpscl_emp = pd.crpscl_emp_na.wmean(pd.wgt_na);
@@ -257,54 +260,65 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
       }
    }
 
-   // Compute ME and RMSE values
-   fbar = obar = ffbar = oobar = fobar = 0.0;
-   for(i=0; i<pd.n_obs; i++) {
-
-      if(pd.skip_ba[i]) continue;
-
-      // Track running sums
-      w      = pd.wgt_na[i]/w_sum;
-      obar  += w *  pd.o_na[i];
-      oobar += w *  pd.o_na[i] *  pd.o_na[i];
-      fbar  += w * pd.mn_na[i];
-      ffbar += w * pd.mn_na[i] * pd.mn_na[i];
-      fobar += w * pd.mn_na[i] *  pd.o_na[i];
-   }
-
-   // Derive ME and RMSE from partial sums
-   me   = fbar - obar;
-   rmse = sqrt(ffbar + oobar - 2.0*fobar);
-
-   // Compute the square root of the average variance
-   spread = square_root(pd.var_na.wmean(pd.wgt_na));
-
-   // If observation error was specified, compute ME_OERR and RMSE_OERR
-   if(pd.has_obs_error()) {
+   // Compute ensemble mean based statistics
+   if(!skip_mean) {
 
+      // Compute ME and RMSE values
       fbar = obar = ffbar = oobar = fobar = 0.0;
       for(i=0; i<pd.n_obs; i++) {
 
          if(pd.skip_ba[i]) continue;
 
          // Track running sums
          w      = pd.wgt_na[i]/w_sum;
-         obar  += w *       pd.o_na[i];
-         oobar += w *       pd.o_na[i] *       pd.o_na[i];
-         fbar  += w * pd.mn_oerr_na[i];
-         ffbar += w * pd.mn_oerr_na[i] * pd.mn_oerr_na[i];
-         fobar += w * pd.mn_oerr_na[i] *       pd.o_na[i];
+         obar  += w *  pd.o_na[i];
+         oobar += w *  pd.o_na[i] *  pd.o_na[i];
+         fbar  += w * pd.mn_na[i];
+         ffbar += w * pd.mn_na[i] * pd.mn_na[i];
+         fobar += w * pd.mn_na[i] *  pd.o_na[i];
       }
 
-      // Derive ME_OERR and RMSE_OERR from partial sums
-      me_oerr   = fbar - obar;
-      rmse_oerr = sqrt(ffbar + oobar - 2.0*fobar);
+      // Derive ME and RMSE from partial sums
+      me   = fbar - obar;
+      rmse = sqrt(ffbar + oobar - 2.0*fobar);
+
+      // If observation error was specified, compute ME_OERR and RMSE_OERR
+      if(pd.has_obs_error()) {
+
+         fbar = obar = ffbar = oobar = fobar = 0.0;
+         for(i=0; i<pd.n_obs; i++) {
+
+            if(pd.skip_ba[i]) continue;
+
+            // Track running sums
+            w      = pd.wgt_na[i]/w_sum;
+            obar  += w *       pd.o_na[i];
+            oobar += w *       pd.o_na[i] *       pd.o_na[i];
+            fbar  += w * pd.mn_oerr_na[i];
+            ffbar += w * pd.mn_oerr_na[i] * pd.mn_oerr_na[i];
+            fobar += w * pd.mn_oerr_na[i] *       pd.o_na[i];
+         }
+
+         // Derive ME_OERR and RMSE_OERR from partial sums
+         me_oerr   = fbar - obar;
+         rmse_oerr = sqrt(ffbar + oobar - 2.0*fobar);
+      }
+      else {
+         me_oerr   = bad_data_double;
+         rmse_oerr = bad_data_double;
+      }
    }
+   // Set ensemble mean based statistics to bad data
    else {
+      me        = bad_data_double;
+      rmse      = bad_data_double;
       me_oerr   = bad_data_double;
       rmse_oerr = bad_data_double;
    }
 
+   // Compute the square root of the average variance
+   spread = square_root(pd.var_na.wmean(pd.wgt_na));
+
    // Compute the square root of the average perturbed variance
    spread_oerr = square_root(pd.var_oerr_na.wmean(pd.wgt_na));