Per #1575, replace PairDataPoint::subset_pairs_mpr_thresh() with a ut…

…ility function named apply_mpr_thresh_mask(). This is for Grid-Stat to apply the mpr_thresh settings after the DataPlane pairs have been created but prior to applying any smoothing operations.

met/src/libcode/vx_statistics/pair_data_point.cc

@@ -292,75 +292,6 @@ PairDataPoint PairDataPoint::subset_pairs_cnt_thresh(
    return(out_pd);
 }
 
-////////////////////////////////////////////////////////////////////////
-
-PairDataPoint PairDataPoint::subset_pairs_mpr_thresh(
-                 const StringArray &sa, const ThreshArray &ta) const {
-
-   // Check for no work to be done
-   if(sa.n() == 0 && ta.n() == 0) {
-      return(*this);
-   }
-
-   // Check for constant length
-   if(sa.n() != ta.n()) {
-      mlog << Error << "\nPairDataPoint::subset_pairs_mpr_thresh() -> "
-           << "the \"" << conf_key_mpr_column << "\" ("
-           << write_css(sa) << ") and \"" << conf_key_mpr_thresh
-           << "\" (" << write_css(ta)
-           << ") config file entries must have the same length!\n\n";
-      exit(1);
-   }
-
-   int i;
-   PairDataPoint out_pd;
-
-   // Allocate memory for output pairs
-   out_pd.extend(n_obs);
-   out_pd.set_climo_cdf_info(cdf_info);
-
-   bool cmn_flag = set_climo_flag(f_na, cmn_na);
-   bool csd_flag = set_climo_flag(f_na, csd_na);
-   bool wgt_flag = set_climo_flag(f_na, wgt_na);
-
-   // Loop over the pairs
-   for(i=0; i<n_obs; i++) {
-
-      // Check for bad data
-      if(is_bad_data(f_na[i])                 ||
-         is_bad_data(o_na[i])                 ||
-         (cmn_flag && is_bad_data(cmn_na[i])) ||
-         (csd_flag && is_bad_data(csd_na[i])) ||
-         (wgt_flag && is_bad_data(wgt_na[i]))) continue;
-
-      // Keep pairs which meet the threshold criteria
-      if(check_mpr_thresh(f_na[i], o_na[i], cmn_na[i], csd_na[i],
-                          sa, ta)) {
-
-         // Handle point data
-         if(is_point_vx()) {
-            out_pd.add_point_pair(sid_sa[i].c_str(), lat_na[i],
-                      lon_na[i], x_na[i], y_na[i],
-                      vld_ta[i], lvl_na[i], elv_na[i],
-                      f_na[i], o_na[i], o_qc_sa[i].c_str(),
-                      cmn_na[i], csd_na[i], wgt_na[i]);
-         }
-         // Handle gridded data
-         else {
-            out_pd.add_grid_pair(f_na[i], o_na[i], cmn_na[i],
-                      csd_na[i], wgt_na[i]);
-         }
-      }
-   } // end for
-
-   mlog << Debug(3)
-        << "Using " << out_pd.n_obs << " of " << n_obs
-        << " pairs for matched pair filtering columns (" << write_css(sa)
-        << ") and thresholds (" << ta.get_str() << ").\n";
-
-   return(out_pd);
-}
-
 ////////////////////////////////////////////////////////////////////////
 //
 // Code for class VxPairDataPoint
