Feature 1450 hersbach (#1662)

* Per #1450, add new ECNT columns for Hersback CRPS. Still need to actually compute the stats though. * Per #1450, update NumArray functions to only sort if the data is not yet sorted. And check for bad data when computing the standard deviation. * Per #1450, add code to compute the empirical CRPS value. * Per #1450, large change to the new output for the empirical CRPS. In order to aggregate decomposed empirical CRPS reliability and potential correctly, we'd need to write (n+1)*2 additional columns. While the empirical crps can be aggregated as a weighted mean, the decomposition cannot. It just isn't feasible to do this in the ECNT line type. If this reliability and potential really are required, recommend that we add an entirely new CRPS line type instead of tacking onto ECNT. These changes simply remove reliabilit and potential from the output. * Per #1450 and #1451, replacing single CRPS_CLIMO column with CRPSCL and CRPSCL_EMP which will be needed for #1451. * Per #1450, delete temp files I'd accidentally committed. * Per #1450, update the user's guide with CRPS updates. * Fix bug replacing crpss_emp with crpss_gaus.
dtcenter · Feb 17, 2021 · 11a5419 · 11a5419
1 parent 4c9ecde
commit 11a5419
Show file tree

Hide file tree

Showing 20 changed files with 247 additions and 547 deletions.
diff --git a/met/.cproject b/met/.cproject
diff --git a/met/.project b/met/.project
diff --git a/met/data/table_files/met_header_columns_V10.0.txt b/met/data/table_files/met_header_columns_V10.0.txt
@@ -17,7 +17,7 @@ V10.0 : STAT : PJC      : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID
 V10.0 : STAT : PRC      : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_THRESH) THRESH_[0-9]* PODY_[0-9]* POFD_[0-9]*
 V10.0 : STAT : PSTD     : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_THRESH) BASER BASER_NCL BASER_NCU RELIABILITY RESOLUTION UNCERTAINTY ROC_AUC BRIER BRIER_NCL BRIER_NCU BRIERCL BRIERCL_NCL BRIERCL_NCU BSS BSS_SMPL THRESH_[0-9]*
 V10.0 : STAT : ECLV     : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BASER VALUE_BASER (N_PTS) CL_[0-9]* VALUE_[0-9]*
-V10.0 : STAT : ECNT     : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR
+V10.0 : STAT : ECNT     : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP
 V10.0 : STAT : RPS      : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_PROB RPS_REL RPS_RES RPS_UNC RPS RPSS RPSS_SMPL RPS_COMP
 V10.0 : STAT : RHIST    : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_RANK) RANK_[0-9]*
 V10.0 : STAT : PHIST    : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BIN_SIZE (N_BIN) BIN_[0-9]*

diff --git a/met/docs/Users_Guide/appendixC.rst b/met/docs/Users_Guide/appendixC.rst
@@ -887,9 +887,9 @@ ________________________________________________
 CRPS
 ~~~~
 
-Called "CRPS" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`
+Called "CRPS", "CRPSCL", "CRPS_EMP", and "CRPSCL_EMP" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`
 
-The continuous ranked probability score (CRPS) is the integral, over all possible thresholds, of the Brier scores (:ref:`Gneiting et al., 2004 <Gneiting-2004>`). In MET, the CRPS calculation uses a normal distribution fit to the ensemble forecasts. In many cases, use of other distributions would be better.
+The continuous ranked probability score (CRPS) is the integral, over all possible thresholds, of the Brier scores (:ref:`Gneiting et al., 2004 <Gneiting-2004>`). In MET, the CRPS is calculated two ways: using a normal distribution fit to the ensemble forecasts (CRPS and CRPSCL), and using the empirical ensemble distribution (CRPS_EMP and CRPSCL_EMP). In some cases, use of other distributions would be better.
 
 WARNING: The normal distribution is probably a good fit for temperature and pressure, and possibly a not horrible fit for winds. However, the normal approximation will not work on most precipitation forecasts and may fail for many other atmospheric variables.
 
