Feature 1764 point stat

met/docs/Users_Guide/point-stat.rst

@@ -121,7 +121,7 @@ The forecast value at P is chosen as the grid point inside the interpolation are
 HiRA framework
 ~~~~~~~~~~~~~~
 
-The Point-Stat tool has been enhanced to include the High Resolution Assessment (HiRA) verification logic (:ref:`Mittermaier, 2014 <Mittermaier-2014>`). HiRA is analogous to neighborhood verification but for point observations. The HiRA logic interprets the forecast values surrounding each point observation as an ensemble forecast. These ensemble values are processed in two ways. First, the ensemble continuous statistics (ECNT) and the ranked probability score (RPS) line types are computed directly from the ensemble values. Second, for each categorical threshold specified, a fractional coverage value is computed as the ratio of the nearby forecast values that meet the threshold criteria. Point-Stat evaluates those fractional coverage values as if they were a probability forecast. When applying HiRA, users should enable the matched pair (MPR), probabilistic (PCT, PSTD, PJC, or PRC), continuous ensemble statistics (ECNT), or ranked probability score (RPS) line types in the **output_flag** dictionary. The number of probabilistic HiRA output lines is determined by the number of categorical forecast thresholds and HiRA neighborhood widths chosen.
+The Point-Stat tool has been enhanced to include the High Resolution Assessment (HiRA) verification logic (:ref:`Mittermaier, 2014 <Mittermaier-2014>`). HiRA is analogous to neighborhood verification but for point observations. The HiRA logic interprets the forecast values surrounding each point observation as an ensemble forecast. These ensemble values are processed in three ways. First, the ensemble continuous statistics (ECNT), the observation rank statistics (ORANK) and the ranked probability score (RPS) line types are computed directly from the ensemble values. Second, for each categorical threshold specified, a fractional coverage value is computed as the ratio of the nearby forecast values that meet the threshold criteria. Point-Stat evaluates those fractional coverage values as if they were a probability forecast. When applying HiRA, users should enable the matched pair (MPR), probabilistic (PCT, PSTD, PJC, or PRC), continuous ensemble statistics (ECNT), observation rank statistics (ORANK) or ranked probability score (RPS) line types in the **output_flag** dictionary. The number of probabilistic HiRA output lines is determined by the number of categorical forecast thresholds and HiRA neighborhood widths chosen.
 
 The HiRA framework provides a unique method for evaluating models in the neighborhood of point observations, allowing for some spatial and temporal uncertainty in the forecast and/or the observations. Additionally, the HiRA framework can be used to compare deterministic forecasts to ensemble forecasts. In MET, the neighborhood is a circle or square centered on the grid point closest to the observation location. An event is defined, then the proportion of points with events in the neighborhood is calculated. This proportion is treated as an ensemble probability, though it is likely to be uncalibrated.
 
@@ -425,6 +425,7 @@ ________________________
      pjc    = BOTH;
      prc    = BOTH;
      ecnt   = BOTH;  // Only for HiRA
+     orank  = BOTH;  // Only for HiRA
      rps    = BOTH;  // Only for HiRA
      eclv   = BOTH;
      mpr    = BOTH;
@@ -450,26 +451,29 @@ The **output_flag** array controls the type of output that the Point-Stat tool g
 
 9. **VL1L2** for Vector L1L2 Partial Sums
 
-10. **VCNT** for Vector Continuous Statistics (Note that bootstrap confidence intervals are not currently calculated for this line type.)
+10. **VAL1L2** for Vector Anomaly L1L2 Partial Sums when climatological data is supplied
 
-11. **VAL1L2** for Vector Anomaly L1L2 Partial Sums when climatological data is supplied
+11. **PCT** for Contingency Table counts for Probabilistic forecasts
 
-12. **PCT** for Contingency Table counts for Probabilistic forecasts
+12. **PSTD** for contingency table Statistics for Probabilistic forecasts with Dichotomous outcomes
 
-13. **PSTD** for contingency table Statistics for Probabilistic forecasts with Dichotomous outcomes
+13. **PJC** for Joint and Conditional factorization for Probabilistic forecasts
 
-14. **PJC** for Joint and Conditional factorization for Probabilistic forecasts
+14. **PRC** for Receiver Operating Characteristic for Probabilistic forecasts
 
-15. **PRC** for Receiver Operating Characteristic for Probabilistic forecasts
+15. **ECNT** for Ensemble Continuous Statistics is only computed for the HiRA methodology
 
-16. **ECNT** for Ensemble Continuous Statistics is only computed for the HiRA methodology
+16. **ORANK** for Ensemble Matched Pair Information when point observations are supplied for the HiRA methodology
 
 17. **RPS** for Ranked Probability Score is only computed for the HiRA methodology
 
 18. **ECLV** for Economic Cost/Loss Relative Value
 
 19. **MPR** for Matched Pair data
 
