From 66f7fb23226b48c44f1c4d11675348d8aafd1d82 Mon Sep 17 00:00:00 2001 From: johnhg Date: Mon, 1 Mar 2021 20:54:34 -0700 Subject: [PATCH] Update develop-ref after #1688 and #1689 (#1690) * Per #1429, enhance error message from DataLine::get_item(). (#1682) * Feature 1429 tc_log second try (#1686) * Per #1429, enhance error message from DataLine::get_item(). * Per #1429, I realize that the line number actually is readily available in the DataLine class... so include it in the error message. * Feature 1588 ps_log (#1687) * Per #1588, updated pair_data_point.h/.cc to add detailed Debug(4) log messages, as specified in the GitHub issue. Do still need to test each of these cases to confirm that the log messages look good. * Per #1588, switch very detailed interpolation details from debug level 4 to 5. * Per #1588, remove the Debug(4) log message about duplicate obs since it's been moved up to a higher level. * Per #1588, add/update detailed log messages when processing point observations for bad data, off the grid, bad topo, big topo diffs, bad fcst value, and duplicate obs. * #1454 Disabled plot_data_plane_CESM_SSMI_microwave and plot_data_plane_CESM_sea_ice_nc becaues of not evenly spaced * #1454 Moved NC attribute name to nc_utils.h * #1454 Corrected sanity checking for lat/lon projection based on the percentage of the delta instead of fixed tolerance * #1454 Corrected sanity checking for lat/lon projection based on the percentage of the delta instead of fixed tolerance * #1454 Corrected data.delta_lon * #1454 Change bact to use diff instead of absolute value of diff * 454 Deleted instea dof commenting out * 454 Deleted instea dof commenting out * 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. Co-authored-by: Howard Soh Co-authored-by: hsoh-u --- met/docs/Users_Guide/ensemble-stat.rst | 6 +- met/src/basic/vx_util/data_line.cc | 16 +- met/src/basic/vx_util/interp_util.cc | 2 +- met/src/libcode/vx_analysis_util/stat_line.cc | 16 +- met/src/libcode/vx_analysis_util/stat_line.h | 9 +- met/src/libcode/vx_data2d_nc_met/met_file.cc | 5 - .../vx_data2d_nc_pinterp/pinterp_file.cc | 3 - met/src/libcode/vx_data2d_nccf/nccf_file.cc | 102 ++++---- met/src/libcode/vx_nc_util/nc_utils.cc | 5 - met/src/libcode/vx_nc_util/nc_utils.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 | 108 ++++++-- 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 | 134 ++++++++-- .../libcode/vx_statistics/pair_data_point.h | 9 + 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 +++ test/xml/unit_plot_data_plane.xml | 41 --- 31 files changed, 723 insertions(+), 281 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/basic/vx_util/data_line.cc b/met/src/basic/vx_util/data_line.cc index ab16353db1..da012b9a0a 100644 --- a/met/src/basic/vx_util/data_line.cc +++ b/met/src/basic/vx_util/data_line.cc @@ -253,7 +253,12 @@ const char * DataLine::get_item(int k) const if ( (k < 0) || (k >= N_items) ) { - mlog << Error << "\nDataLine::get_item(int) -> range check error\n\n"; + ConcatString cs = (File ? File->filename() : ""); + + mlog << Error << "\nDataLine::get_item(int) -> " + << "range check error reading line number " << LineNumber + << ", item number " << k+1 << " of " << N_items + << " from file \"" << cs << "\"\n\n"; exit ( 1 ); @@ -640,7 +645,8 @@ LineDataFile::LineDataFile(const LineDataFile &) { -mlog << Error << "\nLineDataFile::LineDataFile(const LineDataFile &) -> should never be called!\n\n"; +mlog << Error << "\nLineDataFile::LineDataFile(const LineDataFile &) -> " + << "should never be called!\n\n"; exit ( 1 ); @@ -654,7 +660,8 @@ LineDataFile & LineDataFile::operator=(const LineDataFile &) { -mlog << Error << "\nLineDataFile::operator=(const LineDataFile &) -> should never be called!\n\n"; +mlog << Error << "\nLineDataFile::operator=(const LineDataFile &) -> " + << "should never be called!\n\n"; exit ( 1 ); @@ -698,7 +705,8 @@ in = new ifstream; if ( !in ) { - mlog << Error << "\nLineDataFile::open(const char *) -> can't allocate input stream\n\n"; + mlog << Error << "\nLineDataFile::open(const char *) -> " + << "can't allocate input stream\n\n"; exit ( 1 ); diff --git a/met/src/basic/vx_util/interp_util.cc b/met/src/basic/vx_util/interp_util.cc index 90f3ae8c3f..9cb8a3d552 100644 --- a/met/src/basic/vx_util/interp_util.cc +++ b/met/src/basic/vx_util/interp_util.cc @@ -686,7 +686,7 @@ double interp_geog_match(const DataPlane &dp, const GridTemplate >, } if(!is_bad_data(interp_v)) { - mlog << Debug(4) + mlog << Debug(5) << "For observation value " << obs_v << " at grid (x, y) = (" << obs_x << ", " << obs_y << ") found forecast value " << interp_v << " at nearest matching geography point (" 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_data2d_nc_met/met_file.cc b/met/src/libcode/vx_data2d_nc_met/met_file.cc index 6f922b80e6..83c56a5e8e 100644 --- a/met/src/libcode/vx_data2d_nc_met/met_file.cc +++ b/met/src/libcode/vx_data2d_nc_met/met_file.cc @@ -41,11 +41,6 @@ static const string valid_time_ut_att_name = "valid_time_ut"; static const string init_time_ut_att_name = "init_time_ut"; static const string accum_time_att_name = "accum_time_sec"; -static const string name_att_name = "name"; -static const string long_name_att_name = "long_name"; -static const string level_att_name = "level"; -static const string units_att_name = "units"; - static const int max_met_args = 30; //////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_data2d_nc_pinterp/pinterp_file.cc b/met/src/libcode/vx_data2d_nc_pinterp/pinterp_file.cc index 4f9f24ca16..84cec1cf06 100644 --- a/met/src/libcode/vx_data2d_nc_pinterp/pinterp_file.cc +++ b/met/src/libcode/vx_data2d_nc_pinterp/pinterp_file.cc @@ -58,9 +58,6 @@ static const char hpa_units_str [] = "hPa"; static const string init_time_att_name = "START_DATE"; -static const string description_att_name = "description"; -static const string units_att_name = "units"; - static const int max_pinterp_args = 30; static const double pinterp_missing = 1.0e35; diff --git a/met/src/libcode/vx_data2d_nccf/nccf_file.cc b/met/src/libcode/vx_data2d_nccf/nccf_file.cc index b96ea5ed2d..91696d8a91 100644 --- a/met/src/libcode/vx_data2d_nccf/nccf_file.cc +++ b/met/src/libcode/vx_data2d_nccf/nccf_file.cc @@ -32,7 +32,8 @@ using namespace std; //////////////////////////////////////////////////////////////////////// -static const int max_met_args = 30; +static const int max_met_args = 30; +static const double DELTA_TOLERANCE_PERCENT = 0.05; const double NcCfFile::DELTA_TOLERANCE = 15.0; @@ -2870,31 +2871,32 @@ void NcCfFile::get_grid_from_lat_lon_vars(NcVar *lat_var, NcVar *lon_var, // Figure out the dlat/dlon values from the dimension variables - double lat_values[lat_counts]; - double lon_values[lon_counts]; - - bool two_dim_corrd = false; long x_size = get_data_size(lon_var); long y_size = get_data_size(lat_var); long latlon_counts = lon_counts*lat_counts; - if( x_size == latlon_counts || y_size == latlon_counts ) two_dim_corrd = true; - else if( x_size != lon_counts || y_size != lat_counts) + bool two_dim_corrd = (x_size == latlon_counts) && (y_size == latlon_counts ); + if( !two_dim_corrd && (x_size != lon_counts || y_size != lat_counts)) { mlog << Error << "\n" << method_name << " -> " << "Coordinate variables don't match dimension sizes in netCDF file.\n\n"; exit(1); } + double lat_values[lat_counts]; + double lon_values[lon_counts]; + bool lat_first = false; if (two_dim_corrd) { + lat_first = (lat_counts == get_dim_size(lat_var, 0)); long cur[2], length[2]; - for (int i=0; i<2; i++) { - cur[i] = 0; - length[i] = 1; - } - length[0] = lat_counts; + cur[0] = cur[1] = 0; + length[0] = length[1] = 1; + if (lat_first) length[0] = lat_counts; + else length[1] = lat_counts; get_nc_data(lat_var,lat_values, length, cur); - length[1] = lon_counts; - length[0] = 1; + + length[0] = length[1] = 1; + if (lat_first) length[1] = lon_counts; + else length[0] = lon_counts; get_nc_data(lon_var,lon_values, length, cur); } else { @@ -2902,8 +2904,7 @@ void NcCfFile::get_grid_from_lat_lon_vars(NcVar *lat_var, NcVar *lon_var, get_nc_data(lon_var,lon_values); } - // Calculate dlat and dlon assuming they are constant. MET requires that - // dlat be equal to dlon + // Calculate dlat and dlon assuming they are constant. double dlat = lat_values[1] - lat_values[0]; double dlon = rescale_lon(lon_values[1] - lon_values[0]); @@ -2919,53 +2920,58 @@ void NcCfFile::get_grid_from_lat_lon_vars(NcVar *lat_var, NcVar *lon_var, // entire grid. CF compliancy doesn't require this, but MET does. if (!skip_sanity_check) { + double degree_tolerance; float lat_missing_value = bad_data_double; float lon_missing_value = bad_data_double; - NcVarAtt *missing_value_att = (NcVarAtt*) 0; - missing_value_att = get_nc_att(lat_var, (string)"_FillValue"); - if (IS_VALID_NC_P(missing_value_att)) { - lat_missing_value = get_att_value_double(missing_value_att); - } - if( missing_value_att ) delete missing_value_att; - missing_value_att = get_nc_att(lon_var, (string)"_FillValue"); - if (IS_VALID_NC_P(missing_value_att)) { - lon_missing_value = get_att_value_double(missing_value_att); - } - if( missing_value_att ) delete missing_value_att; + bool sanity_check_failed = false; + get_var_att_float(lat_var, fill_value_att_name, lat_missing_value); + get_var_att_float(lon_var, fill_value_att_name, lon_missing_value); - if (fabs(dlat - dlon) > DELTA_TOLERANCE) - { - mlog << Error << "\n" << method_name << " -> " - << "MET can only process Latitude/Longitude files where the delta lat and delta lon are the same\n\n"; - exit(1); - } - + degree_tolerance = fabs(dlat * DELTA_TOLERANCE_PERCENT); for (int i = 1; i < (int)lat_counts; ++i) { - if ((fabs(lat_missing_value - lat_values[i]) < DELTA_TOLERANCE) || - (fabs(lat_missing_value - lat_values[i-1]) < DELTA_TOLERANCE)) continue; + if ((fabs(lat_missing_value - lat_values[i]) < degree_tolerance) || + (fabs(lat_missing_value - lat_values[i-1]) < degree_tolerance)) continue; double curr_delta = lat_values[i] - lat_values[i-1]; - if (fabs(curr_delta - dlat) > DELTA_TOLERANCE) + if (fabs(curr_delta - dlat) > degree_tolerance) { + mlog << Debug(4) << method_name << " -> lat[" + << i-1 << "]=" << lat_values[i-1] << " lat[" + << i << "]=" << lat_values[i] << " " + << fabs(curr_delta - dlat) << " > " << degree_tolerance << "\n"; mlog << Error << "\n" << method_name << " -> " - << "MET can only process Latitude/Longitude files where the lat delta is constant (dlat=" - << dlat <<", dlon=" << dlon << ")\n\n"; - exit(1); + << "MET can only process Latitude/Longitude files where the latitudes are evenly spaced (dlat=" + << dlat <<", delta[" << i << "]=" << curr_delta << ")\n\n"; + sanity_check_failed = true; + break; } } - + + degree_tolerance = fabs(dlon * DELTA_TOLERANCE_PERCENT); for (int i = 1; i < (int)lon_counts; ++i) { - if ((fabs(lon_missing_value - lon_values[i]) < DELTA_TOLERANCE) || - (fabs(lon_missing_value - lon_values[i-1]) < DELTA_TOLERANCE)) continue; + if ((fabs(lon_missing_value - lon_values[i]) < degree_tolerance) || + (fabs(lon_missing_value - lon_values[i-1]) < degree_tolerance)) continue; double curr_delta = rescale_lon(lon_values[i] - lon_values[i-1]); - if (fabs(curr_delta - dlon) > DELTA_TOLERANCE) + if (fabs(curr_delta - dlon) > degree_tolerance) { + mlog << Debug(4) << method_name << " -> lon[" + << i-1 << "]=" << lon_values[i-1] << " lon[" + << i << "]=" << lon_values[i] << " " + << fabs(curr_delta - dlon) << " > " << degree_tolerance << "\n"; mlog << Error << "\n" << method_name << " -> " - << "MET can only process Latitude/Longitude files where the lon delta is constant\n\n"; - exit(1); + << "MET can only process Latitude/Longitude files where the longitudes are evenly spaced (dlon=" + << dlon <<", delta[" << i << "]=" << curr_delta << ")\n\n"; + sanity_check_failed = true; + break; } } + + if (sanity_check_failed) { + mlog << Error << "\n" << method_name << " -> " + << "Please check the input data is the lat/lon projection\n\n"; + exit(1); + } } // Fill in the data structure. Remember to negate the longitude @@ -2991,9 +2997,9 @@ void NcCfFile::get_grid_from_lat_lon_vars(NcVar *lat_var, NcVar *lon_var, data.lat_ll = lat_values[lat_counts-1]; } - grid.set(data); + grid.set(data); // resets swap_to_north to false if (dlat < 0) grid.set_swap_to_north(true); - + } //////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_nc_util/nc_utils.cc b/met/src/libcode/vx_nc_util/nc_utils.cc index a0c302b96c..62b64e233b 100644 --- a/met/src/libcode/vx_nc_util/nc_utils.cc +++ b/met/src/libcode/vx_nc_util/nc_utils.cc @@ -25,11 +25,6 @@ using namespace netCDF::exceptions; //////////////////////////////////////////////////////////////////////// -static const string level_att_name = "level"; -static const string units_att_name = "units"; -static const string missing_value_att_name = "missing_value"; -static const string fill_value_att_name = "_FillValue"; - //////////////////////////////////////////////////////////////////////// void patch_nc_name(string *var_name) { diff --git a/met/src/libcode/vx_nc_util/nc_utils.h b/met/src/libcode/vx_nc_util/nc_utils.h index 9945afc563..79cc241f07 100644 --- a/met/src/libcode/vx_nc_util/nc_utils.h +++ b/met/src/libcode/vx_nc_util/nc_utils.h @@ -128,6 +128,15 @@ static const string nc_att_use_var_id = "use_var_id"; static const char nc_att_obs_version[] = "MET_Obs_version"; static const char nc_att_met_point_nccf[] = "MET_point_NCCF"; +static const string description_att_name = "description"; +static const string fill_value_att_name = "_FillValue"; +static const string level_att_name = "level"; +static const string long_name_att_name = "long_name"; +static const string missing_value_att_name = "missing_value"; +static const string name_att_name = "name"; +static const string units_att_name = "units"; + + static const char nc_time_unit_exp[] = "^[a-z|A-Z]* since [0-9]\\{1,4\\}-[0-9]\\{1,2\\}-[0-9]\\{1,2\\}"; static const char MET_NC_Obs_ver_1_2[] = "1.02"; 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 bfcebaad0c..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; @@ -424,13 +437,7 @@ bool PairBase::add_point_obs(const char *sid, if(check_unique) { vector::iterator o_it = (*it).second.obs.begin(); for(;o_it != (*it).second.obs.end(); o_it++) { - if( (*o_it).ut == ut) { - mlog << Debug(4) - << "Skipping duplicate observation for [lat:lon:level:elevation] = [" - << obs_key << "] valid at " << unix_to_yyyymmdd_hhmmss(ut) - << " with value = " << o << "\n"; - return false; - } + if((*o_it).ut == ut) return false; } } @@ -1016,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); @@ -1066,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; iname() ) { + if(var_name != obs_info->name()) { rej_var++; return; } @@ -805,10 +822,10 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, // Check if the observation quality flag is included in the list if(obs_qty_filt.n() && strcmp(obs_qty, "")) { bool qty_match = false; - for(i=0; imagic_str() << " versus " + << obs_info->magic_str() + << ", skipping observation with bad data value:\n" + << point_obs_to_string(hdr_arr, hdr_typ_str, hdr_sid_str, + hdr_ut, obs_qty, obs_arr, var_name) + << "\n"; rej_obs++; return; } @@ -845,6 +869,15 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, // Check if the observation's lat/lon is on the grid if(x < 0 || x >= gr.nx() || y < 0 || y >= gr.ny()) { + mlog << Debug(4) + << "For " << fcst_info->magic_str() << " versus " + << obs_info->magic_str() + << ", skipping observation off the grid where (x, y) = (" + << x << ", " << y << ") and grid (nx, ny) = (" << gr.nx() + << ", " << gr.ny() << "):\n" + << point_obs_to_string(hdr_arr, hdr_typ_str, hdr_sid_str, + hdr_ut, obs_qty, obs_arr, var_name) + << "\n"; rej_grd++; return; } @@ -861,12 +894,14 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, // Skip bad topography values if(is_bad_data(hdr_elv) || is_bad_data(topo)) { mlog << Debug(4) - << "Skipping observation due to missing topography values for " - << "[msg_typ:sid:lat:lon:elevation] = [" - << hdr_typ_str << ":" << hdr_sid_str << ":" - << hdr_lat << ":" << -1.0*hdr_lon << ":" - << hdr_elv << "] and model topography = " - << topo << ".\n"; + << "For " << fcst_info->magic_str() << " versus " + << obs_info->magic_str() + << ", skipping observation due to bad topography values " + << "where observation elevation = " << hdr_elv + << " and model topography = " << topo << ":\n" + << point_obs_to_string(hdr_arr, hdr_typ_str, hdr_sid_str, + hdr_ut, obs_qty, obs_arr, var_name) + << "\n"; rej_topo++; return; } @@ -874,14 +909,16 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, // Check the topography difference threshold if(!sfc_info.topo_use_obs_thresh.check(topo - hdr_elv)) { mlog << Debug(4) - << "Skipping observation for topography difference since " + << "For " << fcst_info->magic_str() << " versus " + << obs_info->magic_str() + << ", skipping observation due to topography difference " + << "where observation elevation (" << hdr_elv + << ") minus model topography (" << topo << ") = " << topo - hdr_elv << " is not " - << sfc_info.topo_use_obs_thresh.get_str() << " for " - << "[msg_typ:sid:lat:lon:elevation] = [" - << hdr_typ_str << ":" << hdr_sid_str << ":" - << hdr_lat << ":" << -1.0*hdr_lon << ":" - << hdr_elv << "] and model topography = " - << topo << ".\n"; + << sfc_info.topo_use_obs_thresh.get_str() << ":\n" + << point_obs_to_string(hdr_arr, hdr_typ_str, hdr_sid_str, + hdr_ut, obs_qty, obs_arr, var_name) + << "\n"; rej_topo++; return; } @@ -1099,6 +1136,14 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, } if(is_bad_data(fcst_v)) { + mlog << Debug(4) + << "For " << fcst_info->magic_str() << " versus " + << obs_info->magic_str() + << ", skipping observation due to bad data in the interpolated " + << "forecast value:\n" + << point_obs_to_string(hdr_arr, hdr_typ_str, hdr_sid_str, + hdr_ut, obs_qty, obs_arr, var_name) + << "\n"; inc_count(rej_fcst, i, j, k); continue; } @@ -1113,6 +1158,13 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, hdr_lat, hdr_lon, obs_x, obs_y, hdr_ut, obs_lvl, obs_hgt, fcst_v, obs_v, obs_qty, cmn_v, csd_v, wgt_v)) { + mlog << Debug(4) + << "For " << fcst_info->magic_str() << " versus " + << obs_info->magic_str() + << ", skipping observation since it is a duplicate:\n" + << point_obs_to_string(hdr_arr, hdr_typ_str, hdr_sid_str, + hdr_ut, obs_qty, obs_arr, var_name) + << "\n"; inc_count(rej_dup, i, j, k); } @@ -1298,6 +1350,7 @@ PairDataPoint subset_pairs(const PairDataPoint &pd, // Allocate memory for output pairs out_pd.extend(pd.n_obs); + out_pd.set_climo_cdf_info(pd.cdf_info); bool cmn_flag = set_climo_flag(pd.f_na, pd.cmn_na); bool csd_flag = set_climo_flag(pd.f_na, pd.csd_na); @@ -1365,6 +1418,8 @@ void subset_wind_pairs(const PairDataPoint &pd_u, const PairDataPoint &pd_v, out_pd_v.erase(); out_pd_u.extend(pd_u.n_obs); out_pd_v.extend(pd_v.n_obs); + out_pd_u.set_climo_cdf_info(pd_u.cdf_info); + out_pd_v.set_climo_cdf_info(pd_v.cdf_info); bool cmn_flag = set_climo_flag(pd_u.f_na, pd_u.cmn_na) && set_climo_flag(pd_v.f_na, pd_v.cmn_na); @@ -1452,6 +1507,7 @@ PairDataPoint subset_climo_cdf_bin(const PairDataPoint &pd, // Allocate memory for output pairs out_pd.extend(pd.n_obs); + out_pd.set_climo_cdf_info(pd.cdf_info); bool cmn_flag = set_climo_flag(pd.f_na, pd.cmn_na); bool csd_flag = set_climo_flag(pd.f_na, pd.csd_na); @@ -1494,6 +1550,36 @@ PairDataPoint subset_climo_cdf_bin(const PairDataPoint &pd, return(out_pd); } +//////////////////////////////////////////////////////////////////////// + +// Write the point observation in the MET point format for logging +ConcatString point_obs_to_string(float *hdr_arr, const char *hdr_typ_str, + const char *hdr_sid_str, unixtime hdr_ut, + const char *obs_qty, float *obs_arr, + const char *var_name) { + ConcatString obs_cs, name; + + if((var_name != 0) && (0 < strlen(var_name))) name = var_name; + else name = obs_arr[1]; + + // + // Write the 11-column MET point format: + // Message_Type Station_ID Valid_Time(YYYYMMDD_HHMMSS) + // Lat(Deg North) Lon(Deg East) Elevation(msl) + // Var_Name(or GRIB_Code) Level Height(msl or agl) + // QC_String Observation_Value + // + obs_cs << " " + << hdr_typ_str << " " << hdr_sid_str << " " + << unix_to_yyyymmdd_hhmmss(hdr_ut) << " " + << hdr_arr[0] << " " << -1.0*hdr_arr[1] << " " + << hdr_arr[2] << " " << name << " " + << obs_arr[2] << " " << obs_arr[3] << " " + << obs_qty << " " << obs_arr[4]; + + return(obs_cs); +} + //////////////////////////////////////////////////////////////////////// // // End miscellaneous functions diff --git a/met/src/libcode/vx_statistics/pair_data_point.h b/met/src/libcode/vx_statistics/pair_data_point.h index 4c04625bd6..84494cb2d0 100644 --- a/met/src/libcode/vx_statistics/pair_data_point.h +++ b/met/src/libcode/vx_statistics/pair_data_point.h @@ -205,6 +205,8 @@ class VxPairDataPoint { void set_interp(int i_interp, InterpMthd mthd, int width, GridTemplateFactory::GridTemplates shape); + void set_climo_cdf_info(const ClimoCDFInfo &); + void set_msg_typ_sfc(const StringArray &); void set_msg_typ_lnd(const StringArray &); void set_msg_typ_wtr(const StringArray &); @@ -254,6 +256,13 @@ extern void subset_wind_pairs(const PairDataPoint &, extern PairDataPoint subset_climo_cdf_bin(const PairDataPoint &, const ThreshArray &, int i_bin); +// Write the point observation in the MET point format for logging +extern ConcatString point_obs_to_string( + float *hdr_arr, const char *hdr_typ_str, + const char *hdr_sid_str, unixtime hdr_ut, + const char *obs_qty, float *obs_arr, + const char *var_name); + //////////////////////////////////////////////////////////////////////// #endif // __PAIR_DATA_POINT_H__ diff --git a/met/src/libcode/vx_statistics/read_climo.cc b/met/src/libcode/vx_statistics/read_climo.cc index e292a381b3..27d79cbb42 100644 --- a/met/src/libcode/vx_statistics/read_climo.cc +++ b/met/src/libcode/vx_statistics/read_climo.cc @@ -245,7 +245,7 @@ void read_climo_file(const char *climo_file, GrdFileType ctype, // Regrid, if needed if(!(mtddf->grid() == 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 + + + diff --git a/test/xml/unit_plot_data_plane.xml b/test/xml/unit_plot_data_plane.xml index 1eb25e614c..c655ff7116 100644 --- a/test/xml/unit_plot_data_plane.xml +++ b/test/xml/unit_plot_data_plane.xml @@ -343,47 +343,6 @@ - - &MET_BIN;/plot_data_plane - \ - &DATA_DIR_MODEL;/cesm/ssmi.ifrac.gx3v5.2003.nc \ - &OUTPUT_DIR;/plot_data_plane/CESM_ssmi.ifrac.gx3v5.2003.ps \ - 'name="ifrac"; level="(0,*,*)";' \ - -v 4 - - - &OUTPUT_DIR;/plot_data_plane/CESM_ssmi.ifrac.gx3v5.2003.ps - - - - - &MET_BIN;/plot_data_plane - \ - &DATA_DIR_MODEL;/cesm/b.e11.B20TRC5CNBDRD.f09_g16.002.cice.h1.aice_d_nh.19200101-20051231.nc \ - &OUTPUT_DIR;/plot_data_plane/CESM_b.e11.B20TRC5CNBDRD.f09_g16.002.cice.h1.aice_d_nh.19200101-20051231.ps \ - 'name="aice_d"; level="(0,*,*)";' \ - -v 4 - - - &OUTPUT_DIR;/plot_data_plane/CESM_b.e11.B20TRC5CNBDRD.f09_g16.002.cice.h1.aice_d_nh.19200101-20051231.ps - - - - - &MET_BIN;/plot_data_plane \