Update develop-ref after #1955 and #1956 (#1959)

Co-authored-by: Julie Prestopnik <jpresto@seneca.rap.ucar.edu>
Co-authored-by: johnhg <johnhg@ucar.edu>
Co-authored-by: Seth Linden <linden@kiowa.rap.ucar.edu>
Co-authored-by: John Halley Gotway <johnhg@ucar.edu>
Co-authored-by: j-opatz <59586397+j-opatz@users.noreply.github.com>
Co-authored-by: jprestop <jpresto@ucar.edu>
Co-authored-by: Seth Linden <linden@ucar.edu>

met/docs/Users_Guide/config_options.rst

@@ -93,12 +93,14 @@ The configuration file language supports the following data types:
       e.g. "<USP90(2.5)" means less than the 90-th percentile values which
       the user has already determined to be 2.5 outside of MET.
 
-    * "==FBIAS1" to automatically de-bias the data. This option must be
-      used in conjunction with a simple threshold in the other field.
-      For example, when "obs.cat_thresh = >5.0" and
-      "fcst.cat_thresh = ==FBIAS1;", MET applies the >5.0 threshold to the
-      observations and then chooses a forecast threshold which results in a
-      frequency bias of 1.
+    * "==FBIAS" for a user-specified frequency bias value.
+      e.g. "==FBIAS1" to automatically de-bias the data, "==FBIAS0.9" to select a low-bias threshold, or "==FBIAS1.1" to select a high-bias threshold.
+      This option must be used in
+      conjunction with a simple threshold in the other field. For example,
+      when "obs.cat_thresh = >5.0" and "fcst.cat_thresh = ==FBIAS1;",
+      MET applies the >5.0 threshold to the observations and then chooses a
+      forecast threshold which results in a frequency bias of 1.
+      The frequency bias can be any float value > 0.0.
 
     * "CDP" for climatological distribution percentile thresholds.
       These thresholds require that the climatological mean and standard
@@ -119,7 +121,7 @@ The configuration file language supports the following data types:
     For example, "==CDP25" is automatically expanded to 4 percentile bins:
     >=CDP0&&<CDP25,>=CDP25&&<CDP50,>=CDP50&&<CDP75,>=CDP75&&<=CDP100
 
-  * When sample percentile thresholds of type SFP, SOP, SCP, or FBIAS1 are
+  * When sample percentile thresholds of type SFP, SOP, SCP, or FBIAS are
     requested, MET recomputes the actual percentile that the threshold
     represents. If the requested percentile and actual percentile differ by
     more than 5%, a warning message is printed. This may occur when the

met/docs/Users_Guide/ensemble-stat.rst

@@ -69,6 +69,7 @@ The usage statement for the Ensemble Stat tool is shown below:
          [-grid_obs file]
          [-point_obs file]
          [-ens_mean file]
+         [-ctrl file]
          [-obs_valid_beg time]
          [-obs_valid_end time]
          [-outdir path]
@@ -92,23 +93,23 @@ Optional arguments for ensemble_stat
 
 4. To produce ensemble statistics using gridded observations, use the **-grid_obs file** option to specify a gridded observation file. This option may be used multiple times if your observations are in several files.
 
+5. To produce ensemble statistics using point observations, use the **-point_obs file** option to specify a NetCDF point observation file. This option may be used multiple times if your observations are in several files.
 
-5. To produce ensemble statistics using point observations, use the **-point_obs file** to specify a NetCDF point observation file. This option may be used multiple times if your observations are in several files.
+6. To override the simple ensemble mean value of the input ensemble members for the ECNT, SSVAR, and ORANK line types, the **-ens_mean file** option specifies an ensemble mean model data file. This option replaces the **-ssvar_mean file** option from earlier versions of MET.
 
+7. The **-ctrl file** option specifies an ensemble control member data file. The control member is included in the computation of the ensemble mean but excluded from the spread.
 
