From 11cb05d13cfa7f92ded04bdf4e48cf463c936295 Mon Sep 17 00:00:00 2001 From: johnhg Date: Fri, 19 Mar 2021 08:54:41 -0600 Subject: [PATCH] Update develop-ref after #1719 (#1727) * Start on write netcdf pickle alternative. * Write dataplane array. * Start on read of netcdf as pickle alternative. * Create attribute variables. * Use global attributes for met_info attrs. * Add grid structure. * Read metadata back into met_info.attrs. * Convert grid.nx and grid.ny to int. * Rename _name key to name. * Removed pickle write. * Fixed write_pickle_dataplane to work for both numpy and xarray. * Use items() to iterate of key, value attrs. * Write temporary text file. * Renamed scripts. * Changed script names in Makefile.am. * Replaced pickle with tmp_nc. * Fixed wrapper script names. * Test for attrs in met_in.met_data. * Initial version of read_tmp_point module. * Added read_tmp_point.py to install list. * Start on Python3_Script::read_tmp_point. * Write MPR tmp ascii file. * Renamed to read_tmp_ascii to use for point point and MPR. * Renamed to read_tmp_ascii to use for point point and MPR. * Define Python3_Script::import_read_tmp_ascii_py. * Call Python3_Script::import_read_tmp_ascii_py. * Append MET_BASE/wrappers to sys.path. * Finished implementation of Python3_Script::import_read_tmp_ascii_py. * Call Python3_Script::read_tmp_ascii in python_handler. * Revised python3_script::read_tmp_ascii with call to run, PyRun_String. * Return PyObject* from Python3_Script::run. * Restored call to run_python_string for now. * Per #1429, enhance error message from DataLine::get_item(). (#1682) * Feature 1429 tc_log second try (#1686) * Per #1429, enhance error message from DataLine::get_item(). * Per #1429, I realize that the line number actually is readily available in the DataLine class... so include it in the error message. * Feature 1588 ps_log (#1687) * Per #1588, updated pair_data_point.h/.cc to add detailed Debug(4) log messages, as specified in the GitHub issue. Do still need to test each of these cases to confirm that the log messages look good. * Per #1588, switch very detailed interpolation details from debug level 4 to 5. * Per #1588, remove the Debug(4) log message about duplicate obs since it's been moved up to a higher level. * Per #1588, add/update detailed log messages when processing point observations for bad data, off the grid, bad topo, big topo diffs, bad fcst value, and duplicate obs. * #1454 Disabled plot_data_plane_CESM_SSMI_microwave and plot_data_plane_CESM_sea_ice_nc becaues of not evenly spaced * #1454 Moved NC attribute name to nc_utils.h * #1454 Corrected sanity checking for lat/lon projection based on the percentage of the delta instead of fixed tolerance * #1454 Corrected sanity checking for lat/lon projection based on the percentage of the delta instead of fixed tolerance * #1454 Corrected data.delta_lon * #1454 Change bact to use diff instead of absolute value of diff * 454 Deleted instea dof commenting out * 454 Deleted instea dof commenting out * Feature 1684 bss and 1685 single reference model (#1689) * Per #1684, move an instance of the ClimoCDFInfo class into PairBase. Also define derive_climo_vals() and derive_climo_prob() utility functions. * Add to VxPairDataPoint and VxPairDataEnsemble functions to set the ClimoCDFInfo class. * Per #1684, update ensemble_stat and point_stat to set the ClimoCDFInfo object based on the contents of the config file. * Per #1684, update the vx_statistics library and stat_analysis to make calls to the new derive_climo_vals() and derive_climo_prob() functions. * Per #1684, since cdf_info is a member of PairBase class, need to handle it in the PairDataPoint and PairDataEnsemble assignment and subsetting logic. * Per #1684, during development, I ran across and then updated this log message. * Per #1684, working on log messages and figured that the regridding climo data should be moved from Debug(1) to at least Debug(2). * Per #1684 and #1685, update the logic for the derive_climo_vals() utility function. If only a single climo bin is requested, just return the climo mean. Otherwise, sample the requested number of values. * Per #1684, just fixing the format of this log message. * Per #1684, add a STATLine::get_offset() member function. * Per #1684, update parse_orank_line() logic. Rather than calling NumArray::clear() call NumArray::erase() to preserve allocated memory. Also, instead of parsing ensemble member values by column name, parse them by offset number. * Per #1684, call EnsemblePairData::extend() when parsing ORANK data to allocate one block of memory instead of bunches of litte ones. * Per #1684 and #1685, add another call to Ensemble-Stat to test computing the CRPSCL_EMP from a single climo mean instead of using the full climo distribution. * Per #1684 and #1685, update ensemble-stat docs about computing CRPSS_EMP relative to a single reference model. * Per #1684, need to update Grid-Stat to store the climo cdf info in the PairDataPoint objects. * Per #1684, remove debug print statements. * Per #1684, need to set cdf_info when aggregating MPR lines in Stat-Analysis. * Per #1684 and #1685, update PairDataEnsemble::compute_pair_vals() to print a log message indicating the climo data being used as reference: For a climo distribution defined by mean and stdev: DEBUG 3: Computing ensemble statistics relative to a 9-member climatological ensemble. For a single deterministic reference: DEBUG 3: Computing ensemble statistics relative to the climatological mean. * Per #1691, add met-10.0.0-beta4 release notes. (#1692) * Updated Python documentation * Per #1694, add VarInfo::magic_str_attr() to construct a field summary string from the name_attr() and level_attr() functions. * Per #1694, fixing 2 issues here. There was a bug in the computation of the max value. Had a less-than sign that should have been greater-than. Also, switch from tracking data by it's magic_str() to simply using VAR_i and VAR_j strings. We *could* have just used the i, j integers directly, but constructing the ij joint histogram integer could have been tricky since we start numbering with 0 instead of 1. i=0, j=1 would result in 01 which is the same as integer of 1. If we do want to switch to integers, we just need to make them 1-based and add +1 all over the place. * Per #1694, just switching to consistent variable name. * Just consistent spacing. * Added python3_script::import_read_tmp_ascii. * Restored read_tmp_ascii call. * Added lookup into ascii module. * Adding files for ReadTheDocs * Adding .yaml file for ReadTheDocs * Updated path to requirements.txt file * Updated path to conf.py file * Removing ReadTheDocs files and working in separate branch * Return PyObject* from read_tmp_ascii. * Put point_data in global namespace. * Remove temporary ascii file. * Added tmp_ascii_path. * Removed read_obs_from_pickle. * Trying different options for formats (#1702) * Per #1706, add bugfix to the develop branch. Also add a new job to unit_stat_analysis.xml to test out the aggregation of the ECNT line type. This will add new unit test output and cause the NB to fail. (#1708) * Feature 1471 python_grid (#1704) * Per #1471, defined a parse_grid_string() function in the vx_statistics library and then updated vx_data2d_python to call that function. However, this creates a circular dependency because vx_data2d_python now depends on vx_statistics. * Per #1471, because of the change in dependencies, I had to modify many, many Makefile.am files to link to the -lvx_statistics after -lvx_data2d_python. This is not great, but I didn't find a better solution. * Per #1471, add a sanity check to make sure the grid and data dimensions actually match. * Per #1471, add 3 new unit tests to demonstrate setting the python grid as a named grid, grid specification string, or a gridded data file. * Per #1471, document python grid changes in appendix F. * Per #1471, just spacing. * Per #1471, lots of Makefile.am changes to get this code to compile on kiowa. Worringly, it compiled and linked fine on my Mac laptop but not on kiowa. Must be some large differences in the linker logic. Co-authored-by: John Halley Gotway * Committing a fix for unit_python.xml directly to the develop branch. We referenced in a place where it's not defined. * Add *.dSYM to the .gitignore files in the src and internal_tests directories. * Replaced tmp netcdf _name attribute with name_str. * Append user script path to system path. * Revert "Feature 1319 no pickle" (#1717) * Fixed typos, added content, and modified release date format * #1715 Initial release * #1715 Do not combined if there are no overlapping beteewn TQZ and UV records * #1715 Added pb2nc_compute_pbl_cape * #1715 Added pb2nc_compute_pbl_cape * #1715 Reduced obs_bufr_var. Removed pb_report_type * #1715 Added a blank line for Error/Warning * Per #1725, return good status from TrackInfoArray::add() when using an ATCF line to create a new track. (#1726) * Per #1705, update the threshold node heirarchy by adding a climo_prob() function to determine the climatological probability of a CDP-type threshold. Also update derive_climo_prob() in pair_base.cc to call the new climo_prob() function. (#1724) * Bugfix 1716 develop perc_thresh (#1722) * Per #1716, committing changes from Randy Bullock to support floating point percentile thresholds. * Per #1716, no code changes, just consistent formatting. * Per #1716, change SFP50 example to SFP33.3 to show an example of using floating point percentile values. Co-authored-by: David Fillmore Co-authored-by: Howard Soh Co-authored-by: hsoh-u Co-authored-by: Julie.Prestopnik Co-authored-by: David Fillmore Co-authored-by: John Halley Gotway Co-authored-by: MET Tools Test Account --- met/docs/Users_Guide/config_options.rst | 2 +- met/docs/conf.py | 2 +- met/docs/index.rst | 85 +++++++++- met/src/basic/vx_config/my_config_scanner.cc | 78 ++++----- met/src/basic/vx_config/threshold.cc | 150 +++++++++++++++++ met/src/basic/vx_config/threshold.h | 29 ++-- met/src/libcode/vx_statistics/pair_base.cc | 39 +---- met/src/libcode/vx_tc_util/track_info.cc | 1 + met/src/tools/other/pb2nc/pb2nc.cc | 89 ++++++---- test/config/PB2NCConfig_pbl | 162 +++++++++++++++++++ test/xml/unit_pb2nc.xml | 19 +++ 11 files changed, 544 insertions(+), 112 deletions(-) create mode 100644 test/config/PB2NCConfig_pbl diff --git a/met/docs/Users_Guide/config_options.rst b/met/docs/Users_Guide/config_options.rst index 564e222e95..d8acbab286 100644 --- a/met/docs/Users_Guide/config_options.rst +++ b/met/docs/Users_Guide/config_options.rst @@ -81,7 +81,7 @@ The configuration file language supports the following data types: * The following percentile threshold types are supported: * "SFP" for a percentile of the sample forecast values. - e.g. ">SFP50" means greater than the 50-th forecast percentile. + e.g. ">SFP33.3" means greater than the 33.3-rd forecast percentile. * "SOP" for a percentile of the sample observation values. e.g. ">SOP75" means greater than the 75-th observation percentile. diff --git a/met/docs/conf.py b/met/docs/conf.py index efa6f948c9..13f65b1b9d 100644 --- a/met/docs/conf.py +++ b/met/docs/conf.py @@ -24,7 +24,7 @@ verinfo = version release = f'{version}' release_year = '2021' -release_date = f'{release_year}0302' +release_date = f'{release_year}-03-31' copyright = f'{release_year}, {author}' # -- General configuration --------------------------------------------------- diff --git a/met/docs/index.rst b/met/docs/index.rst index 4e00f3168d..a66a482aa1 100644 --- a/met/docs/index.rst +++ b/met/docs/index.rst @@ -1,16 +1,95 @@ ===================== MET version |version| ===================== -Developed by the `Developmental Testbed Center `_, Boulder, CO +Developed by the `Developmental Testbed Center `_, +Boulder, CO .. image:: _static/METplus_banner_photo_web.png History ------- -The Model Evaluation Tools (MET) were developed by the Developmental Testbed Center (DTC) and released in January 2008. The goal of the tools was to provide the community with a platform independent and extensible framework for reproducible verification. The DTC partners, including NCAR, NOAA, and the USAF, decided to start by replicating the NOAA EMC (see list of acronyms below) Mesoscale Branch verification package, called VSDB. In the first release, MET included several pre-processing, statistical, and analysis tools to provided the primary functionality as the EMC VSDB system, and also included a spatial verification package called MODE. +The Model Evaluation Tools (MET) were developed by the Developmental Testbed +Center (DTC) and released in January 2008. The goal of the tools was to +provide the community with a platform-independent and extensible framework +for reproducible verification. +The DTC partners, including NCAR, NOAA, and the USAF, decided to start by +replicating the NOAA EMC (see list of acronyms below) Mesoscale Branch +verification package, called VSDB. +In the first release, MET included several pre-processing, statistical, +and analysis tools to provide the same primary functionality as the EMC VSDB +system, and also included a spatial verification package called MODE. -Over the years, MET and VSDB packages grew in complexity. Verification capability at other NOAA laboratories, such as ESRL, were also under heavy development. An effort to unify verification capability was first started under the HIWPP project and led by NOAA ESRL. In 2015, the NGGPS Program Office started working groups to focus on several aspects of the next gen system, including the Verification and Validation Working Group. This group made the recommendation to use MET as the foundation for a unified verification capability. In 2016, NCAR and GSD leads visited EMC to gather requirements. At that time, the concept of METplus was developed as it extends beyond the original code base. It was originally called METplus but several constraints have driven the transition to the use of METplus. METplus is now the unified verification, validation, and diagnostics capability for NOAA's UFS and a component of NCAR's SIMA modeling frameworks. It being actively developed by NCAR, ESRL, EMC and is open to community contributions. +Over the years, MET and VSDB packages grew in complexity. Verification +capability at other NOAA laboratories, such as ESRL, were also under heavy +development. An effort to unify verification capability was first started +under the HIWPP project and led by NOAA ESRL. In 2015, the NGGPS +Program Office started working groups to focus on several aspects of the +next gen system, including the Verification and Validation Working Group. +This group made the recommendation to use MET as the foundation for a +unified verification capability. In 2016, NCAR and GSD leads visited EMC +to gather requirements. At that time, the concept of METplus was developed +as it extends beyond the original code base. It was originally MET+ but +several constraints have driven the transition to the use of METplus. +METplus is now the unified verification, validation, and +diagnostics capability for NOAA's UFS and a component of NCAR's SIMA +modeling frameworks. It is being actively developed by NCAR, ESRL, EMC +and is open to community contributions. +METplus Concept +--------------- +METplus is the overarching, or umbrella, repository and hence framework for +the Unified Forecast System verification capability. It is intended to be +extensible through adding additional capability developed by the community. +The core components of the framework include MET, the associated database and +display systems called METviewer and METexpress, and a suite of Python +wrappers to provide low-level automation and examples, also called use-cases. +A description of each tool along with some ancillary repositories are as +follows: + +* **MET** - core statistical tool that matches up grids with either gridded + analyses or point observations and applies configurable methods to compute + statistics and diagnostics +* **METviewer** - core database and display system intended for deep analysis + of MET output +* **METexpress** - core database and display system intended for quick + analysis via pre-defined queries of MET output +* **METplus wrappers** - suite of Python-based wrappers that provide + low-level automation of MET tools and newly developed plotting capability +* **METplus use-cases** - configuration files and sample data to show how to + invoke METplus wrappers to make using MET tools easier and reproducible +* **METcalcpy** - suite of Python-based scripts to be used by other + components of METplus tools for statistical aggregation, event + equalization, and other analysis needs +* **METplotpy** - suite of Python-based scripts to plot MET output, + and in come cases provide additional post-processing of output prior + to plotting +* **METdatadb** - database to store MET output and to be used by both + METviewer and METexpress + +The umbrella repository will be brought together by using a software package +called `manage_externals `_ +developed by the Community Earth System Modeling (CESM) team, hosted at NCAR +and NOAA Earth System's Research Laboratory. The manage_externals paackage +was developed because CESM is comprised of a number of different components +that are developed and managed independently. Each component also may have +additional "external" dependencies that need to be maintained independently. + +Acronyms +-------- + +* **MET** - Model Evaluation Tools +* **DTC** - Developmental Testbed Center +* **NCAR** - National Center for Atmospheric Research +* **NOAA** - National Oceanic and Atmospheric Administration +* **EMC** - Environmental Modeling Center +* **VSDB** - Verification Statistics Data Base +* **MODE** - Method for Object-Based Diagnostic Evaluation +* **UFS** - Unified Forecast System +* **SIMA** -System for Integrated Modeling of the Atmosphere +* **ESRL** - Earth Systems Research Laboratory +* **HIWPP** - High Impact Weather Predication Project +* **NGGPS** - Next Generation Global Predicatio System +* **GSD** - Global Systems Division .. toctree:: :hidden: diff --git a/met/src/basic/vx_config/my_config_scanner.cc b/met/src/basic/vx_config/my_config_scanner.cc index 57246913cf..1acae0582b 100644 --- a/met/src/basic/vx_config/my_config_scanner.cc +++ b/met/src/basic/vx_config/my_config_scanner.cc @@ -169,6 +169,8 @@ static bool replace_env(ConcatString &); static bool is_fort_thresh_no_spaces(); +static bool is_simple_perc_thresh(); + static int do_simple_perc_thresh(); @@ -370,6 +372,8 @@ if ( is_float_v2() ) { if ( do_float() ) return ( token(FLOAT) ); } if ( is_fort_thresh_no_spaces() ) { return ( do_fort_thresh() ); } +if ( is_simple_perc_thresh() ) { return ( do_simple_perc_thresh() ); } + int t; if ( is_id() ) { t = do_id(); return ( token(t) ); } @@ -533,7 +537,6 @@ if ( is_lhs ) { strncpy(configlval.text, configtext, max_id_length); return ( if ( strcmp(configtext, "print" ) == 0 ) { return ( PRINT ); } - // // boolean? // @@ -554,17 +557,13 @@ for (j=0; jlookup(configtext); if ( e && (e->is_number()) && (! is_lhs) ) { - // cout << "=================== id = \"" << configtext << "\" is_lhs = " << (is_lhs ? "true" : "false") << "\n"; - - // cout << "do_id() -> \n"; - // e->dump(cout); - if ( e->type() == IntegerType ) { set_int(configlval.nval, e->i_value()); @@ -613,28 +607,20 @@ if ( e && (! is_lhs) && (e->type() == UserFunctionType) ) { } - /////////////////////////////////////////////////////////////////////// - - - - // // fortran threshold without spaces? (example: "le150") // -if ( (strncmp(configtext, "lt", 2) == 0) && is_number(configtext + 2, max_id_length - 2) ) { return ( do_fort_thresh() ); } - for (j=0; j " @@ -1482,11 +1493,8 @@ if ( index < 0 ) { } - configlval.pc_info.perc_index = index; -// configlval.pc_info.is_simple = true; -configlval.pc_info.value = value; -// configlval.pc_info.value2 = bad_data_double;; +configlval.pc_info.value = value; return ( SIMPLE_PERC_THRESH ); @@ -1495,9 +1503,3 @@ return ( SIMPLE_PERC_THRESH ); //////////////////////////////////////////////////////////////////////// - - - - - - diff --git a/met/src/basic/vx_config/threshold.cc b/met/src/basic/vx_config/threshold.cc index 7879f7090a..b75ec9784a 100644 --- a/met/src/basic/vx_config/threshold.cc +++ b/met/src/basic/vx_config/threshold.cc @@ -166,6 +166,37 @@ return ( n ); //////////////////////////////////////////////////////////////////////// +double Or_Node::climo_prob() const + +{ + +if ( !left_child || !right_child ) { + + mlog << Error << "\nOr_Node::climo_prob() -> " + << "node not populated!\n\n"; + + exit ( 1 ); + +} + +double prob = bad_data_double; +double prob_left = left_child->climo_prob(); +double prob_right = right_child->climo_prob(); + +if ( !is_bad_data(prob_left) && !is_bad_data(prob_right) ) { + + prob = min(prob_left + prob_right, 1.0); + +} + +return ( prob ); + +} + + +//////////////////////////////////////////////////////////////////////// + + bool Or_Node::need_perc() const { @@ -356,6 +387,55 @@ return ( n ); //////////////////////////////////////////////////////////////////////// +double And_Node::climo_prob() const + +{ + +if ( !left_child || !right_child ) { + + mlog << Error << "\nAnd_Node::climo_prob() -> " + << "node not populated!\n\n"; + + exit ( 1 ); + +} + +double prob = bad_data_double; +double prob_left = left_child->climo_prob(); +double prob_right = right_child->climo_prob(); + + // + // For opposing inequalities, compute the difference in percentiles + // + +if ( !is_bad_data(prob_left) && !is_bad_data(prob_right) ) { + + // + // Support complex threshold types >a&&b + // + + if ( ( left_child->type() == thresh_gt || left_child->type() == thresh_ge ) && + ( right_child->type() == thresh_lt || right_child->type() == thresh_le ) ) { + + prob = max( 0.0, prob_right - ( 1.0 - prob_left ) ); + + } + else if ( ( left_child->type() == thresh_lt || left_child->type() == thresh_le ) && + ( right_child->type() == thresh_gt || right_child->type() == thresh_ge ) ) { + + prob = max( 0.0, prob_left - ( 1.0 - prob_right ) ); + + } +} + +return ( prob ); + +} + + +//////////////////////////////////////////////////////////////////////// + + bool And_Node::need_perc() const { @@ -540,6 +620,23 @@ return ( n ); //////////////////////////////////////////////////////////////////////// +double Not_Node::climo_prob() const + +{ + +double prob = bad_data_double; +double prob_child = child->climo_prob(); + +if ( !is_bad_data(prob_child) ) prob = 1.0 - prob_child; + +return ( prob ); + +} + + +//////////////////////////////////////////////////////////////////////// + + bool Not_Node::need_perc() const { @@ -1065,6 +1162,59 @@ return; //////////////////////////////////////////////////////////////////////// +double Simple_Node::climo_prob() const + +{ + +double prob = bad_data_double; + +if ( Ptype == perc_thresh_climo_dist ) { + + // Climo probability varies based on the threshold type + switch ( op ) { + + case thresh_lt: + case thresh_le: + + prob = PT/100.0; + break; + + case thresh_eq: + + prob = 0.0; + break; + + case thresh_ne: + + prob = 1.0; + break; + + case thresh_gt: + case thresh_ge: + + prob = 1.0 - PT/100.0; + break; + + default: + + mlog << Error << "\nSimple_Node::climo_prob() -> " + << "cannot convert climatological distribution percentile " + << "threshold to a probability!\n\n"; + + exit ( 1 ); + break; + + } // switch +} + +return ( prob ); + +} + + +//////////////////////////////////////////////////////////////////////// + + bool Simple_Node::need_perc() const { diff --git a/met/src/basic/vx_config/threshold.h b/met/src/basic/vx_config/threshold.h index ebca96a81c..493173e58d 100644 --- a/met/src/basic/vx_config/threshold.h +++ b/met/src/basic/vx_config/threshold.h @@ -157,6 +157,8 @@ class ThreshNode { virtual double pvalue() const = 0; + virtual double climo_prob() const = 0; + virtual bool need_perc() const = 0; virtual void set_perc(const NumArray *, const NumArray *, const NumArray *) = 0; @@ -197,6 +199,8 @@ class Or_Node : public ThreshNode { double pvalue() const; + double climo_prob() const; + bool need_perc() const; void set_perc(const NumArray *, const NumArray *, const NumArray *); @@ -217,10 +221,10 @@ class Or_Node : public ThreshNode { //////////////////////////////////////////////////////////////////////// -inline ThreshType Or_Node::type() const { return ( thresh_complex ); } -inline double Or_Node::value() const { return ( bad_data_double ); } -inline PercThreshType Or_Node::ptype() const { return ( no_perc_thresh_type ); } -inline double Or_Node::pvalue() const { return ( bad_data_double ); } +inline ThreshType Or_Node::type() const { return ( thresh_complex ); } +inline double Or_Node::value() const { return ( bad_data_double ); } +inline PercThreshType Or_Node::ptype() const { return ( no_perc_thresh_type ); } +inline double Or_Node::pvalue() const { return ( bad_data_double ); } //////////////////////////////////////////////////////////////////////// @@ -244,6 +248,8 @@ class And_Node : public ThreshNode { double pvalue() const; + double climo_prob() const; + bool need_perc() const; void set_perc(const NumArray *, const NumArray *, const NumArray *); @@ -293,6 +299,8 @@ class Not_Node : public ThreshNode { double pvalue() const; + double climo_prob() const; + bool need_perc() const; void set_perc(const NumArray *, const NumArray *, const NumArray *); @@ -363,6 +371,8 @@ class Simple_Node : public ThreshNode { double pvalue() const; + double climo_prob() const; + bool need_perc() const; void get_simple_nodes(vector &) const; @@ -435,6 +445,7 @@ class SingleThresh { double get_value() const; PercThreshType get_ptype() const; double get_pvalue() const; + double get_climo_prob() const; void get_simple_nodes(vector &) const; void multiply_by(const double); @@ -451,11 +462,11 @@ class SingleThresh { //////////////////////////////////////////////////////////////////////// -inline ThreshType SingleThresh::get_type() const { return ( node ? node->type() : thresh_na ); } -inline double SingleThresh::get_value() const { return ( node ? node->value() : bad_data_double ); } -inline PercThreshType SingleThresh::get_ptype() const { return ( node ? node->ptype() : no_perc_thresh_type ); } -inline double SingleThresh::get_pvalue() const { return ( node ? node->pvalue() : bad_data_double ); } - +inline ThreshType SingleThresh::get_type() const { return ( node ? node->type() : thresh_na ); } +inline double SingleThresh::get_value() const { return ( node ? node->value() : bad_data_double ); } +inline PercThreshType SingleThresh::get_ptype() const { return ( node ? node->ptype() : no_perc_thresh_type ); } +inline double SingleThresh::get_pvalue() const { return ( node ? node->pvalue() : bad_data_double ); } +inline double SingleThresh::get_climo_prob() const { return ( node ? node->climo_prob() : bad_data_double ); } //////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_statistics/pair_base.cc b/met/src/libcode/vx_statistics/pair_base.cc index 8066ed262f..0fe6a1b006 100644 --- a/met/src/libcode/vx_statistics/pair_base.cc +++ b/met/src/libcode/vx_statistics/pair_base.cc @@ -1064,46 +1064,21 @@ NumArray derive_climo_prob(const ClimoCDFInfo &cdf_info, const NumArray &mn_na, const NumArray &sd_na, const SingleThresh &othresh) { int i, n_mn, n_sd; - double prob; NumArray climo_prob, climo_vals; + double prob; // Number of valid climo mean and standard deviation n_mn = mn_na.n_valid(); n_sd = sd_na.n_valid(); - // For CDP threshold types, the climo probability is constant - if(othresh.get_ptype() == perc_thresh_climo_dist) { - - // Climo probability varies based on the threshold type - switch(othresh.get_type()) { - - case thresh_lt: - case thresh_le: - prob = othresh.get_pvalue()/100.0; - break; - - case thresh_eq: - prob = 0.0; - break; - - case thresh_ne: - prob = 1.0; - break; + // Check for constant climo probability + if(!is_bad_data(prob = othresh.get_climo_prob())) { - case thresh_gt: - case thresh_ge: - prob = 1.0 - othresh.get_pvalue()/100.0; - break; - - default: - mlog << Error << "\nderive_climo_prob() -> " - << "climatological threshold \"" << othresh.get_str() - << "\" cannot be converted to a probability!\n\n"; - exit(1); - break; - } + mlog << Debug(4) + << "For threshold " << othresh.get_str() + << ", using a constant climatological probability value of " + << prob << ".\n"; - // Add constant climo probability value climo_prob.add_const(prob, n_mn); } // If both mean and standard deviation were provided, use them to diff --git a/met/src/libcode/vx_tc_util/track_info.cc b/met/src/libcode/vx_tc_util/track_info.cc index b7c443c0f3..312d9aa620 100644 --- a/met/src/libcode/vx_tc_util/track_info.cc +++ b/met/src/libcode/vx_tc_util/track_info.cc @@ -801,6 +801,7 @@ bool TrackInfoArray::add(const ATCFTrackLine &l, bool check_dup, bool check_anly TrackInfo t; t.add(l, check_dup, check_anly); Track.push_back(t); + status = true; } return(status); diff --git a/met/src/tools/other/pb2nc/pb2nc.cc b/met/src/tools/other/pb2nc/pb2nc.cc index b2b6eba1fe..b36095ccf5 100644 --- a/met/src/tools/other/pb2nc/pb2nc.cc +++ b/met/src/tools/other/pb2nc/pb2nc.cc @@ -2428,7 +2428,7 @@ void write_netcdf_hdr_data() { // Check for no messages retained if(dim_count <= 0) { - mlog << Error << method_name << " -> " + mlog << Error << "\n" << method_name << " -> " << "No PrepBufr messages retained. Nothing to write.\n\n"; // Delete the NetCDF file remove_temp_file(ncfile); @@ -2920,16 +2920,27 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, float *pqtzuv_tq, *pqtzuv_uv; float *pqtzuv_merged = (float *) 0; float *next_pqtzuv, *prev_pqtzuv; + float tq_pres_max, tq_pres_min, uv_pres_max, uv_pres_min; std::map::iterator it, it_tq, it_uv; // Gets pressure levels for TQZ records - for (it=pqtzuv_map_tq.begin(); it!=pqtzuv_map_tq.end(); ++it) { - tq_levels.add(int(it->first)); + it = pqtzuv_map_tq.begin(); + tq_pres_min = tq_pres_max = it->first; + for (; it!=pqtzuv_map_tq.end(); ++it) { + float pres_v = it->first; + if (tq_pres_min > pres_v) tq_pres_min = pres_v; + if (tq_pres_max < pres_v) tq_pres_max = pres_v; + tq_levels.add(nint(pres_v)); } // Gets pressure levels for common records - for (it=pqtzuv_map_uv.begin(); it!=pqtzuv_map_uv.end(); ++it) { - if (tq_levels.has(int(it->first))) { - common_levels.add(int(it->first)); + it = pqtzuv_map_uv.begin(); + uv_pres_min = uv_pres_max = it->first; + for (; it!=pqtzuv_map_uv.end(); ++it) { + float pres_v = it->first; + if (uv_pres_min > pres_v) uv_pres_min = pres_v; + if (uv_pres_max < pres_v) uv_pres_max = pres_v; + if (tq_levels.has(nint(pres_v))) { + common_levels.add(nint(pres_v)); } } @@ -2937,22 +2948,36 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, log_tqz_and_uv(pqtzuv_map_tq, pqtzuv_map_uv, method_name); } + bool no_overlap = (tq_pres_max < uv_pres_min) || (tq_pres_min > uv_pres_max); + mlog << Debug(6) << method_name << "TQZ pressures: " << tq_pres_max + << " to " << tq_pres_min << " UV pressures: " << uv_pres_max + << " to " << uv_pres_min << (no_overlap ? " no overlap!" : " overlapping") << "\n"; + if( no_overlap ) { + mlog << Warning << "\n" << method_name + << "Can not combine TQ and UV records because of no overlapping." + << " TQZ count: " << tq_count << ", UV count: " << uv_count + << " common_levels: " << common_levels.n() << "\n\n"; + return pqtzuv_map_merged.size(); + } + // Select first record by 1) merging two records with the same pressure // level or 2) interpolate + int tq_pres, uv_pres; next_pqtzuv = (float *)0; it_tq = pqtzuv_map_tq.begin(); it_uv = pqtzuv_map_uv.begin(); pqtzuv_tq = (float *)it_tq->second; pqtzuv_uv = (float *)it_uv->second;; pqtzuv_merged = new float[mxr8vt]; - if (common_levels.has(int(it_tq->first)) - || common_levels.has(int(it_uv->first))) { + tq_pres = nint(it_tq->first); + uv_pres = nint(it_uv->first); + if (common_levels.has(tq_pres) || common_levels.has(uv_pres)) { // Found the records with the same precsure level - if (it_tq->first != it_uv->first) { - if (common_levels.has(int(it_uv->first))) { + if (tq_pres != uv_pres) { + if (common_levels.has(uv_pres)) { pqtzuv_uv = pqtzuv_map_uv[it_uv->first]; } - else if (common_levels.has(int(it_tq->first))) { + else if (common_levels.has(tq_pres)) { pqtzuv_tq = pqtzuv_map_tq[it_tq->first]; } } @@ -2968,7 +2993,7 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, prev_pqtzuv = (float *)it_uv->second; ++it_uv; } - next_pqtzuv = it_uv->second; + next_pqtzuv = (float *)it_uv->second; } else { //Interpolate TQZ into UV @@ -2978,7 +3003,7 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, prev_pqtzuv = (float *)it_tq->second; ++it_tq; } - next_pqtzuv = it_tq->second; + next_pqtzuv = (float *)it_tq->second; } interpolate_pqtzuv(prev_pqtzuv, pqtzuv_merged, next_pqtzuv); } @@ -2996,6 +3021,7 @@ int combine_tqz_and_uv(map pqtzuv_map_tq, if(mlog.verbosity_level() >= PBL_DEBUG_LEVEL) { log_merged_tqz_uv(pqtzuv_map_tq, pqtzuv_map_uv, pqtzuv_map_merged, method_name); } + delete [] pqtzuv_merged; } return pqtzuv_map_merged.size(); @@ -3034,7 +3060,7 @@ float compute_pbl(map pqtzuv_map_tq, hgt_cnt = spfh_cnt = 0; for (it=pqtzuv_map_merged.begin(); it!=pqtzuv_map_merged.end(); ++it) { if (index < 0) { - mlog << Error << method_name << "negative index: " << index << "\n"; + mlog << Error << "\n" << method_name << "negative index: " << index << "\n\n"; break; } @@ -3048,13 +3074,13 @@ float compute_pbl(map pqtzuv_map_tq, pbl_data_vgrd[index] = pqtzuv[5]; if (!is_eq(pbl_data_spfh[index], bad_data_float)) spfh_cnt++; if (!is_eq(pbl_data_hgt[index], bad_data_float)) hgt_cnt++; - selected_levels.add(int(it->first)); + selected_levels.add(nint(it->first)); } index--; } if (index != -1) { - mlog << Error << method_name << "Missing some levels (" << index << ")\n"; + mlog << Error << "\n" << method_name << "Missing some levels (" << index << ")\n"; } if (pbl_level > MAX_PBL_LEVEL) { @@ -3070,7 +3096,7 @@ float compute_pbl(map pqtzuv_map_tq, if (!is_eq(highest_pressure, bad_data_float)) { index = MAX_PBL_LEVEL - 1; for (; it!=pqtzuv_map_tq.end(); ++it) { - int pres_level = int(it->first); + int pres_level = nint(it->first); if (selected_levels.has(pres_level)) break; float *pqtzuv = pqtzuv_map_merged[it->first]; @@ -3139,10 +3165,10 @@ void insert_pbl(float *obs_arr, const float pbl_value, const int pbl_code, hdr_info << unix_to_yyyymmdd_hhmmss(hdr_vld_ut) << " " << hdr_typ << " " << hdr_sid; if (is_eq(pbl_value, bad_data_float)) { - mlog << Warning << "Failed to compute PBL " << hdr_info << "\n\n"; + mlog << Warning << "\nFailed to compute PBL " << hdr_info << "\n\n"; } else if (pbl_value < hdr_elv) { - mlog << Warning << "Not saved because the computed PBL (" << pbl_value + mlog << Warning << "\nNot saved because the computed PBL (" << pbl_value << ") is less than the station elevation (" << hdr_elv << "). " << hdr_info << "\n\n"; obs_arr[4] = 0; @@ -3156,7 +3182,7 @@ void insert_pbl(float *obs_arr, const float pbl_value, const int pbl_code, << " lat: " << hdr_lat << ", lon: " << hdr_lon << ", elv: " << hdr_elv << " " << hdr_info << "\n\n"; if (obs_arr[4] > MAX_PBL) { - mlog << Warning << " Computed PBL (" << obs_arr[4] << " from " + mlog << Warning << "\nComputed PBL (" << obs_arr[4] << " from " << pbl_value << ") is too high, Reset to " << MAX_PBL << " " << hdr_info<< "\n\n"; obs_arr[4] = MAX_PBL; @@ -3192,9 +3218,14 @@ int interpolate_by_pressure(int length, float *pres_data, float *var_data) { << var_data[idx_start] << " and " << var_data[idx_end] << "\n"; float data_diff = var_data[idx_end] - var_data[idx_start]; for (idx2 = idx_start+1; idx2 pqtzuv_map_pivot, if (first_pres < it_pivot->first) break; } mlog << Debug(8) << method_name << "pivot->first: " << it_pivot->first - << " aux->first " << it_aux->first << " first_pres: " << first_pres - << " prev_pqtzuv[0]" << prev_pqtzuv[0] << "\n"; + << " aux->first: " << it_aux->first << " first_pres: " << first_pres + << " prev_pqtzuv[0]: " << prev_pqtzuv[0] << "\n"; // Find next UV level for (; it_aux!=pqtzuv_map_aux.end(); ++it_aux) { // Skip the records below the first mathcing/interpolated level diff --git a/test/config/PB2NCConfig_pbl b/test/config/PB2NCConfig_pbl new file mode 100644 index 0000000000..eeedd7a3e4 --- /dev/null +++ b/test/config/PB2NCConfig_pbl @@ -0,0 +1,162 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// PB2NC configuration file. +// +// For additional information, see the MET_BASE/config/README file. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// PrepBufr message type +// +message_type = ["ONLYSF", "ADPUPA"]; + +// +// Mapping of message type group name to comma-separated list of values +// Derive PRMSL only for SURFACE message types +// +message_type_group_map = [ + { key = "SURFACE"; val = "ADPSFC,SFCSHP,MSONET"; }, + { key = "ANYAIR"; val = "AIRCAR,AIRCFT"; }, + { key = "ANYSFC"; val = "ADPSFC,SFCSHP,ADPUPA,PROFLR,MSONET"; }, + { key = "ONLYSF"; val = "ADPSFC,SFCSHP"; } +]; + +// +// Mapping of input PrepBufr message types to output message types +// +message_type_map = []; + +// +// PrepBufr station ID +// +station_id = []; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Observation time window +// +obs_window = { + beg = -2700; + end = 2700; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Observation retention regions +// +mask = { + grid = ""; + poly = ""; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Observing location elevation +// +elevation_range = { + beg = -1000; + end = 100000; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Observation types +// +pb_report_type = []; + +in_report_type = []; + +instrument_type = []; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Vertical levels to retain +// +level_range = { + beg = 1; + end = 511; +} + +level_category = [0, 1, 4, 5, 6]; + +/////////////////////////////////////////////////////////////////////////////// + +// +// BUFR variable names to retain or derive. +// Use obs_bufr_map to rename variables in the output. +// If empty, process all available variables. +// +obs_bufr_var = ["D_CAPE", "D_PBL"]; +//////////////////////////////////////////////////////////////////////////////// + +// +// Mapping of input BUFR variable names to output variables names. +// The default PREPBUFR map, obs_prepbufr_map, is appended to this map. +// +obs_bufr_map = []; + +// +// Default mapping for PREPBUFR. Replace input BUFR variable names with GRIB +// abbreviations in the output. This default map is appended to obs_bufr_map. +// This should not typically be overridden. +// +obs_prefbufr_map = [ + { key = "POB"; val = "PRES"; }, + { key = "QOB"; val = "SPFH"; }, + { key = "TOB"; val = "TMP"; }, + { key = "UOB"; val = "UGRD"; }, + { key = "VOB"; val = "VGRD"; }, + { key = "D_DPT"; val = "DPT"; }, + { key = "D_WDIR"; val = "WDIR"; }, + { key = "D_WIND"; val = "WIND"; }, + { key = "D_RH"; val = "RH"; }, + { key = "D_MIXR"; val = "MIXR"; }, + { key = "D_PBL"; val = "HPBL"; }, + { key = "D_PRMSL"; val = "PRMSL"; }, + { key = "D_CAPE"; val = "CAPE"; }, + { key = "TDO"; val = "DPT"; }, + { key = "PMO"; val = "PRMSL"; }, + { key = "TOCC"; val = "TCDC"; }, + { key = "HOVI"; val = "VIS"; }, + { key = "CEILING"; val = "HGT"; }, + { key = "MXGS"; val = "GUST"; } +]; + +//////////////////////////////////////////////////////////////////////////////// + +quality_mark_thresh = 9; +event_stack_flag = TOP; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Time periods for the summarization +// obs_var (string array) is added and works like grib_code (int array) +// when use_var_id is enabled and variable names are saved. +// +time_summary = { + flag = FALSE; + raw_data = FALSE; + beg = "000000"; + end = "235959"; + step = 300; + width = 600; + grib_code = []; + obs_var = [ "TMP", "WDIR", "RH" ]; + type = [ "min", "max", "range", "mean", "stdev", "median", "p80" ]; + vld_freq = 0; + vld_thresh = 0.0; +} + +//////////////////////////////////////////////////////////////////////////////// + +tmp_dir = "/tmp"; +version = "V10.0"; + +//////////////////////////////////////////////////////////////////////////////// diff --git a/test/xml/unit_pb2nc.xml b/test/xml/unit_pb2nc.xml index a65f52110d..0366d902f9 100644 --- a/test/xml/unit_pb2nc.xml +++ b/test/xml/unit_pb2nc.xml @@ -131,6 +131,25 @@ + + &MET_BIN;/pb2nc + + STATION_ID + MASK_GRID + MASK_POLY + QUALITY_MARK_THRESH 2 + + \ + &DATA_DIR_OBS;/prepbufr/nam.20210311.t00z.prepbufr.tm00 \ + &OUTPUT_DIR;/pb2nc/nam.20210311.t00z.prepbufr.tm00.pbl.nc \ + &CONFIG_DIR;/PB2NCConfig_pbl \ + -v 1 + + + &OUTPUT_DIR;/pb2nc/nam.20210311.t00z.prepbufr.tm00.pbl.nc + + + &MET_BIN;/pb2nc