From 996197cff8c7272906602de8245ae26cd91549d5 Mon Sep 17 00:00:00 2001 From: johnhg Date: Mon, 1 Mar 2021 18:03:30 -0700 Subject: [PATCH] Feature 1684 bss and 1685 single reference model (#1689) * Per #1684, move an instance of the ClimoCDFInfo class into PairBase. Also define derive_climo_vals() and derive_climo_prob() utility functions. * Add to VxPairDataPoint and VxPairDataEnsemble functions to set the ClimoCDFInfo class. * Per #1684, update ensemble_stat and point_stat to set the ClimoCDFInfo object based on the contents of the config file. * Per #1684, update the vx_statistics library and stat_analysis to make calls to the new derive_climo_vals() and derive_climo_prob() functions. * Per #1684, since cdf_info is a member of PairBase class, need to handle it in the PairDataPoint and PairDataEnsemble assignment and subsetting logic. * Per #1684, during development, I ran across and then updated this log message. * Per #1684, working on log messages and figured that the regridding climo data should be moved from Debug(1) to at least Debug(2). * Per #1684 and #1685, update the logic for the derive_climo_vals() utility function. If only a single climo bin is requested, just return the climo mean. Otherwise, sample the requested number of values. * Per #1684, just fixing the format of this log message. * Per #1684, add a STATLine::get_offset() member function. * Per #1684, update parse_orank_line() logic. Rather than calling NumArray::clear() call NumArray::erase() to preserve allocated memory. Also, instead of parsing ensemble member values by column name, parse them by offset number. * Per #1684, call EnsemblePairData::extend() when parsing ORANK data to allocate one block of memory instead of bunches of litte ones. * Per #1684 and #1685, add another call to Ensemble-Stat to test computing the CRPSCL_EMP from a single climo mean instead of using the full climo distribution. * Per #1684 and #1685, update ensemble-stat docs about computing CRPSS_EMP relative to a single reference model. * Per #1684, need to update Grid-Stat to store the climo cdf info in the PairDataPoint objects. * Per #1684, remove debug print statements. * Per #1684, need to set cdf_info when aggregating MPR lines in Stat-Analysis. * Per #1684 and #1685, update PairDataEnsemble::compute_pair_vals() to print a log message indicating the climo data being used as reference: For a climo distribution defined by mean and stdev: DEBUG 3: Computing ensemble statistics relative to a 9-member climatological ensemble. For a single deterministic reference: DEBUG 3: Computing ensemble statistics relative to the climatological mean. --- met/docs/Users_Guide/ensemble-stat.rst | 6 +- met/src/libcode/vx_analysis_util/stat_line.cc | 16 +- met/src/libcode/vx_analysis_util/stat_line.h | 9 +- .../libcode/vx_statistics/compute_stats.cc | 3 +- met/src/libcode/vx_statistics/ens_stats.cc | 11 +- met/src/libcode/vx_statistics/pair_base.cc | 100 ++++++-- met/src/libcode/vx_statistics/pair_base.h | 48 ++-- .../vx_statistics/pair_data_ensemble.cc | 58 ++--- .../vx_statistics/pair_data_ensemble.h | 8 +- .../libcode/vx_statistics/pair_data_point.cc | 21 ++ .../libcode/vx_statistics/pair_data_point.h | 2 + met/src/libcode/vx_statistics/read_climo.cc | 2 +- .../tools/core/ensemble_stat/ensemble_stat.cc | 2 +- .../ensemble_stat/ensemble_stat_conf_info.cc | 4 +- met/src/tools/core/grid_stat/grid_stat.cc | 73 +++--- met/src/tools/core/point_stat/point_stat.cc | 2 +- .../core/point_stat/point_stat_conf_info.cc | 3 + .../core/stat_analysis/aggr_stat_line.cc | 45 ++-- .../tools/core/stat_analysis/aggr_stat_line.h | 1 + .../core/stat_analysis/parse_stat_line.cc | 9 +- .../tools/core/stat_analysis/stat_analysis.cc | 4 +- test/config/EnsembleStatConfig_one_cdf_bin | 237 ++++++++++++++++++ test/xml/unit_climatology.xml | 29 +++ 23 files changed, 550 insertions(+), 143 deletions(-) create mode 100644 test/config/EnsembleStatConfig_one_cdf_bin diff --git a/met/docs/Users_Guide/ensemble-stat.rst b/met/docs/Users_Guide/ensemble-stat.rst index 9a4b8f476e..cfe920c0e6 100644 --- a/met/docs/Users_Guide/ensemble-stat.rst +++ b/met/docs/Users_Guide/ensemble-stat.rst @@ -36,7 +36,11 @@ The ranked probability score (RPS) is included in the Ranked Probability Score ( Climatology data ~~~~~~~~~~~~~~~~ -The Ensemble-Stat output includes at least three statistics computed relative to external climatology data. The climatology is defined by mean and standard deviation fields, and both are required in the computation of ensemble skill score statistics. MET assumes that the climatology follows a normal distribution, defined by the mean and standard deviation at each point. When computing the CRPS skill score for (:ref:`Gneiting et al., 2004 `) the reference CRPS statistic is computed using the climatological mean and standard deviation directly. When computing the CRPS skill score for (:ref:`Hersbach, 2000 `) the reference CRPS statistic is computed by selecting equal-area-spaced values from the assumed normal climatological distribution. The number of points selected is determined by the *cdf_bins* setting in the *climo_cdf* dictionary. The reference CRPS is computed empirically from this ensemble of climatology values. The climatological distribution is also used for the RPSS. The forecast RPS statistic is computed from a probabilistic contingency table in which the probabilities are derived from the ensemble member values. In a simliar fashion, the climatogical probability for each observed value is derived from the climatological distribution. The area of the distribution to the left of the observed value is interpreted as the climatological probability. These climatological probabilities are also evaluated using a probabilistic contingency table from which the reference RPS score is computed. The skill scores are derived by comparing the forecast statistic to the reference climatology statistic. +The Ensemble-Stat output includes at least three statistics computed relative to external climatology data. The climatology is defined by mean and standard deviation fields, and typically both are required in the computation of ensemble skill score statistics. MET assumes that the climatology follows a normal distribution, defined by the mean and standard deviation at each point. + +When computing the CRPS skill score for (:ref:`Gneiting et al., 2004 `) the reference CRPS statistic is computed using the climatological mean and standard deviation directly. When computing the CRPS skill score for (:ref:`Hersbach, 2000 `) the reference CRPS statistic is computed by selecting equal-area-spaced values from the assumed normal climatological distribution. The number of points selected is determined by the *cdf_bins* setting in the *climo_cdf* dictionary. The reference CRPS is computed empirically from this ensemble of climatology values. If the number bins is set to 1, the climatological CRPS is computed using only the climatological mean value. In this way, the empirical CRPSS may be computed relative to a single model rather than a climatological distribution. + +The climatological distribution is also used for the RPSS. The forecast RPS statistic is computed from a probabilistic contingency table in which the probabilities are derived from the ensemble member values. In a simliar fashion, the climatogical probability for each observed value is derived from the climatological distribution. The area of the distribution to the left of the observed value is interpreted as the climatological probability. These climatological probabilities are also evaluated using a probabilistic contingency table from which the reference RPS score is computed. The skill scores are derived by comparing the forecast statistic to the reference climatology statistic. Ensemble observation error ~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/met/src/libcode/vx_analysis_util/stat_line.cc b/met/src/libcode/vx_analysis_util/stat_line.cc index b1d66871e0..9c2dadb7e8 100644 --- a/met/src/libcode/vx_analysis_util/stat_line.cc +++ b/met/src/libcode/vx_analysis_util/stat_line.cc @@ -327,6 +327,18 @@ bool STATLine::has(const char *col_str) const { +return ( !is_bad_data(get_offset(col_str)) ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +int STATLine::get_offset(const char *col_str) const + +{ + int offset = bad_data_int; int dim = bad_data_int; @@ -353,10 +365,10 @@ if ( is_bad_data(offset) ) { } // - // Return whether a valid offset value was found + // Return the offset value // -return ( !is_bad_data(offset) ); +return ( offset ); } diff --git a/met/src/libcode/vx_analysis_util/stat_line.h b/met/src/libcode/vx_analysis_util/stat_line.h index ed52829950..6362031d58 100644 --- a/met/src/libcode/vx_analysis_util/stat_line.h +++ b/met/src/libcode/vx_analysis_util/stat_line.h @@ -61,10 +61,11 @@ class STATLine : public DataLine { // retrieve values of the header columns // - bool has (const char *) const; - ConcatString get (const char *, bool check_na = true) const; - const char * get_item (const char *, bool check_na = true) const; - const char * get_item (int, bool check_na = true) const; + bool has (const char *) const; + int get_offset(const char *) const; + ConcatString get (const char *, bool check_na = true) const; + const char * get_item (const char *, bool check_na = true) const; + const char * get_item (int, bool check_na = true) const; const char * version () const; const char * model () const; diff --git a/met/src/libcode/vx_statistics/compute_stats.cc b/met/src/libcode/vx_statistics/compute_stats.cc index 7e6e74f66f..211fe9860e 100644 --- a/met/src/libcode/vx_statistics/compute_stats.cc +++ b/met/src/libcode/vx_statistics/compute_stats.cc @@ -743,7 +743,8 @@ void compute_pctinfo(const PairDataPoint &pd, bool pstd_flag, // Use input climatological probabilities or derive them if(cmn_flag) { if(cprob_in) climo_prob = *cprob_in; - else climo_prob = derive_climo_prob(pd.cmn_na, pd.csd_na, + else climo_prob = derive_climo_prob(pd.cdf_info, + pd.cmn_na, pd.csd_na, pct_info.othresh); } diff --git a/met/src/libcode/vx_statistics/ens_stats.cc b/met/src/libcode/vx_statistics/ens_stats.cc index 4748b186ac..ebfe4b5e28 100644 --- a/met/src/libcode/vx_statistics/ens_stats.cc +++ b/met/src/libcode/vx_statistics/ens_stats.cc @@ -484,10 +484,10 @@ void RPSInfo::set(const PairDataEnsemble &pd) { // Check that thresholds are actually defined if(fthresh.n() == 0) { mlog << Error << "\nRPSInfo::set(const PairDataEnsemble &) -> " - << "no thresholds provided to compute the RPS line type! " - << "Specify thresholds using the \"" - << conf_key_prob_cat_thresh - << "\" configuration file option.\n\n"; + << "no thresholds provided to compute the RPS line type!\n" + << "Specify thresholds using the \"" << conf_key_prob_cat_thresh + << "\" configuration file option or by providing climatological " + << "mean and standard deviation data.\n\n"; exit(1); } @@ -522,7 +522,8 @@ void RPSInfo::set(const PairDataEnsemble &pd) { climo_pct.zero_out(); // Derive climatological probabilities - if(cmn_flag) climo_prob = derive_climo_prob(pd.cmn_na, pd.csd_na, + if(cmn_flag) climo_prob = derive_climo_prob(pd.cdf_info, + pd.cmn_na, pd.csd_na, fthresh[i]); // Loop over the observations diff --git a/met/src/libcode/vx_statistics/pair_base.cc b/met/src/libcode/vx_statistics/pair_base.cc index 490037c78b..8066ed262f 100644 --- a/met/src/libcode/vx_statistics/pair_base.cc +++ b/met/src/libcode/vx_statistics/pair_base.cc @@ -74,6 +74,8 @@ void PairBase::clear() { interp_mthd = InterpMthd_None; interp_shape = GridTemplateFactory::GridTemplate_None; + cdf_info.clear(); + o_na.clear(); x_na.clear(); y_na.clear(); @@ -121,6 +123,8 @@ void PairBase::erase() { interp_mthd = InterpMthd_None; interp_shape = GridTemplateFactory::GridTemplate_None; + cdf_info.clear(); + o_na.erase(); x_na.erase(); y_na.erase(); @@ -267,6 +271,15 @@ void PairBase::set_interp_shape(GridTemplateFactory::GridTemplates shape) { //////////////////////////////////////////////////////////////////////// +void PairBase::set_climo_cdf_info(const ClimoCDFInfo &info) { + + cdf_info = info; + + return; +} + +//////////////////////////////////////////////////////////////////////// + void PairBase::set_fcst_ut(unixtime ut){ fcst_ut = ut; @@ -1010,11 +1023,49 @@ bool set_climo_flag(const NumArray &f_na, const NumArray &c_na) { //////////////////////////////////////////////////////////////////////// -NumArray derive_climo_prob(const NumArray &mn_na, const NumArray &sd_na, +void derive_climo_vals(const ClimoCDFInfo &cdf_info, + double m, double s, + NumArray &climo_vals) { + + // Initialize + climo_vals.erase(); + + // cdf_info.cdf_ta starts with >=0.0 and ends with >=1.0. + // The number of bins is the number of thresholds minus 1. + + // Check for bad mean value + if(is_bad_data(m) || cdf_info.cdf_ta.n() < 2) { + return; + } + // Single climo bin + else if(cdf_info.cdf_ta.n() == 2) { + climo_vals.add(m); + } + // Check for bad standard deviation value + else if(is_bad_data(s)) { + return; + } + // Extract climo distribution values + else { + + // Skip the first and last thresholds + for(int i=1; i " + mlog << Error << "\nderive_climo_prob() -> " << "climatological threshold \"" << othresh.get_str() << "\" cannot be converted to a probability!\n\n"; exit(1); @@ -1060,23 +1111,17 @@ NumArray derive_climo_prob(const NumArray &mn_na, const NumArray &sd_na, // threshold else if(n_mn > 0 && n_sd > 0) { - mlog << Debug(2) - << "Deriving normal approximation of climatological " - << "probabilities for threshold " << othresh.get_str() - << ".\n"; + // The first (>=0.0) and last (>=1.0) climo thresholds are omitted + mlog << Debug(4) + << "Deriving climatological probabilities for threshold " + << othresh.get_str() << " by sampling " << cdf_info.cdf_ta.n()-2 + << " values from the normal climatological distribution.\n"; - // Compute probability value for each point + // Compute the probability by sampling from the climo distribution + // and deriving the event frequency for(i=0; i 2) { + mlog << Debug(3) + << "Computing ensemble statistics relative to a " + << cdf_info.cdf_ta.n() - 2 + << "-member climatological ensemble.\n"; + } + else { + mlog << Debug(3) + << "No reference climatology data provided.\n"; + } + // Compute the rank for each observation for(i=0, n_pair=0, n_skip_const=0, n_skip_vld=0; i=0.0) and last (>=1.0) climo CDF thresholds - for(int i=1; igrid() == vx_grid)) { - mlog << Debug(1) + mlog << Debug(2) << "Regridding the " << cur_ut_cs << " \"" << info->magic_str() << "\" climatology field to the verification grid.\n"; diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc index 31ab7b8b56..f3affe5895 100644 --- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc +++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc @@ -1735,7 +1735,7 @@ void process_grid_vx() { // Initialize pd_all.clear(); pd_all.set_ens_size(n_vx_vld[i]); - pd_all.set_climo_cdf(conf_info.vx_opt[i].cdf_info); + pd_all.set_climo_cdf_info(conf_info.vx_opt[i].cdf_info); pd_all.skip_const = conf_info.vx_opt[i].vx_pd.pd[0][0][0].skip_const; // Apply the current mask to the fields and compute the pairs 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 index 0461bc1b7f..e4361e041a 100644 --- a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc +++ b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc @@ -877,8 +877,8 @@ void EnsembleStatVxOpt::set_vx_pd(EnsembleStatConfInfo *conf_info) { // Define the dimensions vx_pd.set_pd_size(n_msg_typ, n_mask, n_interp); - // Store climo CDF - vx_pd.set_climo_cdf(cdf_info); + // Store the climo CDF info + vx_pd.set_climo_cdf_info(cdf_info); // Store the list of surface message types vx_pd.set_msg_typ_sfc(conf_info->msg_typ_sfc); diff --git a/met/src/tools/core/grid_stat/grid_stat.cc b/met/src/tools/core/grid_stat/grid_stat.cc index 5a90a6464f..14b09c6d5a 100644 --- a/met/src/tools/core/grid_stat/grid_stat.cc +++ b/met/src/tools/core/grid_stat/grid_stat.cc @@ -145,7 +145,8 @@ static void build_outfile_name(unixtime, int, const char *, static void process_scores(); -static void get_mask_points(const MaskPlane &, const DataPlane *, +static void get_mask_points(const GridStatVxOpt &, + const MaskPlane &, const DataPlane *, const DataPlane *, const DataPlane *, const DataPlane *, const DataPlane *, PairDataPoint &); @@ -797,7 +798,8 @@ void process_scores() { } // Apply the current mask to the current fields - get_mask_points(mask_mp, &fcst_dp_smooth, &obs_dp_smooth, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp_smooth, &obs_dp_smooth, &cmn_dp, &csd_dp, &wgt_dp, pd); // Set the mask name @@ -981,7 +983,8 @@ void process_scores() { } // Apply the current mask to the U-wind fields - get_mask_points(mask_mp, &fu_dp_smooth, &ou_dp_smooth, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fu_dp_smooth, &ou_dp_smooth, &cmnu_dp, &csdu_dp, &wgt_dp, pd_u); // Compute VL1L2 @@ -1136,9 +1139,11 @@ void process_scores() { mask_bad_data(mask_mp, ogy_dp); // Apply the current mask to the current fields - get_mask_points(mask_mp, &fgx_dp, &ogx_dp, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fgx_dp, &ogx_dp, 0, 0, &wgt_dp, pd_gx); - get_mask_points(mask_mp, &fgy_dp, &ogy_dp, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fgy_dp, &ogy_dp, 0, 0, &wgt_dp, pd_gy); // Set the mask name @@ -1217,7 +1222,8 @@ void process_scores() { conf_info.vx_opt[i].ocat_ta.need_perc()) { // Apply the current mask - get_mask_points(mask_mp, &fcst_dp, &obs_dp, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp, &obs_dp, &cmn_dp, 0, 0, pd); // Process percentile thresholds @@ -1271,9 +1277,11 @@ void process_scores() { // Apply the current mask to the distance map and // thresholded fields - get_mask_points(mask_mp, &fcst_dp_dmap, &obs_dp_dmap, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp_dmap, &obs_dp_dmap, 0, 0, 0, pd); - get_mask_points(mask_mp, &fcst_dp_thresh, &obs_dp_thresh, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp_thresh, &obs_dp_thresh, 0, 0, 0, pd_thr); dmap_info.set_options( @@ -1346,7 +1354,8 @@ void process_scores() { conf_info.vx_opt[i].ocat_ta.need_perc()) { // Apply the current mask - get_mask_points(mask_mp, &fcst_dp, &obs_dp, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp, &obs_dp, &cmn_dp, 0, 0, pd); // Process percentile thresholds @@ -1445,9 +1454,11 @@ void process_scores() { // Apply the current mask to the fractional coverage // and thresholded fields - get_mask_points(mask_mp, &fcst_dp_smooth, &obs_dp_smooth, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp_smooth, &obs_dp_smooth, 0, 0, &wgt_dp, pd); - get_mask_points(mask_mp, &fcst_dp_thresh, &obs_dp_thresh, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp_thresh, &obs_dp_thresh, 0, 0, 0, pd_thr); // Store climatology values as bad data @@ -1618,7 +1629,8 @@ void process_scores() { } // Apply the current mask to the current fields - get_mask_points(mask_mp, &fcst_dp_smooth, &obs_dp_smooth, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fcst_dp_smooth, &obs_dp_smooth, &cmn_dp_smooth, &csd_dp, &wgt_dp, pd); // Set the mask name @@ -1706,7 +1718,8 @@ void process_scores() { } // Apply the current mask to the U-wind fields - get_mask_points(mask_mp, &fu_dp_smooth, &ou_dp_smooth, + get_mask_points(conf_info.vx_opt[i], mask_mp, + &fu_dp_smooth, &ou_dp_smooth, &cmnu_dp_smooth, 0, &wgt_dp, pd_u); // Compute VL1L2 @@ -1790,29 +1803,33 @@ void process_scores() { //////////////////////////////////////////////////////////////////////// -void get_mask_points(const MaskPlane &mask_mp, +void get_mask_points(const GridStatVxOpt &vx_opt, + const MaskPlane &mask_mp, const DataPlane *fcst_ptr, const DataPlane *obs_ptr, const DataPlane *cmn_ptr, const DataPlane *csd_ptr, const DataPlane *wgt_ptr, PairDataPoint &pd) { - // Initialize - pd.erase(); + // Initialize + pd.erase(); - // Apply the mask the data fields or fill with default values - apply_mask(*fcst_ptr, mask_mp, pd.f_na); - apply_mask(*obs_ptr, mask_mp, pd.o_na); - pd.n_obs = pd.o_na.n(); + // Store the climo CDF info + pd.set_climo_cdf_info(vx_opt.cdf_info); + + // Apply the mask the data fields or fill with default values + apply_mask(*fcst_ptr, mask_mp, pd.f_na); + apply_mask(*obs_ptr, mask_mp, pd.o_na); + pd.n_obs = pd.o_na.n(); - if(cmn_ptr) apply_mask(*cmn_ptr, mask_mp, pd.cmn_na); - else pd.cmn_na.add_const(bad_data_double, pd.n_obs); - if(csd_ptr) apply_mask(*csd_ptr, mask_mp, pd.csd_na); - else pd.csd_na.add_const(bad_data_double, pd.n_obs); - if(wgt_ptr) apply_mask(*wgt_ptr, mask_mp, pd.wgt_na); - else pd.wgt_na.add_const(default_grid_weight, pd.n_obs); + if(cmn_ptr) apply_mask(*cmn_ptr, mask_mp, pd.cmn_na); + else pd.cmn_na.add_const(bad_data_double, pd.n_obs); + if(csd_ptr) apply_mask(*csd_ptr, mask_mp, pd.csd_na); + else pd.csd_na.add_const(bad_data_double, pd.n_obs); + if(wgt_ptr) apply_mask(*wgt_ptr, mask_mp, pd.wgt_na); + else pd.wgt_na.add_const(default_grid_weight, pd.n_obs); - if(cmn_ptr && csd_ptr) pd.add_climo_cdf(); + if(cmn_ptr && csd_ptr) pd.add_climo_cdf(); - return; + return; } //////////////////////////////////////////////////////////////////////// diff --git a/met/src/tools/core/point_stat/point_stat.cc b/met/src/tools/core/point_stat/point_stat.cc index f61583f3ee..d1edfdb935 100644 --- a/met/src/tools/core/point_stat/point_stat.cc +++ b/met/src/tools/core/point_stat/point_stat.cc @@ -1762,7 +1762,7 @@ void do_hira_ens(int i_vx, const PairDataPoint *pd_ptr) { hira_pd.clear(); hira_pd.extend(pd_ptr->n_obs); hira_pd.set_ens_size(gt->size()); - hira_pd.set_climo_cdf(conf_info.vx_opt[i_vx].cdf_info); + hira_pd.set_climo_cdf_info(conf_info.vx_opt[i_vx].cdf_info); f_ens.extend(gt->size()); // Process each observation point diff --git a/met/src/tools/core/point_stat/point_stat_conf_info.cc b/met/src/tools/core/point_stat/point_stat_conf_info.cc index 0d8dd3fa43..ecd6b8b3dc 100644 --- a/met/src/tools/core/point_stat/point_stat_conf_info.cc +++ b/met/src/tools/core/point_stat/point_stat_conf_info.cc @@ -932,6 +932,9 @@ void PointStatVxOpt::set_vx_pd(PointStatConfInfo *conf_info) { // Define the dimensions vx_pd.set_pd_size(n_msg_typ, n_mask, n_interp); + // Store the climo CDF info + vx_pd.set_climo_cdf_info(cdf_info); + // Store the surface message type group cs = surface_msg_typ_group_str; if(conf_info->msg_typ_group_map.count(cs) == 0) { 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 3bd72ae766..db7cd98ba9 100644 --- a/met/src/tools/core/stat_analysis/aggr_stat_line.cc +++ b/met/src/tools/core/stat_analysis/aggr_stat_line.cc @@ -575,7 +575,22 @@ ConcatString StatHdrInfo::get_shc_str(const ConcatString &cur_case, //////////////////////////////////////////////////////////////////////// // -// Code for AggrTimeSeriesInfo structure. +// Code for AggrENSInfo structure +// +//////////////////////////////////////////////////////////////////////// + +void AggrENSInfo::clear() { + hdr.clear(); + ens_pd.clear(); + me_na.clear(); + mse_na.clear(); + me_oerr_na.clear(); + mse_oerr_na.clear(); +} + +//////////////////////////////////////////////////////////////////////// +// +// Code for AggrTimeSeriesInfo structure // //////////////////////////////////////////////////////////////////////// @@ -2190,6 +2205,9 @@ void aggr_mpr_lines(LineDataFile &f, STATAnalysisJob &job, // if(m.count(key) == 0) { + bool center = false; + aggr.pd.cdf_info.set_cdf_ta(nint(1.0/job.out_bin_size), center); + aggr.pd.f_na.clear(); aggr.pd.o_na.clear(); aggr.pd.cmn_na.clear(); @@ -2208,6 +2226,7 @@ void aggr_mpr_lines(LineDataFile &f, STATAnalysisJob &job, aggr.fcst_var = cur.fcst_var; aggr.obs_var = cur.obs_var; aggr.hdr.clear(); + m[key] = aggr; } // @@ -2552,8 +2571,7 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, // Add a new map entry, if necessary // if(m.count(key) == 0) { - aggr.ens_pd.clear(); - aggr.hdr.clear(); + aggr.clear(); m[key] = aggr; } @@ -2775,8 +2793,7 @@ void aggr_rhist_lines(LineDataFile &f, STATAnalysisJob &job, // Add a new map entry, if necessary // if(m.count(key) == 0) { - aggr.ens_pd.clear(); - aggr.hdr.clear(); + aggr.clear(); for(i=0; i::iterator it; @@ -3045,17 +3062,17 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job, // Add a new map entry, if necessary // if(m.count(key) == 0) { - aggr.ens_pd.clear(); + aggr.clear(); + bool center = false; + aggr.ens_pd.cdf_info.set_cdf_ta(nint(1.0/job.out_bin_size), center); aggr.ens_pd.obs_error_flag = !is_bad_data(cur.ens_mean_oerr); aggr.ens_pd.set_ens_size(cur.n_ens); + aggr.ens_pd.extend(cur.total); for(i=0; i " - << "the \"N_ENS\" column must remain constant. " + << "the \"N_ENS\" column must remain constant. " << "Try setting \"-column_eq N_ENS n\".\n\n"; throw(1); } @@ -3099,12 +3116,12 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job, m[key].ens_pd.v_na.add(n_valid); // 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); + derive_climo_vals(m[key].ens_pd.cdf_info, + cur.climo_mean, cur.climo_stdev, climo_vals); // Store empirical CRPS stats m[key].ens_pd.crps_emp_na.add(compute_crps_emp(cur.obs, cur.ens_na)); - m[key].ens_pd.crpscl_emp_na.add(compute_crps_emp(cur.obs, cur_clm)); + m[key].ens_pd.crpscl_emp_na.add(compute_crps_emp(cur.obs, climo_vals)); // Store Gaussian CRPS stats m[key].ens_pd.crps_gaus_na.add(compute_crps_gaus(cur.obs, cur.ens_mean, cur.spread)); 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 2e6ab5b72e..e8ed7809f9 100644 --- a/met/src/tools/core/stat_analysis/aggr_stat_line.h +++ b/met/src/tools/core/stat_analysis/aggr_stat_line.h @@ -152,6 +152,7 @@ struct AggrENSInfo { StatHdrInfo hdr; PairDataEnsemble ens_pd; NumArray me_na, mse_na, me_oerr_na, mse_oerr_na; + void clear(); }; struct AggrRPSInfo { 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 2f79ae4c3a..14bf981ea6 100644 --- a/met/src/tools/core/stat_analysis/parse_stat_line.cc +++ b/met/src/tools/core/stat_analysis/parse_stat_line.cc @@ -461,8 +461,7 @@ void parse_relp_line(STATLine &l, RELPData &r_data) { //////////////////////////////////////////////////////////////////////// void parse_orank_line(STATLine &l, ORANKData &o_data) { - int i; - char col_str[max_str_len]; + int i, ens1; o_data.total = atoi(l.get_item("TOTAL")); o_data.index = atoi(l.get_item("INDEX")); @@ -480,10 +479,10 @@ void parse_orank_line(STATLine &l, ORANKData &o_data) { o_data.n_ens = atoi(l.get_item("N_ENS")); // Parse out ENS_i - o_data.ens_na.clear(); + o_data.ens_na.erase(); + ens1 = l.get_offset("ENS_1"); for(i=0; i + + + + echo "&DATA_DIR_MODEL;/grib1/arw-fer-gep1/arw-fer-gep1_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-fer-gep5/arw-fer-gep5_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep2/arw-sch-gep2_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep6/arw-sch-gep6_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep3/arw-tom-gep3_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep7/arw-tom-gep7_2012040912_F024.grib" \ + > &OUTPUT_DIR;/climatology/ensemble_stat_input_file_list; \ + &MET_BIN;/ensemble_stat + + OUTPUT_PREFIX ONE_CDF_BIN + CLIMO_MEAN_FILE_LIST "&DATA_DIR_CLIMO;/NCEP_NCAR_40YR_1.0deg/cmean_1d.19590410" + + \ + &OUTPUT_DIR;/climatology/ensemble_stat_input_file_list \ + &CONFIG_DIR;/EnsembleStatConfig_one_cdf_bin \ + -point_obs &OUTPUT_DIR;/pb2nc/ndas.20120410.t12z.prepbufr.tm00.nc \ + -grid_obs &DATA_DIR_OBS;/laps/laps_2012041012_F000.grib \ + -outdir &OUTPUT_DIR;/climatology + + + &OUTPUT_DIR;/climatology/ensemble_stat_ONE_CDF_BIN_20120410_120000V.stat + &OUTPUT_DIR;/climatology/ensemble_stat_ONE_CDF_BIN_20120410_120000V_ecnt.txt + &OUTPUT_DIR;/climatology/ensemble_stat_ONE_CDF_BIN_20120410_120000V_ens.nc + + +