Feature 1746 wavelet stat

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;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

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 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,6 +230,21 @@ _______________________
 
 .. code-block:: none
 
+
+  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.  
+
+_______________________
+
+.. code-block:: none
+
   grid_decomp_flag = AUTO;
 
   tile = {

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;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

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,12 +934,15 @@ 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;
 
+   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()) {
@@ -996,17 +1000,33 @@ void do_intensity_scale(const NumArray &f_na, const NumArray &o_na,
       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 (missing), skipping applying threshold.\n";
+         apply_fcst_thresh = 0;
+      }
+
+      // 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 (missing), skipping applying threshold.\n";
+         apply_obs_thresh = 0;
+      }
+
       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
@@ -1023,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
@@ -1102,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,
@@ -1114,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
@@ -2265,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;
@@ -2286,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.
@@ -2302,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
@@ -2331,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;

met/src/tools/core/wavelet_stat/wavelet_stat_conf_info.cc

@@ -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);
 
@@ -250,7 +253,35 @@ 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);
+      }
+
+      // 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()) {
 

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;
 }
 
 ////////////////////////////////////////////////////////////////////////////////