-6. To override the simple ensemble mean value of the input ensemble members for the ECNT, SSVAR, and ORANK line types, the **-ens_mean file** specifies an ensemble mean model data file. This option replaces the **-ssvar_mean file** from earlier versions of MET.
+8. To filter point observations by time, use **-obs_valid_beg time** in YYYYMMDD[_HH[MMSS]] format to set the beginning of the matching observation time window.
 
-7. To filter point observations by time, use **-obs_valid_beg time** in YYYYMMDD[_HH[MMSS]] format to set the beginning of the matching observation time window.
+9. As above, use **-obs_valid_end time** in YYYYMMDD[_HH[MMSS]] format to set the end of the matching observation time window.
 
-8. As above, use **-obs_valid_end time** in YYYYMMDD[_HH[MMSS]] format to set the end of the matching observation time window.
+10. Specify the **-outdir path** option to override the default output directory (./).
 
-9. Specify the **-outdir path** option to override the default output directory (./).
+11. The **-log** file outputs log messages to the specified file.
 
-10. The **-log** file outputs log messages to the specified file.
+12. The **-v level** option indicates the desired level of verbosity. The value of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity will increase the amount of logging.
 
-11. The **-v level** option indicates the desired level of verbosity. The value of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity will increase the amount of logging.
-
-12. The **-compress level** option indicates the desired level of compression (deflate level) for NetCDF variables. The valid level is between 0 and 9. The value of "level" will override the default setting of 0 from the configuration file or the environment variable MET_NC_COMPRESS. Setting the compression level to 0 will make no compression for the NetCDF output. Lower number is for fast compression and higher number is for better compression.
+13. The **-compress level** option indicates the desired level of compression (deflate level) for NetCDF variables. The valid level is between 0 and 9. The value of "level" will override the default setting of 0 from the configuration file or the environment variable MET_NC_COMPRESS. Setting the compression level to 0 will make no compression for the NetCDF output. Lower number is for fast compression and higher number is for better compression.
 
 An example of the ensemble_stat calling sequence is shown below:
 

met/docs/Users_Guide/reformat_grid.rst

@@ -74,7 +74,7 @@ Required arguments for the pcp_combine
 Optional arguments for pcp_combine
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. This option may be used multiple times and output will be created for each.
+3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each. In general, the field string should include the **name** and **level** of the requested data and be enclosed in single quotes. It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options. For example, use **-field ‘name=”ACPCP”; level=”A6”; convert(x)=x/25.4;’** to read 6-hourly accumulated convective precipitation from a GRIB file and convert from millimeters to inches. 
 
 4. The **-name list** option is a comma-separated list of output variable names which override the default choices. If specified, the number of names must match the number of variables to be written to the output file.
 

met/src/basic/vx_config/config.tab.yy

@@ -1661,7 +1661,7 @@ if ( s->PT < 0 || s->PT > 100 )  {
 
 }
 