@@ -1651,6 +1582,66 @@ double get_mpr_column_value(double f, double o, double cmn, double csd,
 
 ////////////////////////////////////////////////////////////////////////
 
+void apply_mpr_thresh_mask(DataPlane &fcst_dp, DataPlane &obs_dp,
+                           DataPlane &cmn_dp, DataPlane &csd_dp,
+                           const StringArray &col_sa, const ThreshArray &col_ta) {
+
+   // Check for no work to be done
+   if(col_sa.n() == 0 && col_ta.n() == 0) return;
+
+   // Check for constant length
+   if(col_sa.n() != col_ta.n()) {
+      mlog << Error << "\napply_mpr_thresh_mask() -> "
+           << "the \"" << conf_key_mpr_column << "\" ("
+           << write_css(col_sa) << ") and \"" << conf_key_mpr_thresh
+           << "\" (" << write_css(col_ta)
+           << ") config file entries must have the same length!\n\n";
+      exit(1);
+   }
+
+   int  nxy = fcst_dp.nx() * fcst_dp.ny();
+   int  n_skip = 0;
+   bool cmn_flag = !(cmn_dp.is_empty());
+   bool csd_flag = !(csd_dp.is_empty());
+
+   // Loop over the pairs
+   for(int i=0; i<nxy; i++) {
+
+      double cmn = (cmn_flag ? cmn_dp.buf()[i] : bad_data_double);
+      double csd = (csd_flag ? csd_dp.buf()[i] : bad_data_double);
+
+      // Check for bad data
+      if(is_bad_data(fcst_dp.buf()[i])  ||
+         is_bad_data(obs_dp.buf()[i])   ||
+         (cmn_flag && is_bad_data(cmn)) ||
+         (csd_flag && is_bad_data(csd))) continue;
+
+      // Discard pairs which do not meet the threshold criteria
+      if(!check_mpr_thresh(fcst_dp.buf()[i], obs_dp.buf()[i], cmn, csd,
+                           col_sa, col_ta)) {
+
+         // Increment skip counter
+         n_skip++;
+
+         // Set point to bad data
+         fcst_dp.buf()[i] = bad_data_double;
+         obs_dp.buf()[i]  = bad_data_double;
+         if(cmn_flag) cmn_dp.buf()[i] = bad_data_double;
+         if(csd_flag) csd_dp.buf()[i] = bad_data_double;
+      }
+   } // end for i
+
+   mlog << Debug(3)
+        << "Discarded " << n_skip << " of " << nxy
+        << " pairs for matched pair filtering columns ("
+        << write_css(col_sa) << ") and thresholds ("
+        << col_ta.get_str() << ").\n";
+
+   return;
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void subset_wind_pairs(const PairDataPoint &pd_u, const PairDataPoint &pd_v,
                        const SingleThresh &ft, const SingleThresh &ot,
                        const SetLogic type,

met/src/libcode/vx_statistics/pair_data_point.h

@@ -68,9 +68,6 @@ class PairDataPoint : public PairBase {
       PairDataPoint subset_pairs_cnt_thresh(const SingleThresh &ft,
                                             const SingleThresh &ot,
                                             const SetLogic type) const;
-
-      PairDataPoint subset_pairs_mpr_thresh(const StringArray &sa,
-                                            const ThreshArray &ta) const;
 };
 
 ////////////////////////////////////////////////////////////////////////
@@ -267,6 +264,10 @@ extern bool check_mpr_thresh(double, double, double, double,
 extern double get_mpr_column_value(double, double, double, double,
                         const char *);
 
+extern void apply_mpr_thresh_mask(DataPlane &, DataPlane &,
+                        DataPlane &, DataPlane &,
+                        const StringArray &, const ThreshArray &);
+
 // Apply conditional thresholds to subset the wind pairs
 extern void subset_wind_pairs(const PairDataPoint &,
                         const PairDataPoint &, const SingleThresh &,

met/src/tools/core/grid_stat/grid_stat.cc

@@ -713,6 +713,13 @@ void process_scores() {
            << " climatology standard deviation field(s) for forecast "
            << conf_info.vx_opt[i].fcst_info->magic_str() << ".\n";
 
+      // Apply MPR threshold filters
+      if(conf_info.vx_opt[i].mpr_sa.n() > 0) {
+         apply_mpr_thresh_mask(fcst_dp, obs_dp, cmn_dp, csd_dp,
+                               conf_info.vx_opt[i].mpr_sa,
+                               conf_info.vx_opt[i].mpr_ta);
+      }
+
       // Setup the first pass through the data
       if(is_first_pass) setup_first_pass(fcst_dp);
 
@@ -1831,11 +1838,6 @@ void get_mask_points(const GridStatVxOpt &vx_opt,
 
    if(cmn_ptr && csd_ptr) pd.add_climo_cdf();
 
-   // Apply any MPR filters
-   if(vx_opt.mpr_sa.n() > 0) {
-      pd = pd.subset_pairs_mpr_thresh(vx_opt.mpr_sa, vx_opt.mpr_ta);
-   }
-
    return;
 }