+20. **VCNT** for Vector Continuous Statistics (Note that bootstrap confidence intervals are not currently calculated for this line type.)
+
+
 Note that the first two line types are easily derived from each other. Users are free to choose which measures are most desired. The output line types are described in more detail in :numref:`point_stat-output`.
 
 Note that writing out matched pair data (MPR lines) for a large number of cases is generally not recommended. The MPR lines create very large output files and are only intended for use on a small set of cases.
@@ -489,9 +493,9 @@ point_stat_PREFIX_HHMMSSL_YYYYMMDD_HHMMSSV.stat where PREFIX indicates the user-
 
 The output ASCII files are named similarly:
 
-point_stat_PREFIX_HHMMSSL_YYYYMMDD_HHMMSSV_TYPE.txt where TYPE is one of mpr, fho, ctc, cts, cnt, mctc, mcts, pct, pstd, pjc, prc, ecnt, rps, eclv, sl1l2, sal1l2, vl1l2, vcnt or val1l2 to indicate the line type it contains.
+point_stat_PREFIX_HHMMSSL_YYYYMMDD_HHMMSSV_TYPE.txt where TYPE is one of mpr, fho, ctc, cts, cnt, mctc, mcts, pct, pstd, pjc, prc, ecnt, orank, rps, eclv, sl1l2, sal1l2, vl1l2, vcnt or val1l2 to indicate the line type it contains.
 
-The first set of header columns are common to all of the output files generated by the Point-Stat tool. Tables describing the contents of the header columns and the contents of the additional columns for each line type are listed in the following tables. The ECNT line type is described in :numref:`table_ES_header_info_es_out_ECNT`. The RPS line type is described in :numref:`table_ES_header_info_es_out_RPS`.
+The first set of header columns are common to all of the output files generated by the Point-Stat tool. Tables describing the contents of the header columns and the contents of the additional columns for each line type are listed in the following tables. The ECNT line type is described in :numref:`table_ES_header_info_es_out_ECNT`. The ORANK line type is described in :numref:`table_ES_header_info_es_out_ORANK`. The RPS line type is described in :numref:`table_ES_header_info_es_out_RPS`.
 
 .. _table_PS_header_info_point-stat_out:
 

met/src/tools/core/point_stat/point_stat.cc