@@ -906,12 +906,14 @@ The score can be interpreted as a continuous version of the mean absolute error
 CRPS Skill Score
 ~~~~~~~~~~~~~~~~
 
-Called "CRPSS" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`
+Called "CRPSS" and "CRPSS_EMP" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`
 
 The continuous ranked probability skill score (CRPSS) is similar to the MSESS and the BSS, in that it compares its namesake score to that of a reference forecast to produce a positively oriented score between 0 and 1.
 
 .. math:: \text{CRPSS} = 1 - \frac{\text{CRPS}_{fcst}}{ \text{CRPS}_{ref}}
 
+For the normal distribution fit (CRPSS), the reference CRPS is computed using the climatological mean and standard deviation. For the empirical distribution (CRPSS_EMP), the reference CRPS is computed by sampling from the assumed normal climatological distribution defined by the mean and standard deviation.  
+
 IGN
 ~~~
 

diff --git a/met/docs/Users_Guide/ensemble-stat.rst b/met/docs/Users_Guide/ensemble-stat.rst
@@ -31,7 +31,7 @@ Often, the goal of ensemble forecasting is to reproduce the distribution of obse
 
 The relative position (RELP) is a count of the number of times each ensemble member is closest to the observation. For stochastic or randomly derived ensembles, this statistic is meaningless. For specified ensemble members, however, it can assist users in determining if any ensemble member is performing consistently better or worse than the others.
 
-The ranked probability score (RPS) is included in the Ranked Probability Score (RPS) line type. It is the mean of the Brier scores computed from ensemble probabilities derived for each probability category threshold (prob_cat_thresh) specified in the configuration file. The continuous ranked probability score (CRPS) is the average the distance between the forecast (ensemble) cumulative distribution function and the observation cumulative distribution function. It is an analog of the Brier score, but for continuous forecast and observation fields. (:ref:`Gneiting et al., 2004 <Gneiting-2004>`). The CRPS statistic is included in the Ensemble Continuous Statistics (ECNT) line type, along with other statistics quantifying the ensemble spread and ensemble mean skill.
+The ranked probability score (RPS) is included in the Ranked Probability Score (RPS) line type. It is the mean of the Brier scores computed from ensemble probabilities derived for each probability category threshold (prob_cat_thresh) specified in the configuration file. The continuous ranked probability score (CRPS) is the average the distance between the forecast (ensemble) cumulative distribution function and the observation cumulative distribution function. It is an analog of the Brier score, but for continuous forecast and observation fields. The CRPS statistic is computed using two methods: assuming a normal distribution defined by the ensemble mean and spread (:ref:`Gneiting et al., 2004 <Gneiting-2004>`) and using the empirical ensemble distribution (:ref:`Hersbach, 2000 <Hersbach-2000>`). The CRPS statistic is included in the Ensemble Continuous Statistics (ECNT) line type, along with other statistics quantifying the ensemble spread and ensemble mean skill.
 
 Ensemble observation error
 ~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -587,10 +587,10 @@ The format of the STAT and ASCII output of the Ensemble-Stat tool are described
     - Number of ensemble values
   * - 27
     - CRPS
-    - The Continuous Ranked Probability Score
+    - The Continuous Ranked Probability Score (normal distribution)
   * - 28
     - CRPSS
-    - The Continuous Ranked Probability Skill Score
+    - The Continuous Ranked Probability Skill Score (normal distribution)
   * - 29
     - IGN
     - The Ignorance Score
@@ -615,6 +615,18 @@ The format of the STAT and ASCII output of the Ensemble-Stat tool are described
   * - 36
     - SPREAD_PLUS_OERR
     - The square root of the sum of unperturbed ensemble variance and the observation error variance
