From 11a5419931e3488e26e5e9359fc43efc46444478 Mon Sep 17 00:00:00 2001 From: johnhg Date: Wed, 17 Feb 2021 12:08:59 -0700 Subject: [PATCH] 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. --- met/.cproject | 362 ------------------ met/.project | 78 ---- .../table_files/met_header_columns_V10.0.txt | 2 +- met/docs/Users_Guide/appendixC.rst | 8 +- met/docs/Users_Guide/ensemble-stat.rst | 18 +- met/docs/Users_Guide/point-stat.rst | 2 +- met/docs/Users_Guide/refs.rst | 6 + met/src/basic/vx_util/num_array.cc | 4 +- met/src/basic/vx_util/stat_column_defs.h | 4 +- met/src/libcode/vx_stat_out/stat_columns.cc | 27 +- .../libcode/vx_statistics/compute_stats.cc | 19 +- met/src/libcode/vx_statistics/ens_stats.cc | 50 ++- met/src/libcode/vx_statistics/ens_stats.h | 5 +- .../vx_statistics/pair_data_ensemble.cc | 116 ++++-- .../vx_statistics/pair_data_ensemble.h | 19 +- .../core/stat_analysis/aggr_stat_line.cc | 54 ++- .../tools/core/stat_analysis/aggr_stat_line.h | 2 +- .../core/stat_analysis/parse_stat_line.cc | 11 +- .../core/stat_analysis/parse_stat_line.h | 5 +- test/hdr/met_10_0.hdr | 2 +- 20 files changed, 247 insertions(+), 547 deletions(-) delete mode 100644 met/.cproject delete mode 100644 met/.project diff --git a/met/.cproject b/met/.cproject deleted file mode 100644 index a4fe10144f..0000000000 --- a/met/.cproject +++ /dev/null @@ -1,362 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/met/.project b/met/.project deleted file mode 100644 index fcf9456036..0000000000 --- a/met/.project +++ /dev/null @@ -1,78 +0,0 @@ - - - MET_svn - - - - - - org.eclipse.cdt.managedbuilder.core.genmakebuilder - clean,full,incremental, - - - ?name? - - - - org.eclipse.cdt.make.core.append_environment - true - - - org.eclipse.cdt.make.core.autoBuildTarget - all - - - org.eclipse.cdt.make.core.buildArguments - - - - org.eclipse.cdt.make.core.buildCommand - make - - - org.eclipse.cdt.make.core.cleanBuildTarget - clean - - - org.eclipse.cdt.make.core.contents - org.eclipse.cdt.make.core.activeConfigSettings - - - org.eclipse.cdt.make.core.enableAutoBuild - false - - - org.eclipse.cdt.make.core.enableCleanBuild - true - - - org.eclipse.cdt.make.core.enableFullBuild - true - - - org.eclipse.cdt.make.core.fullBuildTarget - all - - - org.eclipse.cdt.make.core.stopOnError - true - - - org.eclipse.cdt.make.core.useDefaultBuildCmd - true - - - - - org.eclipse.cdt.managedbuilder.core.ScannerConfigBuilder - - - - - - org.eclipse.cdt.core.cnature - org.eclipse.cdt.core.ccnature - org.eclipse.cdt.managedbuilder.core.managedBuildNature - org.eclipse.cdt.managedbuilder.core.ScannerConfigNature - - diff --git a/met/data/table_files/met_header_columns_V10.0.txt b/met/data/table_files/met_header_columns_V10.0.txt index 851ae5e3ee..d2c8ae5c4f 100644 --- 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 index d5ad94b01a..05905e9c67 100644 --- 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 `). 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 `). 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 index 60ef284b97..99532dcaf6 100644 --- 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 `). 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 `) and using the empirical ensemble distribution (:ref:`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 index 5756e40cd5..1ec6fdacb9 100644 --- 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 `; :ref:`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 `; :ref:`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 index 269242adb3..7e02dcb6d1 100644 --- 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 index 0481f4516c..bf80325736 100644 --- 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 index 33d52d7156..d5a8cf0c1a 100644 --- 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 index 09905b4f08..570977fb17 100644 --- 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 index 674c3c59d3..7e6e74f66f 100644 --- 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= all ensemble members + if(is_eq(obs_cdf, 0.0)) integral += (obs - fcst); + + return(integral); +} + +//////////////////////////////////////////////////////////////////////// +// +// Compute the Gaussian continuous ranked probability score (CRPS), +// the ignorance score (IGN), and probability integral transform (PIT) +// +//////////////////////////////////////////////////////////////////////// + +void compute_crps_gaus_ign_pit(double obs, const NumArray &ens_na, + double &crps, double &ign, double &pit) { double m, s, z; // Mean and standard deviation of the ensemble values @@ -1755,7 +1820,7 @@ void compute_crps_ign_pit(double obs, const NumArray &ens_na, z = (obs - m)/s; - // Compute CRPS + // Compute Gaussian CRPS crps = s*(z*(2.0*znorm(z) - 1) + 2.0*dnorm(z) - 1.0/sqrt(pi)); // Compute IGN @@ -1763,6 +1828,7 @@ void compute_crps_ign_pit(double obs, const NumArray &ens_na, // Compute PIT pit = normal_cdf(obs, m, s); + } return; @@ -1805,9 +1871,10 @@ PairDataEnsemble subset_climo_cdf_bin(const PairDataEnsemble &pd, // required for ensemble output line types. // // Include in subset: - // wgt_na, o_na, cmn_na, csd_na, cdf_na, v_na, r_na, crps_na, - // ign_na, pit_na, var_na, var_oerr_na, - // var_plus_oerr_na, mn_na, mn_oerr_na, e_na + // wgt_na, o_na, cmn_na, csd_na, cdf_na, v_na, r_na, + // crps_emp_na, crps_gaus_na, ign_na, pit_na, + // var_na, var_oerr_na, var_plus_oerr_na, + // mn_na, mn_oerr_na, e_na // // Exclude from subset: // sid_sa, lat_na, lon_na, x_na, y_na, vld_ta, lvl_ta, elv_ta, @@ -1820,7 +1887,8 @@ PairDataEnsemble subset_climo_cdf_bin(const PairDataEnsemble &pd, out_pd.cdf_na.add(pd.cdf_na[i]); out_pd.v_na.add(pd.v_na[i]); out_pd.r_na.add(pd.r_na[i]); - out_pd.crps_na.add(pd.crps_na[i]); + out_pd.crps_emp_na.add(pd.crps_emp_na[i]); + out_pd.crps_gaus_na.add(pd.crps_gaus_na[i]); out_pd.ign_na.add(pd.ign_na[i]); out_pd.pit_na.add(pd.pit_na[i]); out_pd.skip_ba.add(false); diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.h b/met/src/libcode/vx_statistics/pair_data_ensemble.h index 7fa90737bf..49a848dd95 100644 --- a/met/src/libcode/vx_statistics/pair_data_ensemble.h +++ b/met/src/libcode/vx_statistics/pair_data_ensemble.h @@ -78,11 +78,14 @@ class PairDataEnsemble : public PairBase { NumArray *e_na; // Ensemble values [n_ens][n_obs] NumArray v_na; // Number of valid ensemble values [n_obs] NumArray r_na; // Observation ranks [n_obs] - NumArray crps_na; // Continuous Ranked Probability Score [n_obs] + + NumArray crps_emp_na; // Empirical Continuous Ranked Probability Score [n_obs] + NumArray crps_gaus_na; // Gaussian CRPS [n_obs] + NumArray ign_na; // Ignorance Score [n_obs] NumArray pit_na; // Probability Integral Transform [n_obs] - int n_ens; // Number of ensemble members + int n_ens; // Number of ensemble members int n_pair; // Number of valid pairs, n_obs - sum(skip_ba) bool skip_const; // Skip cases where the observation and // all ensemble members are constant @@ -107,7 +110,11 @@ class PairDataEnsemble : public PairBase { double ssvar_bin_size; // Variance bin size for spread/skill SSVARInfo *ssvar_bins; // Ensemble spread/skill bin information [n_ssvar_bin] - double crpss; // Continuous ranked probability skill score + double crpscl_emp; // Empirical climatological CRPS score + double crpss_emp; // Empirical CRPS skill score + double crpscl_gaus; // Guassian climatological CRPS score + double crpss_gaus; // Guassian CRPS skill score + double me; // ME for ensemble mean double rmse; // RMSE for ensemble mean double me_oerr; // ME for mean of perturbed members @@ -290,8 +297,10 @@ class VxPairDataEnsemble { // //////////////////////////////////////////////////////////////////////// -extern void compute_crps_ign_pit(double, const NumArray &, double &, - double &, double &); +extern double compute_crps_emp(double, const NumArray &); + +extern void compute_crps_gaus_ign_pit(double, const NumArray &, + double &, double &, double &); // Subset pairs for a specific climatology CDF bin extern PairDataEnsemble subset_climo_cdf_bin(const PairDataEnsemble &, diff --git a/met/src/tools/core/stat_analysis/aggr_stat_line.cc b/met/src/tools/core/stat_analysis/aggr_stat_line.cc index 686a7cd87a..034afdb0cb 100644 --- a/met/src/tools/core/stat_analysis/aggr_stat_line.cc +++ b/met/src/tools/core/stat_analysis/aggr_stat_line.cc @@ -2514,7 +2514,7 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, AggrENSInfo aggr; ECNTData cur; ConcatString key; - double crps_fcst, crps_climo, v; + double crps_emp, crpscl_emp, crps_gaus, crpscl_gaus, v; map::iterator it; // @@ -2553,7 +2553,8 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, // if(m.count(key) == 0) { aggr.ens_pd.clear(); - aggr.crps_climo_na.clear(); + aggr.crpscl_emp_na.clear(); + aggr.crpscl_gaus_na.clear(); aggr.hdr.clear(); m[key] = aggr; } @@ -2571,7 +2572,10 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, // // Store the current statistics and weight (TOTAL column) // - m[key].ens_pd.crps_na.add(cur.crps); + m[key].ens_pd.crps_emp_na.add(cur.crps_emp); + m[key].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.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)); @@ -2590,18 +2594,6 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, bad_data_double : cur.rmse_oerr * cur.rmse_oerr)); - // - // Compute and store climatological CRPS - // - if(!is_bad_data(cur.crps) && !is_bad_data(cur.crpss) && - !is_eq(cur.crpss, 1.0)) { - crps_climo = cur.crps / (1.0 - cur.crpss); - m[key].crps_climo_na.add(crps_climo); - } - else { - m[key].crps_climo_na.add(bad_data_double); - } - // // Keep track of the unique header column entries // @@ -2627,17 +2619,20 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, v = it->second.mse_oerr_na.wmean(it->second.ens_pd.wgt_na); it->second.ens_pd.rmse_oerr = (is_bad_data(v) ? bad_data_double : sqrt(v)); - crps_fcst = it->second.ens_pd.crps_na.wmean(it->second.ens_pd.wgt_na); - crps_climo = it->second.crps_climo_na.wmean(it->second.ens_pd.wgt_na); + 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); + 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); - // Compute aggregated CRPSS - if(!is_bad_data(crps_fcst) && !is_bad_data(crps_climo) && - !is_eq(crps_climo, 0.0)) { - it->second.ens_pd.crpss = (crps_climo - crps_fcst)/crps_climo; - } - else { - it->second.ens_pd.crpss = bad_data_double; - } + // Compute aggregated empirical CRPSS + it->second.ens_pd.crpss_emp = + (is_bad_data(crps_emp) || is_bad_data(crpscl_emp) && is_eq(crpscl_emp, 0.0) ? + bad_data_double : (crpscl_emp - crps_emp)/crpscl_emp); + + // Compute aggregated Gaussian CRPSS + it->second.ens_pd.crpss_gaus = + (is_bad_data(crps_gaus) || is_bad_data(crpscl_gaus) && is_eq(crpscl_gaus, 0.0) ? + bad_data_double : (crpscl_gaus - crps_gaus)/crpscl_gaus); } // end for it @@ -3008,7 +3003,7 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job, ORANKData cur; ConcatString key; int i, n_valid, n_bin; - double esum, esumsq, crps, ign, pit; + double esum, esumsq, crps_gaus, ign, pit; map::iterator it; // @@ -3102,9 +3097,10 @@ 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 CRPS, IGN, and PIT for the current point - compute_crps_ign_pit(cur.obs, cur.ens_na, crps, ign, pit); - m[key].ens_pd.crps_na.add(crps); + // Compute Empirical and Gaussian CRPS, IGN, and PIT for the current point + 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); diff --git a/met/src/tools/core/stat_analysis/aggr_stat_line.h b/met/src/tools/core/stat_analysis/aggr_stat_line.h index 05f404cf6e..071645ade9 100644 --- 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,7 @@ struct AggrISCInfo { struct AggrENSInfo { StatHdrInfo hdr; PairDataEnsemble ens_pd; - NumArray crps_climo_na; + 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 index e2d10201a0..365be05c90 100644 --- a/met/src/tools/core/stat_analysis/parse_stat_line.cc +++ b/met/src/tools/core/stat_analysis/parse_stat_line.cc @@ -347,10 +347,15 @@ void parse_ecnt_line(STATLine &l, ECNTData &e_data) { e_data.total = atoi(l.get_item("TOTAL")); e_data.n_ens = atof(l.get_item("N_ENS")); - e_data.crps = atof(l.get_item("CRPS")); - e_data.crpss = atof(l.get_item("CRPSS")); - e_data.ign = atof(l.get_item("IGN")); + e_data.crps_emp = atof(l.get_item("CRPS_EMP")); + e_data.crpscl_emp = atof(l.get_item("CRPSCL_EMP")); + e_data.crpss_emp = atof(l.get_item("CRPSS_EMP")); + + e_data.crps_gaus = atof(l.get_item("CRPS")); + e_data.crpscl_gaus = atof(l.get_item("CRPSCL")); + e_data.crpss_gaus = atof(l.get_item("CRPSS")); + e_data.ign = atof(l.get_item("IGN")); e_data.me = atof(l.get_item("ME")); e_data.rmse = atof(l.get_item("RMSE")); e_data.spread = atof(l.get_item("SPREAD")); diff --git a/met/src/tools/core/stat_analysis/parse_stat_line.h b/met/src/tools/core/stat_analysis/parse_stat_line.h index d3c6cbdf91..4f6dafa43b 100644 --- a/met/src/tools/core/stat_analysis/parse_stat_line.h +++ b/met/src/tools/core/stat_analysis/parse_stat_line.h @@ -62,8 +62,9 @@ struct MPRData { // Ensemble continuous statistics (ECNT) data structure struct ECNTData { int total, n_ens; - double crps, crpss, ign; - double me, rmse, spread; + double crps_emp, crpscl_emp, crpss_emp; + double crps_gaus, crpscl_gaus, crpss_gaus; + double ign, me, rmse, spread; double me_oerr, rmse_oerr, spread_oerr; double spread_plus_oerr; }; diff --git a/test/hdr/met_10_0.hdr b/test/hdr/met_10_0.hdr index 1c3004e680..8355e557a3 100644 --- a/test/hdr/met_10_0.hdr +++ b/test/hdr/met_10_0.hdr @@ -17,7 +17,7 @@ PJC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_L 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 _VAR_ 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 _VAR_ 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 BASE N_PTS _VAR_ -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 +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 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 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 CRPS IGN N_RANK CRPSS SPREAD _VAR_ 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 _VAR_