diff --git a/src/libcode/vx_data2d/Makefile.am b/src/libcode/vx_data2d/Makefile.am index c705ee2486..0203547703 100644 --- a/src/libcode/vx_data2d/Makefile.am +++ b/src/libcode/vx_data2d/Makefile.am @@ -15,6 +15,7 @@ libvx_data2d_a_SOURCES = \ leveltype_to_string.cc leveltype_to_string.h \ level_info.cc level_info.h \ var_info.cc var_info.h \ + var_info_pairs.cc var_info_pairs.h \ data_class.cc data_class.h \ data2d_utils.cc data2d_utils.h \ table_lookup.cc table_lookup.h \ diff --git a/src/libcode/vx_data2d/Makefile.in b/src/libcode/vx_data2d/Makefile.in index ead0ee6915..0129891261 100644 --- a/src/libcode/vx_data2d/Makefile.in +++ b/src/libcode/vx_data2d/Makefile.in @@ -111,6 +111,7 @@ am_libvx_data2d_a_OBJECTS = \ libvx_data2d_a-leveltype_to_string.$(OBJEXT) \ libvx_data2d_a-level_info.$(OBJEXT) \ libvx_data2d_a-var_info.$(OBJEXT) \ + libvx_data2d_a-var_info_pairs.$(OBJEXT) \ libvx_data2d_a-data_class.$(OBJEXT) \ libvx_data2d_a-data2d_utils.$(OBJEXT) \ libvx_data2d_a-table_lookup.$(OBJEXT) \ @@ -137,7 +138,8 @@ am__depfiles_remade = ./$(DEPDIR)/libvx_data2d_a-data2d_utils.Po \ ./$(DEPDIR)/libvx_data2d_a-leveltype_to_string.Po \ ./$(DEPDIR)/libvx_data2d_a-mask_filters.Po \ ./$(DEPDIR)/libvx_data2d_a-table_lookup.Po \ - ./$(DEPDIR)/libvx_data2d_a-var_info.Po + ./$(DEPDIR)/libvx_data2d_a-var_info.Po \ + ./$(DEPDIR)/libvx_data2d_a-var_info_pairs.Po am__mv = mv -f AM_V_lt = $(am__v_lt_@AM_V@) am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@) @@ -361,6 +363,7 @@ libvx_data2d_a_SOURCES = \ leveltype_to_string.cc leveltype_to_string.h \ level_info.cc level_info.h \ var_info.cc var_info.h \ + var_info_pairs.cc var_info_pairs.h \ data_class.cc data_class.h \ data2d_utils.cc data2d_utils.h \ table_lookup.cc table_lookup.h \ @@ -423,6 +426,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_data2d_a-mask_filters.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_data2d_a-table_lookup.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_data2d_a-var_info.Po@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_data2d_a-var_info_pairs.Po@am__quote@ # am--include-marker $(am__depfiles_remade): @$(MKDIR_P) $(@D) @@ -486,6 +490,20 @@ libvx_data2d_a-var_info.obj: var_info.cc @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_data2d_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_data2d_a-var_info.obj `if test -f 'var_info.cc'; then $(CYGPATH_W) 'var_info.cc'; else $(CYGPATH_W) '$(srcdir)/var_info.cc'; fi` +libvx_data2d_a-var_info_pairs.o: var_info_pairs.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_data2d_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_data2d_a-var_info_pairs.o -MD -MP -MF $(DEPDIR)/libvx_data2d_a-var_info_pairs.Tpo -c -o libvx_data2d_a-var_info_pairs.o `test -f 'var_info_pairs.cc' || echo '$(srcdir)/'`var_info_pairs.cc +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_data2d_a-var_info_pairs.Tpo $(DEPDIR)/libvx_data2d_a-var_info_pairs.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='var_info_pairs.cc' object='libvx_data2d_a-var_info_pairs.o' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_data2d_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_data2d_a-var_info_pairs.o `test -f 'var_info_pairs.cc' || echo '$(srcdir)/'`var_info_pairs.cc + +libvx_data2d_a-var_info_pairs.obj: var_info_pairs.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_data2d_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_data2d_a-var_info_pairs.obj -MD -MP -MF $(DEPDIR)/libvx_data2d_a-var_info_pairs.Tpo -c -o libvx_data2d_a-var_info_pairs.obj `if test -f 'var_info_pairs.cc'; then $(CYGPATH_W) 'var_info_pairs.cc'; else $(CYGPATH_W) '$(srcdir)/var_info_pairs.cc'; fi` +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_data2d_a-var_info_pairs.Tpo $(DEPDIR)/libvx_data2d_a-var_info_pairs.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='var_info_pairs.cc' object='libvx_data2d_a-var_info_pairs.obj' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_data2d_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_data2d_a-var_info_pairs.obj `if test -f 'var_info_pairs.cc'; then $(CYGPATH_W) 'var_info_pairs.cc'; else $(CYGPATH_W) '$(srcdir)/var_info_pairs.cc'; fi` + libvx_data2d_a-data_class.o: data_class.cc @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_data2d_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_data2d_a-data_class.o -MD -MP -MF $(DEPDIR)/libvx_data2d_a-data_class.Tpo -c -o libvx_data2d_a-data_class.o `test -f 'data_class.cc' || echo '$(srcdir)/'`data_class.cc @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_data2d_a-data_class.Tpo $(DEPDIR)/libvx_data2d_a-data_class.Po @@ -673,6 +691,7 @@ distclean: distclean-am -rm -f ./$(DEPDIR)/libvx_data2d_a-mask_filters.Po -rm -f ./$(DEPDIR)/libvx_data2d_a-table_lookup.Po -rm -f ./$(DEPDIR)/libvx_data2d_a-var_info.Po + -rm -f ./$(DEPDIR)/libvx_data2d_a-var_info_pairs.Po -rm -f Makefile distclean-am: clean-am distclean-compile distclean-generic \ distclean-tags @@ -725,6 +744,7 @@ maintainer-clean: maintainer-clean-am -rm -f ./$(DEPDIR)/libvx_data2d_a-mask_filters.Po -rm -f ./$(DEPDIR)/libvx_data2d_a-table_lookup.Po -rm -f ./$(DEPDIR)/libvx_data2d_a-var_info.Po + -rm -f ./$(DEPDIR)/libvx_data2d_a-var_info_pairs.Po -rm -f Makefile maintainer-clean-am: distclean-am maintainer-clean-generic diff --git a/src/libcode/vx_data2d/var_info_pairs.cc b/src/libcode/vx_data2d/var_info_pairs.cc new file mode 100644 index 0000000000..3cbf5dcddc --- /dev/null +++ b/src/libcode/vx_data2d/var_info_pairs.cc @@ -0,0 +1,224 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2024 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: var_info_pairs.cc +// +// Description: +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include + +#include "var_info.h" +#include "var_info_pairs.h" + +#include "util_constants.h" + +#include "vx_math.h" +#include "vx_util.h" +#include "vx_log.h" +#include "vx_data2d.h" + +using namespace std; + +/////////////////////////////////////////////////////////////////////////////// +// +// Code for class VarInfoPairs +// +/////////////////////////////////////////////////////////////////////////////// + +VarInfoPairs::VarInfoPairs() { + + init_from_scratch(); +} + +/////////////////////////////////////////////////////////////////////////////// + +VarInfoPairs::~VarInfoPairs() { + + clear(); +} + +/////////////////////////////////////////////////////////////////////////////// + +VarInfoPairs::VarInfoPairs(const VarInfoPairs &f) { + + init_from_scratch(); + + assign(f); +} + +/////////////////////////////////////////////////////////////////////////////// + +VarInfoPairs & VarInfoPairs::operator=(const VarInfoPairs &f) { + + if ( this == &f ) return *this; + + assign(f); + + return *this; +} + +/////////////////////////////////////////////////////////////////////////////// + +VarInfo *VarInfoPairs::clone() const { + + VarInfoPairs *ret = new VarInfoPairs(*this); + + return (VarInfo *)ret; +} + +/////////////////////////////////////////////////////////////////////////////// + +void VarInfoPairs::init_from_scratch() { + + // First call the parent's initialization + VarInfo::init_from_scratch(); + + clear(); + + return; +} + +/////////////////////////////////////////////////////////////////////////////// + +void VarInfoPairs::assign(const VarInfoPairs &v) { + + // First call the parent's assign + VarInfo::assign(v); + + return; +} + +/////////////////////////////////////////////////////////////////////////////// + +void VarInfoPairs::clear() { + + // First call the parent's clear + VarInfo::clear(); + + return; +} + +/////////////////////////////////////////////////////////////////////////////// + +void VarInfoPairs::dump(ostream &out) const { + + // Dump out the contents + out << "VarInfoPairs::dump():\n"; + + VarInfo::dump(out); + + return; +} + +/////////////////////////////////////////////////////////////////////////////// + +void VarInfoPairs::set_dict(Dictionary & dict) { + VarInfo::set_dict(dict); +} + +/////////////////////////////////////////////////////////////////////////////// + +bool VarInfoPairs::is_precipitation() const { + bool status = false; + + // + // Check set_attrs entry + // + if(!is_bad_data(SetAttrIsPrecipitation)) { + return(SetAttrIsPrecipitation != 0); + } + + return status; +} + +/////////////////////////////////////////////////////////////////////////////// + +bool VarInfoPairs::is_specific_humidity() const { + bool status = false; + + // + // Check set_attrs entry + // + if(!is_bad_data(SetAttrIsSpecificHumidity)) { + return(SetAttrIsSpecificHumidity != 0); + } + + return status; +} + +/////////////////////////////////////////////////////////////////////////////// + +bool VarInfoPairs::is_u_wind() const { + bool status = false; + + // + // Check set_attrs entry + // + if(!is_bad_data(SetAttrIsUWind)) { + return(SetAttrIsUWind != 0); + } + + return status; +} + +/////////////////////////////////////////////////////////////////////////////// + +bool VarInfoPairs::is_v_wind() const { + bool status = false; + + // + // Check set_attrs entry + // + if(!is_bad_data(SetAttrIsVWind)) { + return(SetAttrIsVWind != 0); + } + + return status; +} + +/////////////////////////////////////////////////////////////////////////////// + +bool VarInfoPairs::is_wind_speed() const { + bool status = false; + + // + // Check set_attrs entry + // + if(!is_bad_data(SetAttrIsWindSpeed)) { + return(SetAttrIsWindSpeed != 0); + } + + return status; +} + +/////////////////////////////////////////////////////////////////////////////// + +bool VarInfoPairs::is_wind_direction() const { + bool status = false; + + // + // Check set_attrs entry + // + if(!is_bad_data(SetAttrIsWindDirection)) { + return(SetAttrIsWindDirection != 0); + } + + return status; +} + +/////////////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_data2d/var_info_pairs.h b/src/libcode/vx_data2d/var_info_pairs.h new file mode 100644 index 0000000000..9d5644a6fc --- /dev/null +++ b/src/libcode/vx_data2d/var_info_pairs.h @@ -0,0 +1,73 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2024 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + + +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __VAR_INFO_PAIRS_H__ +#define __VAR_INFO_PAIRS_H__ + +/////////////////////////////////////////////////////////////////////////////// + +#include "var_info.h" +#include "vx_config.h" + +#include "data_file_type.h" + +/////////////////////////////////////////////////////////////////////////////// + +class VarInfoPairs : public VarInfo { + + private: + + void init_from_scratch(); + void assign(const VarInfoPairs &); + + public: + VarInfoPairs(); + ~VarInfoPairs(); + VarInfoPairs(const VarInfoPairs &); + VarInfoPairs & operator=(const VarInfoPairs &); + VarInfo *clone() const; + + void dump(std::ostream &) const; + void clear(); + + // + // get stuff + // + + GrdFileType file_type() const; + + // + // set stuff + // + + void set_dict(Dictionary &); + + // + // do stuff + // + + bool is_precipitation() const; + bool is_specific_humidity() const; + bool is_u_wind() const; + bool is_v_wind() const; + bool is_wind_speed() const; + bool is_wind_direction() const; +}; + +/////////////////////////////////////////////////////////////////////////////// + +inline GrdFileType VarInfoPairs::file_type() const { return FileType_Pairs; } + +/////////////////////////////////////////////////////////////////////////////// + +#endif // __VAR_INFO_PAIRS_H__ + +/////////////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_data2d_factory/var_info_factory.cc b/src/libcode/vx_data2d_factory/var_info_factory.cc index 32a398fb39..968d68dfd4 100644 --- a/src/libcode/vx_data2d_factory/var_info_factory.cc +++ b/src/libcode/vx_data2d_factory/var_info_factory.cc @@ -27,6 +27,7 @@ #include "var_info_nc_cf.h" #include "var_info_nc_met.h" #include "var_info_nc_wrf.h" +#include "var_info_pairs.h" #include "var_info_ugrid.h" #ifdef WITH_PYTHON @@ -115,6 +116,10 @@ VarInfo * VarInfoFactory::new_var_info(GrdFileType type) ugrid_compile_error(method_name); #endif + case FileType_Pairs: + vi = new VarInfoPairs; + break; + case FileType_HdfEos: mlog << Error << "\n" << method_name << "Support for GrdFileType = " << grdfiletype_to_string(type) diff --git a/src/libcode/vx_statistics/pair_base.cc b/src/libcode/vx_statistics/pair_base.cc index 88fdca66bd..4f07b87362 100644 --- a/src/libcode/vx_statistics/pair_base.cc +++ b/src/libcode/vx_statistics/pair_base.cc @@ -989,8 +989,9 @@ void VxPairBase::clear() { obs_qty_inc_filt.clear(); obs_qty_exc_filt.clear(); - mpr_column.clear(); - mpr_thresh.clear(); + mpr_thr_inc_map.clear(); + mpr_str_inc_map.clear(); + mpr_str_exc_map.clear(); fcst_ut = (unixtime) 0; beg_ut = (unixtime) 0; @@ -1061,8 +1062,9 @@ void VxPairBase::assign(const VxPairBase &vx_pb) { obs_qty_inc_filt = vx_pb.obs_qty_inc_filt; obs_qty_exc_filt = vx_pb.obs_qty_exc_filt; - mpr_column = vx_pb.mpr_column; - mpr_thresh = vx_pb.mpr_thresh; + mpr_thr_inc_map = vx_pb.mpr_thr_inc_map; + mpr_str_inc_map = vx_pb.mpr_str_inc_map; + mpr_str_exc_map = vx_pb.mpr_str_exc_map; msg_typ_sfc = vx_pb.msg_typ_sfc; msg_typ_lnd = vx_pb.msg_typ_lnd; @@ -1434,20 +1436,27 @@ void VxPairBase::set_interp(int i_interp, //////////////////////////////////////////////////////////////////////// -void VxPairBase::set_mpr_thresh(const StringArray &sa, const ThreshArray &ta) { +void VxPairBase::set_mpr_thr_inc_map(const map &m) { - // Check for constant length - if(sa.n() != ta.n()) { - mlog << Error << "\nVxPairBase::set_mpr_thresh() -> " - << "the \"" << conf_key_mpr_column << "\" (" - << write_css(sa) << ") and \"" << conf_key_mpr_thresh - << "\" (" << write_css(ta) - << ") config file entries must have the same length!\n\n"; - exit(1); - } + mpr_thr_inc_map = m; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void VxPairBase::set_mpr_str_inc_map(const map &m) { + + mpr_str_inc_map = m; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void VxPairBase::set_mpr_str_exc_map(const map &m) { - mpr_column = sa; - mpr_thresh = ta; + mpr_str_exc_map = m; return; } diff --git a/src/libcode/vx_statistics/pair_base.h b/src/libcode/vx_statistics/pair_base.h index c829a8efe7..4da5c29a51 100644 --- a/src/libcode/vx_statistics/pair_base.h +++ b/src/libcode/vx_statistics/pair_base.h @@ -286,8 +286,15 @@ class VxPairBase { ////////////////////////////////////////////////////////////////// - StringArray mpr_column; // Names of MPR columns or diffs of columns - ThreshArray mpr_thresh; // Filtering thresholds for the MPR columns + // Mapping of numeric MPR columns or diffs of colums to + // inclusion thresholds + std::map mpr_thr_inc_map; + + // Mapping of string MPR columns to list of inclusion + // and exclusion strings + std::map mpr_str_inc_map; + std::map mpr_str_exc_map; + ////////////////////////////////////////////////////////////////// @@ -378,7 +385,9 @@ class VxPairBase { void set_interp(int i_interp, InterpMthd mthd, int width, GridTemplateFactory::GridTemplates shape); - void set_mpr_thresh(const StringArray &, const ThreshArray &); + void set_mpr_thr_inc_map(const std::map &); + void set_mpr_str_inc_map(const std::map &); + void set_mpr_str_exc_map(const std::map &); void set_climo_cdf_info_ptr(const ClimoCDFInfo *); diff --git a/src/libcode/vx_statistics/pair_data_point.cc b/src/libcode/vx_statistics/pair_data_point.cc index 85595d69d3..d34b65798f 100644 --- a/src/libcode/vx_statistics/pair_data_point.cc +++ b/src/libcode/vx_statistics/pair_data_point.cc @@ -606,7 +606,7 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, // Check matched pair filtering options ConcatString reason_cs; if(!check_mpr_thresh(fcst_v, obs_v, cpi, - mpr_column, mpr_thresh, &reason_cs)) { + mpr_thr_inc_map, &reason_cs)) { if(mlog.verbosity_level() >= REJECT_DEBUG_LEVEL) { mlog << Debug(REJECT_DEBUG_LEVEL) @@ -742,13 +742,13 @@ bool check_fo_thresh(double f, double o, const ClimoPntInfo &cpi, //////////////////////////////////////////////////////////////////////// bool check_mpr_thresh(double f, double o, const ClimoPntInfo &cpi, - const StringArray &col_sa, const ThreshArray &col_ta, + const map &m, ConcatString *reason_ptr) { // Initialize if(reason_ptr) reason_ptr->erase(); // Check arrays - if(col_sa.n() == 0 || col_ta.n() == 0) return true; + if(m.size() == 0) return true; bool keep = true; bool absv = false; @@ -756,18 +756,18 @@ bool check_mpr_thresh(double f, double o, const ClimoPntInfo &cpi, ConcatString cs; double v, v_cur; - // Loop over all the column filter names - for(int i=0; i &m) { // Check for no work to be done - if(col_sa.n() == 0 && col_ta.n() == 0) return; - - // Check for constant length - if(col_sa.n() != col_ta.n()) { - mlog << Error << "\napply_mpr_thresh_mask() -> " - << "the \"" << conf_key_mpr_column << "\" (" - << write_css(col_sa) << ") and \"" << conf_key_mpr_thresh - << "\" (" << write_css(col_ta) - << ") config file entries must have the same length!\n\n"; - exit(1); - } + if(m.size() == 0) return; int nxy = fcst_dp.nx() * fcst_dp.ny(); int n_skip = 0; @@ -884,8 +877,7 @@ void apply_mpr_thresh_mask(DataPlane &fcst_dp, DataPlane &obs_dp, (ocsd_flag && is_bad_data(cpi.ocsd))) continue; // Discard pairs which do not meet the threshold criteria - if(!check_mpr_thresh(fcst_dp.buf()[i], obs_dp.buf()[i], cpi, - col_sa, col_ta)) { + if(!check_mpr_thresh(fcst_dp.buf()[i], obs_dp.buf()[i], cpi, m)) { // Increment skip counter n_skip++; @@ -902,9 +894,7 @@ void apply_mpr_thresh_mask(DataPlane &fcst_dp, DataPlane &obs_dp, mlog << Debug(3) << "Discarded " << n_skip << " of " << nxy - << " pairs for matched pair filtering columns (" - << write_css(col_sa) << ") and thresholds (" - << col_ta.get_str() << ").\n"; + << " pairs when applying matched pair filtering thresholds.\n"; return; } diff --git a/src/libcode/vx_statistics/pair_data_point.h b/src/libcode/vx_statistics/pair_data_point.h index 5f43a5e185..57084675d1 100644 --- a/src/libcode/vx_statistics/pair_data_point.h +++ b/src/libcode/vx_statistics/pair_data_point.h @@ -141,7 +141,7 @@ extern bool check_fo_thresh(double, double, const ClimoPntInfo &, const SetLogic); extern bool check_mpr_thresh(double, double, const ClimoPntInfo &, - const StringArray &, const ThreshArray &, + const std::map &, ConcatString * = 0); extern double get_mpr_column_value(double, double, const ClimoPntInfo &, @@ -150,7 +150,7 @@ extern double get_mpr_column_value(double, double, const ClimoPntInfo &, extern void apply_mpr_thresh_mask(DataPlane &, DataPlane &, DataPlane &, DataPlane &, DataPlane &, DataPlane &, - const StringArray &, const ThreshArray &); + const std::map &); extern bool check_seeps_thresh(double, double, const StringArray &, const ThreshArray &, diff --git a/src/tools/core/grid_stat/grid_stat.cc b/src/tools/core/grid_stat/grid_stat.cc index a1c43f0824..acd788993a 100644 --- a/src/tools/core/grid_stat/grid_stat.cc +++ b/src/tools/core/grid_stat/grid_stat.cc @@ -823,12 +823,11 @@ void process_scores() { << (ocsd_dp.is_empty() ? 0 : 1) << " standard deviation field(s).\n"; // Apply MPR threshold filters - if(conf_info.vx_opt[i].mpr_sa.n() > 0) { + if(conf_info.vx_opt[i].mpr_thr_inc_map.size() > 0) { apply_mpr_thresh_mask(fcst_dp, obs_dp, fcmn_dp, fcsd_dp, ocmn_dp, ocsd_dp, - conf_info.vx_opt[i].mpr_sa, - conf_info.vx_opt[i].mpr_ta); + conf_info.vx_opt[i].mpr_thr_inc_map); } // Setup the first pass through the data diff --git a/src/tools/core/grid_stat/grid_stat_conf_info.cc b/src/tools/core/grid_stat/grid_stat_conf_info.cc index f9b704e455..a049ceeab2 100644 --- a/src/tools/core/grid_stat/grid_stat_conf_info.cc +++ b/src/tools/core/grid_stat/grid_stat_conf_info.cc @@ -599,8 +599,7 @@ void GridStatVxOpt::clear() { var_name.clear(); var_suffix.clear(); - mpr_sa.clear(); - mpr_ta.clear(); + mpr_thr_inc_map.clear(); fcat_ta.clear(); ocat_ta.clear(); @@ -734,8 +733,26 @@ void GridStatVxOpt::process_config( for(i=0; i " + << "The length of \"" << conf_key_mpr_column << "\" and \"" + << conf_key_mpr_thresh << "\" must match (" << mpr_sa.n() + << " != " << mpr_ta.n() << ")!\n\n"; + exit(1); + } + + // Store in map + for(int i=0; mpr_sa.n(); i++) { + if(mpr_thr_inc_map.count(mpr_sa[i]) == 0) { + ThreshArray ta; + mpr_thr_inc_map[(mpr_sa[i])] = ta; + } + mpr_thr_inc_map[(mpr_sa[i])].add(mpr_ta[i]); + } // Conf: cat_thresh fcat_ta = fdict.lookup_thresh_array(conf_key_cat_thresh); @@ -1016,7 +1033,7 @@ bool GridStatVxOpt::is_uv_match(const GridStatVxOpt &v) const { // // The following do not impact matched pairs: // desc, var_name, var_suffix, - // mpr_sa, mpr_ta, + // mpr_thr_inc_map, // fcat_ta, ocat_ta, // fcnt_ta, ocnt_ta, cnt_logic, // fwind_ta, owind_ta, wind_logic, diff --git a/src/tools/core/grid_stat/grid_stat_conf_info.h b/src/tools/core/grid_stat/grid_stat_conf_info.h index e24326c7cd..c46eef12a1 100644 --- a/src/tools/core/grid_stat/grid_stat_conf_info.h +++ b/src/tools/core/grid_stat/grid_stat_conf_info.h @@ -149,8 +149,8 @@ class GridStatVxOpt { ConcatString var_suffix; // nc_pairs_var_suffix string // nc_pairs_var_str is deprecated - StringArray mpr_sa; // MPR filtering columns - ThreshArray mpr_ta; // MPR filtering thresholds + // Matched pair inclusion thresholds + std::map mpr_thr_inc_map; ThreshArray fcat_ta; // fcst categorical thresholds ThreshArray ocat_ta; // obs categorical thresholds diff --git a/src/tools/core/pair_stat/pair_stat.cc b/src/tools/core/pair_stat/pair_stat.cc index 9749aa6635..508883ef14 100644 --- a/src/tools/core/pair_stat/pair_stat.cc +++ b/src/tools/core/pair_stat/pair_stat.cc @@ -75,8 +75,6 @@ static void do_mcts (MCTSInfo &, int, const PairDataPoint *); static void do_cnt_sl1l2 (const PairStatVxOpt &, const PairDataPoint *); static void do_vl1l2 (VL1L2Info *&, int, const PairDataPoint *, const PairDataPoint *); static void do_pct (const PairStatVxOpt &, const PairDataPoint *); -static void do_hira_ens ( int, const PairDataPoint *); -static void do_hira_prob ( int, const PairDataPoint *); static void finish_txt_files(); @@ -165,6 +163,20 @@ void process_command_line(int argc, char **argv) { // Check for error. There should be no arguments left if(cline.n() != 0) usage(); + // Expand the input pairs file lists + StringArray tmp_src(pairs_files); + pairs_files.clear(); + + for(int i=0; i 0) pairs_files.add(sa); + else pairs_files.add(tmp_src[i]); + } + // Check for required argruments if(pairs_files.n() == 0) { mlog << Error << "\n" << method_name @@ -205,9 +217,11 @@ void process_command_line(int argc, char **argv) { // Process the configuration conf_info.process_config(pairs_format); - // List the pairs files + // List the input pair files mlog << Debug(1) - << "Pairs File(s): " << write_css(pairs_files) << "\n"; + << "Reading " << pairs_files.n() << " " + << pairsformat_to_string(pairs_format) << " Pairs File(s): " + << write_css(pairs_files) << "\n"; // Set the model name if(conf_info.model.empty()) { @@ -257,7 +271,7 @@ void setup_first_pass(const DataPlane &dp, const Grid &data_grid) { void setup_txt_files() { int max_col, max_prob_col, max_mctc_col, max_orank_col; - int n_prob, n_cat, n_eclv, n_ens; + int n_prob, n_cat, n_eclv; ConcatString base_name; // Create output file names for the stat file and optional text files @@ -270,15 +284,12 @@ void setup_txt_files() { ///////////////////////////////////////////////////////////////////// // Get the maximum number of data columns - n_prob = max(conf_info.get_max_n_fprob_thresh(), - conf_info.get_max_n_hira_prob()); + n_prob = conf_info.get_max_n_fprob_thresh(); n_cat = conf_info.get_max_n_cat_thresh() + 1; n_eclv = conf_info.get_max_n_eclv_points(); - n_ens = conf_info.get_max_n_hira_ens(); max_prob_col = get_n_pjc_columns(n_prob); max_mctc_col = get_n_mctc_columns(n_cat); - max_orank_col = get_n_orank_columns(n_ens); // Determine the maximum number of data columns max_col = (max_prob_col > max_stat_col ? max_prob_col : max_stat_col); @@ -356,10 +367,6 @@ void setup_txt_files() { max_col = get_n_eclv_columns(n_eclv) + n_header_columns + 1; break; - case i_orank: - max_col = get_n_orank_columns(n_ens) + n_header_columns + 1; - break; - default: max_col = n_txt_columns[i] + n_header_columns + 1; break; @@ -396,10 +403,6 @@ void setup_txt_files() { write_eclv_header_row(1, n_eclv, txt_at[i], 0, 0); break; - case i_orank: - write_orank_header_row(1, n_ens, txt_at[i], 0, 0); - break; - default: write_header_row(txt_columns[i], n_txt_columns[i], 1, txt_at[i], 0, 0); @@ -848,7 +851,6 @@ void process_scores() { ConcatString cs; // Initialize pointers - PairDataPoint *pd_ptr = (PairDataPoint *) nullptr; CTSInfo *cts_info = (CTSInfo *) nullptr; MCTSInfo mcts_info; VL1L2Info *vl1l2_info = (VL1L2Info *) nullptr; @@ -870,389 +872,320 @@ void process_scores() { // Compute scores for each PairData object and write output for(int i_vx=0; i_vxvx_pd.fcst_dpa.n_planes() == 0) continue; // Store the description - if(conf_info.vx_opt[i_vx].vx_pd.desc.empty()) { + if(vx_ptr->vx_pd.desc.empty()) { shc.set_desc(na_str); } else { - shc.set_desc(conf_info.vx_opt[i_vx].vx_pd.desc.c_str()); + shc.set_desc(vx_ptr->vx_pd.desc.c_str()); } // Store the forecast variable name - shc.set_fcst_var(conf_info.vx_opt[i_vx].vx_pd.fcst_info->name_attr()); + shc.set_fcst_var(vx_ptr->vx_pd.fcst_info->name_attr()); // Store the forecast variable units - shc.set_fcst_units(conf_info.vx_opt[i_vx].vx_pd.fcst_info->units_attr()); + shc.set_fcst_units(vx_ptr->vx_pd.fcst_info->units_attr()); // Set the forecast level name - shc.set_fcst_lev(conf_info.vx_opt[i_vx].vx_pd.fcst_info->level_attr().c_str()); + shc.set_fcst_lev(vx_ptr->vx_pd.fcst_info->level_attr().c_str()); // Store the observation variable name - shc.set_obs_var(conf_info.vx_opt[i_vx].vx_pd.obs_info->name_attr()); + shc.set_obs_var(vx_ptr->vx_pd.obs_info->name_attr()); // Store the observation variable units - cs = conf_info.vx_opt[i_vx].vx_pd.obs_info->units_attr(); + cs = vx_ptr->vx_pd.obs_info->units_attr(); if(cs.empty()) cs = na_string; shc.set_obs_units(cs); // Set the observation level name - shc.set_obs_lev(conf_info.vx_opt[i_vx].vx_pd.obs_info->level_attr().c_str()); + shc.set_obs_lev(vx_ptr->vx_pd.obs_info->level_attr().c_str()); // Set the forecast lead time - shc.set_fcst_lead_sec(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa[0].lead()); + shc.set_fcst_lead_sec(vx_ptr->vx_pd.fcst_dpa[0].lead()); // Set the forecast valid time - shc.set_fcst_valid_beg(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa[0].valid()); - shc.set_fcst_valid_end(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa[0].valid()); + shc.set_fcst_valid_beg(vx_ptr->vx_pd.fcst_dpa[0].valid()); + shc.set_fcst_valid_end(vx_ptr->vx_pd.fcst_dpa[0].valid()); // Set the observation lead time shc.set_obs_lead_sec(0); // Set the observation valid time - shc.set_obs_valid_beg(conf_info.vx_opt[i_vx].vx_pd.beg_ut); - shc.set_obs_valid_end(conf_info.vx_opt[i_vx].vx_pd.end_ut); - - // Loop through the message types - for(int i_msg_typ=0; i_msg_typmagic_str() - << " versus " - << conf_info.vx_opt[i_vx].vx_pd.obs_info->magic_str() - << ", for observation type " << pd_ptr->msg_typ - << ", over region " << pd_ptr->mask_name - << ", for interpolation method " - << shc.get_interp_mthd() << "(" - << shc.get_interp_pnts_str() - << "), using " << pd_ptr->n_obs << " matched pairs.\n"; - - // List counts for reasons why observations were rejected - cs << cs_erase - << "Number of matched pairs = " << pd_ptr->n_obs << "\n" - << "Observations processed = " << conf_info.vx_opt[i_vx].vx_pd.n_try << "\n" - << "Rejected: station id = " << conf_info.vx_opt[i_vx].vx_pd.rej_sid << "\n" - << "Rejected: obs var name = " << conf_info.vx_opt[i_vx].vx_pd.rej_var << "\n" - << "Rejected: valid time = " << conf_info.vx_opt[i_vx].vx_pd.rej_vld << "\n" - << "Rejected: bad obs value = " << conf_info.vx_opt[i_vx].vx_pd.rej_obs << "\n" - << "Rejected: off the grid = " << conf_info.vx_opt[i_vx].vx_pd.rej_grd << "\n" - << "Rejected: topography = " << conf_info.vx_opt[i_vx].vx_pd.rej_topo << "\n" - << "Rejected: level mismatch = " << conf_info.vx_opt[i_vx].vx_pd.rej_lvl << "\n" - << "Rejected: quality marker = " << conf_info.vx_opt[i_vx].vx_pd.rej_qty << "\n" - << "Rejected: message type = " << conf_info.vx_opt[i_vx].vx_pd.rej_typ[n] << "\n" - << "Rejected: masking region = " << conf_info.vx_opt[i_vx].vx_pd.rej_mask[n] << "\n" - << "Rejected: bad fcst value = " << conf_info.vx_opt[i_vx].vx_pd.rej_fcst[n] << "\n" - << "Rejected: bad climo mean = " << conf_info.vx_opt[i_vx].vx_pd.rej_cmn[n] << "\n" - << "Rejected: bad climo stdev = " << conf_info.vx_opt[i_vx].vx_pd.rej_csd[n] << "\n" - << "Rejected: mpr filter = " << conf_info.vx_opt[i_vx].vx_pd.rej_mpr[n] << "\n" - << "Rejected: duplicates = " << conf_info.vx_opt[i_vx].vx_pd.rej_dup[n] << "\n"; - - // Print report based on the number of matched pairs - if(pd_ptr->n_obs > 0) { - mlog << Debug(3) << cs; - } - // Continue for zero matched pairs - else { - mlog << Debug(2) << cs; - continue; - } + shc.set_obs_valid_beg(vx_ptr->vx_pd.beg_ut); + shc.set_obs_valid_end(vx_ptr->vx_pd.end_ut); - // Process percentile thresholds - conf_info.vx_opt[i_vx].set_perc_thresh(pd_ptr); + // Store the message type in the obtype column + shc.set_obtype(na_str); - // Write out the MPR lines - if(conf_info.vx_opt[i_vx].output_flag[i_mpr] != STATOutputType::None) { - write_mpr_row(shc, pd_ptr, - conf_info.vx_opt[i_vx].output_flag[i_mpr], - stat_at, i_stat_row, - txt_at[i_mpr], i_txt_row[i_mpr], false); + // Loop through the verification masking regions + for(int i_mask=0; i_maskget_n_mask(); i_mask++) { - // Reset the obtype column - shc.set_obtype(conf_info.vx_opt[i_vx].msg_typ[i_msg_typ].c_str()); + // Store the verification masking region + shc.set_mask(vx_ptr->mask_name[i_mask].c_str()); - // Reset the observation valid time - shc.set_obs_valid_beg(conf_info.vx_opt[i_vx].vx_pd.beg_ut); - shc.set_obs_valid_end(conf_info.vx_opt[i_vx].vx_pd.end_ut); - } + // Store the interpolation method as nearest + shc.set_interp_mthd(InterpMthd::Nearest); + shc.set_interp_wdth(1); - // Write out the SEEPS MPR lines - if(conf_info.vx_opt[i_vx].output_flag[i_seeps_mpr] != STATOutputType::None) { - write_seeps_mpr_row(shc, pd_ptr, - conf_info.vx_opt[i_vx].output_flag[i_seeps_mpr], - stat_at, i_stat_row, - txt_at[i_seeps_mpr], i_txt_row[i_seeps_mpr], false); + PairDataPoint *pd_ptr = &vx_ptr->vx_pd.pd[i_mask]; + + mlog << Debug(2) + << "Processing " << vx_ptr->vx_pd.fcst_info->magic_str() + << " versus " << vx_ptr->vx_pd.obs_info->magic_str() + << ", over region " << pd_ptr->mask_name + << "), using " << pd_ptr->n_obs << " matched pairs.\n"; + + // List counts for reasons why observations were rejected + cs << cs_erase + << "Number of matched pairs = " << pd_ptr->n_obs << "\n" + << "Observations processed = " << vx_ptr->vx_pd.n_try << "\n" + << "Rejected: station id = " << vx_ptr->vx_pd.rej_sid << "\n" + << "Rejected: obs var name = " << vx_ptr->vx_pd.rej_var << "\n" + << "Rejected: valid time = " << vx_ptr->vx_pd.rej_vld << "\n" + << "Rejected: bad obs value = " << vx_ptr->vx_pd.rej_obs << "\n" + << "Rejected: off the grid = " << vx_ptr->vx_pd.rej_grd << "\n" + << "Rejected: topography = " << vx_ptr->vx_pd.rej_topo << "\n" + << "Rejected: level mismatch = " << vx_ptr->vx_pd.rej_lvl << "\n" + << "Rejected: quality marker = " << vx_ptr->vx_pd.rej_qty << "\n" + << "Rejected: message type = " << vx_ptr->vx_pd.rej_typ[i_mask] << "\n" + << "Rejected: masking region = " << vx_ptr->vx_pd.rej_mask[i_mask] << "\n" + << "Rejected: bad fcst value = " << vx_ptr->vx_pd.rej_fcst[i_mask] << "\n" + << "Rejected: bad climo mean = " << vx_ptr->vx_pd.rej_cmn[i_mask] << "\n" + << "Rejected: bad climo stdev = " << vx_ptr->vx_pd.rej_csd[i_mask] << "\n" + << "Rejected: mpr filter = " << vx_ptr->vx_pd.rej_mpr[i_mask] << "\n" + << "Rejected: duplicates = " << vx_ptr->vx_pd.rej_dup[i_mask] << "\n"; + + // Print report based on the number of matched pairs + if(pd_ptr->n_obs > 0) mlog << Debug(3) << cs; + else mlog << Debug(2) << cs; + + // Process percentile thresholds + vx_ptr->set_perc_thresh(pd_ptr); + + // Write out the MPR lines + if(vx_ptr->output_flag[i_mpr] != STATOutputType::None) { + write_mpr_row(shc, pd_ptr, + vx_ptr->output_flag[i_mpr], + stat_at, i_stat_row, + txt_at[i_mpr], i_txt_row[i_mpr], false); + + // Reset the obtype column + shc.set_obtype(na_str); + + // Reset the observation valid time + shc.set_obs_valid_beg(vx_ptr->vx_pd.beg_ut); + shc.set_obs_valid_end(vx_ptr->vx_pd.end_ut); + } + + // Write out the SEEPS MPR lines + if(vx_ptr->output_flag[i_seeps_mpr] != STATOutputType::None) { + write_seeps_mpr_row(shc, pd_ptr, + vx_ptr->output_flag[i_seeps_mpr], + stat_at, i_stat_row, + txt_at[i_seeps_mpr], i_txt_row[i_seeps_mpr], false); + + // Reset the obtype column + shc.set_obtype(na_str); + + // Reset the observation valid time + shc.set_obs_valid_beg(vx_ptr->vx_pd.beg_ut); + shc.set_obs_valid_end(vx_ptr->vx_pd.end_ut); + } + + // Write out the SEEPS lines + if(vx_ptr->output_flag[i_seeps] != STATOutputType::None) { + compute_aggregated_seeps(pd_ptr, &pd_ptr->seeps_agg); + write_seeps_row(shc, &pd_ptr->seeps_agg, + vx_ptr->output_flag[i_seeps], + stat_at, i_stat_row, + txt_at[i_seeps], i_txt_row[i_seeps]); + } + + // Compute CTS scores + if(!vx_ptr->vx_pd.fcst_info->is_prob() && + vx_ptr->fcat_ta.n() > 0 && + (vx_ptr->output_flag[i_fho] != STATOutputType::None || + vx_ptr->output_flag[i_ctc] != STATOutputType::None || + vx_ptr->output_flag[i_cts] != STATOutputType::None || + vx_ptr->output_flag[i_eclv] != STATOutputType::None)) { + + // Initialize + for(int i_cat=0; i_catfcat_ta.n(); i_cat++) { + + if(cts_info[i_cat].cts.n_pairs() == 0) continue; + + // Write out FHO + if(vx_ptr->output_flag[i_fho] != STATOutputType::None) { + write_fho_row(shc, cts_info[i_cat], + vx_ptr->output_flag[i_fho], + stat_at, i_stat_row, + txt_at[i_fho], i_txt_row[i_fho]); } - // Write out the SEEPS lines - if(conf_info.vx_opt[i_vx].output_flag[i_seeps] != STATOutputType::None) { - compute_aggregated_seeps(pd_ptr, &pd_ptr->seeps_agg); - write_seeps_row(shc, &pd_ptr->seeps_agg, - conf_info.vx_opt[i_vx].output_flag[i_seeps], + // Write out CTC + if(vx_ptr->output_flag[i_ctc] != STATOutputType::None) { + write_ctc_row(shc, cts_info[i_cat], + vx_ptr->output_flag[i_ctc], stat_at, i_stat_row, - txt_at[i_seeps], i_txt_row[i_seeps]); + txt_at[i_ctc], i_txt_row[i_ctc]); } - // Compute CTS scores - if(!conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_prob() && - conf_info.vx_opt[i_vx].fcat_ta.n() > 0 && - (conf_info.vx_opt[i_vx].output_flag[i_fho] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_ctc] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_cts] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_eclv] != STATOutputType::None)) { - - // Initialize - for(int i_cat=0; i_catis_prob() && - conf_info.vx_opt[i_vx].fcat_ta.n() > 1 && - (conf_info.vx_opt[i_vx].output_flag[i_mctc] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_mcts] != STATOutputType::None)) { - - // Initialize - mcts_info.clear(); - - // Compute MCTS Info - do_mcts(mcts_info, i_vx, pd_ptr); - - if(mcts_info.cts.n_pairs() == 0) continue; - - // Write out MCTC - if(conf_info.vx_opt[i_vx].output_flag[i_mctc] != STATOutputType::None) { - write_mctc_row(shc, mcts_info, - conf_info.vx_opt[i_vx].output_flag[i_mctc], - stat_at, i_stat_row, - txt_at[i_mctc], i_txt_row[i_mctc]); - } - - // Write out MCTS - if(conf_info.vx_opt[i_vx].output_flag[i_mcts] != STATOutputType::None) { - write_mcts_row(shc, mcts_info, - conf_info.vx_opt[i_vx].output_flag[i_mcts], - stat_at, i_stat_row, - txt_at[i_mcts], i_txt_row[i_mcts]); - } - } // end Compute MCTS scores - - // Compute CNT, SL1L2, and SAL1L2 scores - if(!conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_prob() && - (conf_info.vx_opt[i_vx].output_flag[i_cnt] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_sl1l2] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_sal1l2] != STATOutputType::None)) { - do_cnt_sl1l2(conf_info.vx_opt[i_vx], pd_ptr); + // Write out CTS + if(vx_ptr->output_flag[i_cts] != STATOutputType::None) { + write_cts_row(shc, cts_info[i_cat], + vx_ptr->output_flag[i_cts], + stat_at, i_stat_row, + txt_at[i_cts], i_txt_row[i_cts]); } - // Compute VL1L2 and VAL1L2 partial sums for UGRD and VGRD - if(!conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_prob() && - conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_v_wind() && - conf_info.vx_opt[i_vx].vx_pd.fcst_info->uv_index() >= 0 && - (conf_info.vx_opt[i_vx].output_flag[i_vl1l2] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_val1l2] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_vcnt] != STATOutputType::None)) { - - // Store the forecast variable name - shc.set_fcst_var(ugrd_vgrd_abbr_str); - - // Store the observation variable name - shc.set_obs_var(ugrd_vgrd_abbr_str); - - // Initialize - for(int i_wind=0; i_winduv_index(); - - // Check to make sure message types, masking regions, - // and interpolation methods match - if(conf_info.vx_opt[i_vx].get_n_msg_typ() != - conf_info.vx_opt[u_vx].get_n_msg_typ() || - conf_info.vx_opt[i_vx].get_n_mask() != - conf_info.vx_opt[u_vx].get_n_mask() || - conf_info.vx_opt[i_vx].get_n_interp() != - conf_info.vx_opt[u_vx].get_n_interp()) { - mlog << Warning << "\nprocess_scores() -> " - << "when computing VL1L2 and/or VAL1L2 vector " - << "partial sums, the U and V components must " - << "be processed using the same set of message " - << "types, masking regions, and interpolation " - << "methods. Failing to do so will cause " - << "unexpected results!\n\n"; - } - - // Compute VL1L2 and VAL1L2 - do_vl1l2(vl1l2_info, i_vx, - &conf_info.vx_opt[u_vx].vx_pd.pd[n], - &conf_info.vx_opt[i_vx].vx_pd.pd[n]); - - // Loop through all of the wind speed thresholds - for(int i_wind=0; i_wind 0) { - write_vl1l2_row(shc, vl1l2_info[i_wind], - conf_info.vx_opt[i_vx].output_flag[i_vl1l2], - stat_at, i_stat_row, - txt_at[i_vl1l2], i_txt_row[i_vl1l2]); - } - - // Write out VAL1L2 - if(conf_info.vx_opt[i_vx].output_flag[i_val1l2] != STATOutputType::None && - vl1l2_info[i_wind].vacount > 0) { - write_val1l2_row(shc, vl1l2_info[i_wind], - conf_info.vx_opt[i_vx].output_flag[i_val1l2], - stat_at, i_stat_row, - txt_at[i_val1l2], i_txt_row[i_val1l2]); - } - - // Write out VCNT - if(conf_info.vx_opt[i_vx].output_flag[i_vcnt] != STATOutputType::None && - vl1l2_info[i_wind].vcount > 0) { - write_vcnt_row(shc, vl1l2_info[i_wind], - conf_info.vx_opt[i_vx].output_flag[i_vcnt], - stat_at, i_stat_row, - txt_at[i_vcnt], i_txt_row[i_vcnt]); - } - - } // end for i - - // Reset the forecast variable name - shc.set_fcst_var(conf_info.vx_opt[i_vx].vx_pd.fcst_info->name_attr()); - - // Reset the observation variable name - shc.set_obs_var(conf_info.vx_opt[i_vx].vx_pd.obs_info->name_attr()); - - } // end Compute VL1L2 and VAL1L2 - - // Compute PCT counts and scores - if(conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_prob() && - (conf_info.vx_opt[i_vx].output_flag[i_pct] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_pstd] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_pjc] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_prc] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_eclv] != STATOutputType::None)) { - do_pct(conf_info.vx_opt[i_vx], pd_ptr); + // Write out ECLV + if(vx_ptr->output_flag[i_eclv] != STATOutputType::None) { + write_eclv_row(shc, cts_info[i_cat], vx_ptr->eclv_points, + vx_ptr->output_flag[i_eclv], + stat_at, i_stat_row, + txt_at[i_eclv], i_txt_row[i_eclv]); } + } // end for i_cat + } // end Compute CTS scores - // Reset the verification masking region - shc.set_mask(conf_info.vx_opt[i_vx].mask_name[i_mask].c_str()); + // Compute MCTS scores + if(!vx_ptr->vx_pd.fcst_info->is_prob() && + vx_ptr->fcat_ta.n() > 1 && + (vx_ptr->output_flag[i_mctc] != STATOutputType::None || + vx_ptr->output_flag[i_mcts] != STATOutputType::None)) { - } // end for i_interp + // Initialize + mcts_info.clear(); - // Apply HiRA ensemble verification logic - if(!conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_prob() && - conf_info.vx_opt[i_vx].hira_info.flag && - (conf_info.vx_opt[i_vx].output_flag[i_ecnt] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_rps] != STATOutputType::None)) { + // Compute MCTS Info + do_mcts(mcts_info, i_vx, pd_ptr); - int n = conf_info.vx_opt[i_vx].vx_pd.three_to_one(i_msg_typ, i_mask, 0); + if(mcts_info.cts.n_pairs() == 0) continue; - pd_ptr = &conf_info.vx_opt[i_vx].vx_pd.pd[n]; + // Write out MCTC + if(vx_ptr->output_flag[i_mctc] != STATOutputType::None) { + write_mctc_row(shc, mcts_info, + vx_ptr->output_flag[i_mctc], + stat_at, i_stat_row, + txt_at[i_mctc], i_txt_row[i_mctc]); + } - // Process percentile thresholds - conf_info.vx_opt[i_vx].set_perc_thresh(pd_ptr); + // Write out MCTS + if(vx_ptr->output_flag[i_mcts] != STATOutputType::None) { + write_mcts_row(shc, mcts_info, + vx_ptr->output_flag[i_mcts], + stat_at, i_stat_row, + txt_at[i_mcts], i_txt_row[i_mcts]); + } + } // end Compute MCTS scores + + // Compute CNT, SL1L2, and SAL1L2 scores + if(!vx_ptr->vx_pd.fcst_info->is_prob() && + (vx_ptr->output_flag[i_cnt] != STATOutputType::None || + vx_ptr->output_flag[i_sl1l2] != STATOutputType::None || + vx_ptr->output_flag[i_sal1l2] != STATOutputType::None)) { + do_cnt_sl1l2(*vx_ptr, pd_ptr); + } - // Appy HiRA verification and write ensemble output - do_hira_ens(i_vx, pd_ptr); + // Compute VL1L2 and VAL1L2 partial sums for UGRD and VGRD + if(!vx_ptr->vx_pd.fcst_info->is_prob() && + vx_ptr->vx_pd.fcst_info->is_v_wind() && + vx_ptr->vx_pd.fcst_info->uv_index() >= 0 && + (vx_ptr->output_flag[i_vl1l2] != STATOutputType::None || + vx_ptr->output_flag[i_val1l2] != STATOutputType::None || + vx_ptr->output_flag[i_vcnt] != STATOutputType::None)) { + + // Store the forecast variable name + shc.set_fcst_var(ugrd_vgrd_abbr_str); + + // Store the observation variable name + shc.set_obs_var(ugrd_vgrd_abbr_str); + + // Initialize + for(int i_wind=0; i_windvx_pd.fcst_info->uv_index(); + + // Check to make sure the masking regions match + if(conf_info.vx_opt[i_vx].get_n_mask() != + conf_info.vx_opt[u_vx].get_n_mask()) { + mlog << Warning << "\nprocess_scores() -> " + << "when computing VL1L2 and/or VAL1L2 vector " + << "partial sums, the U and V components must " + << "be processed using the same set of mask regions. " + << "Failing to do so will cause unexpected results!\n\n"; + } - } // end HiRA for ensembles + // Compute VL1L2 and VAL1L2 + do_vl1l2(vl1l2_info, i_vx, + &conf_info.vx_opt[u_vx].vx_pd.pd[i_mask], + &conf_info.vx_opt[i_vx].vx_pd.pd[i_mask]); - // Apply HiRA probabilistic verification logic - if(!conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_prob() && - conf_info.vx_opt[i_vx].hira_info.flag && - (conf_info.vx_opt[i_vx].output_flag[i_mpr] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_pct] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_pstd] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_pjc] != STATOutputType::None || - conf_info.vx_opt[i_vx].output_flag[i_prc] != STATOutputType::None)) { + // Loop through all of the wind speed thresholds + for(int i_wind=0; i_windfwind_ta.n(); i_wind++) { - int n = conf_info.vx_opt[i_vx].vx_pd.three_to_one(i_msg_typ, i_mask, 0); + // Write out VL1L2 + if(vx_ptr->output_flag[i_vl1l2] != STATOutputType::None && + vl1l2_info[i_wind].vcount > 0) { + write_vl1l2_row(shc, vl1l2_info[i_wind], + vx_ptr->output_flag[i_vl1l2], + stat_at, i_stat_row, + txt_at[i_vl1l2], i_txt_row[i_vl1l2]); + } - pd_ptr = &conf_info.vx_opt[i_vx].vx_pd.pd[n]; + // Write out VAL1L2 + if(vx_ptr->output_flag[i_val1l2] != STATOutputType::None && + vl1l2_info[i_wind].vacount > 0) { + write_val1l2_row(shc, vl1l2_info[i_wind], + vx_ptr->output_flag[i_val1l2], + stat_at, i_stat_row, + txt_at[i_val1l2], i_txt_row[i_val1l2]); + } - // Process percentile thresholds - conf_info.vx_opt[i_vx].set_perc_thresh(pd_ptr); + // Write out VCNT + if(vx_ptr->output_flag[i_vcnt] != STATOutputType::None && + vl1l2_info[i_wind].vcount > 0) { + write_vcnt_row(shc, vl1l2_info[i_wind], + vx_ptr->output_flag[i_vcnt], + stat_at, i_stat_row, + txt_at[i_vcnt], i_txt_row[i_vcnt]); + } + } // end for i_wind - // Apply HiRA verification and write probabilistic output - do_hira_prob(i_vx, pd_ptr); + // Reset the forecast variable name + shc.set_fcst_var(vx_ptr->vx_pd.fcst_info->name_attr()); - } // end HiRA for probabilities + // Reset the observation variable name + shc.set_obs_var(vx_ptr->vx_pd.obs_info->name_attr()); - } // end for i_mask - } // end for i_msg_typ + } // end Compute VL1L2 and VAL1L2 + + // Compute PCT counts and scores + if(vx_ptr->vx_pd.fcst_info->is_prob() && + (vx_ptr->output_flag[i_pct] != STATOutputType::None || + vx_ptr->output_flag[i_pstd] != STATOutputType::None || + vx_ptr->output_flag[i_pjc] != STATOutputType::None || + vx_ptr->output_flag[i_prc] != STATOutputType::None || + vx_ptr->output_flag[i_eclv] != STATOutputType::None)) { + do_pct(conf_info.vx_opt[i_vx], pd_ptr); + } + + // Reset the verification masking region + shc.set_mask(vx_ptr->mask_name[i_mask].c_str()); + + } // end for i_mask mlog << Debug(2) << "\n" << sep_str << "\n\n"; @@ -1730,377 +1663,6 @@ void do_pct(const PairStatVxOpt &vx_opt, const PairDataPoint *pd_ptr) { //////////////////////////////////////////////////////////////////////// -void do_hira_ens(int i_vx, const PairDataPoint *pd_ptr) { - PairDataEnsemble hira_pd; - int i, j, k, lvl_blw, lvl_abv; - NumArray f_ens; - - // Set flag for specific humidity - bool spfh_flag = conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_specific_humidity() && - conf_info.vx_opt[i_vx].vx_pd.obs_info->is_specific_humidity(); - - shc.set_interp_mthd(InterpMthd::Nbrhd, - conf_info.vx_opt[i_vx].hira_info.shape); - - // Loop over the HiRA widths - for(i=0; i " - << "failed to get GridTemplate for " << i << "-th width.\n\n"; - continue; - } - - // Initialize - hira_pd.clear(); - hira_pd.extend(pd_ptr->n_obs); - hira_pd.set_ens_size(gt->size()); - hira_pd.set_climo_cdf_info_ptr(&conf_info.vx_opt[i_vx].cdf_info); - f_ens.extend(gt->size()); - - // Process each observation point - for(j=0; jn_obs; j++) { - - // Determine the forecast level values - find_vert_lvl(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa, - pd_ptr->lvl_na[j], lvl_blw, lvl_abv); - - // Get the nearby forecast values - get_interp_points(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa, - pd_ptr->x_na[j], pd_ptr->y_na[j], - InterpMthd::Nbrhd, conf_info.vx_opt[i_vx].hira_info.width[i], - conf_info.vx_opt[i_vx].hira_info.shape, grid.wrap_lon(), - conf_info.vx_opt[i_vx].hira_info.vld_thresh, spfh_flag, - conf_info.vx_opt[i_vx].vx_pd.fcst_info->level().type(), - pd_ptr->lvl_na[j], lvl_blw, lvl_abv, f_ens); - - // Check for values - if(f_ens.n() == 0) continue; - - // TODO: Add has_climo member function instead - - // Skip points where climatology has been specified but is bad data - if((conf_info.vx_opt[i_vx].vx_pd.fcmn_dpa.n_planes() > 0 && - is_bad_data(pd_ptr->fcmn_na[j])) || - (conf_info.vx_opt[i_vx].vx_pd.ocmn_dpa.n_planes() > 0 && - is_bad_data(pd_ptr->ocmn_na[j]))) continue; - - // Store climo data - ClimoPntInfo cpi(pd_ptr->fcmn_na[j], pd_ptr->fcsd_na[j], - pd_ptr->ocmn_na[j], pd_ptr->ocsd_na[j]); - - // Store the observation value - hira_pd.add_point_obs( - pd_ptr->typ_sa[j].c_str(), pd_ptr->sid_sa[j].c_str(), - pd_ptr->lat_na[j], pd_ptr->lon_na[j], - pd_ptr->x_na[j], pd_ptr->y_na[j], pd_ptr->vld_ta[j], - pd_ptr->lvl_na[j], pd_ptr->elv_na[j], - pd_ptr->o_na[j], pd_ptr->o_qc_sa[j].c_str(), - cpi, pd_ptr->wgt_na[j]); - - // Store the ensemble mean and member values - hira_pd.mn_na.add(f_ens.mean()); - for(k=0; kmagic_str() - << " versus " - << conf_info.vx_opt[i_vx].vx_pd.obs_info->magic_str() - << ", for observation type " << pd_ptr->msg_typ - << ", over region " << pd_ptr->mask_name - << ", for interpolation method HiRA Ensemble NBRHD(" - << shc.get_interp_pnts_str() - << "), using " << hira_pd.n_obs << " matched pairs.\n"; - - // Check for zero matched pairs - if(hira_pd.o_na.n() == 0) { - if(gt) { delete gt; gt = nullptr; } - continue; - } - - // Compute the pair values - hira_pd.compute_pair_vals(rng_ptr); - - // Write out the ECNT line - if(conf_info.vx_opt[i_vx].output_flag[i_ecnt] != STATOutputType::None) { - - // Compute ensemble statistics - ECNTInfo ecnt_info; - ecnt_info.set(hira_pd); - - write_ecnt_row(shc, ecnt_info, - conf_info.vx_opt[i_vx].output_flag[i_ecnt], - stat_at, i_stat_row, - txt_at[i_ecnt], i_txt_row[i_ecnt]); - } // end if ECNT - - // Write out the ORANK line - if(conf_info.vx_opt[i_vx].output_flag[i_orank] != STATOutputType::None) { - - write_orank_row(shc, &hira_pd, - conf_info.vx_opt[i_vx].output_flag[i_orank], - stat_at, i_stat_row, - txt_at[i_orank], i_txt_row[i_orank], false); - - // Reset the obtype column - shc.set_obtype(pd_ptr->msg_typ.c_str()); - - // Reset the observation valid time - shc.set_obs_valid_beg(conf_info.vx_opt[i_vx].vx_pd.beg_ut); - shc.set_obs_valid_end(conf_info.vx_opt[i_vx].vx_pd.end_ut); - } // end if ORANK - - // Write out the RPS line - if(conf_info.vx_opt[i_vx].output_flag[i_rps] != STATOutputType::None) { - - // Store ensemble RPS thresholds - RPSInfo rps_info; - rps_info.set_prob_cat_thresh(conf_info.vx_opt[i_vx].hira_info.prob_cat_ta); - - // If prob_cat_thresh is empty, try to select other thresholds - if(rps_info.fthresh.n() == 0) { - - // Use observation climo data, if avaiable - if(hira_pd.ocmn_na.n_valid() > 0 && - hira_pd.ocsd_na.n_valid() > 0 && - conf_info.vx_opt[i_vx].cdf_info.cdf_ta.n() > 0) { - mlog << Debug(3) << "Resetting the empty HiRA \"" - << conf_key_prob_cat_thresh << "\" thresholds to " - << "climatological distribution thresholds.\n"; - rps_info.set_cdp_thresh(conf_info.vx_opt[i_vx].cdf_info.cdf_ta); - } - // Otherwise, use categorical observation thresholds - else { - mlog << Debug(3) << "Resetting the empty HiRA \"" - << conf_key_prob_cat_thresh << "\" thresholds to the " - << "observed categorical thresholds.\n"; - rps_info.set_prob_cat_thresh(conf_info.vx_opt[i_vx].ocat_ta); - } - } - - // Check for no thresholds - if(rps_info.fthresh.n() == 0) { - mlog << Debug(3) << "Skipping HiRA RPS output since no " - << "\"" << conf_key_prob_cat_thresh << "\" thresholds are " - << "defined in the \"" << conf_key_hira - << "\" dictionary.\n"; - if(gt) { delete gt; gt = nullptr; } - break; - } - - // Compute ensemble RPS statistics - rps_info.set(hira_pd); - - write_rps_row(shc, rps_info, - conf_info.vx_opt[i_vx].output_flag[i_rps], - stat_at, i_stat_row, - txt_at[i_rps], i_txt_row[i_rps]); - } // end if RPS - - if(gt) { delete gt; gt = nullptr; } - - } // end for i - - return; -} - -//////////////////////////////////////////////////////////////////////// - -void do_hira_prob(int i_vx, const PairDataPoint *pd_ptr) { - PairDataPoint hira_pd; - int i, j, k, lvl_blw, lvl_abv; - double f_cov, ocmn_cov; - NumArray ocmn_cov_na; - SingleThresh cat_thresh; - PCTInfo pct_info; - - // Set flag for specific humidity - bool spfh_flag = conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_specific_humidity() && - conf_info.vx_opt[i_vx].vx_pd.obs_info->is_specific_humidity(); - bool precip_flag = conf_info.vx_opt[i_vx].vx_pd.fcst_info->is_precipitation() && - conf_info.vx_opt[i_vx].vx_pd.obs_info->is_precipitation(); - - shc.set_interp_mthd(InterpMthd::Nbrhd, - conf_info.vx_opt[i_vx].hira_info.shape); - - // Loop over categorical thresholds and HiRA widths - for(i=0; in_obs; k++) { - - // Store climo data - ClimoPntInfo cpi(pd_ptr->fcmn_na[k], pd_ptr->fcsd_na[k], - pd_ptr->ocmn_na[k], pd_ptr->ocsd_na[k]); - - // Compute the fractional coverage forecast value using the - // observation level value - find_vert_lvl(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa, - pd_ptr->lvl_na[k], lvl_blw, lvl_abv); - - f_cov = compute_interp(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa, - pd_ptr->x_na[k], pd_ptr->y_na[k], pd_ptr->o_na[k], &cpi, - InterpMthd::Nbrhd, conf_info.vx_opt[i_vx].hira_info.width[j], - conf_info.vx_opt[i_vx].hira_info.shape, grid.wrap_lon(), - conf_info.vx_opt[i_vx].hira_info.vld_thresh, spfh_flag, - conf_info.vx_opt[i_vx].vx_pd.fcst_info->level().type(), - pd_ptr->lvl_na[k], lvl_blw, lvl_abv, &cat_thresh); - - // Check for bad data - if(is_bad_data(f_cov)) continue; - - // Compute the climatological event probability as the fractional - // coverage of the observation climatology mean field - if(conf_info.vx_opt[i_vx].vx_pd.ocmn_dpa.n_planes() > 0) { - - // Interpolate to the observation level - find_vert_lvl(conf_info.vx_opt[i_vx].vx_pd.ocmn_dpa, - pd_ptr->lvl_na[k], lvl_blw, lvl_abv); - - ocmn_cov = compute_interp(conf_info.vx_opt[i_vx].vx_pd.ocmn_dpa, - pd_ptr->x_na[k], pd_ptr->y_na[k], pd_ptr->o_na[k], &cpi, - InterpMthd::Nbrhd, conf_info.vx_opt[i_vx].hira_info.width[j], - conf_info.vx_opt[i_vx].hira_info.shape, grid.wrap_lon(), - conf_info.vx_opt[i_vx].hira_info.vld_thresh, spfh_flag, - conf_info.vx_opt[i_vx].vx_pd.fcst_info->level().type(), - pd_ptr->lvl_na[k], lvl_blw, lvl_abv, &cat_thresh); - - // Check for bad data - if(is_bad_data(ocmn_cov)) continue; - else ocmn_cov_na.add(ocmn_cov); - } - - // Store the fractional coverage pair - hira_pd.add_point_pair( - pd_ptr->typ_sa[k].c_str(), - pd_ptr->sid_sa[k].c_str(), - pd_ptr->lat_na[k], pd_ptr->lon_na[k], - pd_ptr->x_na[k], pd_ptr->y_na[k], pd_ptr->vld_ta[k], - pd_ptr->lvl_na[k], pd_ptr->elv_na[k], - f_cov, pd_ptr->o_na[k], pd_ptr->o_qc_sa[k].c_str(), - cpi, pd_ptr->wgt_na[k]); - } // end for k - - mlog << Debug(2) - << "Processing " - << conf_info.vx_opt[i_vx].vx_pd.fcst_info->magic_str() - << conf_info.vx_opt[i_vx].fcat_ta[i].get_str() - << " versus " - << conf_info.vx_opt[i_vx].vx_pd.obs_info->magic_str() - << conf_info.vx_opt[i_vx].ocat_ta[i].get_str() - << ", for observation type " << pd_ptr->msg_typ - << ", over region " << pd_ptr->mask_name - << ", for interpolation method HiRA Probability NBRHD(" - << shc.get_interp_pnts_str() - << "), using " << hira_pd.n_obs << " matched pairs.\n"; - - // Check for zero matched pairs - if(hira_pd.f_na.n() == 0 || hira_pd.o_na.n() == 0) continue; - - // Set up the PCTInfo thresholds and alpha values - pct_info.fthresh = conf_info.vx_opt[i_vx].hira_info.cov_ta; - pct_info.othresh = conf_info.vx_opt[i_vx].ocat_ta[i]; - pct_info.allocate_n_alpha(conf_info.vx_opt[i_vx].get_n_ci_alpha()); - - for(k=0; kmsg_typ.c_str()); - - // Reset the observation valid time - shc.set_obs_valid_beg(conf_info.vx_opt[i_vx].vx_pd.beg_ut); - shc.set_obs_valid_end(conf_info.vx_opt[i_vx].vx_pd.end_ut); - } - - // Set cov_thresh column using the HiRA coverage thresholds - shc.set_cov_thresh(conf_info.vx_opt[i_vx].hira_info.cov_ta); - - // Write out PCT - if(conf_info.vx_opt[i_vx].output_flag[i_pct] != STATOutputType::None) { - write_pct_row(shc, pct_info, - conf_info.vx_opt[i_vx].output_flag[i_pct],1, 1, - stat_at, i_stat_row, - txt_at[i_pct], i_txt_row[i_pct], false); - } - - // Write out PSTD - if(conf_info.vx_opt[i_vx].output_flag[i_pstd] != STATOutputType::None) { - write_pstd_row(shc, pct_info, - conf_info.vx_opt[i_vx].output_flag[i_pstd], 1, 1, - stat_at, i_stat_row, - txt_at[i_pstd], i_txt_row[i_pstd], false); - } - - // Write out PJC - if(conf_info.vx_opt[i_vx].output_flag[i_pjc] != STATOutputType::None) { - write_pjc_row(shc, pct_info, - conf_info.vx_opt[i_vx].output_flag[i_pjc], 1, 1, - stat_at, i_stat_row, - txt_at[i_pjc], i_txt_row[i_pjc], false); - } - - // Write out PRC - if(conf_info.vx_opt[i_vx].output_flag[i_prc] != STATOutputType::None) { - write_prc_row(shc, pct_info, - conf_info.vx_opt[i_vx].output_flag[i_prc], 1, 1, - stat_at, i_stat_row, - txt_at[i_prc], i_txt_row[i_prc], false); - } - - } // end for j - } // end for i - - return; -} - -//////////////////////////////////////////////////////////////////////// - void finish_txt_files() { int i; @@ -2152,18 +1714,17 @@ void usage() { << ") ***\n\n" << "Usage: " << program_name << "\n" - << "\t-pairs file\n" + << "\t-pairs file_1 ... file_n | file_list\n" << "\t-format type\n" << "\t-config config_file\n" << "\t[-outdir path]\n" << "\t[-log file]\n" << "\t[-v level]\n\n" - -// JHG change -pairs file to -pairs file_list to support a long list of inputs - << "\twhere\t\"-pairs file\" is one or more files containing " - << "forecast/observation pairs. May be used multiple times " - << "(required).\n" + << "\twhere\t\"-pairs\" defines one or more input files containing " + << "forecast/observation pairs. May be set as a list of file names " + << "(file_1 ... file_n) or as an ASCII file containing a list of " + << "file names (file_list). May be used multiple times (required)." << "\t\t\"-format type\" defines the input pairs file format " << "and may be set to \"mpr\" or \"ioda\" (required).\n" diff --git a/src/tools/core/pair_stat/pair_stat.h b/src/tools/core/pair_stat/pair_stat.h index 7d35492d0b..4186811a11 100644 --- a/src/tools/core/pair_stat/pair_stat.h +++ b/src/tools/core/pair_stat/pair_stat.h @@ -60,38 +60,35 @@ static const char * ioda_data_config_filename = // Header columns static const char * const * txt_columns[n_txt] = { - fho_columns, ctc_columns, cts_columns, - mctc_columns, mcts_columns, cnt_columns, - sl1l2_columns, sal1l2_columns, vl1l2_columns, - val1l2_columns, pct_columns, pstd_columns, - pjc_columns, prc_columns, ecnt_columns, - orank_columns, rps_columns, eclv_columns, - mpr_columns, vcnt_columns, seeps_mpr_columns, - seeps_columns + fho_columns, ctc_columns, cts_columns, + mctc_columns, mcts_columns, + cnt_columns, sl1l2_columns, sal1l2_columns, + vcnt_columns, vl1l2_columns, val1l2_columns, + pct_columns, pstd_columns, pjc_columns, + prc_columns, eclv_columns, + mpr_columns, seeps_mpr_columns, seeps_columns }; // Length of header columns static const int n_txt_columns[n_txt] = { - n_fho_columns, n_ctc_columns, n_cts_columns, - n_mctc_columns, n_mcts_columns, n_cnt_columns, - n_sl1l2_columns, n_sal1l2_columns, n_vl1l2_columns, - n_val1l2_columns, n_pct_columns, n_pstd_columns, - n_pjc_columns, n_prc_columns, n_ecnt_columns, - n_orank_columns, n_rps_columns, n_eclv_columns, - n_mpr_columns, n_vcnt_columns, n_seeps_mpr_columns, - n_seeps_columns + n_fho_columns, n_ctc_columns, n_cts_columns, + n_mctc_columns, n_mcts_columns, + n_cnt_columns, n_sl1l2_columns, n_sal1l2_columns, + n_vcnt_columns, n_vl1l2_columns, n_val1l2_columns, + n_pct_columns, n_pstd_columns, n_pjc_columns, + n_prc_columns, n_eclv_columns, + n_mpr_columns, n_seeps_mpr_columns, n_seeps_columns }; // Text file abbreviations static const char * const txt_file_abbr[n_txt] = { - "fho", "ctc", "cts", - "mctc", "mcts", "cnt", - "sl1l2", "sal1l2", "vl1l2", - "val1l2", "pct", "pstd", - "pjc", "prc", "ecnt", - "orank", "rps", "eclv", - "mpr", "vcnt", "seeps_mpr", - "seeps" + "fho", "ctc", "cts", + "mctc", "mcts", + "cnt", "sl1l2", "sal1l2", + "vcnt", "vl1l2", "val1l2", + "pct", "pstd", "pjc", + "prc", "eclv", + "mpr", "seeps_mpr", "seeps" }; /////////////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/pair_stat/pair_stat_conf_info.cc b/src/tools/core/pair_stat/pair_stat_conf_info.cc index 957c526056..0d093dac66 100644 --- a/src/tools/core/pair_stat/pair_stat_conf_info.cc +++ b/src/tools/core/pair_stat/pair_stat_conf_info.cc @@ -610,26 +610,6 @@ int PairStatConfInfo::get_max_n_eclv_points() const { //////////////////////////////////////////////////////////////////////// -int PairStatConfInfo::get_max_n_hira_ens() const { - int n = 0; - - for(int i=0; i " + << "The length of \"" << conf_key_mpr_column << "\" and \"" + << conf_key_mpr_thresh << "\" must match (" << mpr_sa.n() + << " != " << mpr_ta.n() << ")!\n\n"; + exit(1); + } + + // Store in map + for(int i=0; mpr_sa.n(); i++) { + if(mpr_thr_inc_map.count(mpr_sa[i]) == 0) { + ThreshArray ta; + mpr_thr_inc_map[(mpr_sa[i])] = ta; + } + mpr_thr_inc_map[(mpr_sa[i])].add(mpr_ta[i]); + } // Dump the contents of the current thresholds if(mlog.verbosity_level() >= 5) { @@ -953,32 +939,14 @@ void PairStatVxOpt::process_config(PairsFormat ftype, // Conf: boot boot_info = parse_conf_boot(&odict); - // Conf: interp - interp_info = parse_conf_interp(&odict, conf_key_interp); - - // Conf: hira - hira_info = parse_conf_hira(&odict); - // Conf: hss_ec_value hss_ec_value = odict.lookup_double(conf_key_hss_ec_value); // Conf: rank_corr_flag rank_corr_flag = odict.lookup_bool(conf_key_rank_corr_flag); - // Conf: message_type - msg_typ = parse_conf_message_type(&odict); - - // Conf: duplicate_flag - duplicate_flag = parse_conf_duplicate_flag(&odict); - - // Conf: obs_summary - obs_summary = parse_conf_obs_summary(&odict); - - // Conf: obs_perc_value - obs_perc = parse_conf_percentile(&odict); - // Conf: desc - vx_pd.set_desc(parse_conf_string(&odict, conf_key_desc).c_str()); + vx_pd.set_desc(parse_conf_string(&odict, conf_key_desc, false).c_str()); // Conf: sid_inc vx_pd.set_sid_inc_filt(parse_conf_sid_list(&odict, conf_key_sid_inc)); @@ -998,23 +966,11 @@ void PairStatVxOpt::process_config(PairsFormat ftype, //////////////////////////////////////////////////////////////////////// void PairStatVxOpt::set_vx_pd(PairStatConfInfo *conf_info) { - int i, n; - int n_msg_typ = msg_typ.n(); - int n_mask = mask_name.n(); - int n_interp = interp_info.n_interp; + int n_mask = mask_name.n(); ConcatString cs; StringArray sa; - // Setup the VxPairDataPoint object with these dimensions: - // [n_msg_typ][n_mask][n_interp] - - // Check for at least one message type - if(n_msg_typ == 0) { - mlog << Error << "\nPairStatVxOpt::set_vx_pd() -> " - << "At least one output message type must be requested in \"" - << conf_key_message_type << "\".\n\n"; - exit(1); - } + // Setup the VxPairDataPoint object for each mask // Check for at least one masking region if(n_mask == 0) { @@ -1027,19 +983,11 @@ void PairStatVxOpt::set_vx_pd(PairStatConfInfo *conf_info) { exit(1); } - // Check for at least one interpolation method - if(n_interp == 0) { - mlog << Error << "\nPairStatVxOpt::set_vx_pd() -> " - << "At least one interpolation method must be requested in \"" - << conf_key_interp << "\".\n\n"; - exit(1); - } - - // Define the dimensions - vx_pd.set_size(n_msg_typ, n_mask, n_interp); + // Define the dimensions with n_msg_typ = n_interp = 1 + vx_pd.set_size(1, n_mask, 1); - // Store the MPR filter threshold - vx_pd.set_mpr_thresh(mpr_sa, mpr_ta); + // Store the MPR filtering thresholds + vx_pd.set_mpr_thr_inc_map(mpr_thr_inc_map); // Store the climo CDF info vx_pd.set_climo_cdf_info_ptr(&cdf_info); @@ -1074,59 +1022,43 @@ void PairStatVxOpt::set_vx_pd(PairStatConfInfo *conf_info) { vx_pd.set_msg_typ_wtr(sa); } - // Define the verifying message type name and values - for(i=0; imsg_typ_group_map[msg_typ[i]]; - if(sa.n() == 0) sa.add(msg_typ[i]); - vx_pd.set_msg_typ_vals(i, sa); - } - // Define the masking information: grid, poly, sid, point + int n; // Define the grid masks - for(i=0; imask_area_map[mask_name[n]])); } // Define the poly masks - for(i=0; imask_area_map[mask_name[n]])); } // Define the station ID masks - for(i=0; imask_sid_map[mask_name[n]])); } // Define the Lat/Lon point masks - for(i=0; i<(int) mask_llpnt.size(); i++) { + for(int i=0; i<(int) mask_llpnt.size(); i++) { n = i + mask_grid.n() + mask_poly.n() + mask_sid.n(); vx_pd.set_mask_llpnt(n, mask_name[n].c_str(), &mask_llpnt[i]); } - // Define the interpolation methods - for(i=0; iseeps_climo_name); vx_pd.set_seeps_thresh(conf_info->seeps_p1_thresh); } + return; } @@ -1185,7 +1117,7 @@ int PairStatVxOpt::n_txt_row(int i_txt_row) const { bool vect_flag = vx_pd.fcst_info->is_v_wind() && vx_pd.fcst_info->uv_index() >= 0; - int n_pd = get_n_msg_typ() * get_n_mask() * get_n_interp(); + int n_pd = get_n_mask(); // Determine row multiplier for climatology bins if(cdf_info.write_bins) { @@ -1242,14 +1174,6 @@ int PairStatVxOpt::n_txt_row(int i_txt_row) const { n = (prob_flag ? 0 : n_pd * get_n_cnt_thresh() * n_bin); break; - case i_vl1l2: - case i_val1l2: - // Number of VL1L2 or VAL1L2 lines = - // Message Types * Masks * Interpolations * Thresholds - n = (!vect_flag ? 0 : n_pd * - get_n_wind_thresh()); - break; - case i_vcnt: // Number of VCNT lines = // Message Types * Masks * Interpolations * Thresholds * @@ -1258,6 +1182,14 @@ int PairStatVxOpt::n_txt_row(int i_txt_row) const { get_n_wind_thresh() * get_n_ci_alpha()); break; + case i_vl1l2: + case i_val1l2: + // Number of VL1L2 or VAL1L2 lines = + // Message Types * Masks * Interpolations * Thresholds + n = (!vect_flag ? 0 : n_pd * + get_n_wind_thresh()); + break; + case i_pct: case i_pjc: case i_prc: @@ -1265,14 +1197,6 @@ int PairStatVxOpt::n_txt_row(int i_txt_row) const { // Message Types * Masks * Interpolations * Thresholds * // Climo Bins n = (!prob_flag ? 0 : n_pd * get_n_oprob_thresh() * n_bin); - - // Number of HiRA PCT, PJC, or PRC lines = - // Message Types * Masks * HiRA widths * Thresholds - if(hira_info.flag) { - n += (prob_flag ? 0 : n_pd * get_n_cat_thresh() * - hira_info.width.n()); - } - break; case i_pstd: @@ -1281,44 +1205,6 @@ int PairStatVxOpt::n_txt_row(int i_txt_row) const { // Alphas * Climo Bins n = (!prob_flag ? 0 : n_pd * get_n_oprob_thresh() * get_n_ci_alpha() * n_bin); - - // Number of HiRA PSTD lines = - // Message Types * Masks * HiRA widths * Thresholds * - // Alphas - if(hira_info.flag) { - n += (prob_flag ? 0 : n_pd * - get_n_cat_thresh() * hira_info.width.n() * - get_n_ci_alpha()); - } - - break; - - case i_ecnt: - case i_rps: - // Number of HiRA ECNT and RPS lines = - // Message Types * Masks * HiRA widths * - // Alphas - if(hira_info.flag) { - n = n_pd * hira_info.width.n() * get_n_ci_alpha(); - } - else { - n = 0; - } - - break; - - case i_orank: - // Number of HiRA ORANK lines possible = - // Number of pairs * Categorical Thresholds * - // HiRA widths - if(hira_info.flag) { - n = vx_pd.get_n_pair() * get_n_cat_thresh() * - hira_info.width.n(); - } - else { - n = 0; - } - break; case i_eclv: @@ -1333,34 +1219,21 @@ int PairStatVxOpt::n_txt_row(int i_txt_row) const { // Forecast Probability Thresholds * Climo Bins n += (!prob_flag ? 0 : n_pd * get_n_oprob_thresh() * get_n_fprob_thresh() * n_bin); - break; case i_mpr: // Compute the number of matched pairs to be written n = vx_pd.get_n_pair(); - - // Maximum number of HiRA MPR lines possible = - // Number of pairs * Max Scalar Categorical Thresholds * - // HiRA widths - if(hira_info.flag) { - n += (prob_flag ? 0 : - vx_pd.get_n_pair() * get_n_cat_thresh() * - hira_info.width.n()); - } - break; case i_seeps_mpr: // Compute the number of matched pairs to be written n = vx_pd.get_n_pair(); - break; case i_seeps: // Compute the number of matched pairs to be written n = vx_pd.get_n_pair(); - break; default: @@ -1408,16 +1281,3 @@ int PairStatVxOpt::get_n_oprob_thresh() const { } //////////////////////////////////////////////////////////////////////// - -int PairStatVxOpt::get_n_hira_ens() const { - int n = (hira_info.flag ? hira_info.width.max() : 0); - return n*n; -} - -//////////////////////////////////////////////////////////////////////// - -int PairStatVxOpt::get_n_hira_prob() const { - return hira_info.flag ? hira_info.cov_ta.n() : 0; -} - -//////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/pair_stat/pair_stat_conf_info.h b/src/tools/core/pair_stat/pair_stat_conf_info.h index 1ebea15d53..dd5dd72da4 100644 --- a/src/tools/core/pair_stat/pair_stat_conf_info.h +++ b/src/tools/core/pair_stat/pair_stat_conf_info.h @@ -37,24 +37,22 @@ static const int i_mcts = 4; static const int i_cnt = 5; static const int i_sl1l2 = 6; static const int i_sal1l2 = 7; -static const int i_vl1l2 = 8; -static const int i_val1l2 = 9; -static const int i_pct = 10; -static const int i_pstd = 11; -static const int i_pjc = 12; -static const int i_prc = 13; -static const int i_ecnt = 14; +static const int i_vcnt = 8; +static const int i_vl1l2 = 9; +static const int i_val1l2 = 10; -static const int i_orank = 15; -static const int i_rps = 16; -static const int i_eclv = 17; -static const int i_mpr = 18; -static const int i_vcnt = 19; -static const int i_seeps_mpr = 20; -static const int i_seeps = 21; +static const int i_pct = 11; +static const int i_pstd = 12; +static const int i_pjc = 13; +static const int i_prc = 14; +static const int i_eclv = 15; -static const int n_txt = 22; +static const int i_mpr = 16; +static const int i_seeps_mpr = 17; +static const int i_seeps = 18; + +static const int n_txt = 19; // Text file type static const STATLineType txt_file_type[n_txt] = { @@ -68,23 +66,20 @@ static const STATLineType txt_file_type[n_txt] = { STATLineType::cnt, // 5 STATLineType::sl1l2, // 6 STATLineType::sal1l2, // 7 - STATLineType::vl1l2, // 8 - STATLineType::val1l2, // 9 - - STATLineType::pct, // 10 - STATLineType::pstd, // 11 - STATLineType::pjc, // 12 - STATLineType::prc, // 13 - STATLineType::ecnt, // 14 - - STATLineType::orank, // 15 - STATLineType::rps, // 16 - STATLineType::eclv, // 17 - STATLineType::mpr, // 18 - STATLineType::vcnt, // 19 - - STATLineType::seeps_mpr, // 20 - STATLineType::seeps // 21 + + STATLineType::vcnt, // 8 + STATLineType::vl1l2, // 9 + STATLineType::val1l2, // 10 + + STATLineType::pct, // 11 + STATLineType::pstd, // 12 + STATLineType::pjc, // 13 + STATLineType::prc, // 14 + STATLineType::eclv, // 15 + + STATLineType::mpr, // 16 + STATLineType::seeps_mpr, // 17 + STATLineType::seeps // 18 }; /////////////////////////////////////////////////////////////////////////////// @@ -124,7 +119,7 @@ class PairStatVxOpt { ////////////////////////////////////////////////////////////////// - VxPairDataPoint vx_pd; // Matched pair data [n_msg_typ][n_mask][n_interp] + VxPairDataPoint vx_pd; // Matched pair data [n_mask] int beg_ds; // Begin observation time window offset int end_ds; // End observation time window offset @@ -147,8 +142,8 @@ class PairStatVxOpt { StringArray mask_poly; // Masking polyline strings StringArray mask_sid; // Masking station ID's - StringArray mpr_sa; // MPR column names - ThreshArray mpr_ta; // MPR column thresholds + // Matched pair inclusion thresholds + std::map mpr_thr_inc_map; // Vector of MaskLatLon objects defining Lat/Lon Point masks std::vector mask_llpnt; @@ -162,18 +157,12 @@ class PairStatVxOpt { NumArray ci_alpha; // Alpha value for confidence intervals BootInfo boot_info; // Bootstrapping information - InterpInfo interp_info; // Interpolation information - HiRAInfo hira_info; // HiRA verification logic double hss_ec_value; // HSS expected correct value bool rank_corr_flag; // Flag for computing rank correlations StringArray msg_typ; // Array of message types - DuplicateType duplicate_flag; // Duplicate observations - ObsSummary obs_summary; // Summarize observations - int obs_perc; // Summary percentile value - // Output file options STATOutputType output_flag[n_txt]; // Flag for each output line type @@ -190,9 +179,7 @@ class PairStatVxOpt { // Compute the number of output lines for this task int n_txt_row(int i) const; - int get_n_msg_typ() const; int get_n_mask() const; - int get_n_interp() const; int get_n_cnt_thresh() const; int get_n_cat_thresh() const; @@ -202,17 +189,13 @@ class PairStatVxOpt { int get_n_oprob_thresh() const; int get_n_eclv_points() const; - int get_n_hira_ens() const; - int get_n_hira_prob() const; int get_n_cdf_bin() const; int get_n_ci_alpha() const; }; //////////////////////////////////////////////////////////////////////// -inline int PairStatVxOpt::get_n_msg_typ() const { return msg_typ.n(); } inline int PairStatVxOpt::get_n_mask() const { return mask_name.n(); } -inline int PairStatVxOpt::get_n_interp() const { return interp_info.n_interp; } inline int PairStatVxOpt::get_n_eclv_points() const { return eclv_points.n(); } inline int PairStatVxOpt::get_n_cdf_bin() const { return cdf_info.n_bin; } @@ -298,8 +281,6 @@ class PairStatConfInfo { int get_max_n_fprob_thresh() const; int get_max_n_oprob_thresh() const; int get_max_n_eclv_points() const; - int get_max_n_hira_ens() const; - int get_max_n_hira_prob() const; // Check for any verification of vectors bool get_vflag() const; diff --git a/src/tools/core/point_stat/point_stat_conf_info.cc b/src/tools/core/point_stat/point_stat_conf_info.cc index cfd5e10432..c083d2fe4d 100644 --- a/src/tools/core/point_stat/point_stat_conf_info.cc +++ b/src/tools/core/point_stat/point_stat_conf_info.cc @@ -762,8 +762,7 @@ void PointStatVxOpt::clear() { mask_sid.clear(); mask_llpnt.clear(); - mpr_sa.clear(); - mpr_ta.clear(); + mpr_thr_inc_map.clear(); mask_name.clear(); @@ -919,8 +918,26 @@ void PointStatVxOpt::process_config(GrdFileType ftype, int_to_setlogic(odict.lookup_int(conf_key_wind_logic))); // Conf: mpr_column and mpr_thresh - mpr_sa = odict.lookup_string_array(conf_key_mpr_column); - mpr_ta = odict.lookup_thresh_array(conf_key_mpr_thresh); + StringArray mpr_sa(odict.lookup_string_array(conf_key_mpr_column)); + ThreshArray mpr_ta(odict.lookup_thresh_array(conf_key_mpr_thresh)); + + // Check for the same length + if(mpr_sa.n() != mpr_ta.n()) { + mlog << Error << "\nPointStatVxOpt::process_config() -> " + << "The length of \"" << conf_key_mpr_column << "\" and \"" + << conf_key_mpr_thresh << "\" must match (" << mpr_sa.n() + << " != " << mpr_ta.n() << ")!\n\n"; + exit(1); + } + + // Store in map + for(int i=0; mpr_sa.n(); i++) { + if(mpr_thr_inc_map.count(mpr_sa[i]) == 0) { + ThreshArray ta; + mpr_thr_inc_map[(mpr_sa[i])] = ta; + } + mpr_thr_inc_map[(mpr_sa[i])].add(mpr_ta[i]); + } // Dump the contents of the current thresholds if(mlog.verbosity_level() >= 5) { @@ -1089,7 +1106,7 @@ void PointStatVxOpt::set_vx_pd(PointStatConfInfo *conf_info) { vx_pd.set_size(n_msg_typ, n_mask, n_interp); // Store the MPR filter threshold - vx_pd.set_mpr_thresh(mpr_sa, mpr_ta); + vx_pd.set_mpr_thr_inc_map(mpr_thr_inc_map); // Store the climo CDF info vx_pd.set_climo_cdf_info_ptr(&cdf_info); diff --git a/src/tools/core/point_stat/point_stat_conf_info.h b/src/tools/core/point_stat/point_stat_conf_info.h index 9db3081dd7..3660cb78d8 100644 --- a/src/tools/core/point_stat/point_stat_conf_info.h +++ b/src/tools/core/point_stat/point_stat_conf_info.h @@ -129,8 +129,8 @@ class PointStatVxOpt { StringArray mask_poly; // Masking polyline strings StringArray mask_sid; // Masking station ID's - StringArray mpr_sa; // MPR column names - ThreshArray mpr_ta; // MPR column thresholds + // Matched pair inclusion thresholds + std::map mpr_thr_inc_map; // Vector of MaskLatLon objects defining Lat/Lon Point masks std::vector mask_llpnt;