Feature #2339 stat_analysis_seeps

docs/Users_Guide/grid-stat.rst

@@ -412,7 +412,7 @@ The **output_flag** array controls the type of output that the Grid-Stat tool ge
 
 21. **DMAP** for Distance Map Statistics
 
-22. **SEEPS** for SEEPS (Stable Equitable Error in Probability Space) score. It's described in :numref:`_table_PS_format_info_SEEPS`. The SEEPS score of matched pair data is saved into the NetCDF.
+22. **SEEPS** for SEEPS (Stable Equitable Error in Probability Space) score. It's described in :numref:`table_PS_format_info_SEEPS`. The SEEPS score of matched pair data is saved into the NetCDF.
 
 
 Note that the first two line types are easily derived from one another. The user is free to choose which measure is most desired. The output line types are described in more detail in :numref:`grid_stat-output`.

docs/Users_Guide/stat-analysis.rst

@@ -40,12 +40,12 @@ The **-derive** job command option can be used to perform the derivation of stat
 Aggregated values from multiple STAT lines
 ------------------------------------------
 
-The Stat-Analysis "aggregate" job aggregates values from multiple STAT lines of the same type. The user may specify the specific line type of interest and any other relevant search criteria. The Stat-Analysis tool then creates sums of each of the values in all lines matching the search criteria. The aggregated data are output as the same line type as the user specified. The STAT line types which may be aggregated in this way are the contingency table (FHO, CTC, PCT, MCTC, NBRCTC), partial sums (SL1L2, SAL1L2, VL1L2, and VAL1L2), and other (ISC, ECNT, RPS, RHIST, PHIST, RELP, NBRCNT, SSVAR, and GRAD) line types.
+The Stat-Analysis "aggregate" job aggregates values from multiple STAT lines of the same type. The user may specify the specific line type of interest and any other relevant search criteria. The Stat-Analysis tool then creates sums of each of the values in all lines matching the search criteria. The aggregated data are output as the same line type as the user specified. The STAT line types which may be aggregated in this way are the contingency table (FHO, CTC, PCT, MCTC, NBRCTC), partial sums (SL1L2, SAL1L2, VL1L2, and VAL1L2), and other (ISC, ECNT, RPS, RHIST, PHIST, RELP, NBRCNT, SSVAR, GRAD, and SEEPS) line types.
 
 Aggregate STAT lines and produce aggregated statistics
 ------------------------------------------------------
 
-The Stat-Analysis "aggregate-stat" job aggregates multiple STAT lines of the same type together and produces relevant statistics from the aggregated line. This may be done in the same manner listed above in :numref:`StA_Aggregated-values-from`. However, rather than writing out the aggregated STAT line itself, the relevant statistics generated from that aggregated line are provided in the output. Specifically, if a contingency table line type (FHO, CTC, PCT, MCTC, or NBRCTC) has been aggregated, contingency table statistics (CTS, ECLV, PSTD, MCTS, or NBRCTS) line types can be computed. If a partial sums line type (SL1L2 or SAL1L2) has been aggregated, the continuous statistics (CNT) line type can be computed. If a vector partial sums line type (VL1L2 or VAL1L2) has been aggregated, the vector continuous statistics (VCNT) line type can be computed. For ensembles, the ORANK line type can be accumulated into ECNT, RPS, RHIST, PHIST, RELP, or SSVAR output. If the matched pair line type (MPR) has been aggregated, may output line types (FHO, CTC, CTS, CNT, MCTC, MCTS, SL1L2, SAL1L2, VL1L2, VCNT, WDIR, PCT, PSTD, PJC, PRC, or ECLV) can be computed. Multiple output line types may be specified for each "aggregate-stat" job, as long as each output is derivable from the input.
+The Stat-Analysis "aggregate-stat" job aggregates multiple STAT lines of the same type together and produces relevant statistics from the aggregated line. This may be done in the same manner listed above in :numref:`StA_Aggregated-values-from`. However, rather than writing out the aggregated STAT line itself, the relevant statistics generated from that aggregated line are provided in the output. Specifically, if a contingency table line type (FHO, CTC, PCT, MCTC, or NBRCTC) has been aggregated, contingency table statistics (CTS, ECLV, PSTD, MCTS, or NBRCTS) line types can be computed. If a partial sums line type (SL1L2 or SAL1L2) has been aggregated, the continuous statistics (CNT) line type can be computed. If a vector partial sums line type (VL1L2 or VAL1L2) has been aggregated, the vector continuous statistics (VCNT) line type can be computed. For ensembles, the ORANK line type can be accumulated into ECNT, RPS, RHIST, PHIST, RELP, or SSVAR output. If a SEEPS matched pair line type (SEEPS_MPR) has been aggregated, the aggregated SEEPS line type (SEEPS) can be computed. If the matched pair line type (MPR) has been aggregated, may output line types (FHO, CTC, CTS, CNT, MCTC, MCTS, SL1L2, SAL1L2, VL1L2, VCNT, WDIR, PCT, PSTD, PJC, PRC, or ECLV) can be computed. Multiple output line types may be specified for each "aggregate-stat" job, as long as each output is derivable from the input.
 
 When aggregating the matched pair line type (MPR), additional required job command options are determined by the requested output line type(s). For example, the "-out_thresh" (or "-out_fcst_thresh" and  "-out_obs_thresh" options) are required to compute contingnecy table counts (FHO, CTC) or statistics (CTS). Those same job command options can also specify filtering thresholds when computing continuous partial sums (SL1L2, SAL1L2) or statistics (CNT). Output is written for each threshold specified.
 