+  * - 37 
+    - CRPSCL
+    - Climatological Continuous Ranked Probability Score (normal distribution)
+  * - 38
+    - CRPS_EMP 
+    - The Continuous Ranked Probability Score (empirical distribution)
+  * - 39
+    - CRPSCL_EMP 
+    - Climatological Continuous Ranked Probability Score (empirical distribution)
+  * - 40 
+    - CRPSS_EMP
+    - The Continuous Ranked Probability Skill Score (empirical distribution)
 
 .. _table_ES_header_info_es_out_RPS:
 

diff --git a/met/docs/Users_Guide/point-stat.rst b/met/docs/Users_Guide/point-stat.rst
@@ -170,7 +170,7 @@ When the "prob" entry is set as a dictionary to define the field of interest, se
 Measures for comparison against climatology
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-For each of the types of statistics mentioned above (categorical, continuous, and probabilistic), it is possible to calculate measures of skill relative to climatology. MET will accept a climatology file provided by the user, and will evaluate it as a reference forecast. Further, anomalies, i.e. departures from average conditions, can be calculated. As with all other statistics, the available measures will depend on the nature of the forecast. Common statistics that use a climatological reference include: the mean squared error skill score (MSESS), the Anomaly Correlation (ANOM_CORR and ANOM_CORR_UNCNTR), scalar and vector anomalies (SAL1L2 and VAL1L2), continuous ranked probability skill score (CRPSS), Brier Skill Score (BSS) (:ref:`Wilks, 2011 <Wilks-2011>`; :ref:`Mason, 2004 <Mason-2004>`).
+For each of the types of statistics mentioned above (categorical, continuous, and probabilistic), it is possible to calculate measures of skill relative to climatology. MET will accept a climatology file provided by the user, and will evaluate it as a reference forecast. Further, anomalies, i.e. departures from average conditions, can be calculated. As with all other statistics, the available measures will depend on the nature of the forecast. Common statistics that use a climatological reference include: the mean squared error skill score (MSESS), the Anomaly Correlation (ANOM_CORR and ANOM_CORR_UNCNTR), scalar and vector anomalies (SAL1L2 and VAL1L2), continuous ranked probability skill score (CRPSS and CRPSS_EMP), Brier Skill Score (BSS) (:ref:`Wilks, 2011 <Wilks-2011>`; :ref:`Mason, 2004 <Mason-2004>`).
 
 Often, the sample climatology is used as a reference by a skill score. The sample climatology is the average over all included observations and may be transparent to the user. This is the case in most categorical skill scores. The sample climatology will probably prove more difficult to improve upon than a long term climatology, since it will be from the same locations and time periods as the forecasts. This may mask legitimate forecast skill. However, a more general climatology, perhaps covering many years, is often easier to improve upon and is less likely to mask real forecast skill.
 

diff --git a/met/docs/Users_Guide/refs.rst b/met/docs/Users_Guide/refs.rst
@@ -132,6 +132,12 @@ References
 |   forecasts. *Monthly Weather Review*, 129, 550-560.
 | 
 
+.. _Hersbach-2000:
+
+| Hersbach, H., 2000: Decomposition of the Continuous Ranked Probability Score
+| for Ensemble Prediction Systems. *Weather and Forecasting*, 15, 559-570.
+| 
+
 .. _Jolliffe-2012:
 
 | Jolliffe, I.T., and D.B. Stephenson, 2012: *Forecast verification. A*

