From 428a9c9a7da9535fe95ce6c57608c1885923e145 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Mon, 23 Jan 2023 09:58:53 -0700 Subject: [PATCH 1/4] Per #2412, port changes over from main_v11.0 to the develop branch. --- src/basic/vx_cal/doyhms_to_unix.cc | 65 +++++++++++++ src/basic/vx_cal/vx_cal.h | 10 ++ src/basic/vx_util/interp_util.cc | 10 +- src/libcode/vx_statistics/read_climo.cc | 119 +++++++++--------------- 4 files changed, 125 insertions(+), 79 deletions(-) diff --git a/src/basic/vx_cal/doyhms_to_unix.cc b/src/basic/vx_cal/doyhms_to_unix.cc index 8c03137841..2c130131bb 100644 --- a/src/basic/vx_cal/doyhms_to_unix.cc +++ b/src/basic/vx_cal/doyhms_to_unix.cc @@ -97,8 +97,28 @@ return ( s ); } + +//////////////////////////////////////////////////////////////////////// + + +int unix_to_day_of_year(unixtime u) { + +int mon, day, yr, hr, min, sec; + +unix_to_mdyhms(u, mon, day, yr, hr, min, sec); + +int sec_of_year = mdyhms_to_unix(mon, day, 1970, 0, 0, 0); + +int sec_per_day = 60 * 60 * 24; + +return ( sec_of_year / sec_per_day ); + +} + + //////////////////////////////////////////////////////////////////////// + long unix_to_long_yyyymmddhh(unixtime u) { int mon, day, yr, hr, min, sec; @@ -112,4 +132,49 @@ return ( yyyymmddhh ); } + +//////////////////////////////////////////////////////////////////////// + + +int sec_of_day_diff(unixtime ut1, unixtime ut2) { + +int sec_per_day = 60 * 60 * 24; + +int s1 = unix_to_sec_of_day(ut1); + +int s2 = unix_to_sec_of_day(ut2); + +int dt = s2 - s1; + + if ( dt < -1 * sec_per_day/2 ) dt += sec_per_day; +else if ( dt > sec_per_day/2 ) dt -= sec_per_day; + +return ( dt ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +int day_of_year_diff(unixtime ut1, unixtime ut2) { + + +int sec_per_day = 60 * 60 * 24; +int day_per_year = 365; + +int d1 = unix_to_day_of_year(ut1); + +int d2 = unix_to_day_of_year(ut2); + +int dt = d2 - d1; + + if ( dt < -1 * day_per_year/2.0 ) dt += day_per_year; +else if ( dt > day_per_year/2.0 ) dt -= day_per_year; + +return ( dt ); + +} + + //////////////////////////////////////////////////////////////////////// diff --git a/src/basic/vx_cal/vx_cal.h b/src/basic/vx_cal/vx_cal.h index b06b6fc7fe..74353e1e91 100644 --- a/src/basic/vx_cal/vx_cal.h +++ b/src/basic/vx_cal/vx_cal.h @@ -81,6 +81,8 @@ extern int hms_to_sec (int hour, int min, int sec); extern int unix_to_sec_of_day (unixtime u); +extern int unix_to_day_of_year (unixtime u); + extern long unix_to_long_yyyymmddhh (unixtime u); // Parse time strings @@ -128,6 +130,14 @@ extern ConcatString sec_to_timestring(int); //////////////////////////////////////////////////////////////////////// +extern int sec_of_day_diff (unixtime ut1, unixtime ut2); + +extern int day_of_year_diff (unixtime ut1, unixtime ut2); + + +//////////////////////////////////////////////////////////////////////// + + extern bool is_datestring(const char * text); extern bool is_yyyymmdd(const char * text); diff --git a/src/basic/vx_util/interp_util.cc b/src/basic/vx_util/interp_util.cc index bdaf7a573d..b123e972ce 100644 --- a/src/basic/vx_util/interp_util.cc +++ b/src/basic/vx_util/interp_util.cc @@ -1280,14 +1280,16 @@ DataPlane valid_time_interp(const DataPlane &in1, const DataPlane &in2, double w1 = bad_data_double; double w2 = bad_data_double; - // Store min and max valid times. + // Store min and max valid times dp1 = (in1.valid() <= in2.valid() ? in1 : in2); dp2 = (in1.valid() > in2.valid() ? in1 : in2); + // Check for matching valid times + if(dp1.valid() == dp2.valid()) return(dp1); + // Range check the times - if(dp1.valid() > to_ut || dp2.valid() < to_ut || - dp1.valid() == dp2.valid()) { - mlog << Error << "\time_interp() -> " + if(dp1.valid() > to_ut || dp2.valid() < to_ut) { + mlog << Error << "\ntime_interp() -> " << "the interpolation time " << unix_to_yyyymmdd_hhmmss(to_ut) << " must fall between the input times: " << unix_to_yyyymmdd_hhmmss(dp1.valid()) << " and " diff --git a/src/libcode/vx_statistics/read_climo.cc b/src/libcode/vx_statistics/read_climo.cc index cb6a0efb17..142833337f 100644 --- a/src/libcode/vx_statistics/read_climo.cc +++ b/src/libcode/vx_statistics/read_climo.cc @@ -150,21 +150,19 @@ void read_climo_file(const char *climo_file, GrdFileType ctype, int day_ts, int hour_ts, const Grid &vx_grid, const RegridInfo ®rid_info, DataPlaneArray &dpa) { + Met2dDataFileFactory mtddf_factory; Met2dDataFile *mtddf = (Met2dDataFile *) 0; VarInfoFactory info_factory; VarInfo *info = (VarInfo *) 0; - DataPlaneArray cur_dpa; + DataPlaneArray clm_dpa; DataPlane dp; - int n_climo, i; - int vld_mon, vld_day, vld_yr, vld_hr, vld_min, vld_sec; - int clm_mon, clm_day, clm_yr, clm_hr, clm_min, clm_sec; - unixtime ut; + int i, n_clm, day_diff_sec, hms_diff_sec; - ConcatString cur_ut_cs; + ConcatString clm_ut_cs; // Allocate memory for data file if(!(mtddf = mtddf_factory.new_met_2d_data_file(climo_file, ctype))) { @@ -179,94 +177,65 @@ void read_climo_file(const char *climo_file, GrdFileType ctype, info->set_dict(*dict); // Read data planes - n_climo = mtddf->data_plane_array(*info, cur_dpa); - - // Compute the valid month and day - unix_to_mdyhms(vld_ut, vld_mon, vld_day, vld_yr, - vld_hr, vld_min, vld_sec); + n_clm = mtddf->data_plane_array(*info, clm_dpa); // Loop through matching records - for(i=0; i= hour_ts) { - mlog << Debug(3) << "Skipping the " << cur_ut_cs << " \"" - << info->magic_str() - << "\" climatology field since the time offset (" - << labs(ut - vld_ut) - << " seconds) >= the \"" << conf_key_hour_interval - << "\" entry (" << hour_ts << " seconds) from file: " - << climo_file << "\n"; + for(i=0; i= day_ts) { + mlog << Debug(3) << "Skipping " << clm_ut_cs << " \"" << info->magic_str() + << "\" climatology field with " << day_diff_sec / sec_per_day + << " day offset (" << conf_key_day_interval << " = " + << day_ts / sec_per_day << ") from file \"" + << climo_file << "\".\n"; continue; } - // Recompute the unixtime to check the day of the year. - // Use the valid YYYY and climo MMDD_HHMMSS. - ut = mdyhms_to_unix(clm_mon, clm_day, vld_yr, - clm_hr, clm_min, clm_sec); - - // Check the day time step. - if(!is_bad_data(day_ts)) { - - // For daily climatology, check the hour timestep. - if(day_ts <= 3600*24 && labs(ut - vld_ut) >= hour_ts) { - mlog << Debug(3) << "Skipping the " << cur_ut_cs << " \"" - << info->magic_str() - << "\" climatology field since the time offset (" - << labs(ut - vld_ut) << " seconds) >= the \"" - << conf_key_hour_interval << "\" entry (" << hour_ts - << " seconds) from daily climatology file: " - << climo_file << "\n"; - continue; - } - // For non-daily climatology, check the day timestep. - else if(labs(ut - vld_ut) >= day_ts) { - mlog << Debug(3) << "Skipping the " << cur_ut_cs << " \"" - << info->magic_str() - << "\" climatology field since the time offset (" - << labs(ut - vld_ut) << " seconds) >= the \"" - << conf_key_day_interval << "\" entry (" << day_ts - << " seconds) from file: " << climo_file << "\n"; - continue; - } + // Check the hour time step + if(!is_bad_data(hour_ts) && abs(hms_diff_sec) >= hour_ts) { + mlog << Debug(3) << "Skipping " << clm_ut_cs << " \"" << info->magic_str() + << "\" climatology field with " << (double) hms_diff_sec / sec_per_hour + << " hour offset (" << conf_key_hour_interval << " = " + << hour_ts / sec_per_hour << ") from file \"" + << climo_file << "\".\n"; + continue; } + // Compute climo timestamp relative to the valid time + unixtime clm_vld_ut = vld_ut + day_diff_sec + hms_diff_sec; + // Print log message for matching record - mlog << Debug(4) - << "Found matching " << cur_ut_cs << " \"" - << info->magic_str() << "\" climatology field in file \"" - << climo_file << "\".\n"; + mlog << Debug(3) << "Storing " << clm_ut_cs << " \"" << info->magic_str() + << "\" climatology field with " << day_diff_sec / sec_per_day + << " day, " << (double) hms_diff_sec / sec_per_hour << " hour offset as time " + << unix_to_yyyymmdd_hhmmss(clm_vld_ut) << " from file \"" + << climo_file << "\".\n"; // Regrid, if needed if(!(mtddf->grid() == vx_grid)) { - mlog << Debug(2) - << "Regridding the " << cur_ut_cs << " \"" + mlog << Debug(2) << "Regridding " << clm_ut_cs << " \"" << info->magic_str() << "\" climatology field to the verification grid.\n"; - dp = met_regrid(cur_dpa[i], mtddf->grid(), vx_grid, + dp = met_regrid(clm_dpa[i], mtddf->grid(), vx_grid, regrid_info); } else { - dp = cur_dpa[i]; + dp = clm_dpa[i]; } - // Set the climo time as valid YYYY and climo MMDD_HHMMSS. - dp.set_valid(ut); + // Set the climo time relative to the valid time + dp.set_valid(clm_vld_ut); // Store the match - dpa.add(dp, cur_dpa.lower(i), cur_dpa.upper(i)); + dpa.add(dp, clm_dpa.lower(i), clm_dpa.upper(i)); } // end for i @@ -376,7 +345,7 @@ DataPlaneArray climo_time_interp(const DataPlaneArray &dpa, int day_ts, // For equality, do a single time interpolation. if(prv_hms == nxt_hms) { - ut = (vld_ut / sec_per_day) + prv_hms; + ut = (vld_ut / sec_per_day)*sec_per_day + prv_hms; interp_dpa.add(climo_hms_interp( dpa, it->second, ut, mthd), dpa.lower(it->second[0]), From 51c2757e69da247430e925d0ab083b72c1155864 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Tue, 24 Jan 2023 16:35:26 -0700 Subject: [PATCH 2/4] Per #2412 add unit test. --- .../config/GridStatConfig_climo_wrap_year | 266 ++++++++++++++++++ .../test_unit/xml/unit_climatology_2.5deg.xml | 31 ++ 2 files changed, 297 insertions(+) create mode 100644 internal/test_unit/config/GridStatConfig_climo_wrap_year diff --git a/internal/test_unit/config/GridStatConfig_climo_wrap_year b/internal/test_unit/config/GridStatConfig_climo_wrap_year new file mode 100644 index 0000000000..0c048db507 --- /dev/null +++ b/internal/test_unit/config/GridStatConfig_climo_wrap_year @@ -0,0 +1,266 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Grid-Stat configuration file. +// +// For additional information, please see the MET User's Guide. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Output model name to be written +// +model = "GFSANL"; + +// +// Output description to be written +// May be set separately in each "obs.field" entry +// +desc = "NA"; + +// +// Output observation type to be written +// +obtype = "GFSANL"; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Verification grid +// +regrid = { + to_grid = NONE; + method = NEAREST; + width = 1; + vld_thresh = 0.5; + shape = SQUARE; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// May be set separately in each "field" entry +// +censor_thresh = []; +censor_val = []; +mpr_column = []; +mpr_thresh = []; +cat_thresh = []; +cnt_thresh = [ NA ]; +cnt_logic = UNION; +wind_thresh = [ NA ]; +wind_logic = UNION; +eclv_points = 0.05; +nc_pairs_var_name = ""; +nc_pairs_var_suffix = ""; +hss_ec_value = NA; +rank_corr_flag = FALSE; + +// +// Forecast and observation fields to be verified +// +fcst = { + + name = "TMP"; + level = "P500"; + + field = [ + { set_attr_init = "20201225_12"; set_attr_valid = "20201225_12"; nc_pairs_var_suffix = "20201225_12"; }, + { set_attr_init = "20210105_12"; set_attr_valid = "20210105_12"; nc_pairs_var_suffix = "20210105_12"; }, + { set_attr_init = "20201225_03"; set_attr_valid = "20201225_03"; nc_pairs_var_suffix = "20201225_03"; }, + { set_attr_init = "20201225_21"; set_attr_valid = "20201225_21"; nc_pairs_var_suffix = "20201225_21"; } + ]; +} +obs = fcst; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Climatology data +// +climo_mean = { + + field = [ + { name = "TMP"; level = "P500"; }, + { name = "TMP"; level = "P500"; }, + { name = "TMP"; level = "P500"; }, + { name = "TMP"; level = "P500"; } + ]; + + file_name = [ ${CLIMO_MEAN_FILE_LIST} ]; + + regrid = { + method = BILIN; + width = 2; + vld_thresh = 0.5; + } + + time_interp_method = DW_MEAN; + day_interval = ${DAY_INTERVAL}; + hour_interval = ${HOUR_INTERVAL}; +} + +climo_stdev = climo_mean; +climo_stdev = { + file_name = [ ${CLIMO_STDEV_FILE_LIST} ]; +} + +// +// May be set separately in each "obs.field" entry +// +climo_cdf = { + cdf_bins = 1; + center_bins = TRUE; + write_bins = FALSE; + direct_prob = FALSE; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Verification masking regions +// +mask = { + grid = [ "FULL" ]; + poly = []; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Confidence interval settings +// +ci_alpha = [ 0.05 ]; + +boot = { + interval = PCTILE; + rep_prop = 1.0; + n_rep = 0; + rng = "mt19937"; + seed = ""; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Data smoothing methods +// +interp = { + field = BOTH; + vld_thresh = 1.0; + shape = SQUARE; + + type = [ + { + method = NEAREST; + width = 1; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Neighborhood methods +// +nbrhd = { + width = [ 1 ]; + cov_thresh = [ >=0.5 ]; + vld_thresh = 1.0; + shape = SQUARE; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Fourier decomposition +// +fourier = { + wave_1d_beg = []; + wave_1d_end = []; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1, 2, 5 ]; + dy = [ 1, 3, 5 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Distance Map statistics +// May be set separately in each "obs.field" entry +// +distance_map = { + baddeley_p = 2; + baddeley_max_dist = NA; + fom_alpha = 0.1; + zhu_weight = 0.5; + beta_value(n) = n * n / 2.0; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Statistical output types +// +output_flag = { + fho = NONE; + ctc = NONE; + cts = NONE; + mctc = NONE; + mcts = NONE; + cnt = NONE; + sl1l2 = STAT; + sal1l2 = STAT; + vl1l2 = NONE; + val1l2 = NONE; + vcnt = NONE; + pct = NONE; + pstd = NONE; + pjc = NONE; + prc = NONE; + eclv = NONE; + nbrctc = NONE; + nbrcts = NONE; + nbrcnt = NONE; + grad = NONE; + dmap = NONE; + seeps = NONE; +} + +// +// NetCDF matched pairs output file +// +nc_pairs_flag = { + latlon = TRUE; + raw = FALSE; + diff = FALSE; + climo = TRUE; + climo_cdp = FALSE; + weight = FALSE; + nbrhd = FALSE; + fourier = FALSE; + gradient = FALSE; + distance_map = FALSE; + apply_mask = FALSE; +} + +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + +//////////////////////////////////////////////////////////////////////////////// + +grid_weight_flag = NONE; +tmp_dir = "/tmp"; +output_prefix = "${OUTPUT_PREFIX}"; +version = "V11.0.0"; + +//////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/xml/unit_climatology_2.5deg.xml b/internal/test_unit/xml/unit_climatology_2.5deg.xml index 7c8c0b300a..148fbf3ea8 100644 --- a/internal/test_unit/xml/unit_climatology_2.5deg.xml +++ b/internal/test_unit/xml/unit_climatology_2.5deg.xml @@ -50,4 +50,35 @@ + + + + &MET_BIN;/grid_stat + + OUTPUT_PREFIX YEAR_WRAP_CLIMO_2.5DEG + DAY_INTERVAL 31 + HOUR_INTERVAL 6 + CLIMO_MEAN_FILE_LIST + "&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_mean.19591215", + "&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_mean.19590115" + + + CLIMO_STDEV_FILE_LIST + "&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_stdv.19591215", + "&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_stdv.19590115" + + + + \ + &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_1200_000.grb2 \ + &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_1200_000.grb2 \ + &CONFIG_DIR;/GridStatConfig_climo_wrap_year \ + -outdir &OUTPUT_DIR;/climatology_2.5deg -v 3 + + + &OUTPUT_DIR;/climatology_2.5deg/grid_stat_WRAP_YEAR_CLIMO_2.5DEG_000000L_20201225_120000V.stat + &OUTPUT_DIR;/climatology_2.5deg/grid_stat_WRAP_YEAR_CLIMO_2.5DEG_000000L_20201225_120000V_pairs.nc + + + From abecaf081393349fc3af70905a26d91611e863bb Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Wed, 25 Jan 2023 09:55:39 -0700 Subject: [PATCH 3/4] Per #2412, refine the new climo unit test to only write SAL1L2 output. --- docs/Users_Guide/config_options.rst | 5 +++-- internal/test_unit/config/GridStatConfig_climo_wrap_year | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/Users_Guide/config_options.rst b/docs/Users_Guide/config_options.rst index 5cc2432fd6..66c2926af0 100644 --- a/docs/Users_Guide/config_options.rst +++ b/docs/Users_Guide/config_options.rst @@ -1318,8 +1318,9 @@ of several entires defining the climatology file names and fields to be used. * The "hour_interval" entry is an integer specifying the spacing in hours of the climatology data for each day. This should be set between 0 and 24, - with 6 and 12 being common choices. Use "NA" if the timing of the - climatology data should not be checked. + with 6 and 12 being common choices. For example, use 6 for climatology data + with 4 times per day, such as 00Z, 06Z, 12Z, and 18Z. Use "NA" if the timing + of the climatology data should not be checked. * The "day_interval" and "hour_interval" entries replace the deprecated entries "match_month", "match_day", and "time_step". diff --git a/internal/test_unit/config/GridStatConfig_climo_wrap_year b/internal/test_unit/config/GridStatConfig_climo_wrap_year index 0c048db507..38e8ffcd0c 100644 --- a/internal/test_unit/config/GridStatConfig_climo_wrap_year +++ b/internal/test_unit/config/GridStatConfig_climo_wrap_year @@ -216,7 +216,7 @@ output_flag = { mctc = NONE; mcts = NONE; cnt = NONE; - sl1l2 = STAT; + sl1l2 = NONE; sal1l2 = STAT; vl1l2 = NONE; val1l2 = NONE; From 90308564581c1cf8f87d4c854e62a85ae0043b2c Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Wed, 25 Jan 2023 11:46:28 -0700 Subject: [PATCH 4/4] Per #2412, fix OUTPUT_PREFIX naming inconsistency in the new unit test. --- internal/test_unit/xml/unit_climatology_2.5deg.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/internal/test_unit/xml/unit_climatology_2.5deg.xml b/internal/test_unit/xml/unit_climatology_2.5deg.xml index 148fbf3ea8..ed534d0191 100644 --- a/internal/test_unit/xml/unit_climatology_2.5deg.xml +++ b/internal/test_unit/xml/unit_climatology_2.5deg.xml @@ -55,7 +55,7 @@ &MET_BIN;/grid_stat - OUTPUT_PREFIX YEAR_WRAP_CLIMO_2.5DEG + OUTPUT_PREFIX WRAP_YEAR_CLIMO_2.5DEG DAY_INTERVAL 31 HOUR_INTERVAL 6 CLIMO_MEAN_FILE_LIST