From 653a3fe24dab34426174f40859cdb25f158ee6a5 Mon Sep 17 00:00:00 2001 From: Seth Linden Date: Thu, 8 Jul 2021 11:05:14 -0600 Subject: [PATCH 01/14] 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 --- .../tools/core/wavelet_stat/wavelet_stat.cc | 30 +++++++++++++++---- .../wavelet_stat/wavelet_stat_conf_info.cc | 19 +++++++++++- 2 files changed, 42 insertions(+), 7 deletions(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat.cc b/met/src/tools/core/wavelet_stat/wavelet_stat.cc index e557d039a0..2ee3bd5413 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -454,7 +454,7 @@ void process_scores() { shc.set_fcst_thresh(conf_info.fcat_ta[i][k]); shc.set_obs_thresh(conf_info.ocat_ta[i][k]); - write_isc_row(shc, isc_info[j][k], + write_isc_row(shc, isc_info[j][k], conf_info.output_flag[i_isc], stat_at, i_stat_row, isc_at, i_isc_row); } // end for k @@ -939,6 +939,9 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, int i, j, k; ConcatString fcst_thresh_str, obs_thresh_str; + int apply_fcst_thresh = 1; + int apply_obs_thresh = 1; + // Check the NumArray lengths n = f_na.n_elements(); if(n != o_na.n_elements()) { @@ -992,10 +995,22 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, // Apply each threshold for(i=0; imagic_str() << " " << fcst_thresh_str << " versus " @@ -1004,11 +1019,14 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, // Apply the threshold to each point to create 0/1 mask fields for(j=0; j= 5) conf.dump(cout); @@ -250,7 +253,21 @@ 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_elements() == 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_elements() == 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()) { From 1578a4e577a79f51c346c07f29751838d19338cf Mon Sep 17 00:00:00 2001 From: Seth Linden Date: Thu, 8 Jul 2021 13:37:36 -0600 Subject: [PATCH 02/14] For issue #1746, added new unit test that uses a config file that has empty lists for the forecast and observation thresholds. SL --- test/xml/unit_wavelet_stat.xml | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/test/xml/unit_wavelet_stat.xml b/test/xml/unit_wavelet_stat.xml index 0b4d33510d..9e5f6efdd9 100644 --- a/test/xml/unit_wavelet_stat.xml +++ b/test/xml/unit_wavelet_stat.xml @@ -15,7 +15,10 @@ &TEST_DIR; true - + + + + &MET_BIN;/wavelet_stat @@ -35,5 +38,27 @@ + + + + + + &MET_BIN;/wavelet_stat + + OUTPUT_PREFIX GRIB1_NAM_STAGE4_NO_THRESH + + \ + &DATA_DIR_MODEL;/grib1/nam_st4/nam_2012040900_F012_gSt4.grib \ + &OUTPUT_DIR;/pcp_combine/stage4_2012040912_F012_APCP12.nc \ + &CONFIG_DIR;/WaveletStatConfig_no_thresholds \ + -outdir &OUTPUT_DIR;/wavelet_stat -v 1 + + + &OUTPUT_DIR;/wavelet_stat/wavelet_stat_GRIB1_NAM_STAGE4_NO_THRESH_120000L_20120409_120000V.stat + &OUTPUT_DIR;/wavelet_stat/wavelet_stat_GRIB1_NAM_STAGE4_NO_THRESH_120000L_20120409_120000V_isc.txt + &OUTPUT_DIR;/wavelet_stat/wavelet_stat_GRIB1_NAM_STAGE4_NO_THRESH_120000L_20120409_120000V.nc + &OUTPUT_DIR;/wavelet_stat/wavelet_stat_GRIB1_NAM_STAGE4_NO_THRESH_120000L_20120409_120000V.ps + + From f95389afc25b3ec709c623bc829b0d762eee8f45 Mon Sep 17 00:00:00 2001 From: Seth Linden Date: Thu, 8 Jul 2021 15:56:32 -0600 Subject: [PATCH 03/14] 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 --- met/docs/Users_Guide/wavelet-stat.rst | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/met/docs/Users_Guide/wavelet-stat.rst b/met/docs/Users_Guide/wavelet-stat.rst index ae26221f34..47f3722cd4 100644 --- a/met/docs/Users_Guide/wavelet-stat.rst +++ b/met/docs/Users_Guide/wavelet-stat.rst @@ -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 all assigned to 1. :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. @@ -230,6 +230,22 @@ _______________________ .. code-block:: none + fcst = { + field = [ + { + name = "APCP"; + level = [ "A12" ]; + cat_thresh = []; + } + ]; + } + +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 grid-point values will be set to 1 and used for statistics + +_______________________ + +.. code-block:: none + grid_decomp_flag = AUTO; tile = { From da0b030a242f213fb60280023d1c29bcbf63f847 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Fri, 9 Jul 2021 12:38:00 -0600 Subject: [PATCH 04/14] Per #1746, cleaning up for consistent indentation. --- .../tools/core/wavelet_stat/wavelet_stat.cc | 23 ++++++++----------- .../wavelet_stat/wavelet_stat_conf_info.cc | 8 +++---- 2 files changed, 14 insertions(+), 17 deletions(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat.cc b/met/src/tools/core/wavelet_stat/wavelet_stat.cc index 2ee3bd5413..f35c0eeabe 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -454,7 +454,7 @@ void process_scores() { shc.set_fcst_thresh(conf_info.fcat_ta[i][k]); shc.set_obs_thresh(conf_info.ocat_ta[i][k]); - write_isc_row(shc, isc_info[j][k], + write_isc_row(shc, isc_info[j][k], conf_info.output_flag[i_isc], stat_at, i_stat_row, isc_at, i_isc_row); } // end for k @@ -995,20 +995,20 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, // Apply each threshold for(i=0; i Date: Fri, 9 Jul 2021 12:39:23 -0600 Subject: [PATCH 05/14] Per #1746, cleaning up for consistent indentation. --- met/src/tools/core/wavelet_stat/wavelet_stat.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat.cc b/met/src/tools/core/wavelet_stat/wavelet_stat.cc index f35c0eeabe..daab9bd712 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -1021,7 +1021,7 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, for(j=0; j Date: Fri, 9 Jul 2021 16:41:43 -0600 Subject: [PATCH 06/14] 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. --- .../tools/core/wavelet_stat/wavelet_stat.cc | 53 +++++++++++-------- 1 file changed, 31 insertions(+), 22 deletions(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat.cc b/met/src/tools/core/wavelet_stat/wavelet_stat.cc index daab9bd712..de3848cd8e 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -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. // //////////////////////////////////////////////////////////////////////// @@ -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); @@ -933,7 +934,7 @@ 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; @@ -1017,11 +1018,15 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, << 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 mad) mad = fabs(diff[j]); } // end for j // Compute the contingency table for the binary fields @@ -1038,15 +1043,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 @@ -1117,7 +1122,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, @@ -1129,7 +1134,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 @@ -2280,7 +2285,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; @@ -2301,12 +2307,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. @@ -2317,6 +2317,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 @@ -2346,10 +2355,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; From 06a1c75900e954fbfb119caeb410cca6b71de59b Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Fri, 9 Jul 2021 16:44:47 -0600 Subject: [PATCH 07/14] 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. --- met/data/config/WaveletStatConfig_default | 4 ++-- met/scripts/config/WaveletStatConfig_APCP_12 | 8 ++++---- test/config/WaveletStatConfig | 4 ++-- test/config/WaveletStatConfig_no_thresholds | 4 ++-- test/config/WaveletStatConfig_python | 4 ++-- test/config/WaveletStatConfig_python_mixed | 4 ++-- 6 files changed, 14 insertions(+), 14 deletions(-) diff --git a/met/data/config/WaveletStatConfig_default b/met/data/config/WaveletStatConfig_default index 5a50fe2f0d..6a7356c7ee 100644 --- a/met/data/config/WaveletStatConfig_default +++ b/met/data/config/WaveletStatConfig_default @@ -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; } //////////////////////////////////////////////////////////////////////////////// diff --git a/met/scripts/config/WaveletStatConfig_APCP_12 b/met/scripts/config/WaveletStatConfig_APCP_12 index 8abdaddb57..67079ef1c7 100644 --- a/met/scripts/config/WaveletStatConfig_APCP_12 +++ b/met/scripts/config/WaveletStatConfig_APCP_12 @@ -50,7 +50,7 @@ fcst = { { name = "APCP"; level = [ "A12" ]; - cat_thresh = [ >0.0, >=5.0 ]; + cat_thresh = [ >0.0, >=5.0, NA ]; } ]; } @@ -59,7 +59,7 @@ obs = { { name = "APCP_12"; level = [ "(*,*)" ]; - cat_thresh = [ >0.0, >=5.0 ]; + cat_thresh = [ >0.0, >=5.0, NA ]; } ]; } @@ -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; } //////////////////////////////////////////////////////////////////////////////// diff --git a/test/config/WaveletStatConfig b/test/config/WaveletStatConfig index b9a3983bdf..06963c8f51 100644 --- a/test/config/WaveletStatConfig +++ b/test/config/WaveletStatConfig @@ -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; } //////////////////////////////////////////////////////////////////////////////// diff --git a/test/config/WaveletStatConfig_no_thresholds b/test/config/WaveletStatConfig_no_thresholds index 9b7baef8e0..42e02ca680 100644 --- a/test/config/WaveletStatConfig_no_thresholds +++ b/test/config/WaveletStatConfig_no_thresholds @@ -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; } //////////////////////////////////////////////////////////////////////////////// diff --git a/test/config/WaveletStatConfig_python b/test/config/WaveletStatConfig_python index 28df2065a8..b9e2d62f04 100644 --- a/test/config/WaveletStatConfig_python +++ b/test/config/WaveletStatConfig_python @@ -118,8 +118,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; } //////////////////////////////////////////////////////////////////////////////// diff --git a/test/config/WaveletStatConfig_python_mixed b/test/config/WaveletStatConfig_python_mixed index a5e09298b7..899766b007 100644 --- a/test/config/WaveletStatConfig_python_mixed +++ b/test/config/WaveletStatConfig_python_mixed @@ -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; } //////////////////////////////////////////////////////////////////////////////// From 1fd1c591180a069c14364d6954f8792c06ebade5 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Mon, 12 Jul 2021 09:46:19 -0600 Subject: [PATCH 08/14] Per #1746, used apply_fcst_thresh where it should have been apply_obs_thresh. --- met/src/tools/core/wavelet_stat/wavelet_stat.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat.cc b/met/src/tools/core/wavelet_stat/wavelet_stat.cc index de3848cd8e..b477e5e18c 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -1021,7 +1021,7 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, // Apply the threshold, if specified for(j=0, mad=bad_data_double; j Date: Mon, 12 Jul 2021 15:09:42 -0600 Subject: [PATCH 09/14] 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 --- met/docs/Users_Guide/wavelet-stat.rst | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/met/docs/Users_Guide/wavelet-stat.rst b/met/docs/Users_Guide/wavelet-stat.rst index 47f3722cd4..47007d162d 100644 --- a/met/docs/Users_Guide/wavelet-stat.rst +++ b/met/docs/Users_Guide/wavelet-stat.rst @@ -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. This can also be done with no thresholds indicated at all and in that case the grid-point values are all assigned to 1. :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. @@ -230,17 +230,16 @@ _______________________ .. code-block:: none - fcst = { - field = [ - { - name = "APCP"; - level = [ "A12" ]; - cat_thresh = []; - } - ]; - } - -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 grid-point values will be set to 1 and used for statistics + + cat_thresh = []; + 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 threhsolds is an empy 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 set the threshold to NA internally and skip applying the thresholds. + +Since the application has the ability to loop through multiple thresholds (for multiple fields), a user can specify a list of thresholds and include NA at the end of the list in order to produce statistics without applying any threshold for the given variable. _______________________ From 249f79c6e5d584e60e516b7568d8d23a3ff76dda Mon Sep 17 00:00:00 2001 From: Seth Linden Date: Mon, 12 Jul 2021 15:41:47 -0600 Subject: [PATCH 10/14] 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 --- .../core/wavelet_stat/wavelet_stat_conf_info.cc | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc b/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc index d66ad25789..b02723de0e 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc @@ -257,7 +257,7 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype, // 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_elements() == 0) { - mlog << Debug(2) << "Found empty list for forecast threshold, setting threshold type to NA.\n"; + mlog << Debug(2) << "Found empty list for forecast threshold, setting threshold type to NA.\n"; fcat_ta[i].add(st_NA); } @@ -268,6 +268,20 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype, ocat_ta[i].add(st_NA); } + // Send a warning if forecast threshold is NA but observation threshold is not + if( (strcmp((fcat_ta[i].get_str()).c_str(), "NA") == 0) && + (strcmp((ocat_ta[i].get_str()).c_str(), "NA") != 0) ) { + mlog << Warning << "The forecast threshold is set to NA but the observation threshold is " << ocat_ta[i].get_str() << "...\n"; + mlog << Warning << "if you are getting unexpected results change the forecast threshold to a numeric threshold.\n"; + } + + // Send a warning if observation threshold is NA but forecast threshold is not + if( (strcmp((ocat_ta[i].get_str()).c_str(), "NA") == 0) && + (strcmp((fcat_ta[i].get_str()).c_str(), "NA") != 0) ) { + mlog << Warning << "The observation threshold is set to NA but the forecast threshold is " << fcat_ta[i].get_str() << "...\n"; + mlog << Warning << "if you are getting unexpected results change the observation threshold to a numeric threshold.\n"; + } + // Check for the same number of fcst and obs thresholds if(fcat_ta[i].n_elements() != ocat_ta[i].n_elements()) { From 20a288f2ff1d89a7322993dda6874823f266f6d6 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 15 Jul 2021 10:44:33 -0600 Subject: [PATCH 11/14] Per #1746, fix a couple of typos and tweak wording in the wavelet-stat chapter. --- met/docs/Users_Guide/wavelet-stat.rst | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/met/docs/Users_Guide/wavelet-stat.rst b/met/docs/Users_Guide/wavelet-stat.rst index 47007d162d..19db97d1ed 100644 --- a/met/docs/Users_Guide/wavelet-stat.rst +++ b/met/docs/Users_Guide/wavelet-stat.rst @@ -230,16 +230,18 @@ _______________________ .. 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. +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 threhsolds is an empy 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 set the threshold to NA internally and skip applying the thresholds. +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 specify a list of thresholds and include NA at the end of the list in order to produce statistics without applying any threshold for the given variable. +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. _______________________ From abd2a34df167bbca0bd8c280e34234732dd76da5 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 15 Jul 2021 10:45:17 -0600 Subject: [PATCH 12/14] Per #1746, loop over each pair of fcst/obs thresholds to check for inconsistent use of the NA threshold type. --- .../wavelet_stat/wavelet_stat_conf_info.cc | 28 ++++++++----------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc b/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc index b02723de0e..77bd56616e 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc @@ -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; mapoutput_map; Dictionary *fcst_dict = (Dictionary *) 0; @@ -267,24 +267,9 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype, mlog << Debug(2) << "Found empty list for observation threshold, setting threshold type to NA.\n"; ocat_ta[i].add(st_NA); } - - // Send a warning if forecast threshold is NA but observation threshold is not - if( (strcmp((fcat_ta[i].get_str()).c_str(), "NA") == 0) && - (strcmp((ocat_ta[i].get_str()).c_str(), "NA") != 0) ) { - mlog << Warning << "The forecast threshold is set to NA but the observation threshold is " << ocat_ta[i].get_str() << "...\n"; - mlog << Warning << "if you are getting unexpected results change the forecast threshold to a numeric threshold.\n"; - } - // Send a warning if observation threshold is NA but forecast threshold is not - if( (strcmp((ocat_ta[i].get_str()).c_str(), "NA") == 0) && - (strcmp((fcat_ta[i].get_str()).c_str(), "NA") != 0) ) { - mlog << Warning << "The observation threshold is set to NA but the forecast threshold is " << fcat_ta[i].get_str() << "...\n"; - mlog << Warning << "if you are getting unexpected results change the observation threshold to a numeric threshold.\n"; - } - // Check for the same number of fcst and obs thresholds if(fcat_ta[i].n_elements() != ocat_ta[i].n_elements()) { - mlog << Error << "\nWaveletStatConfInfo::process_config() -> " << "The number of thresholds for each field in \"fcst." << conf_key_cat_thresh @@ -293,6 +278,17 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype, exit(1); } + // Send a warning about inconsistent use of the NA threshold + for(j=0; j " + << "Skipping thresholding for either the forecast (" + << fcat_ta[i][j].get_str() << ") or observation (" << ocat_ta[i][j].get_str() + << ") but not the other. This 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(); From 8650b7403b0cac43e73cfe12204e9d962fa80181 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 15 Jul 2021 10:47:14 -0600 Subject: [PATCH 13/14] Per #1746, a bit of code cleanup replacing calls to n_elements() with n() to make the code more concise. --- met/src/tools/core/wavelet_stat/wavelet_stat.cc | 14 +++++++------- .../core/wavelet_stat/wavelet_stat_conf_info.cc | 14 +++++++------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat.cc b/met/src/tools/core/wavelet_stat/wavelet_stat.cc index b477e5e18c..677a0d7de3 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -420,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 " << "the forecast and observation arrays must have equal " << "length.\n\n"; @@ -966,7 +966,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 " << "The number of thresholds for each field in \"fcst." << conf_key_cat_thresh @@ -290,7 +290,7 @@ void WaveletStatConfInfo::process_config(GrdFileType ftype, } // 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 @@ -532,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 " @@ -668,9 +668,9 @@ int WaveletStatConfInfo::n_isc_row() { // Compute the number of output lines for each verification field for(i=0,n=0; i 1) n += (n_scale + 2) * fcat_ta[i].n_elements(); + if(n_tile > 1) n += (n_scale + 2) * fcat_ta[i].n(); } return(n); From fb1a109b4b7d794c21db3e3e50c48e8897a5f3b3 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 15 Jul 2021 11:07:06 -0600 Subject: [PATCH 14/14] 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. --- met/src/tools/core/wavelet_stat/wavelet_stat.cc | 15 ++++++++------- .../core/wavelet_stat/wavelet_stat_conf_info.cc | 6 +++--- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/met/src/tools/core/wavelet_stat/wavelet_stat.cc b/met/src/tools/core/wavelet_stat/wavelet_stat.cc index 677a0d7de3..c2dbc61635 100644 --- a/met/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/met/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -939,10 +939,8 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, int bnd, row, col; int i, j, k; ConcatString fcst_thresh_str, obs_thresh_str; + bool apply_fcst_thresh, apply_obs_thresh; - int apply_fcst_thresh = 1; - int apply_obs_thresh = 1; - // Check the NumArray lengths n = f_na.n(); if(n != o_na.n()) { @@ -997,19 +995,22 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na, // Apply each threshold for(i=0; i " - << "Skipping thresholding for either the forecast (" - << fcat_ta[i][j].get_str() << ") or observation (" << ocat_ta[i][j].get_str() - << ") but not the other. This may produce unexpected results!\n\n"; + << "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"; } }