-if ( s->Ptype == perc_thresh_freq_bias && !is_eq(s->PT, 1.0) )  {
+if ( s->Ptype == perc_thresh_freq_bias && s->PT <= 0 )  {
 
    mlog << Error << "\ndo_simple_perc_thresh() -> "
         << "unsupported frequency bias percentile threshold!\n\n";

met/src/basic/vx_config/config_file.cc

@@ -312,9 +312,10 @@ if ( empty(name) )  {
 
 DictionaryStack DS(*this);
 ConcatString temp_filename = get_tmp_dir();
-
-temp_filename << "/" << make_temp_file_name("met_config", 0);
-
+
+temp_filename << "/" << "met_config";
+temp_filename = make_temp_file_name(temp_filename.c_str(), 0);
+
 recursive_envs(name, temp_filename.c_str());
 
 bison_input_filename = (const char *) temp_filename.c_str();
@@ -422,9 +423,10 @@ if ( empty(s) )  {
 ofstream out;
 ConcatString temp_filename = get_tmp_dir();
 
-temp_filename << "/" << make_temp_file_name("config", ".temp");
-
- out.open(temp_filename.c_str());
+temp_filename << "/" << "met_config";
+temp_filename = make_temp_file_name(temp_filename.c_str(), 0);
+
+out.open(temp_filename.c_str());
 
 if ( ! out )  {
 

met/src/basic/vx_config/threshold.cc

@@ -934,10 +934,11 @@ void Simple_Node::set_perc(const NumArray *fptr, const NumArray *optr, const Num
 
 {
 
-int i, count;
+int i;
 double ptile, diff;
 NumArray data;
 const NumArray * ptr = 0;
+bool fbias_fcst = false;
 
    //
    //  handle sample percentile types
@@ -972,14 +973,16 @@ else if ( Ptype == perc_thresh_freq_bias )  {
         fthr->get_ptype() == no_perc_thresh_type   &&
         fthr->get_type()  != thresh_complex )  {
 
+      fbias_fcst = false;
+
       ptr = optr;
       op  = fthr->get_type();
       PT  = fptr->compute_percentile(fthr->get_value(),
                                      is_inclusive(fthr->get_type()));
 
       mlog << Debug(3)
-           << "The forecast threshold \"" << fthr->get_str()
-           << "\" includes " << PT * 100.0 << "% of the data.\n";
+           << "The forecast threshold value \"" << fthr->get_str()
+           << "\" represents the " << PT * 100.0 << "-th percentile.\n";
 
    }
 
@@ -991,14 +994,16 @@ else if ( Ptype == perc_thresh_freq_bias )  {
              othr->get_ptype() == no_perc_thresh_type   &&
              othr->get_type()  != thresh_complex )  {
 
+      fbias_fcst = true;
+
       ptr = fptr;
       op  = othr->get_type();
       PT  = optr->compute_percentile(othr->get_value(),
                                      is_inclusive(othr->get_type()));
 
       mlog << Debug(3)
-           << "The observation threshold \"" << othr->get_str()
-           << "\" includes " << PT * 100.0 << "% of the data.\n";
+           << "The observation threshold value \"" << othr->get_str()
+           << "\" represents the " << PT * 100.0 << "-th percentile.\n";
 
    }
 
@@ -1029,7 +1034,8 @@ else if ( Ptype == perc_thresh_freq_bias )  {
 
    PT *= 100.0;
 
-}
+} // end else if PT == perc_thresh_freq_bias
+
    //
    //  nothing to do
    //
@@ -1082,7 +1088,7 @@ if ( data.n() == 0 )  {
         << " threshold \"" << s
         << "\" because no valid data was provided.\n\n";
 
-      exit ( 1 );
+   exit ( 1 );
 
 }
 else  {
@@ -1094,6 +1100,64 @@ else  {
    s.strip_paren();
    abbr_s.strip_paren();
 
+   //
+   //  parse the frequency bias value from the threshold string
+   //
+
+   if ( Ptype == perc_thresh_freq_bias )  {
+
+      ConcatString fs = s;
+
+      fs.replace("==FBIAS", " ", false);
+
+      double fbias_val = atof(fs.c_str());
+
+      //
+      //  range check requested bias value
+      //
+
+      if ( fbias_val <= 0.0 )  {
+
+         mlog << Error << "\nSimple_Node::set_perc() -> "
+              << "the requested frequency bias value (" << fbias_val
+              << ") must be > 0 in threshold \"" << s << "\".\n\n";
+
+      }
+
+      //
+      //  adjust PT by the requested frequency bias amount
+      //
+
+      double PT_new;
+
+      if ( fbias_fcst )  {
+               if ( op == thresh_le || op == thresh_lt )  PT_new = PT * fbias_val;
+          else if ( op == thresh_ge || op == thresh_gt )  PT_new = PT / fbias_val;     
+      }
+      else  {
+               if ( op == thresh_le || op == thresh_lt )  PT_new = PT / fbias_val;
+          else if ( op == thresh_ge || op == thresh_gt )  PT_new = PT * fbias_val;     
+      }
+
+      if ( PT_new > 100.0 )  {
+         mlog << Warning << "\nSimple_Node::set_perc() -> "
+              << "For " << (fbias_fcst ? "forecast" : "observation" )
+              << " threshold \"" << s << "\" the required percentile of "
+              << PT_new << " exceeds the maximum possible value. "
+              << "Resetting to 100.\n\n";
+
+         PT_new = 100.0;
+      }
+
+      mlog << Debug(3)
+           << "For " << (fbias_fcst ? "forecast" : "observation" )
+           << " threshold \"" << s << "\" with type \"" << thresh_type_str[op]
+           << "\" update the requested percentile from " << PT << " to "
+           << PT_new << ".\n";
+
+      PT = PT_new;
+   }
+
    //
    //  compute the percentile and update the strings
    //
@@ -1117,27 +1181,16 @@ else  {
    //  compute the actual percentile and check tolerance
    //
 
-   if ( op == thresh_le || op == thresh_ge || op == thresh_eq )  {
-
-      for ( i=count=0; i<data.n(); i++ )  if ( data[i] <= T ) count++;
-
-   }
-   else  {
-
-      for ( i=count=0; i<data.n(); i++ )  if ( data[i] <  T ) count++;
-
-   }
-
-   ptile = (double) count / data.n();
+   ptile = data.compute_percentile(T, is_inclusive(op));
    diff  = abs(PT / 100.0 - ptile);
 
    if ( !is_eq(PT / 100.0, ptile, perc_thresh_default_tol) )  {
 
       mlog << Warning << "\nSimple_Node::set_perc() -> "
            << "the requested percentile (" << PT
-           << "%) for threshold \"" << s
+           << ") for threshold \"" << s
            << "\" differs from the actual percentile ("
-           << ptile * 100.0 << ") by " << diff * 100.0 << "%.\n"
+           << ptile * 100.0 << ") by " << diff * 100.0 << ".\n"
            << "This is common for small samples or data that contains "
            << "ties.\n\n";
 
@@ -1146,8 +1199,8 @@ else  {
 
       mlog << Debug(3)
            << "The requested percentile (" << PT
-           << "%) for threshold threshold \"" << s
-           << "\" includes " << ptile * 100.0 << "% of the data.\n";
+           << ") for threshold \"" << s << "\" includes "
+           << ptile * 100.0 << "% of the data.\n";
 
    }
 

met/src/basic/vx_util/num_array.cc

@@ -1157,6 +1157,59 @@ double NumArray::wmean_fisher(const NumArray &wgt) const
 
 }
 
+
+////////////////////////////////////////////////////////////////////////
+
+
+double NumArray::variance(int skip_index) const
+
+{
+
+   if(n() == 0)  return ( bad_data_double );
+
+   int j, count;
+   double s, s_sq, var;
+
+   s = s_sq = 0.0;
+   count = 0;
+
+   for(j=0; j<n(); j++) {
+      if(is_bad_data(e[j]) || j == skip_index) continue;
+      s    += e[j];
+      s_sq += e[j]*e[j];
+      count++;
+   }
+
+   // Check for slightly negative precision error
+   if(count > 1) {
+      var = (s_sq - s*s/(double) count)/((double) (count - 1));
+      if(is_eq(var, 0.0)) var = 0.0;
+   }
+   else {
+      var = bad_data_double;
+   }
+
+   return(var);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+
+double NumArray::stdev(int skip_index) const
+
+{
+
+   double v = variance(skip_index);
+
+   if ( !is_bad_data(v) )  v = sqrt(v);
+
+   return(v);
+
+}
+
+
 ////////////////////////////////////////////////////////////////////////
 
 //

met/src/basic/vx_util/num_array.h

@@ -21,6 +21,7 @@
 #include <iostream>
 
 #include "concat_string.h"
+#include "is_bad_data.h"
 
 
 ////////////////////////////////////////////////////////////////////////
@@ -107,6 +108,9 @@ class NumArray {
       double mean_sqrt() const;
       double mean_fisher() const;
 
+      double variance(int skip_index = bad_data_int) const;
+      double stdev(int skip_index = bad_data_int) const;
+
       double wmean(const NumArray &) const;
       double wmean_sqrt(const NumArray &) const;
       double wmean_fisher(const NumArray &) const;