diff --git a/met/src/basic/vx_util/num_array.cc b/met/src/basic/vx_util/num_array.cc
@@ -540,7 +540,7 @@ void NumArray::sort_array()
 
 {
 
-sort(e, Nelements);
+if(!Sorted) sort(e, Nelements);
 
 Sorted = true;
 
@@ -784,7 +784,7 @@ double var;
 
 compute_mean_variance(mn, var);
 
-stdev = square_root(var);
+stdev = (is_bad_data(var) ? bad_data_double : square_root(var));
 
 return;
 

diff --git a/met/src/basic/vx_util/stat_column_defs.h b/met/src/basic/vx_util/stat_column_defs.h
@@ -260,7 +260,9 @@ static const char * ecnt_columns [] = {
    "TOTAL",       "N_ENS",       "CRPS",
    "CRPSS",       "IGN",         "ME",
    "RMSE",        "SPREAD",      "ME_OERR",
-   "RMSE_OERR",   "SPREAD_OERR", "SPREAD_PLUS_OERR"
+   "RMSE_OERR",   "SPREAD_OERR", "SPREAD_PLUS_OERR",
+   "CRPSCL",      "CRPS_EMP",    "CRPSCL_EMP",
+   "CRPSS_EMP"
 };
 
 static const char * rps_columns [] = {

diff --git a/met/src/libcode/vx_stat_out/stat_columns.cc b/met/src/libcode/vx_stat_out/stat_columns.cc
@@ -3659,11 +3659,12 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info,
    //
    // Ensemble Continuous Statistics
    // Dump out the ECNT line:
-   //    TOTAL,        N_ENS,
-   //    CRPS,         CRPSS,        IGN,
-   //    ME,           RMSE,         SPREAD,
-   //    ME_OERR,      RMSE_OERR,    SPREAD_OERR,
-   //    SPREAD_PLUS_OERR
+   //    TOTAL,            N_ENS,        CRPS,
+   //    CRPSS,            IGN,          ME,
+   //    RMSE,             SPREAD,       ME_OERR,
+   //    RMSE_OERR,        SPREAD_OERR,  SPREAD_PLUS_OERR,
+   //    CRPSCL,           CRPS_EMP,     CRPSCL_EMP,
+   //    CRPSS_EMP
    //
    at.set_entry(r, c+0,  // Total Number of Pairs
       ecnt_info.n_pair);
@@ -3672,10 +3673,10 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info,
       ecnt_info.n_ens);
 
    at.set_entry(r, c+2,  // Continuous Ranked Probability Score
-      ecnt_info.crps);
+      ecnt_info.crps_gaus);
 
    at.set_entry(r, c+3,  // Continuous Ranked Probability Skill Score
-      ecnt_info.crpss);
+      ecnt_info.crpss_gaus);
 
    at.set_entry(r, c+4,  // Ignorance Score
       ecnt_info.ign);
@@ -3701,6 +3702,18 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info,
    at.set_entry(r, c+11,  // Mean of unperturbed spread plus observation error
       ecnt_info.spread_plus_oerr);
 
+   at.set_entry(r, c+12,  // Gaussian climatological CRPS
+      ecnt_info.crpscl_gaus);
+
+   at.set_entry(r, c+13,  // Empirical ensemble CRPS
+      ecnt_info.crps_emp);
+
+   at.set_entry(r, c+14,  // Empirical climatological CRPS
+      ecnt_info.crpscl_emp);
+
+   at.set_entry(r, c+15,  // Empirical CRPSS
+      ecnt_info.crpss_emp);
+
    return;
 }
 

diff --git a/met/src/libcode/vx_statistics/compute_stats.cc b/met/src/libcode/vx_statistics/compute_stats.cc
@@ -1306,12 +1306,23 @@ void compute_ecnt_mean(const ECNTInfo *ecnt_info, int n,
    ecnt_mean.n_pair = na.sum();
 
    // Compute unweighted mean for each statistic
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crps_emp);
+   ecnt_mean.crps_emp = na.mean();
 
-   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crps);
-   ecnt_mean.crps = na.mean();
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crpscl_emp);
+   ecnt_mean.crpscl_emp = na.mean();
 
-   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crpss);
-   ecnt_mean.crpss = na.mean();
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crpss_emp);
+   ecnt_mean.crpss_emp = na.mean();
+
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crps_gaus);
+   ecnt_mean.crps_gaus = na.mean();
+
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crpscl_gaus);
+   ecnt_mean.crpscl_gaus = na.mean();
+
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].crpss_gaus);
+   ecnt_mean.crpss_gaus = na.mean();
 
    for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].ign);
    ecnt_mean.ign = na.mean();

