Per #1451, switch from constant pointer to ClimoCDFInfo object to a c…

…opy to make the logic of doing this in Stat-Analysis a little easier.
dtcenter · Feb 19, 2021 · 2505d42 · 2505d42
1 parent c0db18b
commit 2505d42
Show file tree

Hide file tree

Showing 8 changed files with 50 additions and 36 deletions.
diff --git a/met/src/libcode/vx_analysis_util/stat_job.cc b/met/src/libcode/vx_analysis_util/stat_job.cc
@@ -2439,7 +2439,8 @@ ConcatString STATAnalysisJob::get_jobstring() const {
    // Jobs which use out_bin_size
    if(line_type.n_elements() > 0) {
       if(string_to_statlinetype(line_type[0].c_str()) == stat_orank &&
-         out_line_type.has(stat_phist_str)) {
+         (out_line_type.has(stat_phist_str) ||
+          out_line_type.has(stat_ecnt_str))) {
 
          // out_bin_size
         js << "-out_bin_size " << out_bin_size << " ";

diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.cc b/met/src/libcode/vx_statistics/pair_data_ensemble.cc
@@ -92,7 +92,7 @@ void PairDataEnsemble::clear() {
    obs_error_entry.clear();
    obs_error_flag = false;
 
-   cdf_info = (const ClimoCDFInfo *) 0;
+   cdf_info.clear();
 
    for(i=0; i<n_ens; i++) e_na[i].clear();
    if(e_na) { delete [] e_na; e_na = (NumArray *) 0; }
@@ -316,9 +316,9 @@ void PairDataEnsemble::set_ens_size(int n) {
 
 ////////////////////////////////////////////////////////////////////////
 
-void PairDataEnsemble::set_climo_cdf(const ClimoCDFInfo *ptr) {
+void PairDataEnsemble::set_climo_cdf(const ClimoCDFInfo &info) {
 
-   cdf_info = ptr;
+   cdf_info = info;
 
    return;
 }
@@ -452,7 +452,7 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
          }
 
          // Derive ensemble from climo mean and standard deviation
-         cur_clm = derive_normal_cdf_inv(cdf_info, cmn_na[i], csd_na[i]);
+         cur_clm = derive_climo_cdf_inv(cdf_info, cmn_na[i], csd_na[i]);
 
          // Store empirical CRPS stats
          crps_emp_na.add(compute_crps_emp(o_na[i], cur_ens));
@@ -1274,12 +1274,12 @@ void VxPairDataEnsemble::set_ens_size(int n) {
 
 ////////////////////////////////////////////////////////////////////////
 
-void VxPairDataEnsemble::set_climo_cdf(const ClimoCDFInfo *ptr) {
+void VxPairDataEnsemble::set_climo_cdf(const ClimoCDFInfo &info) {
 
    for(int i=0; i<n_msg_typ; i++) {
       for(int j=0; j<n_mask; j++) {
          for(int k=0; k<n_interp; k++) {
-            pd[i][j][k].set_climo_cdf(ptr);
+            pd[i][j][k].set_climo_cdf(info);
          }
       }
    }
@@ -1884,18 +1884,15 @@ double compute_ens_pit(double obs, double m, double s) {
 
 ////////////////////////////////////////////////////////////////////////
 
-NumArray derive_normal_cdf_inv(const ClimoCDFInfo *cdf_info, double m, double s) {
+NumArray derive_climo_cdf_inv(const ClimoCDFInfo &cdf_info, double m, double s) {
    NumArray cdf_inv_na;
 
-   // Check pointer
-   if(!cdf_info) return(cdf_inv_na);
-
    // Check for bad data
    if(is_bad_data(m) || is_bad_data(s)) return(cdf_inv_na);
 
    // Skip the first (>=0.0) and last (>=1.0) climo CDF thresholds
-   for(int i=1; i<cdf_info->cdf_ta.n()-1; i++) {
-      cdf_inv_na.add(normal_cdf_inv(cdf_info->cdf_ta[i].get_value(), m, s));
+   for(int i=1; i<cdf_info.cdf_ta.n()-1; i++) {
+      cdf_inv_na.add(normal_cdf_inv(cdf_info.cdf_ta[i].get_value(), m, s));
    }
 
    return(cdf_inv_na);

diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.h b/met/src/libcode/vx_statistics/pair_data_ensemble.h
@@ -74,8 +74,8 @@ class PairDataEnsemble : public PairBase {
       ObsErrorEntryPtrArray obs_error_entry;
       bool                  obs_error_flag;
 
-      // ClimoCDFInfo pointer
-      const ClimoCDFInfo *cdf_info;
+      // Climatology info
+      ClimoCDFInfo cdf_info;
 
       // Ensemble, valid count, and rank values
       NumArray  *e_na;            // Ensemble values [n_ens][n_obs]
@@ -135,7 +135,7 @@ class PairDataEnsemble : public PairBase {
       void add_ens(int, double);
       void add_ens_var_sums(int, double);
       void set_ens_size(int);
-      void set_climo_cdf(const ClimoCDFInfo *);
+      void set_climo_cdf(const ClimoCDFInfo &);
 
       void add_obs_error_entry(ObsErrorEntry *);
 
@@ -272,7 +272,7 @@ class VxPairDataEnsemble {
       // Call set_ens_size before add_ens
       void set_ens_size(int n);
 
-      void set_climo_cdf(const ClimoCDFInfo *);
+      void set_climo_cdf(const ClimoCDFInfo &);
 
       void set_ssvar_bin_size(double);
       void set_phist_bin_size(double);
@@ -309,7 +309,7 @@ extern double compute_crps_gaus(double, double, double);
 extern double compute_ens_ign(double, double, double);
 extern double compute_ens_pit(double, double, double);
 
-extern NumArray derive_normal_cdf_inv(const ClimoCDFInfo *, double, double);
+extern NumArray derive_climo_cdf_inv(const ClimoCDFInfo &, double, double);
 
 ////////////////////////////////////////////////////////////////////////
 

diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc
@@ -601,7 +601,7 @@ void EnsembleStatConfInfo::set_vx_pd(const IntArray &ens_size) {
       vx_opt[i].vx_pd.set_ens_size(ens_size[i]);
 
       // Store the ClimoCDFInfo
-      vx_opt[i].vx_pd.set_climo_cdf(&vx_opt[i].cdf_info);
+      vx_opt[i].vx_pd.set_climo_cdf(vx_opt[i].cdf_info);
 
    }
 

diff --git a/met/src/tools/core/stat_analysis/aggr_stat_line.cc b/met/src/tools/core/stat_analysis/aggr_stat_line.cc
@@ -2553,8 +2553,6 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job,
          //
          if(m.count(key) == 0) {
             aggr.ens_pd.clear();
-            aggr.crpscl_emp_na.clear();
-            aggr.crpscl_gaus_na.clear();
             aggr.hdr.clear();
             m[key] = aggr;
          }
@@ -2573,9 +2571,9 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job,
          // Store the current statistics and weight (TOTAL column)
          //
          m[key].ens_pd.crps_emp_na.add(cur.crps_emp);
-         m[key].crpscl_emp_na.add(cur.crpscl_emp);
+         m[key].ens_pd.crpscl_emp_na.add(cur.crpscl_emp);
          m[key].ens_pd.crps_gaus_na.add(cur.crps_gaus);
-         m[key].crpscl_gaus_na.add(cur.crpscl_gaus);
+         m[key].ens_pd.crpscl_gaus_na.add(cur.crpscl_gaus);
          m[key].ens_pd.ign_na.add(cur.ign);
          m[key].ens_pd.var_na.add(square(cur.spread));
          m[key].ens_pd.var_oerr_na.add(square(cur.spread_oerr));
@@ -2620,9 +2618,9 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job,
       it->second.ens_pd.rmse_oerr = (is_bad_data(v) ? bad_data_double : sqrt(v));
 
       crps_emp    = it->second.ens_pd.crps_emp_na.wmean(it->second.ens_pd.wgt_na);
-      crpscl_emp  = it->second.crpscl_emp_na.wmean(it->second.ens_pd.wgt_na);
+      crpscl_emp  = it->second.ens_pd.crpscl_emp_na.wmean(it->second.ens_pd.wgt_na);
       crps_gaus   = it->second.ens_pd.crps_gaus_na.wmean(it->second.ens_pd.wgt_na);
-      crpscl_gaus = it->second.crpscl_gaus_na.wmean(it->second.ens_pd.wgt_na);
+      crpscl_gaus = it->second.ens_pd.crpscl_gaus_na.wmean(it->second.ens_pd.wgt_na);
 
       // Compute aggregated empirical CRPSS
       it->second.ens_pd.crpss_emp =
@@ -3002,8 +3000,9 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job,
    AggrENSInfo aggr;
    ORANKData cur;
    ConcatString key;
+   NumArray cur_clm;
    int i, n_valid, n_bin;
-   double esum, esumsq, crps_gaus, ign, pit;
+   double esum, esumsq;
    map<ConcatString, AggrENSInfo>::iterator it;
 
    //
@@ -3054,6 +3053,8 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job,
             n_bin = ceil(1.0/aggr.ens_pd.phist_bin_size);
             for(i=0; i<n_bin; i++) aggr.ens_pd.phist_na.add(0);
             aggr.ens_pd.ssvar_bin_size = job.out_bin_size;
+            bool center = false;
+            aggr.ens_pd.cdf_info.set_cdf_ta(nint(1.0/job.out_bin_size), center);
             aggr.hdr.clear();
             m[key] = aggr;
          }
@@ -3073,8 +3074,8 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job,
          // ensemble spread, ensemble member values, and
          // valid ensemble count
          //
-         m[key].ens_pd.add_grid_obs(cur.obs, cur.climo,
-                                    bad_data_double, default_grid_weight);
+         m[key].ens_pd.add_grid_obs(cur.obs, cur.climo_mean,
+                                    cur.climo_stdev, default_grid_weight);
          m[key].ens_pd.skip_ba.add(false);
          m[key].ens_pd.n_pair++;
          m[key].ens_pd.r_na.add(cur.rank);
@@ -3097,12 +3098,19 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job,
          m[key].ens_pd.esumsq_na.add(esumsq);
          m[key].ens_pd.v_na.add(n_valid);
 
-         // Compute Empirical and Gaussian CRPS, IGN, and PIT for the current point
+         // Derive ensemble from climo mean and standard deviation
+         cur_clm = derive_climo_cdf_inv(m[key].ens_pd.cdf_info,
+                                        cur.climo_mean, cur.climo_stdev);
+
+         // Store empirical CRPS stats
          m[key].ens_pd.crps_emp_na.add(compute_crps_emp(cur.obs, cur.ens_na));
-         compute_crps_gaus_ign_pit(cur.obs, cur.ens_na, crps_gaus, ign, pit);
-         m[key].ens_pd.crps_gaus_na.add(crps_gaus);
-         m[key].ens_pd.ign_na.add(ign);
-         m[key].ens_pd.pit_na.add(pit);
+         m[key].ens_pd.crpscl_emp_na.add(compute_crps_emp(cur.obs, cur_clm));
+
+         // Store Gaussian CRPS stats
+         m[key].ens_pd.crps_gaus_na.add(compute_crps_gaus(cur.obs, cur.ens_mean, cur.spread));
+         m[key].ens_pd.crpscl_gaus_na.add(compute_crps_gaus(cur.obs, cur.climo_mean, cur.climo_stdev));
+         m[key].ens_pd.ign_na.add(compute_ens_ign(cur.obs, cur.ens_mean, cur.spread));
+         m[key].ens_pd.pit_na.add(compute_ens_pit(cur.obs, cur.ens_mean, cur.spread));
 
          //
          // Increment the RHIST counts

diff --git a/met/src/tools/core/stat_analysis/aggr_stat_line.h b/met/src/tools/core/stat_analysis/aggr_stat_line.h
@@ -151,7 +151,6 @@ struct AggrISCInfo {
 struct AggrENSInfo {
    StatHdrInfo hdr;
    PairDataEnsemble ens_pd;
-   NumArray crpscl_emp_na, crpscl_gaus_na;
    NumArray me_na, mse_na, me_oerr_na, mse_oerr_na;
 };
 

diff --git a/met/src/tools/core/stat_analysis/parse_stat_line.cc b/met/src/tools/core/stat_analysis/parse_stat_line.cc
@@ -488,12 +488,21 @@ void parse_orank_line(STATLine &l, ORANKData &o_data) {
 
    o_data.obs_qc           =      l.get_item("OBS_QC", false);
    o_data.ens_mean         = atof(l.get_item("ENS_MEAN"));
-   o_data.climo            = atof(l.get_item("CLIMO"));
    o_data.spread           = atof(l.get_item("SPREAD"));
    o_data.ens_mean_oerr    = atof(l.get_item("ENS_MEAN_OERR"));
    o_data.spread_oerr      = atof(l.get_item("SPREAD_OERR"));
    o_data.spread_plus_oerr = atof(l.get_item("SPREAD_PLUS_OERR"));
 
+   // In met-10.0.0 and later, CLIMO column was replaced by CLIMO_MEAN
+   if(l.has("CLIMO")) {
+      o_data.climo_mean  = atof(l.get_item("CLIMO"));
+      o_data.climo_stdev = bad_data_double;
+   }
+   else {
+      o_data.climo_mean  = atof(l.get_item("CLIMO_MEAN"));
+      o_data.climo_stdev = atof(l.get_item("CLIMO_STDEV"));
+   }
+
    return;
 }
 

diff --git a/met/src/tools/core/stat_analysis/parse_stat_line.h b/met/src/tools/core/stat_analysis/parse_stat_line.h
@@ -93,7 +93,7 @@ struct ORANKData {
    int total, index;
    ConcatString obs_sid, obs_qc;
    double obs_lat, obs_lon, obs_lvl, obs_elv;
-   double obs, pit, climo;
+   double obs, pit, climo_mean, climo_stdev;
    double ens_mean, spread, ens_mean_oerr, spread_oerr;
    double spread_plus_oerr;
    int rank, n_ens_vld, n_ens;