Skip to content

Commit

Permalink
Feature 1746 wavelet stat (#1851)
Browse files Browse the repository at this point in the history 
* For issue #1746 modified code to allow users to pass in an empty list (or NA) for forecast and observation thresholds in order to skip applying the threhsolds, but it will still compute stats with the raw fields. SL

* For issue #1746, added new unit test that uses a config file that has empty lists for the forecast and observation thresholds. SL

* For issue #1746 Added some content related to allowing users to set forecast and observation cat thresholds to an empty list in order to skip the binary masking (and consider all grid-points for stats). SL

* Per #1746, cleaning up for consistent indentation.

* Per #1746, cleaning up for consistent indentation.

* Per #1746, add a revision history note, update the plotting range in the postscript output to be [-n,n] where n is the maximum value of the maximum absolute difference and 1.0, and also fix a bug. When the NA threshold comes AFTER a real threshold, the resulting data and difference values were not being updated.

* Per #1746, change the Wavelet-Stat config file values in the the wvlt_plot dictionary by setting plot_min = plot_max = 0.0. That enables the default logic of the tool to take effect. Choose the plotting range of the wavelet plots as [-n,n], where n is the maximum of 1.0 and the maximum absolute difference.

* Per #1746, used apply_fcst_thresh where it should have been apply_obs_thresh.

* Per issue #1746, modified some content related to users being able to skip applying the categorical threhsolds by putting an empty list or NA in the configuration file. SL

* Per issue #1746 Added some warnings if the forecast threshold is set to NA but the observation threshold is not NA (a numeric threshold) and vice versa. SL

* Per #1746, fix a couple of typos and tweak wording in the wavelet-stat chapter.

* Per #1746, loop over each pair of fcst/obs thresholds to check for inconsistent use of the NA threshold type.

* Per #1746, a bit of code cleanup replacing calls to n_elements() with n() to make the code more concise.

* Per #1746, need to reinitialize apply_fcst_thresh and apply_obs_thresh to true inside the loop since the NA threshold can appear anywhere in the list of thresholds.

Co-authored-by: Seth Linden <linden@kiowa.rap.ucar.edu>
Co-authored-by: John Halley Gotway <johnhg@ucar.edu>
  • Loading branch information
@sethlinden @JohnHalleyGotway
3 people authored Jul 15, 2021
1 parent 8ba6df2 commit 599181a
Show file tree
Hide file tree
Showing 10 changed files with 283 additions and 50 deletions.
4 changes: 2 additions & 2 deletions met/data/config/WaveletStatConfig_default
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,8 +131,8 @@ obs_raw_plot = {

wvlt_plot = {
color_table = "MET_BASE/colortables/NCL_colortables/BlWhRe.ctable";
plot_min = -1.0;
plot_max = 1.0;
plot_min = 0.0;
plot_max = 0.0;
}

////////////////////////////////////////////////////////////////////////////////
Expand Down
19 changes: 18 additions & 1 deletion met/docs/Users_Guide/wavelet-stat.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ The method

The Intensity Scale approach can be summarized in the following 5 steps:

1. For each threshold, the forecast and observation fields are transformed into binary fields: where the grid-point precipitation value meets the threshold criteria it is assigned 1, where the threshold criteria are not met it is assigned 0. :numref:`wavelet-stat_NIMROD_3h_fcst` illustrates an example of a forecast and observation fields, and their corresponding binary fields for a threshold of 1mm/h. This case shows an intense storm of the scale of 160 km displaced almost its entire length. The displacement error is clearly visible from the binary field difference and the contingency table image obtained for the same threshold :numref:`contingency_table_counts`.
1. For each threshold, the forecast and observation fields are transformed into binary fields: where the grid-point precipitation value meets the threshold criteria it is assigned 1, where the threshold criteria are not met it is assigned 0. This can also be done with no thresholds indicated at all and in that case the grid-point values are not transformed to binary fields and instead the raw data is used as is for statistics. :numref:`wavelet-stat_NIMROD_3h_fcst` illustrates an example of a forecast and observation fields, and their corresponding binary fields for a threshold of 1mm/h. This case shows an intense storm of the scale of 160 km displaced almost its entire length. The displacement error is clearly visible from the binary field difference and the contingency table image obtained for the same threshold :numref:`contingency_table_counts`.

2. The binary forecast and observation fields obtained from the thresholding are then decomposed into the sum of components on different scales, by using a 2D Haar wavelet filter (:numref:`wavelet-stat_NIMROD_diff`). Note that the scale components are fields, and their sum adds up to the original binary field. For a forecast defined over square domain of **2ⁿ x 2ⁿ** grid-points, the scale components are **n+1: n** mother wavelet components + the largest father wavelet (or scale-function) component. The **n** mother wavelet components have resolution equal to **1, 2, 4, ...** :math:`\mathbf{2^{n-1}}` grid-points. The largest father wavelet component is a constant field over the **2ⁿ x 2ⁿ** grid-point domain with value equal to the field mean.

Expand Down Expand Up @@ -230,6 +230,23 @@ _______________________

.. code-block:: none
// Empty list of thresholds
cat_thresh = [];
// Or explicitly set the NA threshold type
cat_thresh = [>0.0, >=5.0, NA];
The **cat_thresh** option defines an array of thresholds for each field defined in the fcst and obs dictionaries. The number of forecast and observation categorical thresholds must match. If set to an empty list, the thresholds will not be applied (no binary masking) and all the raw grid-point values will be used for downstream statistics.

If the array of thresholds is an empty list, the application will set the threshold to NA internally and skip applying the thresholds. If the threshold is set to NA explicitly in the list, the application will also skip applying the threshold.

Since the application has the ability to loop through multiple thresholds (for multiple fields), a user can include NA in the list of thresholds to produce statistics for the raw data values for the given field.

_______________________

.. code-block:: none
grid_decomp_flag = AUTO;
tile = {
Expand Down
8 changes: 4 additions & 4 deletions met/scripts/config/WaveletStatConfig_APCP_12
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ fcst = {
{
name = "APCP";
level = [ "A12" ];
cat_thresh = [ >0.0, >=5.0 ];
cat_thresh = [ >0.0, >=5.0, NA ];
}
];
}
Expand All @@ -59,7 +59,7 @@ obs = {
{
name = "APCP_12";
level = [ "(*,*)" ];
cat_thresh = [ >0.0, >=5.0 ];
cat_thresh = [ >0.0, >=5.0, NA ];
}
];
}
Expand Down Expand Up @@ -137,8 +137,8 @@ obs_raw_plot = {

wvlt_plot = {
color_table = "MET_BASE/colortables/NCL_colortables/BlWhRe.ctable";
plot_min = -1.0;
plot_max = 1.0;
plot_min = 0.0;
plot_max = 0.0;
}

////////////////////////////////////////////////////////////////////////////////
Expand Down
83 changes: 54 additions & 29 deletions met/src/tools/core/wavelet_stat/wavelet_stat.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,8 @@
// 010 02/25/15 Halley Gotway Add automated regridding.
// 011 05/15/17 Prestopnik P. Add shape to regrid options.
// 012 04/08/19 Halley Gotway Add percentile thresholds.
// 012 04/01/19 Fillmore Add FCST and OBS units.
// 013 04/01/19 Fillmore Add FCST and OBS units.
// 014 07/09/21 Linden MET #1746 Skip thresholding.
//
////////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -113,8 +114,8 @@ static double mean_array(double *, int);

static void plot_ps_raw(const DataPlane &, const DataPlane &,
const DataPlane &, const DataPlane &, int);
static void plot_ps_wvlt(const double *, int, int, int, ISCInfo &,
int, int);
static void plot_ps_wvlt(const double *, const double, int, int, int,
ISCInfo &, int, int);
static double compute_percentage(double, double);

static void set_plot_dims(int, int);
Expand Down Expand Up @@ -419,7 +420,7 @@ void process_scores() {
// Allocate memory for ISCInfo objects sized as [n_tile][n_thresh]
isc_info = new ISCInfo * [conf_info.get_n_tile()];
for(j=0; j<conf_info.get_n_tile(); j++) {
isc_info[j] = new ISCInfo [conf_info.fcat_ta[i].n_elements()];
isc_info[j] = new ISCInfo [conf_info.fcat_ta[i].n()];
}

// Process percentile thresholds
Expand All @@ -443,7 +444,7 @@ void process_scores() {
do_intensity_scale(f_na, o_na, isc_info[j], i, j);

// Write out the ISC statistics
for(k=0; k<conf_info.fcat_ta[i].n_elements(); k++) {
for(k=0; k<conf_info.fcat_ta[i].n(); k++) {

// Store the tile definition parameters
isc_info[j][k].tile_dim = conf_info.get_tile_dim();
Expand All @@ -467,7 +468,7 @@ void process_scores() {
// Set the mask name
shc.set_mask("TILE_TOT");

for(j=0; j<conf_info.fcat_ta[i].n_elements(); j++) {
for(j=0; j<conf_info.fcat_ta[i].n(); j++) {

// Set the forecast and observation thresholds
shc.set_fcst_thresh(conf_info.fcat_ta[i][j]);
Expand Down Expand Up @@ -933,15 +934,16 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na,
double *f_dwt = (double *) 0, *o_dwt = (double *) 0; // Discrete wavelet transformations
double *f_scl = (double *) 0, *o_scl = (double *) 0; // Binary field decomposed by scale
double *diff = (double *) 0; // Difference field
double mse, fen, oen;
double mse, fen, oen, mad;
int n, ns, n_isc;
int bnd, row, col;
int i, j, k;
ConcatString fcst_thresh_str, obs_thresh_str;
bool apply_fcst_thresh, apply_obs_thresh;

// Check the NumArray lengths
n = f_na.n_elements();
if(n != o_na.n_elements()) {
n = f_na.n();
if(n != o_na.n()) {
mlog << Error << "\ndo_intensity_scale() -> "
<< "the forecast and observation arrays must have equal "
<< "length.\n\n";
Expand All @@ -962,7 +964,7 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na,
ns = conf_info.get_n_scale();

// Set up the ISCInfo thresholds and n_scale
n_isc = conf_info.fcat_ta[i_vx].n_elements();
n_isc = conf_info.fcat_ta[i_vx].n();
for(i=0; i<n_isc; i++) {
isc_info[i].clear();
isc_info[i].fthresh = conf_info.fcat_ta[i_vx][i];
Expand Down Expand Up @@ -991,22 +993,41 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na,
}

// Apply each threshold
for(i=0; i<conf_info.fcat_ta[i_vx].n_elements(); i++) {
for(i=0; i<conf_info.fcat_ta[i_vx].n(); i++) {

// Initialize to true
apply_fcst_thresh = apply_obs_thresh = true;

fcst_thresh_str = isc_info[i].fthresh.get_abbr_str();
obs_thresh_str = isc_info[i].othresh.get_abbr_str();

// If the forecast threshold is set to NA skip masking below
if(isc_info[i].fthresh.get_type() == thresh_na) {
mlog << Debug(2) << "The forecast threshold is NA, skipping applying threshold.\n";
apply_fcst_thresh = false;
}

// If the observation threshold is set to NA skip masking below
if(isc_info[i].othresh.get_type() == thresh_na) {
mlog << Debug(2) << "The observation threshold is NA, skipping applying threshold.\n";
apply_obs_thresh = false;
}

mlog << Debug(2) << "Computing Intensity-Scale decomposition for "
<< conf_info.fcst_info[i_vx]->magic_str() << " "
<< fcst_thresh_str << " versus "
<< conf_info.obs_info[i_vx]->magic_str() << " "
<< obs_thresh_str << ".\n";

// Apply the threshold to each point to create 0/1 mask fields
for(j=0; j<n; j++) {
f_dat[j] = isc_info[i].fthresh.check(f_na[j]);
o_dat[j] = isc_info[i].othresh.check(o_na[j]);
// Apply the threshold, if specified
for(j=0, mad=bad_data_double; j<n; j++) {
f_dat[j] = (apply_fcst_thresh ? isc_info[i].fthresh.check(f_na[j]) : f_na[j]);
o_dat[j] = (apply_obs_thresh ? isc_info[i].othresh.check(o_na[j]) : o_na[j]);
diff[j] = f_dat[j] - o_dat[j];

// Find the maximum absolute difference
if(is_bad_data(mad)) mad = fabs(diff[j]);
else if(fabs(diff[j]) > mad) mad = fabs(diff[j]);
} // end for j

// Compute the contingency table for the binary fields
Expand All @@ -1023,15 +1044,15 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na,
isc_info[i].compute_isc(-1);

// Write the thresholded binary fields to NetCDF
if ( conf_info.nc_info.do_raw || conf_info.nc_info.do_diff ) {
if(conf_info.nc_info.do_raw || conf_info.nc_info.do_diff ) {
write_nc_wav(conf_info.nc_info, f_dat, o_dat, n, i_vx, i_tile, -1,
isc_info[i].fthresh,
isc_info[i].othresh);
}

// Write the thresholded binary difference field to PostScript
if ( ! (conf_info.nc_info.all_false()) ) {
plot_ps_wvlt(diff, n, i_vx, i_tile, isc_info[i], -1, ns);
if(!conf_info.nc_info.all_false()) {
plot_ps_wvlt(diff, mad, n, i_vx, i_tile, isc_info[i], -1, ns);
}

// Initialize the discrete wavelet transforms
Expand Down Expand Up @@ -1102,7 +1123,7 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na,
isc_info[i].compute_isc(j);

// Write the decomposed fields for this scale to NetCDF
if ( ! (conf_info.nc_info.all_false()) ) {
if(!conf_info.nc_info.all_false()) {
write_nc_wav(conf_info.nc_info,
f_scl, o_scl, n, i_vx, i_tile, j,
isc_info[i].fthresh,
Expand All @@ -1114,7 +1135,7 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na,

// Write the decomposed difference field for this scale to PostScript
if(conf_info.ps_plot_flag) {
plot_ps_wvlt(diff, n, i_vx, i_tile, isc_info[i], j, ns);
plot_ps_wvlt(diff, mad, n, i_vx, i_tile, isc_info[i], j, ns);
}

} // end for j
Expand Down Expand Up @@ -2265,7 +2286,8 @@ void plot_ps_raw(const DataPlane &fcst_dp,

////////////////////////////////////////////////////////////////////////

void plot_ps_wvlt(const double *diff, int n, int i_vx, int i_tile,
void plot_ps_wvlt(const double *diff, double mad,
int n, int i_vx, int i_tile,
ISCInfo &isc_info,
int i_scale, int n_scale) {
ConcatString label;
Expand All @@ -2286,12 +2308,6 @@ void plot_ps_wvlt(const double *diff, int n, int i_vx, int i_tile,
v_tab = v_tab_cen;
}

//
// The min and max plotting values should default to [-1.0, 1.0]
// for the decomposed wavelet difference fields.
//
wvlt_ct.rescale(-1.0, 1.0, bad_data_double);

//
// If the wvlt_plot_min or wvlt_plot_max value is set in the
// config file, rescale the colortable to the requested range.
Expand All @@ -2302,6 +2318,15 @@ void plot_ps_wvlt(const double *diff, int n, int i_vx, int i_tile,
conf_info.wvlt_pi.plot_max,
bad_data_double);
}
//
// Otherwise, rescale the colortable to [-d, d] where d is the maximum
// absolute difference and at least 1.0. If thresholds have been applied,
// the plotting range should be [-1.0, 1.0].
//
else {
double max_plot_val = max(1.0, mad);
wvlt_ct.rescale(-1.0*max_plot_val, max_plot_val, bad_data_double);
}

//
// Set the fill color. If a fill value is not specified in the range
Expand Down Expand Up @@ -2331,10 +2356,10 @@ void plot_ps_wvlt(const double *diff, int n, int i_vx, int i_tile,
v_tab -= 1.0*plot_text_sep;
ps_out->write_centered_text(1, 1, h_tab_cen, v_tab, 0.5, 0.5, label.c_str());
if(i_scale == -1)
label.format("Tile %i, Binary, Difference (F-0)",
label.format("Tile %i, Binary, Difference (F-O)",
i_tile+1);
else
label.format("Tile %i, Scale %i, Difference (F-0)",
label.format("Tile %i, Scale %i, Difference (F-O)",
i_tile+1, i_scale+1);

v_tab -= 2.0*plot_text_sep;
Expand Down
41 changes: 34 additions & 7 deletions met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ void WaveletStatConfInfo::read_config(const char *default_file_name,

void WaveletStatConfInfo::process_config(GrdFileType ftype,
GrdFileType otype) {
int i, n;
int i, j, n;
VarInfoFactory info_factory;
map<STATLineType,STATOutputType>output_map;
Dictionary *fcst_dict = (Dictionary *) 0;
Expand All @@ -144,6 +144,9 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype,
Dictionary i_fdict, i_odict;
gsl_wavelet_type type;

SingleThresh st_NA;
st_NA.set_na();

// Dump the contents of the config file
if(mlog.verbosity_level() >= 5) conf.dump(cout);

Expand Down Expand Up @@ -250,10 +253,23 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype,
<< "Forecast categorical thresholds: " << fcat_ta[i].get_str() << "\n"
<< "Observed categorical thresholds: " << ocat_ta[i].get_str() << "\n";
}

// If the forecast threshold array is an empty list (or NA)
// Add the NA threshold type to the list for downstream iteration
if(fcat_ta[i].n() == 0) {
mlog << Debug(2) << "Found empty list for forecast threshold, setting threshold type to NA.\n";
fcat_ta[i].add(st_NA);
}

// If the observation threshold array is an empty list (or NA)
// Add the NA threshold type to the list for downstream iteration
if(ocat_ta[i].n() == 0) {
mlog << Debug(2) << "Found empty list for observation threshold, setting threshold type to NA.\n";
ocat_ta[i].add(st_NA);
}

// Check for the same number of fcst and obs thresholds
if(fcat_ta[i].n_elements() != ocat_ta[i].n_elements()) {

if(fcat_ta[i].n() != ocat_ta[i].n()) {
mlog << Error << "\nWaveletStatConfInfo::process_config() -> "
<< "The number of thresholds for each field in \"fcst."
<< conf_key_cat_thresh
Expand All @@ -262,8 +278,19 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype,
exit(1);
}

// Send a warning about inconsistent use of the NA threshold
for(j=0; j<fcat_ta[i].n(); j++) {
if((fcat_ta[i][j].get_type() == thresh_na) !=
(ocat_ta[i][j].get_type() == thresh_na)) {
mlog << Warning << "\nWaveletStatConfInfo::process_config() -> "
<< "Skipping thresholding for the forecast (" << fcat_ta[i][j].get_str()
<< ") or observation (" << ocat_ta[i][j].get_str()
<< ") but not the other may produce unexpected results!\n\n";
}
}

// Keep track of the maximum number of thresholds
if(fcat_ta[i].n_elements() > max_n_thresh) max_n_thresh = fcat_ta[i].n_elements();
if(fcat_ta[i].n() > max_n_thresh) max_n_thresh = fcat_ta[i].n();

} // end for i

Expand Down Expand Up @@ -505,7 +532,7 @@ void WaveletStatConfInfo::process_tiles(const Grid &grid) {
case(GridDecompType_Tile):

// Number of tiles based on the user-specified locations
n_tile = tile_xll.n_elements();
n_tile = tile_xll.n();

msg << "\nTiling Method: Apply " << n_tile
<< " tile(s) specified in the configuration file "
Expand Down Expand Up @@ -641,9 +668,9 @@ int WaveletStatConfInfo::n_isc_row() {
// Compute the number of output lines for each verification field
for(i=0,n=0; i<n_vx; i++) {

n += (n_scale + 2) * fcat_ta[i].n_elements() * n_tile;
n += (n_scale + 2) * fcat_ta[i].n() * n_tile;

if(n_tile > 1) n += (n_scale + 2) * fcat_ta[i].n_elements();
if(n_tile > 1) n += (n_scale + 2) * fcat_ta[i].n();
}

return(n);
Expand Down
4 changes: 2 additions & 2 deletions test/config/WaveletStatConfig
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,8 @@ obs_raw_plot = {

wvlt_plot = {
color_table = "MET_BASE/colortables/NCL_colortables/BlWhRe.ctable";
plot_min = -1.0;
plot_max = 1.0;
plot_min = 0.0;
plot_max = 0.0;
}

////////////////////////////////////////////////////////////////////////////////
Expand Down
Loading

0 comments on commit 599181a

Please sign in to comment.