internal/test_unit/config/STATAnalysisConfig_seeps

@@ -0,0 +1,118 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// STAT-Analysis configuration file.
+//
+// For additional information, please see the MET User's Guide.
+//
+////////////////////////////////////////////////////////////////////////////////
+
+//
+// Filtering input STAT lines by the contents of each column
+//
+model = [];
+desc  = [];
+
+fcst_lead = [];
+obs_lead  = [];
+
+fcst_valid_beg  = "";
+fcst_valid_end  = "";
+fcst_valid_inc  = [];
+fcst_valid_exc  = [];
+fcst_valid_hour = [];
+
+obs_valid_beg   = "";
+obs_valid_end   = "";
+obs_valid_inc   = [];
+obs_valid_exc   = [];
+obs_valid_hour  = [];
+
+fcst_init_beg   = "";
+fcst_init_end   = "";
+fcst_init_inc   = [];
+fcst_init_exc   = [];
+fcst_init_hour  = [];
+
+obs_init_beg    = "";
+obs_init_end    = "";
+obs_init_inc    = [];
+obs_init_exc    = [];
+obs_init_hour   = [];
+
+fcst_var = [];
+obs_var  = [];
+
+fcst_lev = [];
+obs_lev  = [];
+
+obtype = [];
+
+vx_mask = [];
+
+interp_mthd = [];
+
+interp_pnts = [];
+
+fcst_thresh = [];
+obs_thresh  = [];
+cov_thresh  = [];
+
+alpha = [];
+
+line_type = [];
+
+column = [];
+
+weight = [];
+
+////////////////////////////////////////////////////////////////////////////////
+
+//
+// Array of STAT-Analysis jobs to be performed on the filtered data
+//
+// Job 1 = filter SEEPS lines
+// Job 2 = aggregate SEEPS_MPR lines by interpolation (output equals Job 1)
+// Job 3 = aggregate SEEPS lines
+// Job 4 = aggregate all SEEPS_MPR lines (output equals Job 3)
+// Job 5 = summarize SEEPS_MPR scores
+//
+jobs = [
+   "-job filter         -line_type SEEPS     -dump_row ${OUTPUT_DIR}/CONFIG_POINT_STAT_filter_seeps.stat",
+   "-job aggregate_stat -line_type SEEPS_MPR -out_line_type SEEPS -by INTERP_MTHD,INTERP_PNTS",
+   "-job aggregate      -line_type SEEPS",
+   "-job aggregate_stat -line_type SEEPS_MPR -out_line_type SEEPS",
+   "-job summary        -line_type SEEPS_MPR -column SEEPS -by INTERP_MTHD,INTERP_PNTS"
+];
+
+////////////////////////////////////////////////////////////////////////////////
+
+//
+// Confidence interval settings
+//
+out_alpha = 0.05;
+
+boot = {
+   interval = PCTILE;
+   rep_prop = 1.0;
+   n_rep    = 1000;
+   rng      = "mt19937";
+   seed     = "1";
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+//
+// Skill score index options
+//
+ss_index_name       = "SS_INDEX";
+ss_index_vld_thresh = 1.0;
+
+////////////////////////////////////////////////////////////////////////////////
+
+hss_ec_value   = NA;
+rank_corr_flag = TRUE;
+vif_flag       = FALSE;
+tmp_dir        = "/tmp";
+version        = "V11.0.0";
+
+////////////////////////////////////////////////////////////////////////////////

internal/test_unit/xml/unit_stat_analysis_ps.xml

@@ -100,6 +100,23 @@
     </output>
   </test>
 
+  <test name="stat_analysis_POINT_STAT_SEEPS">
+    <env>
+      <pair><name>OUTPUT_DIR</name> <value>&OUTPUT_DIR;/stat_analysis_ps</value></pair>
+    </env>
+    <exec>&MET_BIN;/stat_analysis</exec>
+    <param> \
+      -lookin &OUTPUT_DIR;/point_stat/point_stat_NCMET_NAM_NDAS_SEEPS_360000L_20120410_120000V.stat \
+      -config &CONFIG_DIR;/STATAnalysisConfig_seeps \
+      -out &OUTPUT_DIR;/stat_analysis_ps/POINT_STAT_SEEPS.out \
+      -v 1
+    </param>
+    <output>
+      <stat>&OUTPUT_DIR;/stat_analysis_ps/CONFIG_POINT_STAT_filter_seeps.stat</stat>
+      <exist>&OUTPUT_DIR;/stat_analysis_ps/POINT_STAT_SEEPS.out</exist>
+    </output>
+  </test>
+
   <test name="stat_analysis_RAMPS">
     <env>
       <pair><name>OUTPUT_DIR</name> <value>&OUTPUT_DIR;/stat_analysis_ps</value></pair>

src/libcode/vx_seeps/seeps.cc

@@ -41,6 +41,8 @@ static const char *var_name_elv       = "elv";
 static const char *dim_name_lat       = "latitude";
 static const char *dim_name_lon       = "longitude";
 
+double weighted_average(double, double, double, double);
+
 ////////////////////////////////////////////////////////////////////////
 
 SeepsClimo *get_seeps_climo() {
@@ -549,6 +551,65 @@ void SeepsAggScore::clear() {
 ////////////////////////////////////////////////////////////////////////
 
 
+SeepsAggScore & SeepsAggScore::operator+=(const SeepsAggScore &c) {
+
+   // Check for degenerate case
+   if(n_obs == 0 && c.n_obs == 0) return(*this);
+
+   // Compute weights
+   double w1 = (double)   n_obs / (n_obs + c.n_obs);
+   double w2 = (double) c.n_obs / (n_obs + c.n_obs);
+
+   // Increment number of obs
+   n_obs += c.n_obs;
+
+   // Increment counts
+   c12 += c.c12;
+   c13 += c.c13;
+   c21 += c.c21;
+   c23 += c.c23;
+   c31 += c.c31;
+   c32 += c.c32;
+
+   // Compute weighted averages
+   s12 = weighted_average(s12, w1, c.s12, w2);
+   s13 = weighted_average(s13, w1, c.s13, w2);
+   s21 = weighted_average(s21, w1, c.s21, w2);
+   s23 = weighted_average(s23, w1, c.s23, w2);
+   s31 = weighted_average(s31, w1, c.s31, w2);
+   s32 = weighted_average(s32, w1, c.s32, w2);
+
+   pv1 = weighted_average(pv1, w1, c.pv1, w2);
+   pv2 = weighted_average(pv2, w1, c.pv2, w2);
+   pv3 = weighted_average(pv3, w1, c.pv3, w2);
+
+   pf1 = weighted_average(pf1, w1, c.pf1, w2);
+   pf2 = weighted_average(pf2, w1, c.pf2, w2);
+   pf3 = weighted_average(pf3, w1, c.pf3, w2);
+
+   mean_fcst = weighted_average(mean_fcst, w1, c.mean_fcst, w2);
+   mean_obs  = weighted_average(mean_obs,  w1, c.mean_obs,  w2);
+
+   score          = weighted_average(score,          w1, c.score,          w2);
+   weighted_score = weighted_average(weighted_score, w1, c.weighted_score, w2);
+
+   return(*this);
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+
+double weighted_average(double v1, double w1, double v2, double w2) {
+   return(is_bad_data(v1) || is_bad_data(v2) ?
+          bad_data_double :
+          v1 * w1 + v2 * w2);
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+
 SeepsClimoGrid::SeepsClimoGrid(int month, int hour) : month{month}, hour{hour}
 {
 

src/libcode/vx_seeps/seeps.h

@@ -73,7 +73,8 @@ struct SeepsScore { // For SEEPS_MPR
 
 struct SeepsAggScore {  // For SEEPS
    void clear();
-
+   SeepsAggScore & operator+=(const SeepsAggScore &);
+
    int   n_obs;
    int   c12;
    int   c13;

src/libcode/vx_statistics/compute_stats.cc

@@ -1700,19 +1700,23 @@ void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &ob
 // ; PV-WAVE prints: 1.00000      5.00000
 // PRINT, SUM(array2, 1)
 // ; PV-WAVE prints: 2.00000      4.00000
-
+////////////////////////////////////////////////////////////////////////
 
 double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps) {
    int seeps_idx;
    SeepsScore *seeps_mpr;
    int seeps_cnt = seeps->n_obs;
-   const double density_radidus_rad = density_radius * rad_per_deg;
-   double rlat[seeps_cnt], rlon[seeps_cnt];
-   double clat[seeps_cnt], clon[seeps_cnt], slat[seeps_cnt], slon[seeps_cnt];
-   double clat_m[seeps_cnt][seeps_cnt], clon_m[seeps_cnt][seeps_cnt];
-   double slat_m[seeps_cnt][seeps_cnt], slon_m[seeps_cnt][seeps_cnt];
-   double clat_slon[seeps_cnt][seeps_cnt], clon_slat[seeps_cnt][seeps_cnt];
-   double density_m[seeps_cnt][seeps_cnt];
+   const double density_radius_rad = density_radius * rad_per_deg;
+   vector<double> rlat(seeps_cnt), rlon(seeps_cnt);
+   vector<double> clat(seeps_cnt), clon(seeps_cnt);
+   vector<double> slat(seeps_cnt), slon(seeps_cnt);
+   vector< vector<double> > clat_m(seeps_cnt, vector<double> (seeps_cnt));
+   vector< vector<double> > clon_m(seeps_cnt, vector<double> (seeps_cnt));
+   vector< vector<double> > slat_m(seeps_cnt, vector<double> (seeps_cnt));
+   vector< vector<double> > slon_m(seeps_cnt, vector<double> (seeps_cnt));
+   vector< vector<double> > clat_slon(seeps_cnt, vector<double> (seeps_cnt));
+   vector< vector<double> > clon_slat(seeps_cnt, vector<double> (seeps_cnt));
+   vector< vector<double> > density_m(seeps_cnt, vector<double> (seeps_cnt));
    static const char *method_name = "compute_seeps_density_vector() -> ";
 
    if (seeps_cnt == 0) {
@@ -1767,11 +1771,11 @@ double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *see
          //IDL: r = acos(r)
          density = acos(density);
          //IDL: if r0 gt 0.0 then r = exp(-(r/r0)^2) * (r le 4. * r0) else r = (r*0.)+1.
-         if (density_radidus_rad <= 0.) density = 1.0;
+         if (density_radius_rad <= 0.) density = 1.0;
          else {
             mask3 = (density <= 4.0) ? 1. : 0.;
-            temp = density / density_radidus_rad;
-            density = exp(-(temp * temp)) * mask3 * density_radidus_rad;
+            temp = density / density_radius_rad;
+            density = exp(-(temp * temp)) * mask3 * density_radius_rad;
          }
          density_vector[j] += density;
       }