Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfix #2412 develop climo #2422

Merged
merged 4 commits into from
Jan 25, 2023
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions docs/Users_Guide/config_options.rst
Original file line number Diff line number Diff line change
@@ -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".
266 changes: 266 additions & 0 deletions internal/test_unit/config/GridStatConfig_climo_wrap_year
Original file line number Diff line number Diff line change
@@ -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 = NONE;
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";

////////////////////////////////////////////////////////////////////////////////
31 changes: 31 additions & 0 deletions internal/test_unit/xml/unit_climatology_2.5deg.xml
Original file line number Diff line number Diff line change
@@ -50,4 +50,35 @@
</output>
</test>

<!-- MET#2412: Test interpolation of monthly climo between December 15 and January 15 -->

<test name="climatology_GRID_STAT_WRAP_YEAR_2.5DEG">
<exec>&MET_BIN;/grid_stat</exec>
<env>
<pair><name>OUTPUT_PREFIX</name> <value>WRAP_YEAR_CLIMO_2.5DEG</value></pair>
<pair><name>DAY_INTERVAL</name> <value>31</value></pair>
<pair><name>HOUR_INTERVAL</name> <value>6</value></pair>
<pair><name>CLIMO_MEAN_FILE_LIST</name>
<value>"&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_mean.19591215",
"&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_mean.19590115"
</value>
</pair>
<pair><name>CLIMO_STDEV_FILE_LIST</name>
<value>"&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_stdv.19591215",
"&DATA_DIR_CLIMO;/NCEP_2.5deg/pgba_stdv.19590115"
</value>
</pair>
</env>
<param> \
&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
</param>
<output>
<stat>&OUTPUT_DIR;/climatology_2.5deg/grid_stat_WRAP_YEAR_CLIMO_2.5DEG_000000L_20201225_120000V.stat</stat>
<grid_nc>&OUTPUT_DIR;/climatology_2.5deg/grid_stat_WRAP_YEAR_CLIMO_2.5DEG_000000L_20201225_120000V_pairs.nc</grid_nc>
</output>
</test>

</met_test>
65 changes: 65 additions & 0 deletions src/basic/vx_cal/doyhms_to_unix.cc
Original file line number Diff line number Diff line change
@@ -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 );

}


////////////////////////////////////////////////////////////////////////
10 changes: 10 additions & 0 deletions src/basic/vx_cal/vx_cal.h
Original file line number Diff line number Diff line change
@@ -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);
10 changes: 6 additions & 4 deletions src/basic/vx_util/interp_util.cc
Original file line number Diff line number Diff line change
@@ -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 "
119 changes: 44 additions & 75 deletions src/libcode/vx_statistics/read_climo.cc
Original file line number Diff line number Diff line change
@@ -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 &regrid_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<n_climo; i++) {

// Store current valid time string
cur_ut_cs = unix_to_yyyymmdd_hhmmss(cur_dpa[i].valid());

// Compute the current month and day
unix_to_mdyhms(cur_dpa[i].valid(), clm_mon, clm_day, clm_yr,
clm_hr, clm_min, clm_sec);

// Recompute the unixtime to check the hour of the day.
// Use the valid YYYYMMDD and climo HHMMSS.
ut = mdyhms_to_unix(vld_mon, vld_day, vld_yr,
clm_hr, clm_min, clm_sec);

// Check the hour time step.
if(!is_bad_data(hour_ts) && 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 file: "
<< climo_file << "\n";
for(i=0; i<n_clm; i++) {

// Store climo time string
clm_ut_cs = unix_to_yyyymmdd_hhmmss(clm_dpa[i].valid());

// Compute day and hour time offsets in seconds
day_diff_sec = day_of_year_diff(vld_ut, clm_dpa[i].valid()) * sec_per_day;
hms_diff_sec = sec_of_day_diff (vld_ut, clm_dpa[i].valid());

// Check the day time step
if(!is_bad_data(day_ts) && abs(day_diff_sec) >= 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]),