@@ -327,9 +327,11 @@ void setup_first_pass(const DataPlane &dp, const Grid &data_grid) {
 ////////////////////////////////////////////////////////////////////////
 
 void setup_txt_files() {
-   int i, max_col, max_prob_col, max_mctc_col, n_prob, n_cat, n_eclv;
+   int i, j;
+   int max_col, max_prob_col, max_mctc_col, max_orank_col;
+   int n_prob, n_cat, n_eclv, n_ens;
    ConcatString base_name;
-
+   
    // Create output file names for the stat file and optional text files
    build_outfile_name(fcst_valid_ut, fcst_lead_sec, "", base_name);
 
@@ -340,23 +342,20 @@ void setup_txt_files() {
    /////////////////////////////////////////////////////////////////////
 
    // Get the maximum number of data columns
-   n_prob = conf_info.get_max_n_fprob_thresh();
+   n_prob = max(conf_info.get_max_n_fprob_thresh(),
+                conf_info.get_max_n_hira_prob());
    n_cat  = conf_info.get_max_n_cat_thresh() + 1;
    n_eclv = conf_info.get_max_n_eclv_points();
+   n_ens  = conf_info.get_max_n_hira_ens();
 
-   // Check for HiRA output
-   for(i=0; i<conf_info.get_n_vx(); i++) {
-      if(conf_info.vx_opt[i].hira_info.flag) {
-         n_prob = max(n_prob, conf_info.vx_opt[i].hira_info.cov_ta.n());
-      }
-   }
-
-   max_prob_col = get_n_pjc_columns(n_prob);
-   max_mctc_col = get_n_mctc_columns(n_cat);
+   max_prob_col  = get_n_pjc_columns(n_prob);
+   max_mctc_col  = get_n_mctc_columns(n_cat);
+   max_orank_col = get_n_orank_columns(n_ens);
 
    // Determine the maximum number of data columns
-   max_col = ( max_prob_col > max_stat_col ? max_prob_col : max_stat_col );
-   max_col = ( max_mctc_col > max_col      ? max_mctc_col : max_col );
+   max_col = (max_prob_col  > max_stat_col ? max_prob_col  : max_stat_col);
+   max_col = (max_mctc_col  > max_col      ? max_mctc_col  : max_col);
+   max_col = (max_orank_col > max_col      ? max_orank_col : max_col);
 
    // Add the header columns
    max_col += n_header_columns + 1;
@@ -429,12 +428,16 @@ void setup_txt_files() {
                max_col = get_n_eclv_columns(n_eclv) + n_header_columns + 1;
                break;
 
+            case(i_orank):
+               max_col = get_n_orank_columns(n_ens) + n_header_columns + 1;
+               break;
+
             default:
                max_col = n_txt_columns[i] + n_header_columns + 1;
                break;
          } // end switch
 
-         // Setup the text AsciiTable
+	 // Setup the text AsciiTable
          txt_at[i].set_size(conf_info.n_txt_row(i) + 1, max_col);
          setup_table(txt_at[i]);
 
@@ -465,6 +468,10 @@ void setup_txt_files() {
                write_eclv_header_row(1, n_eclv, txt_at[i], 0, 0);
                break;
 
+            case(i_orank):
+               write_orank_header_row(1, n_ens, txt_at[i], 0, 0);
+               break;
+
             default:
                write_header_row(txt_columns[i], n_txt_columns[i], 1,
                   txt_at[i], 0, 0);
@@ -1734,11 +1741,11 @@ void do_hira_ens(int i_vx, const PairDataPoint *pd_ptr) {
          if(gt) { delete gt; gt = 0; }
          continue;
       }
-
+      
       // Write out the ECNT line
       if(conf_info.vx_opt[i_vx].output_flag[i_ecnt] != STATOutputType_None) {
-
-         // Compute ensemble statistics
+	 
+	 //Compute ensemble statistics
          hira_pd.compute_pair_vals(rng_ptr);
          ECNTInfo ecnt_info;
          ecnt_info.set(hira_pd);
@@ -1749,6 +1756,18 @@ void do_hira_ens(int i_vx, const PairDataPoint *pd_ptr) {
             txt_at[i_ecnt], i_txt_row[i_ecnt]);
       } // end if ECNT
 
+      // Write out the ORANK line
+      if(conf_info.vx_opt[i_vx].output_flag[i_orank] != STATOutputType_None) {
+
+	 // Compute ensemble statistics
+         hira_pd.compute_pair_vals(rng_ptr);
+
+	 write_orank_row(shc, &hira_pd,
+	    conf_info.vx_opt[i_vx].output_flag[i_orank],
+	    stat_at, i_stat_row,
+	    txt_at[i_orank], i_txt_row[i_orank]);
+      } // end if ORANK
+
       // Write out the RPS line
       if(conf_info.vx_opt[i_vx].output_flag[i_rps] != STATOutputType_None) {
 
@@ -1925,12 +1944,13 @@ void do_hira_prob(int i_vx, const PairDataPoint *pd_ptr) {
 
          // Write out the MPR lines
          if(conf_info.vx_opt[i_vx].output_flag[i_mpr] != STATOutputType_None) {
-            write_mpr_row(shc, &hira_pd,
+
+	    write_mpr_row(shc, &hira_pd,
                conf_info.vx_opt[i_vx].output_flag[i_mpr],
                stat_at, i_stat_row,
                txt_at[i_mpr], i_txt_row[i_mpr], false);
 
-            // Reset the observation valid time
+	    // Reset the observation valid time
             shc.set_obs_valid_beg(conf_info.vx_opt[i_vx].vx_pd.beg_ut);
             shc.set_obs_valid_end(conf_info.vx_opt[i_vx].vx_pd.end_ut);
          }
@@ -1940,35 +1960,39 @@ void do_hira_prob(int i_vx, const PairDataPoint *pd_ptr) {
 
          // Write out PCT
          if(conf_info.vx_opt[i_vx].output_flag[i_pct] != STATOutputType_None) {
-            write_pct_row(shc, pct_info,
+
+	    write_pct_row(shc, pct_info,
                conf_info.vx_opt[i_vx].output_flag[i_pct],1, 1,
                stat_at, i_stat_row,
                txt_at[i_pct], i_txt_row[i_pct], false);
          }
 
          // Write out PSTD
          if(conf_info.vx_opt[i_vx].output_flag[i_pstd] != STATOutputType_None) {
-            write_pstd_row(shc, pct_info,
+
+	    write_pstd_row(shc, pct_info,
                conf_info.vx_opt[i_vx].output_flag[i_pstd], 1, 1,
                stat_at, i_stat_row,
                txt_at[i_pstd], i_txt_row[i_pstd], false);
          }
 
          // Write out PJC
          if(conf_info.vx_opt[i_vx].output_flag[i_pjc] != STATOutputType_None) {
-            write_pjc_row(shc, pct_info,
+
+	    write_pjc_row(shc, pct_info,
                conf_info.vx_opt[i_vx].output_flag[i_pjc], 1, 1,
                stat_at, i_stat_row,
                txt_at[i_pjc], i_txt_row[i_pjc], false);
          }
-
+	 
          // Write out PRC
          if(conf_info.vx_opt[i_vx].output_flag[i_prc] != STATOutputType_None) {
-            write_prc_row(shc, pct_info,
+
+	    write_prc_row(shc, pct_info,
                conf_info.vx_opt[i_vx].output_flag[i_prc], 1, 1,
                stat_at, i_stat_row,
                txt_at[i_prc], i_txt_row[i_prc], false);
-         }
+	 }
 
       } // end for j
    } // end for i

met/src/tools/core/point_stat/point_stat.h

@@ -70,8 +70,8 @@ static const char **txt_columns[n_txt] = {
    sl1l2_columns,  sal1l2_columns, vl1l2_columns,
    val1l2_columns, pct_columns,    pstd_columns,
    pjc_columns,    prc_columns,    ecnt_columns,
-   rps_columns,    eclv_columns,   mpr_columns,
-   vcnt_columns
+   orank_columns,  rps_columns,    eclv_columns,
+   mpr_columns,    vcnt_columns
 };
 
 // Length of header columns
@@ -81,8 +81,8 @@ static const int n_txt_columns[n_txt] = {
    n_sl1l2_columns,  n_sal1l2_columns, n_vl1l2_columns,
    n_val1l2_columns, n_pct_columns,    n_pstd_columns,
    n_pjc_columns,    n_prc_columns,    n_ecnt_columns,
-   n_rps_columns,    n_eclv_columns,   n_mpr_columns,
-   n_vcnt_columns
+   n_orank_columns,  n_rps_columns,    n_eclv_columns,
+   n_mpr_columns,    n_vcnt_columns
 };
 
 // Text file abbreviations
@@ -92,8 +92,8 @@ static const char *txt_file_abbr[n_txt] = {
    "sl1l2",  "sal1l2", "vl1l2",
    "val1l2", "pct",    "pstd",
    "pjc",    "prc",    "ecnt",
-   "rps",    "eclv",   "mpr",
-   "vcnt"
+   "orank",  "rps",    "eclv",
+   "mpr",    "vcnt"
 };
 
 ////////////////////////////////////////////////////////////////////////

met/src/tools/core/point_stat/point_stat_conf_info.cc

@@ -523,6 +523,26 @@ int PointStatConfInfo::get_max_n_eclv_points() const {
 
 ////////////////////////////////////////////////////////////////////////
 
+int PointStatConfInfo::get_max_n_hira_ens() const {
+   int i, n;
+
+   for(i=0,n=0; i<n_vx; i++) n = max(n, vx_opt[i].get_n_hira_ens());
+
+   return(n);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+int PointStatConfInfo::get_max_n_hira_prob() const {
+   int i, n;
+
+   for(i=0,n=0; i<n_vx; i++) n = max(n, vx_opt[i].get_n_hira_prob());
+
+   return(n);
+}
+
+////////////////////////////////////////////////////////////////////////
+
 bool PointStatConfInfo::get_vflag() const {
    int i;
    bool vflag = false;
@@ -744,7 +764,7 @@ void PointStatVxOpt::process_config(GrdFileType ftype,
 
    // Conf: output_flag
    output_map = parse_conf_output_flag(&odict, txt_file_type, n_txt);
-
+   
    // Populate the output_flag array with map values
    for(i=0; i<n_txt; i++) output_flag[i] = output_map[txt_file_type[i]];
 
@@ -1194,7 +1214,7 @@ int PointStatVxOpt::n_txt_row(int i_txt_row) const {
 
       case(i_ecnt):
       case(i_rps):
-         // Number of HiRA ECNT and RPS lines =
+         // Number of HiRA ECNT, ORANK and RPS lines =
          //    Message Types * Masks * Interpolations * HiRA widths *
          //    Alphas
          if(hira_info.flag) {
@@ -1206,6 +1226,20 @@ int PointStatVxOpt::n_txt_row(int i_txt_row) const {
 
          break;
 
+      case(i_orank):
+         // Number HiRA ORANK lines possible =
+         //    Number of pairs * Categorical Thresholds *
+         //    HiRA widths
+         if(hira_info.flag) {
+            n = vx_pd.get_n_pair() * get_n_cat_thresh() *
+                hira_info.width.n();
+         }
+         else {
+            n = 0;
+         }
+
+         break;
+
       case(i_eclv):
          // Number of CTC -> ECLV lines =
          //    Message Types * Masks * Interpolations * Thresholds *
@@ -1282,3 +1316,16 @@ int PointStatVxOpt::get_n_oprob_thresh() const {
 }
 
 ////////////////////////////////////////////////////////////////////////
+
+int PointStatVxOpt::get_n_hira_ens() const {
+   int n = (hira_info.flag ? hira_info.width.max() : 0);
+   return(n*n);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+int PointStatVxOpt::get_n_hira_prob() const {
+   return(hira_info.flag ? hira_info.cov_ta.n() : 0);
+}
+
+////////////////////////////////////////////////////////////////////////