diff --git a/met/src/libcode/vx_statistics/ens_stats.cc b/met/src/libcode/vx_statistics/ens_stats.cc
@@ -176,11 +176,13 @@ void ECNTInfo::init_from_scratch() {
 void ECNTInfo::clear() {
 
    othresh.clear();
-   n_ens   = n_pair    = 0;
-   crps    = crpss     = ign         = bad_data_double;
-   me      = rmse      = spread      = bad_data_double;
-   me_oerr = rmse_oerr = spread_oerr = bad_data_double;
-   spread_plus_oerr    = bad_data_double;
+   n_ens      = n_pair      = 0;
+   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;
+   spread_plus_oerr        = bad_data_double;
 
    return;
 }
@@ -193,11 +195,16 @@ void ECNTInfo::assign(const ECNTInfo &c) {
 
    n_ens            = c.n_ens;
    n_pair           = c.n_pair;
+
+   crps_emp         = c.crps_emp;
+   crpscl_emp       = c.crpscl_emp;
+   crpss_emp        = c.crpss_emp;
 
-   crps             = c.crps;
-   crpss            = c.crpss;
-   ign              = c.ign;
+   crps_gaus        = c.crps_gaus;
+   crpscl_gaus      = c.crpscl_gaus;
+   crpss_gaus       = c.crpss_gaus;
 
+   ign              = c.ign;
    me               = c.me;
    rmse             = c.rmse;
    spread           = c.spread;
@@ -215,15 +222,16 @@ void ECNTInfo::assign(const ECNTInfo &c) {
 void ECNTInfo::set(const PairDataEnsemble &pd) {
    int i;
    double w, w_sum;
-   double crps_climo;
    double fbar, obar, ffbar, oobar, fobar;
    NumArray cur;
 
    // Store the number of ensemble members
    n_ens = pd.n_ens;
-
-   // Get the average CRPS value
-   crps = pd.crps_na.wmean(pd.wgt_na);
+
+   // Get the average empirical and Gaussian CRPS and
+   // components values
+   crps_emp  = pd.crps_emp_na.wmean(pd.wgt_na);
+   crps_gaus = pd.crps_gaus_na.wmean(pd.wgt_na);
 
    // Get the sum of the weights
    for(i=0, n_pair=0, w_sum=0.0; i<pd.wgt_na.n(); i++) {
@@ -234,11 +242,12 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
    }
 
    // Check for bad data
-   if(is_bad_data(crps)            ||
+   if(is_bad_data(crps_emp)        ||
+      is_bad_data(crps_gaus)       ||
       pd.cmn_na.n() != pd.o_na.n() ||
       pd.cmn_na.n() == 0           ||
       pd.cmn_na.has(bad_data_double)) {
-      crpss = bad_data_double;
+      crpss_emp = crpss_gaus = bad_data_double;
    }
    else {
 
@@ -254,11 +263,16 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
          oobar += w * pd.o_na[i]   * pd.o_na[i];
          fobar += w * pd.cmn_na[i] * pd.o_na[i];
       }
-      crps_climo = ffbar + oobar - 2.0*fobar;
 
-      // Compute skill score
-      crpss = (is_eq(crps_climo, 0.0) ?
-               bad_data_double : (crps_climo - crps)/crps_climo);
+      // TODO: MET #1451 correct this computation
+      crpscl_emp  = ffbar + oobar - 2.0*fobar;
+      crpscl_gaus = ffbar + oobar - 2.0*fobar;
+
+      // Compute CRPS skill scores
+      crpss_emp  = (is_eq(crpscl_emp, 0.0) ? bad_data_double :
+                    (crpscl_emp - crps_emp) /crpscl_emp);
+      crpss_gaus = (is_eq(crpscl_gaus, 0.0) ? bad_data_double :
+                    (crpscl_gaus - crps_gaus) /crpscl_gaus);
    }
 
    // Compute the average IGN value