From b70b5e70f3e6d54124d24c37bac474481be1d3aa Mon Sep 17 00:00:00 2001 From: johnhg Date: Sun, 30 Oct 2022 10:38:09 -0600 Subject: [PATCH] Update develop-ref after #2307, #2314, #2315, #2316, and #2317 (#2324) Co-authored-by: Randy Bullock Co-authored-by: davidfillmore Co-authored-by: rgbullock Co-authored-by: John Halley Gotway Co-authored-by: Howard Soh Co-authored-by: johnhg Co-authored-by: John Halley Gotway Co-authored-by: Seth Linden Co-authored-by: George McCabe <23407799+georgemccabe@users.noreply.github.com> Co-authored-by: Julie Prestopnik Co-authored-by: Seth Linden Co-authored-by: j-opatz <59586397+j-opatz@users.noreply.github.com> Co-authored-by: Howard Soh Co-authored-by: jprestop Co-authored-by: Seth Linden Co-authored-by: hsoh-u Co-authored-by: John Halley Gotway Co-authored-by: MET Tools Test Account Co-authored-by: mo-mglover <78152252+mo-mglover@users.noreply.github.com> Co-authored-by: davidalbo Co-authored-by: lisagoodrich <33230218+lisagoodrich@users.noreply.github.com> Co-authored-by: Dan Adriaansen Co-authored-by: Dave Albo Co-authored-by: Lisa Goodrich Co-authored-by: Molly Smith Co-authored-by: Jonathan Vigh Co-authored-by: bikegeek <3753118+bikegeek@users.noreply.github.com> --- .github/jobs/set_job_controls.sh | 2 +- .../build_docker_and_trigger_metplus.yml | 2 +- data/config/IODA2NCConfig_default | 11 +- data/config/TCPairsConfig_default | 41 ++ data/config/TCStatConfig_default | 16 +- data/table_files/met_header_columns_V11.0.txt | 3 +- docs/Users_Guide/appendixC.rst | 14 + docs/Users_Guide/ensemble-stat.rst | 15 + docs/Users_Guide/reformat_point.rst | 13 +- docs/Users_Guide/tc-pairs.rst | 126 +++- docs/Users_Guide/tc-stat.rst | 28 +- internal/scripts/docker/Dockerfile | 2 +- internal/scripts/docker/Dockerfile.copy | 2 +- .../scripts/environment/development.docker | 2 +- .../config/GridDiagConfig_APCP_06_FCST_OBS | 4 +- .../test_unit/config/TCPairsConfig_ALAL2010 | 5 + .../test_unit/config/TCPairsConfig_BASIN_MAP | 5 + .../test_unit/config/TCPairsConfig_CONSENSUS | 5 + .../config/TCPairsConfig_DIAGNOSTICS | 157 +++++ .../test_unit/config/TCPairsConfig_INTERP12 | 6 +- .../test_unit/config/TCPairsConfig_PROBRIRW | 5 + .../test_unit/config/TCStatConfig_ALAL2010 | 14 + .../test_unit/config/TCStatConfig_DIAGNOSTICS | 226 ++++++++ .../test_unit/config/TCStatConfig_PROBRIRW | 14 + internal/test_unit/hdr/met_11_0.hdr | 3 +- internal/test_unit/xml/unit_grid_diag.xml | 40 ++ internal/test_unit/xml/unit_ioda2nc.xml | 38 ++ internal/test_unit/xml/unit_tc_pairs.xml | 23 + internal/test_unit/xml/unit_tc_stat.xml | 13 + src/basic/vx_config/config.tab.cc | 2 +- src/basic/vx_config/config_constants.h | 22 + src/basic/vx_config/config_util.cc | 123 +++- src/basic/vx_config/config_util.h | 12 +- src/basic/vx_util/data_line.cc | 28 + src/basic/vx_util/data_line.h | 2 + src/basic/vx_util/stat_column_defs.h | 14 +- .../vx_data2d_factory/is_netcdf_file.cc | 7 +- src/libcode/vx_nc_obs/nc_obs_util.cc | 35 +- src/libcode/vx_nc_util/nc_utils.cc | 340 +++++++++-- src/libcode/vx_nc_util/nc_utils.h | 38 +- src/libcode/vx_nc_util/nc_utils.hpp | 29 +- src/libcode/vx_stat_out/stat_columns.cc | 29 +- src/libcode/vx_statistics/compute_stats.cc | 15 + src/libcode/vx_statistics/ens_stats.cc | 20 +- src/libcode/vx_statistics/ens_stats.h | 5 + .../vx_statistics/pair_data_ensemble.cc | 103 +++- .../vx_statistics/pair_data_ensemble.h | 10 + src/libcode/vx_tc_util/Makefile.am | 1 + src/libcode/vx_tc_util/Makefile.in | 20 + src/libcode/vx_tc_util/diag_file.cc | 542 ++++++++++++++++++ src/libcode/vx_tc_util/diag_file.h | 147 +++++ src/libcode/vx_tc_util/tc_columns.cc | 102 +++- src/libcode/vx_tc_util/tc_columns.h | 19 +- src/libcode/vx_tc_util/tc_stat_line.cc | 2 + src/libcode/vx_tc_util/tc_stat_line.h | 2 + src/libcode/vx_tc_util/track_info.cc | 130 ++++- src/libcode/vx_tc_util/track_info.h | 50 +- src/libcode/vx_tc_util/track_pair_info.cc | 200 ++++--- src/libcode/vx_tc_util/track_pair_info.h | 23 +- src/libcode/vx_tc_util/track_point.cc | 87 +-- src/libcode/vx_tc_util/track_point.h | 70 +-- src/libcode/vx_tc_util/vx_tc_util.h | 1 + .../core/stat_analysis/aggr_stat_line.cc | 39 +- .../core/stat_analysis/parse_stat_line.cc | 6 + .../core/stat_analysis/parse_stat_line.h | 3 + src/tools/other/ascii2nc/airnow_locations.cc | 12 +- src/tools/other/grid_diag/grid_diag.cc | 133 +++-- src/tools/other/grid_diag/grid_diag.h | 6 + src/tools/other/ioda2nc/ioda2nc.cc | 380 ++++++++---- src/tools/other/ioda2nc/ioda2nc_conf_info.cc | 1 + src/tools/other/ioda2nc/ioda2nc_conf_info.h | 1 + src/tools/tc_utils/tc_gen/Makefile.am | 2 +- src/tools/tc_utils/tc_gen/Makefile.in | 2 +- src/tools/tc_utils/tc_pairs/tc_pairs.cc | 288 +++++++--- src/tools/tc_utils/tc_pairs/tc_pairs.h | 9 +- .../tc_utils/tc_pairs/tc_pairs_conf_info.cc | 83 +++ .../tc_utils/tc_pairs/tc_pairs_conf_info.h | 6 + src/tools/tc_utils/tc_stat/Makefile.am | 2 +- src/tools/tc_utils/tc_stat/Makefile.in | 2 +- src/tools/tc_utils/tc_stat/out.tcst | 0 src/tools/tc_utils/tc_stat/tc_stat.cc | 13 +- .../tc_utils/tc_stat/tc_stat_conf_info.cc | 10 + src/tools/tc_utils/tc_stat/tc_stat_files.cc | 34 +- src/tools/tc_utils/tc_stat/tc_stat_job.cc | 276 +++++++-- src/tools/tc_utils/tc_stat/tc_stat_job.h | 41 +- 85 files changed, 3717 insertions(+), 697 deletions(-) create mode 100644 internal/test_unit/config/TCPairsConfig_DIAGNOSTICS create mode 100644 internal/test_unit/config/TCStatConfig_DIAGNOSTICS create mode 100644 src/libcode/vx_tc_util/diag_file.cc create mode 100644 src/libcode/vx_tc_util/diag_file.h create mode 100644 src/tools/tc_utils/tc_stat/out.tcst diff --git a/.github/jobs/set_job_controls.sh b/.github/jobs/set_job_controls.sh index b26f6d577d..e897aaa912 100755 --- a/.github/jobs/set_job_controls.sh +++ b/.github/jobs/set_job_controls.sh @@ -6,7 +6,7 @@ run_unit_tests=false run_diff=false run_update_truth=false met_base_repo=met-base -met_base_tag=v1.0 +met_base_tag=v1.1 input_data_version=develop truth_data_version=develop diff --git a/.github/workflows/build_docker_and_trigger_metplus.yml b/.github/workflows/build_docker_and_trigger_metplus.yml index e6c8c86a3f..15cd652caf 100644 --- a/.github/workflows/build_docker_and_trigger_metplus.yml +++ b/.github/workflows/build_docker_and_trigger_metplus.yml @@ -29,7 +29,7 @@ jobs: env: SOURCE_BRANCH: ${{ steps.get_branch_name.outputs.branch_name }}-lite MET_BASE_REPO: met-base - MET_BASE_TAG: v1.0 + MET_BASE_TAG: v1.1 - name: Push Docker Image run: .github/jobs/push_docker_image.sh diff --git a/data/config/IODA2NCConfig_default b/data/config/IODA2NCConfig_default index 2ac36c27ce..124854005e 100644 --- a/data/config/IODA2NCConfig_default +++ b/data/config/IODA2NCConfig_default @@ -92,7 +92,16 @@ metadata_map = [ { key = "station_id"; val = "station_id,report_identifier"; }, { key = "pressure"; val = "air_pressure,pressure"; }, { key = "height"; val = "height,height_above_mean_sea_level"; }, - { key = "elevation"; val = ""; } + { key = "datetime"; val = "datetime,dateTime"; }, + { key = "elevation"; val = "elevation,station_elevation"; } +]; + +// +// Default mapping for obs to qc. +// +obs_to_qc_map = [ + { key = "wind_from_direction"; val = "eastward_wind,northward_wind"; }, + { key = "wind_speed"; val = "eastward_wind,northward_wind"; } ]; missing_thresh = [ <=-1e9, >=1e9, ==-9999 ]; diff --git a/data/config/TCPairsConfig_default b/data/config/TCPairsConfig_default index ed3c659d5f..a768ac23e6 100644 --- a/data/config/TCPairsConfig_default +++ b/data/config/TCPairsConfig_default @@ -137,6 +137,47 @@ watch_warn = { time_offset = -14400; } +// +// Diagnostics to be extracted +// +diag_name = []; + +// +// Unit conversions to be applied based on diagnostic names and units +// +diag_convert_map = [ + { source = "TCDIAG"; + key = [ "(10C)", "(10KT)", "(10M/S)" ]; + convert(x) = x / 10; }, + + { source = "LSDIAG_RT"; + key = [ "LAT", "LON", "CSST", "RSST", "DSST", "DSTA", "XDST", "XNST", "NSST", "NSTA", + "NTMX", "NTFR", "U200", "U20C", "V20C", "E000", "EPOS", "ENEG", "EPSS", "ENSS", + "T000", "TLAT", "TLON", "TWAC", "TWXC", "G150", "G200", "G250", "V000", "V850", + "V500", "V300", "SHDC", "SHGC", "T150", "T200", "T250", "SHRD", "SHRS", "SHRG", + "HE07", "HE05", "PW01", "PW02", "PW03", "PW04", "PW05", "PW06", "PW07", "PW08", + "PW09", "PW10", "PW11", "PW12", "PW13", "PW14", "PW15", "PW16", "PW17", "PW18", + "PW20", "PW21" ]; + convert(x) = x / 10; }, + + { source = "LSDIAG_RT"; + key = [ "VVAV", "VMFX", "VVAC" ]; + convert(x) = x / 100; }, + + { source = "LSDIAG_RT"; + key = [ "TADV" ]; + convert(x) = x / 1000000; }, + + { source = "LSDIAG_RT"; + key = [ "Z850", "D200", "TGRD", "DIVC" ]; + convert(x) = x / 10000000; }, + + { source = "LSDIAG_RT"; + key = [ "PENC", "PENV" ]; + convert(x) = x / 10 + 1000; } + +]; + // // Modify basin names to make them consistent across ATCF input files. // diff --git a/data/config/TCStatConfig_default b/data/config/TCStatConfig_default index 5063ad8feb..9925185a7c 100644 --- a/data/config/TCStatConfig_default +++ b/data/config/TCStatConfig_default @@ -93,7 +93,7 @@ line_type = []; // // Stratify by checking the watch/warning status for each track point -// common to both the ADECK and BDECK tracks. If the watch/warning status +// common to both the ADECK and BDECK tracks. If the watch/warning status // of any of the track points appears in the list, retain the entire track. // track_watch_warn = []; @@ -134,6 +134,20 @@ init_str_val = []; init_str_exc_name = []; init_str_exc_val = []; +// +// Stratify track points by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines. +// +diag_thresh_name = []; +diag_thresh_val = []; + +// +// Stratify entire tracks by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines for lead time 0. +// +init_diag_thresh_name = []; +init_diag_thresh_val = []; + // // Stratify by the ADECK and BDECK distances to land. // diff --git a/data/table_files/met_header_columns_V11.0.txt b/data/table_files/met_header_columns_V11.0.txt index 59a97d40a6..fb48b0d16f 100644 --- a/data/table_files/met_header_columns_V11.0.txt +++ b/data/table_files/met_header_columns_V11.0.txt @@ -19,7 +19,7 @@ V11.0 : STAT : PJC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID V11.0 : STAT : PRC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_THRESH) THRESH_[0-9]* PODY_[0-9]* POFD_[0-9]* V11.0 : STAT : PSTD : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_THRESH) BASER BASER_NCL BASER_NCU RELIABILITY RESOLUTION UNCERTAINTY ROC_AUC BRIER BRIER_NCL BRIER_NCU BRIERCL BRIERCL_NCL BRIERCL_NCU BSS BSS_SMPL THRESH_[0-9]* V11.0 : STAT : ECLV : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BASER VALUE_BASER (N_PTS) CL_[0-9]* VALUE_[0-9]* -V11.0 : STAT : ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP CRPS_EMP_FAIR +V11.0 : STAT : ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP CRPS_EMP_FAIR BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS V11.0 : STAT : RPS : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_PROB RPS_REL RPS_RES RPS_UNC RPS RPSS RPSS_SMPL RPS_COMP V11.0 : STAT : RHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_RANK) RANK_[0-9]* V11.0 : STAT : PHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BIN_SIZE (N_BIN) BIN_[0-9]* @@ -37,4 +37,5 @@ V11.0 : MODE : OBJ : VERSION MODEL N_VALID GRID_RES DESC FCST_LEAD FCST_VAL V11.0 : MODE : CTS : VERSION MODEL N_VALID GRID_RES DESC FCST_LEAD FCST_VALID FCST_ACCUM OBS_LEAD OBS_VALID OBS_ACCUM FCST_RAD FCST_THR OBS_RAD OBS_THR FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE FIELD TOTAL FY_OY FY_ON FN_OY FN_ON BASER FMEAN ACC FBIAS PODY PODN POFD FAR CSI GSS HK HSS ODDS V11.0 : TCST : TCMPR : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE TOTAL INDEX LEVEL WATCH_WARN INITIALS ALAT ALON BLAT BLON TK_ERR X_ERR Y_ERR ALTK_ERR CRTK_ERR ADLAND BDLAND AMSLP BMSLP AMAX_WIND BMAX_WIND AAL_WIND_34 BAL_WIND_34 ANE_WIND_34 BNE_WIND_34 ASE_WIND_34 BSE_WIND_34 ASW_WIND_34 BSW_WIND_34 ANW_WIND_34 BNW_WIND_34 AAL_WIND_50 BAL_WIND_50 ANE_WIND_50 BNE_WIND_50 ASE_WIND_50 BSE_WIND_50 ASW_WIND_50 BSW_WIND_50 ANW_WIND_50 BNW_WIND_50 AAL_WIND_64 BAL_WIND_64 ANE_WIND_64 BNE_WIND_64 ASE_WIND_64 BSE_WIND_64 ASW_WIND_64 BSW_WIND_64 ANW_WIND_64 BNW_WIND_64 ARADP BRADP ARRP BRRP AMRD BMRD AGUSTS BGUSTS AEYE BEYE ADIR BDIR ASPEED BSPEED ADEPTH BDEPTH NUM_MEMBERS TRACK_SPREAD DIST_MEAN MSLP_SPREAD MAX_WIND_SPREAD +V11.0 : TCST : TCDIAG : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE TOTAL INDEX SOURCE (N_DIAG) DIAG_[0-9]* VALUE_[0-9]* V11.0 : TCST : PROBRIRW : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE ALAT ALON BLAT BLON INITIALS TK_ERR X_ERR Y_ERR ADLAND BDLAND RIRW_BEG RIRW_END RIRW_WINDOW AWIND_END BWIND_BEG BWIND_END BDELTA BDELTA_MAX BLEVEL_BEG BLEVEL_END (N_THRESH) THRESH_[0-9]* PROB_[0-9]* diff --git a/docs/Users_Guide/appendixC.rst b/docs/Users_Guide/appendixC.rst index 809b3bcb8b..fb6eb4b873 100644 --- a/docs/Users_Guide/appendixC.rst +++ b/docs/Users_Guide/appendixC.rst @@ -470,6 +470,7 @@ Mean Error (ME) --------------- Called "ME" in CNT output :numref:`table_PS_format_info_CNT` +Called "ME_OERR", "ME_GE_OBS", and "ME_LT_OBS" in ECNT output :numref:`table_ES_header_info_es_out_ECNT` The Mean Error, ME, is a measure of overall bias for continuous variables; in particular ME = Bias. It is defined as @@ -1003,6 +1004,19 @@ The continuous ranked probability skill score (CRPSS) is similar to the MSESS an For the normal distribution fit (CRPSS), the reference CRPS is computed using the climatological mean and standard deviation. For the empirical distribution (CRPSS_EMP), the reference CRPS is computed by sampling from the assumed normal climatological distribution defined by the mean and standard deviation. +Bias Ratio +---------- + +Called "BIAS_RATIO" in ECNT output :numref:`table_ES_header_info_es_out_ECNT` + +The bias ratio (BIAS_RATIO) is computed when verifying an ensemble against gridded analyses or point observations. It is defined as the mean error (ME) of ensemble member values greater than or equal to the observation value to which they are matched divided by the absolute value of the mean error (ME) of ensemble member values less than the observation values. + +.. math:: \text{BIAS_RATIO} = \frac{ \text{ME}_{f >= o} }{ |\text{ME}_{f < o}| } + +A perfect forecast has ME = 0. Since BIAS_RATIO is computed as the high bias (ME_GE_OBS) divide by the absolute value of the low bias (ME_LT_OBS), a perfect forecast has BIAS_RATIO = 0/0, which is undefined. In practice, the high and low bias values are unlikely to be 0. + +The range for BIAS_RATIO is 0 to infinity. A score of 1 indicates that the high and low biases are equal. A score greater than 1 indicates that the high bias is larger than the magnitude of the low bias. A score less than 1 indicates the opposite behavior. + IGN --- diff --git a/docs/Users_Guide/ensemble-stat.rst b/docs/Users_Guide/ensemble-stat.rst index 32c1fbd677..db7b393246 100644 --- a/docs/Users_Guide/ensemble-stat.rst +++ b/docs/Users_Guide/ensemble-stat.rst @@ -623,6 +623,21 @@ The format of the STAT and ASCII output of the Ensemble-Stat tool are described * - 41 - CRPS_EMP_FAIR - The Continuous Ranked Probability Skill Score (empirical distribution) adjusted by subtracting 1/2(m) times the mean absolute difference of the ensemble members (m is the ensemble size) + * - 42 + - BIAS_RATIO + - The Bias Ratio + * - 43 + - N_GE_OBS + - The number of ensemble values greater than or equal to their observations + * - 44 + - ME_GE_OBS + - The Mean Error of the ensemble values greater than or equal to their observations + * - 45 + - N_LT_OBS + - The number of ensemble values less than their observations + * - 46 + - ME_LT_OBS + - The Mean Error of the ensemble values less than or equal to their observations .. _table_ES_header_info_es_out_RPS: diff --git a/docs/Users_Guide/reformat_point.rst b/docs/Users_Guide/reformat_point.rst index f4430cd8e1..d1b0cabd6e 100644 --- a/docs/Users_Guide/reformat_point.rst +++ b/docs/Users_Guide/reformat_point.rst @@ -961,13 +961,24 @@ _____________________ { key = "station_id"; val = "station_id,report_identifier"; }, { key = "pressure"; val = "air_pressure,pressure"; }, { key = "height"; val = "height,height_above_mean_sea_level"; }, - { key = "elevation"; val = ""; } + { key = "elevation"; val = "elevation,station_elevation"; } ]; This entry is an array of dictionaries, each containing a **key** string and **val** string which define a mapping of metadata for IODA data files. _____________________ +.. code-block:: none + + obs_to_qc_map = [ + { key = "wind_from_direction"; val = "eastward_wind,northward_wind"; }, + { key = "wind_speed"; val = "eastward_wind,northward_wind"; } + ]; + +This entry is an array of dictionaries, each containing a **key** string and **val** string which define a mapping of QC variable name for IODA data files. + +_____________________ + .. code-block:: none missing_thresh = [ <=-1e9, >=1e9, ==-9999 ]; diff --git a/docs/Users_Guide/tc-pairs.rst b/docs/Users_Guide/tc-pairs.rst index 77ee2af4c7..3e6508446d 100644 --- a/docs/Users_Guide/tc-pairs.rst +++ b/docs/Users_Guide/tc-pairs.rst @@ -9,6 +9,16 @@ Introduction The TC-Pairs tool provides verification for tropical cyclone forecasts in ATCF file format. It matches an ATCF format tropical cyclone (TC) forecast with a second ATCF format reference TC dataset (most commonly the Best Track analysis). The TC-Pairs tool processes both track and intensity adeck data and probabilistic edeck data. The adeck matched pairs contain position errors, as well as wind, sea level pressure, and distance to land values for each TC dataset. The edeck matched pairs contain probabilistic forecast values and the verifying observation values. The pair generation can be subset based on user-defined filtering criteria. Practical aspects of the TC-Pairs tool are described in :numref:`TC-Pairs_Practical-information`. +Scientific and statistical aspects +================================== + +.. _TC-Pairs_Diagnostics: + +Storm Diagnostics +----------------- + +TODO: Add a paragraph about storm diagnostics, describing what they are, why they are important, and how they can be generated. + .. _TC-Pairs_Practical-information: Practical information @@ -24,9 +34,10 @@ The usage statement for tc_pairs is shown below: .. code-block:: none Usage: tc_pairs - -adeck source and/or -edeck source - -bdeck source + -adeck path and/or -edeck path + -bdeck path -config file + [-diag source path] [-out base] [-log file] [-v level] @@ -36,24 +47,30 @@ tc_pairs has required arguments and can accept several optional arguments. Required arguments for tc_pairs ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -1. The **-adeck source** argument indicates the adeck TC-Pairs acceptable format data source containing tropical cyclone model forecast (output from tracker) data to be verified. Acceptable data formats are limited to the standard ATCF format and the one column modified ATCF file, generated by running the tracker in genesis mode. It specifies the name of a TC-Pairs acceptable format file or top-level directory containing TC-Pairs acceptable format files ending in ".dat" to be processed. The **-adeck** or **-edeck** option must be used at least once. +1. The **-adeck path** argument indicates the adeck TC-Pairs acceptable format data containing tropical cyclone model forecast (output from tracker) data to be verified. Acceptable data formats are limited to the standard ATCF format and the one column modified ATCF file, generated by running the tracker in genesis mode. It specifies the name of a TC-Pairs acceptable format file or top-level directory containing TC-Pairs acceptable format files ending in ".dat" to be processed. The **-adeck** or **-edeck** option must be used at least once. -2. The **-edeck source** argument indicates the edeck ATCF format data source containing probabilistic track data to be verified. It specifies the name of an ATCF format file or top-level directory containing ATCF format files ending in ".dat" to be processed. The **-adeck** or **-edeck** option must be used at least once. +2. The **-edeck path** argument indicates the edeck ATCF format data containing probabilistic track data to be verified. It specifies the name of an ATCF format file or top-level directory containing ATCF format files ending in ".dat" to be processed. The **-adeck** or **-edeck** option must be used at least once. -3. The **-bdeck source** argument indicates the TC-Pairs acceptable format data source containing the tropical cyclone reference dataset to be used for verifying the adeck source. This source is typically the NHC Best Track Analysis, but could be any TC-Pairs acceptable formatted reference. The acceptable data formats for bdecks are the same as those for adecks. This argument specifies the name of a TC-Pairs acceptable format file or top-level directory containing TC-Pairs acceptable format files ending in ".dat" to be processed. +3. The **-bdeck path** argument indicates the TC-Pairs acceptable format data containing the tropical cyclone reference dataset to be used for verifying the adeck data. This data is typically the NHC Best Track Analysis, but could be any TC-Pairs acceptable formatted reference. The acceptable data formats for bdecks are the same as those for adecks. This argument specifies the name of a TC-Pairs acceptable format file or top-level directory containing TC-Pairs acceptable format files ending in ".dat" to be processed. 4. The **-config file** argument indicates the name of the configuration file to be used. The contents of the configuration file are discussed below. Optional arguments for tc_pairs ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -5. The -**out base** argument indicates the path of the output file base. This argument overrides the default output file base (**./out_tcmpr**). +5. The **-diag source path** argument indicates the TC-Pairs acceptable format data containing the tropical cyclone diagnostics dataset corresponding to the adeck tracks. The **source** can be set to TCDIAG, LSDIAG_RT, or LSDIAG_DEV to indicate the input diagnostics data source. The **path** argument specifies the name of a TC-Pairs acceptable format file or top-level directory containing TC-Pairs acceptable format files ending in ".dat" to be processed. + +6. The -**out base** argument indicates the path of the output file base. This argument overrides the default output file base (**./out_tcmpr**). + +7. The **-log file** option directs output and errors to the specified log file. All messages will be written to that file as well as standard out and error. Thus, users can save the messages without having to redirect the output on the command line. The default behavior is no log file. -6. The **-log file** option directs output and errors to the specified log file. All messages will be written to that file as well as standard out and error. Thus, users can save the messages without having to redirect the output on the command line. The default behavior is no log file. +8. The **-v level** option indicates the desired level of verbosity. The contents of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity above 1 will increase the amount of logging. -7. The **-v level** option indicates the desired level of verbosity. The contents of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity above 1 will increase the amount of logging. +This tool currently only supports the rapid intensification (**RI**) edeck probability type but support for additional edeck probability types will be added in future releases. -This tool currently only supports the rapid intensification (**RI**) edeck probability type but support for additional edeck probability types will be added in future releases. At least one **-adeck** or **-edeck** option must be specified. The **-adeck, -edeck**, and **-bdeck** options may optionally be followed with **suffix=string** to append that string to all model names found within that data source. This option may be useful when processing track data from two different sources which reuse the same model names. +At least one **-adeck** or **-edeck** option must be specified. The **-adeck, -edeck**, and **-bdeck** options may optionally be followed with **suffix=string** to append that string to all model names found within that data source. This option may be useful when processing track data from two different sources which reuse the same model names. + +The **-diag** option may optionally be followed with **model=string** to override the model name of the tracks to which those diagnostics correspond. The **string** specifies a comma-separated list of one or more ATCF ID's to which these diagnostics should be paired (e.g. **model=OFCL,SHIP**). An example of the tc_pairs calling sequence is shown below: @@ -67,6 +84,8 @@ The TC-Pairs tool implements the following logic: • Parse the adeck, edeck, and bdeck data files and store them as track objects. +• Parse diagnostics data files and add the requested diagnostics to the existing adeck track objects. + • Apply configuration file settings to filter the adeck, edeck, and bdeck track data down to a subset of interest. • Apply configuration file settings to derive additional adeck track data, such as interpolated tracks, consensus tracks, time-lagged tracks, and statistical track and intensity models. @@ -235,6 +254,59 @@ ____________________ The **watch_warn** field specifies the file name and time applied offset to the **watch_warn** flag. The **file_name** string specifies the path of the watch/warning file to be used to determine when a watch or warning is in effect during the forecast initialization and verification times. The default file is named **wwpts_us.txt**, which is found in the installed *share/met/tc_data/* directory within the MET build. The **time_offset** string is the time window (in seconds) assigned to the watch/warning. Due to the non-uniform time watches and warnings are issued, a time window is assigned for which watch/warnings are included in the verification for each valid time. The default watch/warn file is static, and therefore may not include warned storms beyond the current MET code release date; therefore users may wish to create a post in the `METplus GitHub Discussions Forum `_ in order to obtain the most recent watch/warning file if the static file does not contain storms of interest. +____________________ + +.. code-block:: none + + diag_name = []; + +The **diag_name** entry specifies a comma-separated list of strings for the tropical cyclone diagnostics of interest. This applies when the **-tcdiag** and/or **-lsdiag** command line options have been used to provide storm diagnostics data. If a non-zero list of diagnostic names is specified, only those diagnostics appearing in the list are written to the TCDIAG output line type. If defined as an empty list (default), all diagnostics found in the input are written to the TCDIAG output lines. + +A TCMPR line is written to the output for each track point. If diagnostics data is also defined for that track point, a TCDIAG line is written immediately after the corresponding TCMPR line. The contents of that TCDIAG line is determined by diagnostic names requested in the **diag_name** entry. + +____________________ + +.. code-block:: none + + diag_convert_map = [ + { source = "TCDIAG"; + key = [ "(10C)", "(10KT)", "(10M/S)" ]; + convert(x) = x / 10; }, + + { source = "LSDIAG_RT"; + key = [ "LAT", "LON", "CSST", "RSST", "DSST", "DSTA", "XDST", "XNST", "NSST", "NSTA", + "NTMX", "NTFR", "U200", "U20C", "V20C", "E000", "EPOS", "ENEG", "EPSS", "ENSS", + "T000", "TLAT", "TLON", "TWAC", "TWXC", "G150", "G200", "G250", "V000", "V850", + "V500", "V300", "SHDC", "SHGC", "T150", "T200", "T250", "SHRD", "SHRS", "SHRG", + "HE07", "HE05", "PW01", "PW02", "PW03", "PW04", "PW05", "PW06", "PW07", "PW08", + "PW09", "PW10", "PW11", "PW12", "PW13", "PW14", "PW15", "PW16", "PW17", "PW18", + "PW20", "PW21" ]; + convert(x) = x / 10; }, + + { source = "LSDIAG_RT"; + key = [ "VVAV", "VMFX", "VVAC" ]; + convert(x) = x / 100; }, + + { source = "LSDIAG_RT"; + key = [ "TADV" ]; + convert(x) = x / 1000000; }, + + { source = "LSDIAG_RT"; + key = [ "Z850", "D200", "TGRD", "DIVC" ]; + convert(x) = x / 10000000; }, + + { source = "LSDIAG_RT"; + key = [ "PENC", "PENV" ]; + convert(x) = x / 10 + 1000; } + + ]; + +The **diag_convert_map** entries define conversion functions to be applied to diagnostics data read with the **-diag** command line option. Each array element is a dictionary consisting of a **source**, **key**, and **convert(x)** entry. + +The **source** is one of the supported diagnostics data sources. The **key** is an array of strings. The strings can specify diagnostic names or units, although units are only checked for **TCDIAG** sources. If both the name and units are specified, the conversion function for the name takes precedence. **convert(x)** is a function of one variable which defines how the diagnostic data should be converted. The defined function is applied to any diagnostic value whose name or units appears in the **key**. + +____________________ + .. code-block:: none basin_map = [ @@ -487,7 +559,41 @@ TC-Pairs produces output in TCST format. The default output file name can be ove * - 84 - MAX_WIND_SPREAD - consensus variable: the standard deviation of the member's maximum wind speed values - + +.. _TCDIAG Line Type: + +.. list-table:: Format information for TCDIAG (Tropical Cyclone Diagnostics) output line type. + :widths: auto + :header-rows: 2 + + * - + - + - TCDIAG OUTPUT FORMAT + * - Column Number + - Header Column Name + - Description + * - 13 + - TCDIAG + - Tropical Cyclone Diagnostics line type + * - 14 + - TOTAL + - Total number of pairs in track + * - 15 + - INDEX + - Index of the current track pair + * - 16 + - SOURCE + - Diagnostics data source + * - 17 + - N_DIAG + - Number of storm diagnostic name and value columns to follow + * - 18 + - DIAG_i + - Name of the of the ith storm diagnostic (repeated) + * - 19 + - VALUE_i + - Value of the ith storm diagnostic (repeated) + .. _PROBRIRW Line Type: .. list-table:: Format information for PROBRIRW (Probability of Rapid Intensification/Weakening) output line type. diff --git a/docs/Users_Guide/tc-stat.rst b/docs/Users_Guide/tc-stat.rst index 6db00f3160..03013bdf94 100644 --- a/docs/Users_Guide/tc-stat.rst +++ b/docs/Users_Guide/tc-stat.rst @@ -245,7 +245,7 @@ _________________________ column_str_name = []; column_str_val = []; -The **column_str_name** and **column_str_val** fields stratify by performing string matching on non-numeric data columns. Specify a comma-separated list of columns names and values to be **included** in the analysis. The length of the **column_str_val** should match that of the **column_str_name**. Using the **-column_str name val** option within the job command lines may further refine these selections. +The **column_str_name** and **column_str_val** fields stratify by performing string matching on non-numeric data columns. Specify a comma-separated list of columns names and values to be **included** in the analysis. The length of the **column_str_val** should match that of the **column_str_name**. Using the **-column_str name value** option within the job command lines may further refine these selections. _________________________ @@ -254,7 +254,7 @@ _________________________ column_str_exc_name = []; column_str_exc_val = []; -The **column_str_exc_name** and **column_str_exc_val** fields stratify by performing string matching on non-numeric data columns. Specify a comma-separated list of columns names and values to be **excluded** from the analysis. The length of the **column_str_exc_val** should match that of the **column_str_exc_name**. Using the **-column_str_exc name val** option within the job command lines may further refine these selections. +The **column_str_exc_name** and **column_str_exc_val** fields stratify by performing string matching on non-numeric data columns. Specify a comma-separated list of columns names and values to be **excluded** from the analysis. The length of the **column_str_exc_val** should match that of the **column_str_exc_name**. Using the **-column_str_exc name value** option within the job command lines may further refine these selections. _________________________ @@ -263,7 +263,7 @@ _________________________ init_thresh_name = []; init_thresh_val = []; -The **init_thresh_name** and **init_thresh_val** fields stratify by applying thresholds to numeric data columns only when lead = 0. If lead = 0, but the value does not meet the threshold, discard the entire track. The length of the **init_thresh_val** should match that of the **init_thresh_name**. Using the **-init_thresh name val** option within the job command lines may further refine these selections. +The **init_thresh_name** and **init_thresh_val** fields stratify by applying thresholds to numeric data columns only when lead = 0. If lead = 0, but the value does not meet the threshold, discard the entire track. The length of the **init_thresh_val** should match that of the **init_thresh_name**. Using the **-init_thresh name thresh** option within the job command lines may further refine these selections. _________________________ @@ -272,7 +272,7 @@ _________________________ init_str_name = []; init_str_val = []; -The **init_str_name** and **init_str_val** fields stratify by performing string matching on non-numeric data columns only when lead = 0. If lead = 0, but the string **does not** match, discard the entire track. The length of the **init_str_val** should match that of the **init_str_name**. Using the **-init_str name val** option within the job command lines may further refine these selections. +The **init_str_name** and **init_str_val** fields stratify by performing string matching on non-numeric data columns only when lead = 0. If lead = 0, but the string **does not** match, discard the entire track. The length of the **init_str_val** should match that of the **init_str_name**. Using the **-init_str name value** option within the job command lines may further refine these selections. _________________________ @@ -281,7 +281,25 @@ _________________________ init_str_exc_name = []; init_str_exc_val = []; -The **init_str_exc_name** and **init_str_exc_val** fields stratify by performing string matching on non-numeric data columns only when lead = 0. If lead = 0, and the string **does** match, discard the entire track. The length of the **init_str_exc_val** should match that of the **init_str_exc_name**. Using the **-init_str_exc name val** option within the job command lines may further refine these selections. +The **init_str_exc_name** and **init_str_exc_val** fields stratify by performing string matching on non-numeric data columns only when lead = 0. If lead = 0, and the string **does** match, discard the entire track. The length of the **init_str_exc_val** should match that of the **init_str_exc_name**. Using the **-init_str_exc name value** option within the job command lines may further refine these selections. + +_________________________ + +.. code-block:: none + + diag_thresh_name = []; + diag_thresh_val = []; + +The **diag_thresh_name** and **diag_thresh_val** fields stratify individual track points by applying thresholds to numeric data columns from the TCDIAG lines. Specify a comma-separated list of diagnostics names and thresholds to be applied. The length of **diag_thresh_val** should match that of **diag_thresh_name**. If the storm diagnostic does not meet the threshold, discard both the TCMPR and TCDIAG lines for that track point. Using the **-diag_thresh name thresh** option within the job command lines may further refine these selections. + +_________________________ + +.. code-block:: none + + init_diag_thresh_name = []; + init_diag_thresh_val = []; + +The **init_diag_thresh_name** and **init_diag_thresh_val** fields stratify entire tracks by applying thresholds to numeric data columns from the TCDIAG lines, but only when lead = 0. If lead = 0, but the storm diagnostic does not meet the threshold, discard the entire track. The length of the **init_diag_thresh_val** should match that of the **init_diag_thresh_name**. Using the **-init_diag_thresh name thresh** option within the job command lines may further refine these selections. _________________________ diff --git a/internal/scripts/docker/Dockerfile b/internal/scripts/docker/Dockerfile index 4dafd977b3..8ec46e33db 100644 --- a/internal/scripts/docker/Dockerfile +++ b/internal/scripts/docker/Dockerfile @@ -1,5 +1,5 @@ ARG MET_BASE_REPO=met-base -ARG MET_BASE_TAG=v1.0 +ARG MET_BASE_TAG=v1.1 FROM dtcenter/${MET_BASE_REPO}:${MET_BASE_TAG} MAINTAINER John Halley Gotway diff --git a/internal/scripts/docker/Dockerfile.copy b/internal/scripts/docker/Dockerfile.copy index 5c51408fdf..7d0fc6ae6f 100644 --- a/internal/scripts/docker/Dockerfile.copy +++ b/internal/scripts/docker/Dockerfile.copy @@ -1,5 +1,5 @@ ARG MET_BASE_REPO=met-base-unit-test -ARG MET_BASE_TAG=v1.0 +ARG MET_BASE_TAG=v1.1 FROM dtcenter/${MET_BASE_REPO}:${MET_BASE_TAG} MAINTAINER John Halley Gotway diff --git a/internal/scripts/environment/development.docker b/internal/scripts/environment/development.docker index ec4bd202c3..c6674b231b 100644 --- a/internal/scripts/environment/development.docker +++ b/internal/scripts/environment/development.docker @@ -24,7 +24,7 @@ export MET_FREETYPELIB=/usr/lib64 export MET_JASPERLIB=/usr/lib64 export MET_PYTHON=/usr/bin/python3 -export MET_PYTHON_CC="-I/usr/include/python3.6m -I/usr/include/python3.6m" +export MET_PYTHON_CC="-I/usr/include/python3.6m" export MET_PYTHON_LD="-L/usr/lib64 -lpython3.6m -lpthread -ldl -lutil -lm" # -D__64BIT__ is required because we've compiled libgrib2c.a with that flag diff --git a/internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS b/internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS index 9372a6bb0f..5e7c75f5e2 100644 --- a/internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS +++ b/internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS @@ -42,8 +42,8 @@ data = { range = [0,25]; field = [ - { set_attr_name = "FCST_APCP"; }, - { set_attr_name = "OBS_APCP"; } + { ${FIELD1} }, + { ${FIELD2} } ]; } diff --git a/internal/test_unit/config/TCPairsConfig_ALAL2010 b/internal/test_unit/config/TCPairsConfig_ALAL2010 index 0dfc50ee28..a0c2b582c7 100644 --- a/internal/test_unit/config/TCPairsConfig_ALAL2010 +++ b/internal/test_unit/config/TCPairsConfig_ALAL2010 @@ -128,6 +128,11 @@ watch_warn = { time_offset = -14400; } +// +// Diagnostics to be extracted +// +diag_name = []; + // // Modify basin names to make them consistent across ATCF input files. // diff --git a/internal/test_unit/config/TCPairsConfig_BASIN_MAP b/internal/test_unit/config/TCPairsConfig_BASIN_MAP index fb0075ee17..150e4ed4fa 100644 --- a/internal/test_unit/config/TCPairsConfig_BASIN_MAP +++ b/internal/test_unit/config/TCPairsConfig_BASIN_MAP @@ -128,6 +128,11 @@ watch_warn = { time_offset = -14400; } +// +// Diagnostics to be extracted +// +diag_name = []; + // // Modify basin names to make them consistent across ATCF input files. // diff --git a/internal/test_unit/config/TCPairsConfig_CONSENSUS b/internal/test_unit/config/TCPairsConfig_CONSENSUS index 16b2106a2a..b157c71264 100644 --- a/internal/test_unit/config/TCPairsConfig_CONSENSUS +++ b/internal/test_unit/config/TCPairsConfig_CONSENSUS @@ -180,6 +180,11 @@ watch_warn = { time_offset = -14400; } +// +// Diagnostics to be extracted +// +diag_name = []; + // // Modify basin names to make them consistent across ATCF input files. // diff --git a/internal/test_unit/config/TCPairsConfig_DIAGNOSTICS b/internal/test_unit/config/TCPairsConfig_DIAGNOSTICS new file mode 100644 index 0000000000..a2ae20816b --- /dev/null +++ b/internal/test_unit/config/TCPairsConfig_DIAGNOSTICS @@ -0,0 +1,157 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// TC-Pairs configuration file. +// +// For additional information, please see the MET User's Guide. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// ATCF file format reference: +// http://www.nrlmry.navy.mil/atcf_web/docs/database/new/abrdeck.html +// + +// +// Models +// +model = []; + +// +// Description +// +desc = "NA"; + +// +// Storm identifiers +// +storm_id = [ ${STORM_ID} ]; + +// +// Basins +// +basin = []; + +// +// Cyclone numbers +// +cyclone = []; + +// +// Storm names +// +storm_name = []; + +// +// Model initialization time windows to include or exclude +// +init_beg = ${INIT_BEG}; +init_end = ${INIT_END}; +init_inc = []; +init_exc = []; + +// +// Valid model time windows to include or exclude +// +valid_beg = ""; +valid_end = ""; +valid_inc = []; +valid_exc = []; + +// +// Valid times for which output should be written +// +write_valid = []; + +// +// Model initialization hours +// +init_hour = []; + +// +// Required lead time in hours +// +lead_req = []; + +// +// lat/lon polylines defining masking regions +// +init_mask = ""; +valid_mask = ""; + +// +// Specify if the code should check for duplicate ATCF lines +// +check_dup = TRUE; + +// +// Specify special processing to be performed for interpolated models. +// Set to NONE, FILL, or REPLACE. +// +interp12 = NONE; + +// +// Specify how consensus forecasts should be defined +// +consensus = []; + +// +// Forecast lag times +// +lag_time = []; + +// +// CLIPER/SHIFOR baseline forecasts to be derived from the BEST +// and operational (CARQ) tracks. +// +best_baseline = []; +oper_baseline = []; + +// +// Specify if only those track points common to both the ADECK and BDECK +// tracks be written out. +// +match_points = FALSE; + +// +// Specify the NetCDF output of the gen_dland tool containing a gridded +// representation of the minimum distance to land. +// +dland_file = "${MET_TEST_OUTPUT}/tc_dland/tc_dland_half_deg.nc"; + +// +// Specify watch/warning information: +// - Input watch/warning filename +// - Watch/warning time offset in seconds +// +watch_warn = { + file_name = "MET_BASE/tc_data/wwpts_us.txt"; + time_offset = -14400; +} + +// +// Diagnostics to be extracted +// +diag_name = [ ${DIAG_NAME} ]; + +// +// Unit conversions to be applied based on diagnostic names and units +// Commented out to use settings from the default config file +// +// diag_convert_map = []; + +// +// Modify basin names to make them consistent across ATCF input files. +// +basin_map = [ + { key = "SI"; val = "SH"; }, + { key = "SP"; val = "SH"; }, + { key = "AU"; val = "SH"; }, + { key = "AB"; val = "IO"; }, + { key = "BB"; val = "IO"; } +]; + +// +// Indicate a version number for the contents of this configuration file. +// The value should generally not be modified. +// +version = "V11.0.0"; diff --git a/internal/test_unit/config/TCPairsConfig_INTERP12 b/internal/test_unit/config/TCPairsConfig_INTERP12 index 92e00a10b5..2985f596ee 100644 --- a/internal/test_unit/config/TCPairsConfig_INTERP12 +++ b/internal/test_unit/config/TCPairsConfig_INTERP12 @@ -72,7 +72,6 @@ init_hour = []; // lead_req = []; - // // lat/lon polylines defining masking regions // @@ -129,6 +128,11 @@ watch_warn = { time_offset = -14400; } +// +// Diagnostics to be extracted +// +diag_name = []; + // // Modify basin names to make them consistent across ATCF input files. // diff --git a/internal/test_unit/config/TCPairsConfig_PROBRIRW b/internal/test_unit/config/TCPairsConfig_PROBRIRW index 1632603972..c958ec5590 100644 --- a/internal/test_unit/config/TCPairsConfig_PROBRIRW +++ b/internal/test_unit/config/TCPairsConfig_PROBRIRW @@ -128,6 +128,11 @@ watch_warn = { time_offset = -14400; } +// +// Diagnostics to be extracted +// +diag_name = []; + // // Modify basin names to make them consistent across ATCF input files. // diff --git a/internal/test_unit/config/TCStatConfig_ALAL2010 b/internal/test_unit/config/TCStatConfig_ALAL2010 index 85c3141ffe..27d64c0478 100644 --- a/internal/test_unit/config/TCStatConfig_ALAL2010 +++ b/internal/test_unit/config/TCStatConfig_ALAL2010 @@ -136,6 +136,20 @@ init_str_val = []; init_str_exc_name = []; init_str_exc_val = []; +// +// Stratify track points by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines. +// +diag_thresh_name = []; +diag_thresh_val = []; + +// +// Stratify entire tracks by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines for lead time 0. +// +init_diag_thresh_name = []; +init_diag_thresh_val = []; + // // Stratify by the ADECK and BDECK distances to land. // diff --git a/internal/test_unit/config/TCStatConfig_DIAGNOSTICS b/internal/test_unit/config/TCStatConfig_DIAGNOSTICS new file mode 100644 index 0000000000..1801bf5ad9 --- /dev/null +++ b/internal/test_unit/config/TCStatConfig_DIAGNOSTICS @@ -0,0 +1,226 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// TC-Stat configuration file. +// +// For additional information, please see the MET User's Guide. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// The parameters listed below are used to filter the TC-STAT data down to the +// desired subset of lines over which statistics are to be computed. Only those +// lines which meet ALL of the criteria specified will be retained. +// +// The settings that are common to all jobs may be specified once at the top +// level. If no selection is listed for a parameter, that parameter will not +// be used for filtering. If multiple selections are listed for a parameter, +// the analyses will be performed on their union. +// +// Each configuration filtering option may be overridden by a corresponding +// job command option of the same name, as described in the MET User's Guide. +// + +// +// Stratify by the AMODEL or BMODEL columns. +// +amodel = []; +bmodel = []; + +// +// Stratify by the DESC column. +// +desc = []; + +// +// Stratify by the STORM_ID column. +// +storm_id = []; + +// +// Stratify by the BASIN column. +// +basin = []; + +// +// Stratify by the CYCLONE column. +// +cyclone = []; + +// +// Stratify by the STORM_NAME column. +// +storm_name = []; + +// +// Stratify by the INIT times. +// Model initialization time windows to include or exclude +// +init_beg = ""; +init_end = ""; +init_inc = []; +init_exc = []; + +// +// Stratify by the VALID times. +// +valid_beg = ""; +valid_end = ""; +valid_inc = []; +valid_exc = []; + +// +// Stratify by the initialization and valid hours and lead time. +// +init_hour = []; +valid_hour = []; +lead = []; + +// +// Select tracks which contain all required lead times. +// +lead_req = []; + +// +// Stratify by the INIT_MASK and VALID_MASK columns. +// +init_mask = []; +valid_mask = []; + +// +// Stratify by the LINE_TYPE column. +// +line_type = []; + +// +// Stratify by checking the watch/warning status for each track point +// common to both the ADECK and BDECK tracks. If the watch/warning status +// of any of the track points appears in the list, retain the entire track. +// +track_watch_warn = []; + +// +// Stratify by applying thresholds to numeric data columns. +// +column_thresh_name = []; +column_thresh_val = []; + +// +// Stratify by performing string matching on non-numeric data columns. +// +column_str_name = []; +column_str_val = []; + +// +// Stratify by excluding strings in non-numeric data columns. +// +column_str_exc_name = []; +column_str_exc_val = []; + +// +// Similar to the column_thresh options above +// +init_thresh_name = []; +init_thresh_val = []; + +// +// Similar to the column_str options above +// +init_str_name = []; +init_str_val = []; + +// +// Similar to the column_str_exc options above +// +init_str_exc_name = []; +init_str_exc_val = []; + +// +// Stratify track points by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines. +// +diag_thresh_name = []; +diag_thresh_val = []; + +// +// Stratify entire tracks by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines for lead time 0. +// +init_diag_thresh_name = []; +init_diag_thresh_val = []; + +// +// Stratify by the ADECK and BDECK distances to land. +// +water_only = FALSE; + +// +// Specify whether only those track points for which rapid intensification +// or weakening of the maximum wind speed occurred in the previous time +// step should be retained. +// +rirw = { + track = NONE; + adeck = { + time = "24"; + exact = TRUE; + thresh = >=30.0; + } + bdeck = adeck; +} + +// +// Specify whether only those track points occurring near landfall should be +// retained, and define the landfall retention window in HH[MMSS] format +// around the landfall time. +// +landfall = FALSE; +landfall_beg = "-24"; +landfall_end = "00"; + +// +// Specify whether only those track points common to both the ADECK and BDECK +// tracks should be retained. +// +match_points = TRUE; + +// +// Specify whether only those cases common to all models in the dataset should +// be retained. +// +event_equal = FALSE; + +// +// Specify lead times that must be present for a track to be included in the +// event equalization logic. +// +event_equal_lead = []; + +// +// Apply polyline masking logic to the location of the ADECK track at the +// initialization time. +// +out_init_mask = ""; + +// +// Apply polyline masking logic to the location of the ADECK track at the +// valid time. +// +out_valid_mask = ""; + +// +// Array of TCStat analysis jobs to be performed on the filtered data +// +jobs = [ + "-job filter -dump_row ${MET_TEST_OUTPUT}/tc_stat/DIAGNOSTICS_filter_match_points.tcst", + "-job summary -diag_thresh PW01 lt60 -column TK_ERR -by AMODEL", + "-job summary -diag_thresh PW01 ge60 -column TK_ERR -by AMODEL", + "-job summary -diag_thresh TPW lt60 -column TK_ERR -by AMODEL", + "-job summary -diag_thresh TPW ge60 -column TK_ERR -by AMODEL", + "-job summary -init_diag_thresh DTL lt325 -amodel OFCL -column TK_ERR -by AMODEL,LEAD" +]; + +// +// Indicate a version number for the contents of this configuration file. +// The value should generally not be modified. +// +version = "V11.0.0"; diff --git a/internal/test_unit/config/TCStatConfig_PROBRIRW b/internal/test_unit/config/TCStatConfig_PROBRIRW index a48bcc6dba..ba9a79b1d7 100644 --- a/internal/test_unit/config/TCStatConfig_PROBRIRW +++ b/internal/test_unit/config/TCStatConfig_PROBRIRW @@ -136,6 +136,20 @@ init_str_val = []; init_str_exc_name = []; init_str_exc_val = []; +// +// Stratify track points by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines. +// +diag_thresh_name = []; +diag_thresh_val = []; + +// +// Stratify entire tracks by applying thresholds to numeric +// storm diagnostic values in the TCDIAG lines for lead time 0. +// +init_diag_thresh_name = []; +init_diag_thresh_val = []; + // // Stratify by the ADECK and BDECK distances to land. // diff --git a/internal/test_unit/hdr/met_11_0.hdr b/internal/test_unit/hdr/met_11_0.hdr index c1d7e318f0..a559da2695 100644 --- a/internal/test_unit/hdr/met_11_0.hdr +++ b/internal/test_unit/hdr/met_11_0.hdr @@ -17,7 +17,7 @@ PJC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_L PRC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_THRESH _VAR_ PSTD : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_THRESH BASER BASER_NCL BASER_NCU RELIABILITY RESOLUTION UNCERTAINTY ROC_AUC BRIER BRIER_NCL BRIER_NCU BRIERCL BRIERCL_NCL BRIERCL_NCU BSS BSS_SMPL _VAR_ ECLV : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BASE N_PTS _VAR_ -ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP CRPS_EMP_FAIR +ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP CRPS_EMP_FAIR BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS RPS : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_PROB RPS_REL RPS_RES RPS_UNC RPS RPSS RPSS_SMPL RPS_COMP RHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL CRPS IGN N_RANK CRPSS SPREAD _VAR_ PHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BIN_SIZE N_BIN _VAR_ @@ -34,4 +34,5 @@ MODE_SOA : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_L MODE_POA : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE N_VALID GRID_RES OBJECT_ID OBJECT_CAT CENTROID_DIST BOUNDARY_DIST CONVEX_HULL_DIST ANGLE_DIFF ASPECT_DIFF AREA_RATIO INTERSECTION_AREA UNION_AREA SYMMETRIC_DIFF INTERSECTION_OVER_AREA CURVATURE_RATIO COMPLEXITY_RATIO PERCENTILE_INTENSITY_RATIO INTEREST MODE_CTS : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE N_VALID GRID_RES FIELD TOTAL FY_OY FY_ON FN_OY FN_ON BASER FMEAN ACC FBIAS PODY PODN POFD FAR CSI GSS HK HSS ODDS TCST_TCMPR : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE AMODEL BMODEL STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK TOTAL INDEX LEVEL WATCH_WARN INITIALS ALAT ALON BLAT BLON TK_ERR X_ERR Y_ERR ALTK_ERR CRTK_ERR ADLAND BDLAND AMSLP BMSLP AMAX_WIND BMAX_WIND AAL_WIND_34 BAL_WIND_34 ANE_WIND_34 BNE_WIND_34 ASE_WIND_34 BSE_WIND_34 ASW_WIND_34 BSW_WIND_34 ANW_WIND_34 BNW_WIND_34 AAL_WIND_50 BAL_WIND_50 ANE_WIND_50 BNE_WIND_50 ASE_WIND_50 BSE_WIND_50 ASW_WIND_50 BSW_WIND_50 ANW_WIND_50 BNW_WIND_50 AAL_WIND_64 BAL_WIND_64 ANE_WIND_64 BNE_WIND_64 ASE_WIND_64 BSE_WIND_64 ASW_WIND_64 BSW_WIND_64 ANW_WIND_64 BNW_WIND_64 ARADP BRADP ARRP BRRP AMRD BMRD AGUSTS BGUSTS AEYE BEYE ADIR BDIR ASPEED BSPEED ADEPTH BDEPTH NUM_MEMBERS TRACK_SPREAD DIST_MEAN MSLP_SPREAD MAX_WIND_SPREAD +TCST_TCDIAG : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE AMODEL BMODEL STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK TOTAL INDEX SOURCE N_DIAG _VAR_ TCST_PROBRIRW : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE ALAT ALON BLAT BLON INITIALS TK_ERR X_ERR Y_ERR ADLAND BDLAND RI_BEG RI_END RI_WINDOW AWIND_END BWIND_BEG BWIND_END BDELTA BDELTA_MAX BLEVEL_BEG BLEVEL_END N_THRESH _VAR_ diff --git a/internal/test_unit/xml/unit_grid_diag.xml b/internal/test_unit/xml/unit_grid_diag.xml index 10049144ac..96c9cb217f 100644 --- a/internal/test_unit/xml/unit_grid_diag.xml +++ b/internal/test_unit/xml/unit_grid_diag.xml @@ -72,6 +72,10 @@ &DATA_DIR_OBS;/stage4_hmt/stage4_2012041018_06h.grib" \ > &OUTPUT_DIR;/grid_diag/obs_file_list; \ &MET_BIN;/grid_diag + + FIELD1 set_attr_name = "FCST_APCP"; + FIELD2 set_attr_name = "OBS_APCP"; + \ -data &OUTPUT_DIR;/grid_diag/fcst_file_list \ -data &OUTPUT_DIR;/grid_diag/obs_file_list \ @@ -84,4 +88,40 @@ + + + + echo "&DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F006.grib \ + &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F012.grib \ + &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F018.grib \ + &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F024.grib \ + &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F030.grib \ + &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F036.grib \ + &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F042.grib" \ + > &OUTPUT_DIR;/grid_diag/fcst_file_list; \ + echo "&DATA_DIR_OBS;/stage4_hmt/stage4_2012040906_06h.grib \ + &DATA_DIR_OBS;/stage4_hmt/stage4_2012040912_06h.grib \ + &DATA_DIR_OBS;/stage4_hmt/stage4_2012040918_06h.grib \ + &DATA_DIR_OBS;/stage4_hmt/stage4_2012041000_06h.grib \ + &DATA_DIR_OBS;/stage4_hmt/stage4_2012041006_06h.grib \ + &DATA_DIR_OBS;/stage4_hmt/stage4_2012041012_06h.grib \ + &DATA_DIR_OBS;/stage4_hmt/stage4_2012041018_06h.grib" \ + > &OUTPUT_DIR;/grid_diag/obs_file_list; \ + &MET_BIN;/grid_diag + + FIELD1 name = "APCP"; + FIELD2 name = "APCP"; + + \ + -data &OUTPUT_DIR;/grid_diag/fcst_file_list \ + -data &OUTPUT_DIR;/grid_diag/obs_file_list \ + -config &CONFIG_DIR;/GridDiagConfig_APCP_06_FCST_OBS \ + -out &OUTPUT_DIR;/grid_diag/grid_diag_APCP_06_VARN.nc \ + -v 3 + + + &OUTPUT_DIR;/grid_diag/grid_diag_APCP_06_VARN.nc + + + diff --git a/internal/test_unit/xml/unit_ioda2nc.xml b/internal/test_unit/xml/unit_ioda2nc.xml index 91f2a7cbb4..55165568bd 100644 --- a/internal/test_unit/xml/unit_ioda2nc.xml +++ b/internal/test_unit/xml/unit_ioda2nc.xml @@ -83,4 +83,42 @@ + + &MET_BIN;/ioda2nc + + STATION_ID + MASK_GRID + MASK_POLY + MESSAGE_TYPE + + \ + &DATA_DIR_OBS;/ioda/jopa_satwind_20210701T1200Z_out_0000_reduced.nc4 \ + &OUTPUT_DIR;/ioda2nc/jopa_satwind_20210701T1200Z_int_datetime.nc \ + -config &CONFIG_DIR;/IODA2NCConfig_mask \ + -v 2 + + + &OUTPUT_DIR;/ioda2nc/jopa_satwind_20210701T1200Z_int_datetime.nc + + + + + &MET_BIN;/ioda2nc + + STATION_ID + MASK_GRID + MASK_POLY + MESSAGE_TYPE + + \ + &DATA_DIR_OBS;/ioda/2021081612_sonde_small.nc \ + &OUTPUT_DIR;/ioda2nc/2021081612_sonde_small_sid.nc \ + -config &CONFIG_DIR;/IODA2NCConfig_mask \ + -v 2 + + + &OUTPUT_DIR;/ioda2nc/2021081612_sonde_small_sid.nc + + + diff --git a/internal/test_unit/xml/unit_tc_pairs.xml b/internal/test_unit/xml/unit_tc_pairs.xml index 44f986fd10..0080eb53b1 100644 --- a/internal/test_unit/xml/unit_tc_pairs.xml +++ b/internal/test_unit/xml/unit_tc_pairs.xml @@ -182,4 +182,27 @@ + + &MET_BIN;/tc_pairs + + STORM_ID "AL092022" + INIT_BEG "20220926_00" + INIT_END "20220926_18" + DIAG_NAME "DTL", "PW01", "SHRD", "TPW", "LAND", "SHR_MAG", "STM_SPD" + + \ + -adeck &DATA_DIR;/adeck/aal092022_OFCL_SHIP_AVNO.dat \ + -bdeck &DATA_DIR;/bdeck/bal092022.dat \ + -diag TCDIAG &DATA_DIR;/tcdiag/2022/sal092022_avno_doper_20220926*_diag.dat \ + -diag LSDIAG_RT &DATA_DIR;/lsdiag_rt/2022/220926*AL0922_lsdiag.dat model=OFCL,SHIP \ + -config &CONFIG_DIR;/TCPairsConfig_DIAGNOSTICS \ + -out &OUTPUT_DIR;/tc_pairs/al092022_20220926_DIAGNOSTICS \ + -log &OUTPUT_DIR;/tc_pairs/tc_pairs_DIAGNOSTICS.log \ + -v 4 + + + &OUTPUT_DIR;/tc_pairs/al092022_20220926_DIAGNOSTICS.tcst + + + diff --git a/internal/test_unit/xml/unit_tc_stat.xml b/internal/test_unit/xml/unit_tc_stat.xml index 5029a1f9a8..4fd270718d 100644 --- a/internal/test_unit/xml/unit_tc_stat.xml +++ b/internal/test_unit/xml/unit_tc_stat.xml @@ -92,4 +92,17 @@ + + &MET_BIN;/tc_stat + \ + -lookin &OUTPUT_DIR;/tc_pairs/al092022_20220926_DIAGNOSTICS.tcst \ + -config &CONFIG_DIR;/TCStatConfig_DIAGNOSTICS \ + -out &OUTPUT_DIR;/tc_stat/DIAGNOSTICS_stat.out \ + -v 2 + + + &OUTPUT_DIR;/tc_stat/DIAGNOSTICS_stat.out + + + diff --git a/src/basic/vx_config/config.tab.cc b/src/basic/vx_config/config.tab.cc index 5b67009f59..43fe74766f 100644 --- a/src/basic/vx_config/config.tab.cc +++ b/src/basic/vx_config/config.tab.cc @@ -1389,7 +1389,7 @@ yyparse (void) YYDPRINTF ((stderr, "Starting parse\n")); - yystate = 8; + yystate = 0; yyerrstatus = 0; yynerrs = 0; yychar = YYEMPTY; /* Cause a token to be read. */ diff --git a/src/basic/vx_config/config_constants.h b/src/basic/vx_config/config_constants.h index 5cb03579d4..1141e2c75e 100644 --- a/src/basic/vx_config/config_constants.h +++ b/src/basic/vx_config/config_constants.h @@ -85,6 +85,19 @@ enum TrackType { //////////////////////////////////////////////////////////////////////// +// +// Enumeration for tropical cyclone diagnostic types +// + +enum DiagType { + DiagType_None, // Default + TCDiagType, // Tropical Cyclone Diagnostics + LSDiagRTType, // Realtime Large Scale Diagnostics + LSDiagDevType // Development Large Scale Diagnostics +}; + +//////////////////////////////////////////////////////////////////////// + // // Enumeration for 12-hour interpolation logic // @@ -637,6 +650,7 @@ static const char conf_key_trunc_factor[] = "gaussian_trunc_factor"; static const char conf_key_eclv_points[] = "eclv_points"; static const char conf_key_var_name_map[] = "var_name_map"; static const char conf_key_metadata_map[] = "metadata_map"; +static const char conf_key_obs_to_qc_map[] = "obs_to_qc_map"; static const char conf_key_missing_thresh[] = "missing_thresh"; static const char conf_key_control_id[] = "control_id"; static const char conf_key_ens_member_ids[] = "ens_member_ids"; @@ -884,6 +898,7 @@ static const char conf_key_do_polylines_flag [] = "do_polylines"; // PB2NC specific parameter key names // +static const char conf_key_datetime[] = "datetime"; static const char conf_key_station_id[] = "station_id"; static const char conf_key_elevation_range[] = "elevation_range"; static const char conf_key_pb_report_type[] = "pb_report_type"; @@ -1066,6 +1081,9 @@ static const char conf_key_dland_file[] = "dland_file"; static const char conf_key_basin_file[] = "basin_file"; static const char conf_key_track_watch_warn[] = "track_watch_warn"; static const char conf_key_watch_warn[] = "watch_warn"; +static const char conf_key_diag_name[] = "diag_name"; +static const char conf_key_diag_convert_map[] = "diag_convert_map"; +static const char conf_key_source[] = "source"; static const char conf_key_basin_map[] = "basin_map"; static const char conf_key_time_offset[] = "time_offset"; static const char conf_key_amodel[] = "amodel"; @@ -1082,6 +1100,10 @@ static const char conf_key_init_str_name[] = "init_str_name"; static const char conf_key_init_str_val[] = "init_str_val"; static const char conf_key_init_str_exc_name[] = "init_str_exc_name"; static const char conf_key_init_str_exc_val[] = "init_str_exc_val"; +static const char conf_key_diag_thresh_name[] = "diag_thresh_name"; +static const char conf_key_diag_thresh_val[] = "diag_thresh_val"; +static const char conf_key_init_diag_thresh_name[] = "init_diag_thresh_name"; +static const char conf_key_init_diag_thresh_val[] = "init_diag_thresh_val"; static const char conf_key_water_only[] = "water_only"; static const char conf_key_rirw_track[] = "rirw.track"; static const char conf_key_rirw_time_adeck[] = "rirw.adeck.time"; diff --git a/src/basic/vx_config/config_util.cc b/src/basic/vx_config/config_util.cc index 5c1778bf39..80dee7748f 100644 --- a/src/basic/vx_config/config_util.cc +++ b/src/basic/vx_config/config_util.cc @@ -977,7 +977,7 @@ TimeSummaryInfo parse_conf_time_summary(Dictionary *dict) { void parse_add_conf_key_value_map( Dictionary *dict, const char *conf_key_map_name, map *m) { - Dictionary *msg_typ_dict = (Dictionary *) 0; + Dictionary *map_dict = (Dictionary *) 0; ConcatString key, val; int i; @@ -988,14 +988,14 @@ void parse_add_conf_key_value_map( } // Conf: map_name: message_type_map, obs)var_map, etc - msg_typ_dict = dict->lookup_array(conf_key_map_name); + map_dict = dict->lookup_array(conf_key_map_name); // Loop through the array entries - for(i=0; in_entries(); i++) { + for(i=0; in_entries(); i++) { // Lookup the key and value - key = (*msg_typ_dict)[i]->dict_value()->lookup_string(conf_key_key); - val = (*msg_typ_dict)[i]->dict_value()->lookup_string(conf_key_val); + key = (*map_dict)[i]->dict_value()->lookup_string(conf_key_key); + val = (*map_dict)[i]->dict_value()->lookup_string(conf_key_val); if(m->count(key) >= 1) { (*m)[key] = val; @@ -1011,10 +1011,9 @@ void parse_add_conf_key_value_map( /////////////////////////////////////////////////////////////////////////////// - map parse_conf_key_value_map( Dictionary *dict, const char *conf_key_map_name, const char *caller) { - Dictionary *msg_typ_dict = (Dictionary *) 0; + Dictionary *map_dict = (Dictionary *) 0; map m; ConcatString key, val; int i; @@ -1026,14 +1025,14 @@ map parse_conf_key_value_map( } // Conf: map_name: message_type_map, obs_var_map, etc - msg_typ_dict = dict->lookup_array(conf_key_map_name); + map_dict = dict->lookup_array(conf_key_map_name); // Loop through the array entries - for(i=0; in_entries(); i++) { + for(i=0; in_entries(); i++) { // Lookup the key and value - key = (*msg_typ_dict)[i]->dict_value()->lookup_string(conf_key_key); - val = (*msg_typ_dict)[i]->dict_value()->lookup_string(conf_key_val); + key = (*map_dict)[i]->dict_value()->lookup_string(conf_key_key); + val = (*map_dict)[i]->dict_value()->lookup_string(conf_key_val); if(m.count(key) >= 1) { mlog << Warning << "\n" << method_name @@ -1100,6 +1099,69 @@ map parse_conf_obs_name_map(Dictionary *dict) { /////////////////////////////////////////////////////////////////////////////// +map parse_conf_obs_to_qc_map(Dictionary *dict) { + const char *method_name = "parse_conf_obs_to_qc_map() -> "; + return parse_conf_key_values_map(dict, conf_key_obs_to_qc_map, method_name); +} + +/////////////////////////////////////////////////////////////////////////////// + +map parse_conf_key_convert_map( + Dictionary *dict, const char *conf_key_map_name, const char *caller) { + Dictionary *map_dict = (Dictionary *) 0; + int i, j; + StringArray sa; + ConcatString key; + UserFunc_1Arg fx; + map m; + const char *method_name = (0 != caller) ? caller : "parse_conf_key_convert_map() -> "; + + if(!dict) { + mlog << Error << "\n" << method_name << "empty dictionary!\n\n"; + exit(1); + } + + // Conf: tcdiag_convert_map, lsdiag_convert_map, etc + map_dict = dict->lookup_array(conf_key_map_name); + + // Loop through the array entries + for(i=0; in_entries(); i++) { + + // Lookup the key and convert function + sa = (*map_dict)[i]->dict_value()->lookup_string_array(conf_key_key); + fx.clear(); + fx.set((*map_dict)[i]->dict_value()->lookup(conf_key_convert)); + + // Check the function + if(!fx.is_set()) { + mlog << Error << "\n" << method_name + << "lookup for \"" << conf_key_convert << "\" failed in the \"" + << conf_key_map_name << "\" map!\n\n"; + exit(1); + } + + // Add map entry for each string + for(j=0; j= 1) { + mlog << Warning << "\n" << method_name + << "found multiple entries for key \"" << key << "\" in the \"" + << conf_key_map_name << "\" map!\n\n"; + } + + // Add entry to the map + m.insert(pair(key,fx)); + + } // end for j + } // end for i + + return(m); +} + +/////////////////////////////////////////////////////////////////////////////// + void BootInfo::clear() { interval = BootIntervalType_None; rep_prop = bad_data_double; @@ -2685,6 +2747,45 @@ ConcatString tracktype_to_string(TrackType type) { /////////////////////////////////////////////////////////////////////////////// +DiagType string_to_diagtype(const char *s) { + DiagType t = DiagType_None; + + // Convert string to enumerated DiagType + if(strcasecmp(s, conf_val_none) == 0) t = DiagType_None; + else if(strcasecmp(s, "TCDIAG") == 0) t = TCDiagType; + else if(strcasecmp(s, "LSDIAG_RT") == 0) t = LSDiagRTType; + else if(strcasecmp(s, "LSDIAG_DEV") == 0) t = LSDiagDevType; + else { + mlog << Error << "\nstring_to_diagtype() -> " + << "Unexpected DiagType string \"" << s << "\".\n\n"; + exit(1); + } + + return(t); +} + +/////////////////////////////////////////////////////////////////////////////// + +ConcatString diagtype_to_string(DiagType type) { + ConcatString s; + + // Convert enumerated DiagType to string + switch(type) { + case(DiagType_None): s = conf_val_none; break; + case(TCDiagType): s = "TCDIAG"; break; + case(LSDiagRTType): s = "LSDIAG_RT"; break; + case(LSDiagDevType): s = "LSDIAG_DEV"; break; + default: + mlog << Error << "\ndiagtype_to_string() -> " + << "Unexpected DiagType value of " << type << ".\n\n"; + exit(1); + } + + return(s); +} + +/////////////////////////////////////////////////////////////////////////////// + Interp12Type int_to_interp12type(int v) { Interp12Type t = Interp12Type_None; diff --git a/src/basic/vx_config/config_util.h b/src/basic/vx_config/config_util.h index 45de85dd0f..9debcd2d49 100644 --- a/src/basic/vx_config/config_util.h +++ b/src/basic/vx_config/config_util.h @@ -52,9 +52,9 @@ extern NumArray parse_conf_eclv_points(Dictionary *dict); extern ClimoCDFInfo parse_conf_climo_cdf(Dictionary *dict); extern TimeSummaryInfo parse_conf_time_summary(Dictionary *dict); extern std::map parse_conf_key_value_map( - Dictionary *dict, const char *conf_key_map_name, const char *caller=0); + Dictionary *dict, const char *conf_key_map_name, const char *caller=0); extern void parse_add_conf_key_value_map( - Dictionary *dict, const char *conf_key_map_name, std::map *m); + Dictionary *dict, const char *conf_key_map_name, std::map *m); extern std::map parse_conf_message_type_map(Dictionary *dict); extern std::map @@ -62,6 +62,11 @@ extern std::map extern std::map parse_conf_metadata_map(Dictionary *dict); extern std::map parse_conf_obs_name_map(Dictionary *dict); +extern std::map + parse_conf_obs_to_qc_map(Dictionary *dict); +extern std::map + parse_conf_key_convert_map( + Dictionary *dict, const char *conf_key_map_name, const char *caller=0); extern BootInfo parse_conf_boot(Dictionary *dict); extern RegridInfo parse_conf_regrid(Dictionary *dict, bool error_out = default_dictionary_error_out); extern InterpInfo parse_conf_interp(Dictionary *dict, const char *); @@ -109,6 +114,9 @@ extern TrackType int_to_tracktype(int); extern TrackType string_to_tracktype(const char *); extern ConcatString tracktype_to_string(TrackType); +extern DiagType string_to_diagtype(const char *); +extern ConcatString diagtype_to_string(DiagType); + extern Interp12Type int_to_interp12type(int); extern Interp12Type string_to_interp12type(const char *); extern ConcatString interp12type_to_string(Interp12Type); diff --git a/src/basic/vx_util/data_line.cc b/src/basic/vx_util/data_line.cc index f7692fdeb0..dade54254c 100644 --- a/src/basic/vx_util/data_line.cc +++ b/src/basic/vx_util/data_line.cc @@ -815,6 +815,34 @@ return ( 1 ); } +//////////////////////////////////////////////////////////////////////// +// +// Read the next line but do not move the read pointer +// +//////////////////////////////////////////////////////////////////////// + + +int LineDataFile::peek_line(DataLine & a) + +{ + +int status = 0; + +if ( ok() ) { + + int cur_pos = in->tellg(); + + status = a.read_line(this); + + in->seekg(cur_pos); + +} + +return ( 1 ); + +} + + //////////////////////////////////////////////////////////////////////// diff --git a/src/basic/vx_util/data_line.h b/src/basic/vx_util/data_line.h index 4608e8a2bd..fd509d3906 100644 --- a/src/basic/vx_util/data_line.h +++ b/src/basic/vx_util/data_line.h @@ -190,6 +190,8 @@ class LineDataFile { virtual int operator>>(DataLine &); + int peek_line(DataLine &); + int read_fwf_line(DataLine &, const int *wdth, int n_wdth); const char * filename() const; diff --git a/src/basic/vx_util/stat_column_defs.h b/src/basic/vx_util/stat_column_defs.h index 9fca87257c..0101825690 100644 --- a/src/basic/vx_util/stat_column_defs.h +++ b/src/basic/vx_util/stat_column_defs.h @@ -262,12 +262,14 @@ static const char * isc_columns [] = { }; static const char * ecnt_columns [] = { - "TOTAL", "N_ENS", "CRPS", - "CRPSS", "IGN", "ME", - "RMSE", "SPREAD", "ME_OERR", - "RMSE_OERR", "SPREAD_OERR", "SPREAD_PLUS_OERR", - "CRPSCL", "CRPS_EMP", "CRPSCL_EMP", - "CRPSS_EMP", "CRPS_EMP_FAIR" + "TOTAL", "N_ENS", "CRPS", + "CRPSS", "IGN", "ME", + "RMSE", "SPREAD", "ME_OERR", + "RMSE_OERR", "SPREAD_OERR", "SPREAD_PLUS_OERR", + "CRPSCL", "CRPS_EMP", "CRPSCL_EMP", + "CRPSS_EMP", "CRPS_EMP_FAIR", "BIAS_RATIO", + "N_GE_OBS", "ME_GE_OBS", "N_LT_OBS", + "ME_LT_OBS" }; static const char * rps_columns [] = { diff --git a/src/libcode/vx_data2d_factory/is_netcdf_file.cc b/src/libcode/vx_data2d_factory/is_netcdf_file.cc index 3e72029514..b14fa51fa6 100644 --- a/src/libcode/vx_data2d_factory/is_netcdf_file.cc +++ b/src/libcode/vx_data2d_factory/is_netcdf_file.cc @@ -40,9 +40,6 @@ static const char netcdf_magic [] = "CDF"; static const char hdf_magic [] = "HDF"; static const int netcdf_magic_len = m_strlen(netcdf_magic); -static const string nccf_att_name = "Conventions"; -static const string nccf_att_name_l = "conventions"; -static const string nccf_att_name_U = "CONVENTIONS"; static const string nccf_att_value = "CF-"; static const string nccf_att_value2 = "CF "; static const string nccf_att_value3 = "COARDS"; @@ -96,9 +93,7 @@ bool is_nccf_file(const char * filename) NcFile *nc_file = open_ncfile(filename); if (!IS_INVALID_NC_P(nc_file)) { - bool found = get_global_att(nc_file, nccf_att_name, att_val); - if (!found) found = get_global_att(nc_file, nccf_att_name_l, att_val); - if (!found) found = get_global_att(nc_file, nccf_att_name_U, att_val); + bool found = get_cf_conventions(nc_file, att_val); // "Conventions" attrribute if (found) { status = (att_val.compare(0, nccf_att_value.length(), nccf_att_value) == 0 || diff --git a/src/libcode/vx_nc_obs/nc_obs_util.cc b/src/libcode/vx_nc_obs/nc_obs_util.cc index 73db28167a..8f5a019280 100644 --- a/src/libcode/vx_nc_obs/nc_obs_util.cc +++ b/src/libcode/vx_nc_obs/nc_obs_util.cc @@ -480,6 +480,9 @@ void NetcdfObsVars::read_dims_vars(NcFile *f_in) { obs_dim = get_nc_dim(f_in, nc_dim_nobs); // Observation array length hdr_dim = get_nc_dim(f_in, nc_dim_nhdr); // Header array length + use_var_id = false; + get_global_att(f_in, nc_att_use_var_id, use_var_id); + // Get netCDF header variables hdr_typ_var = get_var(f_in, nc_var_hdr_typ); // Message type (String or int) hdr_sid_var = get_var(f_in, nc_var_hdr_sid); // Station ID (String or int) @@ -503,10 +506,14 @@ void NetcdfObsVars::read_dims_vars(NcFile *f_in) { obs_arr_var = get_var(f_in, nc_var_obs_arr); } else { obs_hid_var = ncVar; // Obs. header id array - ncVar = get_var(f_in, nc_var_obs_gc); - if (!IS_INVALID_NC(ncVar)) obs_gc_var = ncVar; // Obs. grib code array - ncVar = get_var(f_in, nc_var_obs_vid); - if (!IS_INVALID_NC(ncVar)) obs_vid_var = ncVar; // Obs. variable id array + if (use_var_id) { + ncVar = get_var(f_in, nc_var_obs_vid); + if (!IS_INVALID_NC(ncVar)) obs_vid_var = ncVar; // Obs. variable id array + } + else { + ncVar = get_var(f_in, nc_var_obs_gc); + if (!IS_INVALID_NC(ncVar)) obs_gc_var = ncVar; // Obs. grib code array + } obs_lvl_var = get_var(f_in, nc_var_obs_lvl); // Obs. pressure level array obs_hgt_var = get_var(f_in, nc_var_obs_hgt); // Obs. highth array obs_val_var = get_var(f_in, nc_var_obs_val); // Obs. value array @@ -516,12 +523,14 @@ void NetcdfObsVars::read_dims_vars(NcFile *f_in) { ncVar = get_var(f_in, nc_var_obs_qty_tbl); if (!IS_INVALID_NC(ncVar)) obs_qty_tbl_var = ncVar; - ncVar = get_var(f_in, nc_var_obs_var); - if (!IS_INVALID_NC(ncVar)) obs_var = ncVar; - ncVar = get_var(f_in, nc_var_unit); - if (!IS_INVALID_NC(ncVar)) unit_var = ncVar; - ncVar = get_var(f_in, nc_var_desc); - if (!IS_INVALID_NC(ncVar)) desc_var = ncVar; + if (use_var_id) { + ncVar = get_var(f_in, nc_var_obs_var); + if (!IS_INVALID_NC(ncVar)) obs_var = ncVar; + ncVar = get_var(f_in, nc_var_unit); + if (!IS_INVALID_NC(ncVar)) unit_var = ncVar; + ncVar = get_var(f_in, nc_var_desc); + if (!IS_INVALID_NC(ncVar)) desc_var = ncVar; + } // PrepBufr only headers ncVar = get_var(f_in, nc_var_hdr_prpt_typ); @@ -531,12 +540,6 @@ void NetcdfObsVars::read_dims_vars(NcFile *f_in) { ncVar = get_var(f_in, nc_var_hdr_inst_typ); if (!IS_INVALID_NC(ncVar)) hdr_inst_typ_var = ncVar; - bool _use_var_id = false; - if (!get_global_att(f_in, nc_att_use_var_id, _use_var_id)) { - _use_var_id = IS_VALID_NC(obs_var); - } - - use_var_id = _use_var_id; } //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_nc_util/nc_utils.cc b/src/libcode/vx_nc_util/nc_utils.cc index 676fa75508..b312655396 100644 --- a/src/libcode/vx_nc_util/nc_utils.cc +++ b/src/libcode/vx_nc_util/nc_utils.cc @@ -165,9 +165,10 @@ bool get_att_value_chars(const NcAtt *att, ConcatString &value) { static const char *method_name = "get_att_value_chars(NcAtt) -> "; if (IS_VALID_NC_P(att)) { nc_type attType = GET_NC_TYPE_ID_P(att); - if (attType == NC_CHAR || attType == NC_STRING) { + if (attType == NC_CHAR) { try { - string att_value; + char att_value[tmp_buf_size]; + memset(att_value, 0, tmp_buf_size); att->getValues(att_value); value = att_value; } @@ -181,6 +182,33 @@ bool get_att_value_chars(const NcAtt *att, ConcatString &value) { << "Please check the encoding of the "<< GET_NC_NAME_P(att) << " attribute.\n\n"; } } + else if (attType == NC_STRING) { + try { + string att_value; + att->getValues(att_value); + value = att_value; + } + catch (exceptions::NcChar ex) { + int num_elements_sub = 8096; + int num_elements = att->getAttLength();; + char *att_value[num_elements]; + for (int i = 0; i < num_elements; i++ ) { + att_value[i] = (char*) calloc(num_elements_sub, sizeof(char)); + } + try { + att->getValues(att_value); + value = att_value[0]; + } + catch (exceptions::NcException ex) { + mlog << Warning << "\n" << method_name + << "Exception: " << ex.what() << "\n" + << "Fail to read " << GET_NC_NAME_P(att) << " attribute (" + << GET_NC_TYPE_NAME_P(att) << " type).\n" + << "Please check the encoding of the "<< GET_NC_NAME_P(att) << " attribute.\n\n"; + } + for (int i = 0; i < num_elements; i++ ) delete att_value[i]; + } + } else { // MET-788: to handle a custom modified NetCDF mlog << Error << "\n" << method_name << "Please convert data type of \"" << GET_NC_NAME_P(att) @@ -336,6 +364,23 @@ bool get_att_no_leap_year(const NcVar *var) { return no_leap_year; } +//////////////////////////////////////////////////////////////////////// + +bool get_cf_conventions(const netCDF::NcFile *nc, ConcatString& conventions_value) { + bool has_attr = false; + multimap::iterator it_att; + multimap map_attrs = nc->getAtts(); + for (it_att = map_attrs.begin(); it_att != map_attrs.end(); it_att++) { + if (to_lower(it_att->first) == to_lower(cf_att_name)) { + NcGroupAtt cf_att = it_att->second; + if (IS_VALID_NC(cf_att)) + has_attr = get_att_value_chars(&cf_att, conventions_value); + } + } + return has_attr; +} + + //////////////////////////////////////////////////////////////////////// ConcatString get_log_msg_for_att(const NcVarAtt *att) { @@ -354,7 +399,7 @@ ConcatString get_log_msg_for_att(const NcVarAtt *att) { //////////////////////////////////////////////////////////////////////// ConcatString get_log_msg_for_att(const NcVarAtt *att, string var_name, - const ConcatString att_name) { + const ConcatString att_name) { ConcatString log_msg; log_msg << "can't read attribute" << " \"" << ((att_name.length() > 0) ? att_name.c_str() : GET_SAFE_NC_NAME_P(att)) @@ -445,7 +490,7 @@ NcGroupAtt *get_nc_att(const NcFile * nc, const ConcatString &att_name, bool exi //////////////////////////////////////////////////////////////////////// bool get_nc_att_value(const NcVar *var, const ConcatString &att_name, - ConcatString &att_val, bool exit_on_error) { + ConcatString &att_val, int grp_id, bool exit_on_error) { bool status = false; NcVarAtt *att = (NcVarAtt *) 0; @@ -833,28 +878,49 @@ void add_att(NcVar *var, const string &att_name, const double att_val) { //////////////////////////////////////////////////////////////////////// -int get_var_names(NcFile *nc, StringArray *varNames) { +int get_var_names(NcFile *nc, StringArray *var_names) { - int i, varCount; NcVar var; + int i = 0; + int var_count = nc->getVarCount(); - varCount = nc->getVarCount(); - - i = 0; - multimap::iterator itVar; multimap mapVar = GET_NC_VARS_P(nc); - for (itVar = mapVar.begin(); itVar != mapVar.end(); ++itVar) { - var = (*itVar).second; - varNames->add(var.getName()); + for (multimap::iterator it_var = mapVar.begin(); + it_var != mapVar.end(); ++it_var) { + var = (*it_var).second; + var_names->add(var.getName()); i++; } - if (i != varCount) { + if (i != var_count) { mlog << Error << "\n\tget_var_names() -> " - << "does not match array, allocated " << varCount << " but assigned " + << "does not match array, allocated " << var_count << " but assigned " << i << ".\n\n"; } - return(varCount); + return(var_count); +} + +//////////////////////////////////////////////////////////////////////// + +int get_var_names(NcFile *nc, StringArray *var_names, StringArray &group_names) { + + NcVar var; + NcGroup nc_group; + int var_count = 0; + multimap var_map; + multimap::iterator it_var; + + for (int idx=0; idxadd(it_var->first); + var_count++; + } + } + } + return(var_count); } //////////////////////////////////////////////////////////////////////// @@ -1657,6 +1723,14 @@ bool get_nc_data(NcVar *var, char *data) { //////////////////////////////////////////////////////////////////////// +bool get_nc_data(NcVar *var, char **data) { + bool return_status = get_nc_data_t(var, data); + + return(return_status); +} + +//////////////////////////////////////////////////////////////////////// + bool get_nc_data(NcVar *var, uchar *data) { bool return_status = false; int data_type = GET_NC_TYPE_ID_P(var); @@ -2114,6 +2188,18 @@ bool args_ok(const LongArray & a) { } //////////////////////////////////////////////////////////////////////// +// Continue even though not exists + +NcGroup get_nc_group(NcFile *nc, const char *group_name) { + NcGroup nc_group; + multimap group_map = nc->getGroups(); + multimap::iterator it = group_map.find(group_name); + if (it != group_map.end()) nc_group = it->second; + return nc_group; +} + +//////////////////////////////////////////////////////////////////////// +// Exit if exists but invalid NcVar get_var(NcFile *nc, const char *var_name) { string new_var_name = var_name; @@ -2123,23 +2209,57 @@ NcVar get_var(NcFile *nc, const char *var_name) { // Retrieve the variable from the NetCDF file. // NcVar var; - multimap varMap = GET_NC_VARS_P(nc); - multimap::iterator it = varMap.find(new_var_name); - if (it != varMap.end()) { - NcVar tmpVar = it->second; - if(IS_INVALID_NC(tmpVar)) { - mlog << Error << "\nget_var() -> " + multimap var_map = GET_NC_VARS_P(nc); + multimap::iterator it = var_map.find(new_var_name); + if (it != var_map.end()) { + var = it->second; + if(IS_INVALID_NC(var)) { + mlog << Error << "\nget_var(var_name) -> " << "can't read \"" << new_var_name << "\" variable.\n\n"; exit(1); } + } - var = tmpVar; + return(var); +} + +//////////////////////////////////////////////////////////////////////// +// Exit if exists but invalid + +NcVar get_var(NcFile *nc, const char *var_name, const char *group_name) { + string nc_var_name; + string new_var_name = var_name; + patch_nc_name(&new_var_name); + + // + // Retrieve the variable from the NetCDF file. + // + NcVar var; + multimap var_map; + NcGroup nc_group = get_nc_group(nc, group_name); + if (IS_VALID_NC(nc_group)) { + nc_var_name = new_var_name; + var_map = GET_NC_VARS(nc_group); + } + else { // This is for IODA data format 1.0 + nc_var_name = new_var_name + "@" + group_name; + var_map = GET_NC_VARS_P(nc); + } + multimap::iterator it = var_map.find(nc_var_name); + if (it != var_map.end()) { + var = it->second; + if(IS_INVALID_NC(var)) { + mlog << Error << "\nget_var(var_name, group_name) -> " + << "can't read \"" << new_var_name << "\" variable.\n\n"; + exit(1); + } } return(var); } //////////////////////////////////////////////////////////////////////// +// Continue even though not exists NcVar get_nc_var(NcFile *nc, const char *var_name, bool log_as_error) { string new_var_name = var_name; @@ -2151,7 +2271,7 @@ NcVar get_nc_var(NcFile *nc, const char *var_name, bool log_as_error) { NcVar var = nc->getVar(new_var_name); if(IS_INVALID_NC(var)) { ConcatString log_message; - log_message << "\nget_nc_var(NcFile) --> The variable \"" + log_message << "\nget_nc_var(var_name) --> The variable \"" << new_var_name << "\" does not exist!\n\n"; if (log_as_error) mlog << Error << log_message; @@ -2162,6 +2282,60 @@ NcVar get_nc_var(NcFile *nc, const char *var_name, bool log_as_error) { return(var); } +//////////////////////////////////////////////////////////////////////// +// Continue even though not exists + +NcVar get_nc_var(NcFile *nc, const ConcatString &var_name, bool log_as_error) { + return get_nc_var(nc, var_name.c_str(), log_as_error); +} + +//////////////////////////////////////////////////////////////////////// +// Continue even though not exists + +NcVar get_nc_var(NcFile *nc, const char *var_name, const char *group_name, + bool log_as_error) { + string nc_var_name; + string new_var_name = var_name; + patch_nc_name(&new_var_name); + + // + // Retrieve the variable from the NetCDF file. + // + NcVar var; + multimap var_map; + NcGroup nc_group = get_nc_group(nc, group_name); + if (IS_VALID_NC(nc_group)) { + nc_var_name = new_var_name; + var_map = GET_NC_VARS(nc_group); + } + else { // This is for IODA data format 1.0 + nc_var_name = new_var_name + "@" + group_name; + var_map = GET_NC_VARS_P(nc); + } + multimap::iterator it = var_map.find(nc_var_name); + if (it != var_map.end()) var = it->second; + + if(IS_INVALID_NC(var)) { + ConcatString log_message; + log_message << "\nget_nc_var(var_name, group_name) --> The variable \"" + << nc_var_name << "\" does not exist!\n\n"; + if (log_as_error) + mlog << Error << log_message; + else + mlog << Warning << log_message; + } + + return(var); +} + +//////////////////////////////////////////////////////////////////////// +// Continue even though not exists + +NcVar get_nc_var(NcFile *nc, const ConcatString &var_name, const char *group_name, + bool log_as_error) { + return get_nc_var(nc, var_name.c_str(), group_name, log_as_error); +} + //////////////////////////////////////////////////////////////////////// void copy_nc_att_byte(NcFile *nc_to, NcGroupAtt *from_att) { @@ -2506,8 +2680,8 @@ void copy_nc_atts(NcFile *nc_from, NcFile *nc_to, const bool all_attrs) { for (multimap::iterator itr = ncAttMap.begin(); itr != ncAttMap.end(); ++itr) { if (all_attrs || - ( (itr->first != "Conventions") - && (itr->first != "missing_value") ) ) { + ( (itr->first != cf_att_name) + && (itr->first != missing_value_att_name) ) ) { NcGroupAtt *from_att = &(itr->second); int dataType = GET_NC_TYPE_ID_P(from_att); switch (dataType) { @@ -2654,6 +2828,7 @@ void copy_nc_data_short(NcVar *var_from, NcVar *var_to, int data_size) { delete[] data; } +//////////////////////////////////////////////////////////////////////// void copy_nc_var_data(NcVar *var_from, NcVar *var_to) { const string method_name = "copy_nc_var_data()"; @@ -2698,13 +2873,60 @@ void copy_nc_var_dims(NcVar *var_from, NcVar *var_to) { //////////////////////////////////////////////////////////////////////// -bool has_var(NcFile *nc, const char * var_name) { - NcVar v = get_var(nc, var_name); +bool has_nc_group(NcFile *nc, const char *group_name) { + multimap group_map = nc->getGroups(); + multimap::iterator it = group_map.find(group_name); + return (it != group_map.end()); +} + +//////////////////////////////////////////////////////////////////////// + +bool has_var(NcFile *nc, const char *var_name) { + string new_var_name = var_name; + patch_nc_name(&new_var_name); + NcVar v = get_var(nc, new_var_name.c_str()); return IS_VALID_NC(v); } //////////////////////////////////////////////////////////////////////// +bool has_var(NcFile *nc, const ConcatString var_name) { + return has_var(nc, var_name.c_str()); +} + +//////////////////////////////////////////////////////////////////////// + +bool has_var(NcFile *nc, const char *var_name, const char *group_name) { + string nc_var_name; + string new_var_name = var_name; + patch_nc_name(&new_var_name); + + // + // Retrieve the variable from the NetCDF file. + // + multimap var_map; + NcGroup nc_group = get_nc_group(nc, group_name); + if (IS_VALID_NC(nc_group)) { + nc_var_name = new_var_name; + var_map = GET_NC_VARS(nc_group); + } + else { // This is for IODA data format 1.0 + nc_var_name = new_var_name + "@" + group_name; + var_map = GET_NC_VARS_P(nc); + } + multimap::iterator it = var_map.find(nc_var_name); + + return (it != var_map.end()); +} + +//////////////////////////////////////////////////////////////////////// + +bool has_var(NcFile *nc, const ConcatString var_name, const char *group_name) { + return has_var(nc, var_name.c_str(), group_name); +} + +//////////////////////////////////////////////////////////////////////// + NcVar add_var(NcFile *nc, const string &var_name, const NcType ncType, const int deflate_level) { vector ncDimVector; string new_var_name = var_name; @@ -3009,26 +3231,24 @@ vector get_dims(const NcVar *var, int *dim_count) { //////////////////////////////////////////////////////////////////////// bool is_nc_name_lat(const ConcatString name) { - bool is_latitude = (name == "lat" || name == "LAT" - || name == "Lat" || name == "Latitude" - || name == "latitude" || name == "LATITUDE"); + ConcatString name_l = to_lower(name); + bool is_latitude = (name_l == "lat" || name_l == "latitude"); return is_latitude; } //////////////////////////////////////////////////////////////////////// bool is_nc_name_lon(const ConcatString name) { - bool is_longitude = (name == "lon" || name == "LON" - || name == "Lon" || name == "Longitude" - || name == "longitude" || name == "LONGITUDE"); + ConcatString name_l = to_lower(name); + bool is_longitude = (name_l == "lon" || name_l == "longitude"); return is_longitude; } //////////////////////////////////////////////////////////////////////// bool is_nc_name_time(const ConcatString name) { - bool is_time = (name == "t" || name == "time" || name == "Time" || name == "TIME" - || name == "datetime" || name == "Datetime" || name == "DATETIME"); + ConcatString name_l = to_lower(name); + bool is_time = (name == "t" || name_l == "time" || name_l == "datetime"); return is_time; } @@ -3061,26 +3281,26 @@ NcVar get_nc_var_lat(const NcFile *nc) { multimap mapVar = GET_NC_VARS_P(nc); static const char *method_name = "get_nc_var_lat() "; - for (multimap::iterator itVar = mapVar.begin(); - itVar != mapVar.end(); ++itVar) { - ConcatString name = (*itVar).first; + for (multimap::iterator it_var = mapVar.begin(); + it_var != mapVar.end(); ++it_var) { + ConcatString name = (*it_var).first; //if (is_nc_name_lat(name)) found = true; - if (get_var_standard_name(&(*itVar).second, name)) { + if (get_var_standard_name(&(*it_var).second, name)) { if (is_nc_name_lat(name)) found = true; } - if (!found && get_var_units(&(*itVar).second, name)) { + if (!found && get_var_units(&(*it_var).second, name)) { if (is_nc_unit_latitude(name.c_str())) { - if (get_nc_att_value(&(*itVar).second, axis_att_name, name)) { + if (get_nc_att_value(&(*it_var).second, axis_att_name, name)) { if (is_nc_attr_lat(name)) found = true; } - else if (get_nc_att_value(&(*itVar).second, + else if (get_nc_att_value(&(*it_var).second, coordinate_axis_type_att_name, name)) { if (is_nc_attr_lat(name)) found = true; } } } if (found) { - var = (*itVar).second; + var = (*it_var).second; break; } } @@ -3103,26 +3323,26 @@ NcVar get_nc_var_lon(const NcFile *nc) { multimap mapVar = GET_NC_VARS_P(nc); static const char *method_name = "get_nc_var_lon() "; - for (multimap::iterator itVar = mapVar.begin(); - itVar != mapVar.end(); ++itVar) { - ConcatString name = (*itVar).first; + for (multimap::iterator it_var = mapVar.begin(); + it_var != mapVar.end(); ++it_var) { + ConcatString name = (*it_var).first; //if (is_nc_name_lon(name)) found = true; - if (get_var_standard_name(&(*itVar).second, name)) { + if (get_var_standard_name(&(*it_var).second, name)) { if (is_nc_name_lon(name)) found = true; } - if (!found && get_var_units(&(*itVar).second, name)) { + if (!found && get_var_units(&(*it_var).second, name)) { if (is_nc_unit_longitude(name.c_str())) { - if (get_nc_att_value(&(*itVar).second, axis_att_name, name)) { + if (get_nc_att_value(&(*it_var).second, axis_att_name, name)) { if (is_nc_attr_lon(name)) found = true; } - else if (get_nc_att_value(&(*itVar).second, + else if (get_nc_att_value(&(*it_var).second, coordinate_axis_type_att_name, name)) { if (is_nc_attr_lon(name)) found = true; } } } if (found) { - var = (*itVar).second; + var = (*it_var).second; break; } } @@ -3145,28 +3365,28 @@ NcVar get_nc_var_time(const NcFile *nc) { multimap mapVar = GET_NC_VARS_P(nc); static const char *method_name = "get_nc_var_time() "; - for (multimap::iterator itVar = mapVar.begin(); - itVar != mapVar.end(); ++itVar) { - ConcatString name = (*itVar).first; + for (multimap::iterator it_var = mapVar.begin(); + it_var != mapVar.end(); ++it_var) { + ConcatString name = (*it_var).first; //if (is_nc_name_time(name)) found = true; - if (get_var_standard_name(&(*itVar).second, name)) { + if (get_var_standard_name(&(*it_var).second, name)) { if (is_nc_name_time(name)) found = true; mlog << Debug(7) << method_name << "checked variable \"" << name << "\" is_time: " << found << "\n"; } - if (!found && get_var_units(&(*itVar).second, name)) { + if (!found && get_var_units(&(*it_var).second, name)) { if (is_nc_unit_time(name.c_str())) { - if (get_nc_att_value(&(*itVar).second, axis_att_name, name)) { + if (get_nc_att_value(&(*it_var).second, axis_att_name, name)) { if (is_nc_attr_time(name)) found = true; } - else if (get_nc_att_value(&(*itVar).second, + else if (get_nc_att_value(&(*it_var).second, coordinate_axis_type_att_name, name)) { if (is_nc_attr_time(name)) found = true; } } } if (found) { - var = (*itVar).second; + var = (*it_var).second; break; } } diff --git a/src/libcode/vx_nc_util/nc_utils.h b/src/libcode/vx_nc_util/nc_utils.h index 199173e263..adb518d402 100644 --- a/src/libcode/vx_nc_util/nc_utils.h +++ b/src/libcode/vx_nc_util/nc_utils.h @@ -133,6 +133,7 @@ static const std::string axis_att_name = "axis"; static const std::string bounds_att_name = "bounds"; static const std::string coordinates_att_name = "coordinates"; static const std::string coordinate_axis_type_att_name = "_CoordinateAxisType"; +static const std::string cf_att_name = "Conventions"; static const std::string description_att_name = "description"; static const std::string fill_value_att_name = "_FillValue"; static const std::string grid_mapping_att_name = "grid_mapping"; @@ -179,14 +180,17 @@ extern long long get_att_value_llong (const netCDF::NcFile *, const ConcatString extern double get_att_value_double(const netCDF::NcFile *, const ConcatString& ); extern bool get_att_no_leap_year(const netCDF::NcVar *); -extern netCDF::NcVarAtt *get_nc_att(const netCDF::NcVar *, const ConcatString &, bool exit_on_error = false); -extern netCDF::NcGroupAtt *get_nc_att(const netCDF::NcFile *, const ConcatString &, bool exit_on_error = false); +extern bool get_cf_conventions(const netCDF::NcFile *, ConcatString&); + +extern netCDF::NcVarAtt *get_nc_att(const netCDF::NcVar *, const ConcatString &, bool exit_on_error = false); +extern netCDF::NcGroupAtt *get_nc_att(const netCDF::NcFile *, const ConcatString &, bool exit_on_error = false); extern bool get_nc_att_value(const netCDF::NcVarAtt *, std::string &); extern bool get_nc_att_value(const netCDF::NcVarAtt *, int &, bool exit_on_error = true); extern bool get_nc_att_value(const netCDF::NcVarAtt *, float &, bool exit_on_error = true); extern bool get_nc_att_value(const netCDF::NcVarAtt *, double &, bool exit_on_error = true); -extern bool get_nc_att_value(const netCDF::NcVar *, const ConcatString &, ConcatString &, bool exit_on_error = false); +extern bool get_nc_att_value(const netCDF::NcVar *, const ConcatString &, ConcatString &, + int nc_id=0, bool exit_on_error = false); extern bool get_nc_att_value(const netCDF::NcVar *, const ConcatString &, int &, bool exit_on_error = false); extern bool get_nc_att_value(const netCDF::NcVar *, const ConcatString &, float &, bool exit_on_error = false); extern bool get_nc_att_value(const netCDF::NcVar *, const ConcatString &, double &, bool exit_on_error = false); @@ -216,7 +220,8 @@ extern void add_att(netCDF::NcVar *, const std::string &, const int ); extern void add_att(netCDF::NcVar *, const std::string &, const float ); extern void add_att(netCDF::NcVar *, const std::string &, const double); -extern int get_var_names(netCDF::NcFile *, StringArray *varNames); +extern int get_var_names(netCDF::NcFile *, StringArray *var_names); +extern int get_var_names(netCDF::NcFile *, StringArray *var_names, StringArray &group_names); extern bool get_var_att_float (const netCDF::NcVar *, const ConcatString &, float &); extern bool get_var_att_double(const netCDF::NcVar *, const ConcatString &, double &); @@ -250,6 +255,7 @@ extern ConcatString* get_string_val(netCDF::NcVar *var, const int index, const i extern bool get_nc_data(netCDF::NcVar *, int *data); extern bool get_nc_data(netCDF::NcVar *, char *data); +extern bool get_nc_data(netCDF::NcVar *, char **data); extern bool get_nc_data(netCDF::NcVar *, uchar *data); extern bool get_nc_data(netCDF::NcVar *, float *data); extern bool get_nc_data(netCDF::NcVar *, double *data); @@ -315,17 +321,31 @@ extern bool put_nc_data_with_dims(netCDF::NcVar *, const double *data, const int extern bool put_nc_data_with_dims(netCDF::NcVar *, const double *data, const long len0, const long len1=0, const long len2=0); -extern netCDF::NcVar get_var(netCDF::NcFile *, const char * var_name); // exit if not exists -extern netCDF::NcVar get_nc_var(netCDF::NcFile *, const char * var_name, bool log_as_error=false); // continue even though not exists - -extern netCDF::NcVar *copy_nc_var(netCDF::NcFile *, netCDF::NcVar *, const int deflate_level=DEF_DEFLATE_LEVEL, const bool all_attrs=true); +extern netCDF::NcGroup get_nc_group(netCDF::NcFile *, const char *group_name); // continue even though not exists + +extern netCDF::NcVar get_var(netCDF::NcFile *, const char *var_name); // exit if exists but invalid +extern netCDF::NcVar get_var(netCDF::NcFile *, const char *var_name, + const char *group_name); // continue even though not exists +extern netCDF::NcVar get_nc_var(netCDF::NcFile *, const char *var_name, + bool log_as_error=false); // continue even though not exists +extern netCDF::NcVar get_nc_var(netCDF::NcFile *, const ConcatString &var_name, + bool log_as_error=false); // continue even though not exists +extern netCDF::NcVar get_nc_var(netCDF::NcFile *, const char *var_name, + const char *group_name, bool log_as_error=false); // continue even though not exists +extern netCDF::NcVar get_nc_var(netCDF::NcFile *, const ConcatString &var_name, + const char *group_name, bool log_as_error=false); // continue even though not exists + +extern netCDF::NcVar *copy_nc_var(netCDF::NcFile *, netCDF::NcVar *, + const int deflate_level=DEF_DEFLATE_LEVEL, const bool all_attrs=true); extern void copy_nc_att(netCDF::NcFile *, netCDF::NcVar *, const ConcatString attr_name); extern void copy_nc_att( netCDF::NcVar *, netCDF::NcVar *, const ConcatString attr_name); extern void copy_nc_atts(netCDF::NcFile *, netCDF::NcFile *, const bool all_attrs=true); extern void copy_nc_atts( netCDF::NcVar *, netCDF::NcVar *, const bool all_attrs=true); extern void copy_nc_var_data(netCDF::NcVar *, netCDF::NcVar *); -extern bool has_var(netCDF::NcFile *, const char * var_name); +extern bool has_nc_group(netCDF::NcFile *, const char *group_name); +extern bool has_var(netCDF::NcFile *, const char *var_name); +extern bool has_var(netCDF::NcFile *, const char *var_name, const char *group_name); extern netCDF::NcVar add_var(netCDF::NcFile *, const std::string &, const netCDF::NcType, const int deflate_level=DEF_DEFLATE_LEVEL); extern netCDF::NcVar add_var(netCDF::NcFile *, const std::string &, const netCDF::NcType, const netCDF::NcDim, const int deflate_level=DEF_DEFLATE_LEVEL); diff --git a/src/libcode/vx_nc_util/nc_utils.hpp b/src/libcode/vx_nc_util/nc_utils.hpp index ceb5bd2807..f7c9c553a6 100644 --- a/src/libcode/vx_nc_util/nc_utils.hpp +++ b/src/libcode/vx_nc_util/nc_utils.hpp @@ -46,12 +46,21 @@ extern void set_def_fill_value(unsigned short *val); template bool get_att_num_value_(const netCDF::NcAtt *att, T &att_val, int matching_type) { - bool status = false; - if (IS_VALID_NC_P(att)) { + bool status = IS_VALID_NC_P(att); + if (status) { int nc_type_id = GET_NC_TYPE_ID_P(att); if (matching_type == nc_type_id) { att->getValues(&att_val); - status = true; + } + else if (NC_FLOAT == nc_type_id) { + float att_value; + att->getValues(&att_value); + att_val = att_value; + } + else if (NC_DOUBLE == nc_type_id) { + double att_value; + att->getValues(&att_value); + att_val = att_value; } else if (NC_CHAR == nc_type_id) { std::string att_value; @@ -62,7 +71,6 @@ bool get_att_num_value_(const netCDF::NcAtt *att, T &att_val, int matching_type) att_val = (double)atof(att_value.c_str()); else // if (matching_type == NC_INT) att_val = atoi(att_value.c_str()); - status = true; } } return(status); @@ -74,8 +82,6 @@ template bool get_nc_att_value_(const netCDF::NcVar *var, const ConcatString &att_name, T &att_val, bool exit_on_error, T bad_data, const char *caller_name) { - bool status = false; - // Initialize att_val = bad_data; @@ -83,7 +89,7 @@ bool get_nc_att_value_(const netCDF::NcVar *var, const ConcatString &att_name, // Retrieve the NetCDF variable attribute. // netCDF::NcVarAtt *att = get_nc_att(var, att_name); - status = get_att_value((netCDF::NcAtt *)att, att_val); + bool status = get_att_value((netCDF::NcAtt *)att, att_val); if (!status) { mlog << Error << "\n" << caller_name << get_log_msg_for_att(att, GET_SAFE_NC_NAME_P(var), att_name); @@ -102,7 +108,6 @@ bool get_nc_att_value_(const netCDF::NcVar *var, const ConcatString &att_name, template bool get_nc_att_value_(const netCDF::NcVarAtt *att, T &att_val, bool exit_on_error, T bad_data, const char *caller_name) { - bool status = true; // Initialize att_val = bad_data; @@ -110,7 +115,7 @@ bool get_nc_att_value_(const netCDF::NcVarAtt *att, T &att_val, bool exit_on_err // // Retrieve the NetCDF variable attribute. // - status = get_att_value((netCDF::NcAtt *)att, att_val); + bool status = get_att_value((netCDF::NcAtt *)att, att_val); if (!status) { mlog << Error << "\n" << caller_name << get_log_msg_for_att(att); @@ -259,12 +264,10 @@ void apply_scale_factor_(T *data, const int cell_count, template bool get_nc_data_t(netCDF::NcVar *var, T *data) { - bool return_status = false; + bool return_status = IS_VALID_NC_P(var); - if (IS_VALID_NC_P(var)) { + if (return_status) { var->getVar(data); - - return_status = true; } return(return_status); } diff --git a/src/libcode/vx_stat_out/stat_columns.cc b/src/libcode/vx_stat_out/stat_columns.cc index 885479bf5c..85f89742bc 100644 --- a/src/libcode/vx_stat_out/stat_columns.cc +++ b/src/libcode/vx_stat_out/stat_columns.cc @@ -4220,12 +4220,14 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info, // // Ensemble Continuous Statistics // Dump out the ECNT line: - // TOTAL, N_ENS, CRPS, - // CRPSS, IGN, ME, - // RMSE, SPREAD, ME_OERR, - // RMSE_OERR, SPREAD_OERR, SPREAD_PLUS_OERR, - // CRPSCL, CRPS_EMP, CRPSCL_EMP, - // CRPSS_EMP CRPS_EMP_FAIR + // TOTAL, N_ENS, CRPS, + // CRPSS, IGN, ME, + // RMSE, SPREAD, ME_OERR, + // RMSE_OERR, SPREAD_OERR, SPREAD_PLUS_OERR, + // CRPSCL, CRPS_EMP, CRPSCL_EMP, + // CRPSS_EMP CRPS_EMP_FAIR, BIAS RATIO, + // N_GE_OBS, ME_GE_OBS, N_LT_OBS, + // ME_LT_OBS // at.set_entry(r, c+0, // Total Number of Pairs ecnt_info.n_pair); @@ -4278,6 +4280,21 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info, at.set_entry(r, c+16, // Empirical ensemble CRPS FAIR ecnt_info.crps_emp_fair); + at.set_entry(r, c+17, // Bias Ratio + ecnt_info.bias_ratio); + + at.set_entry(r, c+18, // Number of ensemble values >= observations + ecnt_info.n_ge_obs); + + at.set_entry(r, c+19, // ME of ensemble values >= observations + ecnt_info.me_ge_obs); + + at.set_entry(r, c+20, // Number of ensemble values < observations + ecnt_info.n_lt_obs); + + at.set_entry(r, c+21, // ME of ensemble values < observations + ecnt_info.me_lt_obs); + return; } diff --git a/src/libcode/vx_statistics/compute_stats.cc b/src/libcode/vx_statistics/compute_stats.cc index 1364fcd11b..8ad12fa36e 100644 --- a/src/libcode/vx_statistics/compute_stats.cc +++ b/src/libcode/vx_statistics/compute_stats.cc @@ -1372,6 +1372,21 @@ void compute_ecnt_mean(const ECNTInfo *ecnt_info, int n, for(i=0,na.erase(); i i && !is_bad_data(e_na[j][i])) { @@ -382,7 +401,7 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) { // Store the current ensemble value cur_ens.add(e_na[j][i]); - // Keep track of the number of ties and the number below + // Track running counts and sums if(is_eq(e_na[j][i], o_na[i])) n_tie++; else if(e_na[j][i] < o_na[i]) n_bel++; } @@ -422,6 +441,10 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) { crpscl_gaus_na.add(bad_data_double); ign_na.add(bad_data_double); pit_na.add(bad_data_double); + n_ge_obs_na.add(bad_data_double); + me_ge_obs_na.add(bad_data_double); + n_lt_obs_na.add(bad_data_double); + me_lt_obs_na.add(bad_data_double); } // Otherwise, compute scores else { @@ -492,6 +515,19 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) { crpscl_gaus_na.add(compute_crps_gaus(o_na[i], cmn_na[i], csd_na[i])); ign_na.add(compute_ens_ign(o_na[i], mean, stdev)); pit_na.add(compute_ens_pit(o_na[i], mean, stdev)); + + // Compute the Bias Ratio terms + int n_ge_obs, n_lt_obs; + double me_ge_obs, me_lt_obs; + compute_bias_ratio_terms(o_na[i], cur_ens, + n_ge_obs, me_ge_obs, + n_lt_obs, me_lt_obs); + + // Store the Bias Ratio terms + n_ge_obs_na.add(n_ge_obs); + me_ge_obs_na.add(me_ge_obs); + n_lt_obs_na.add(n_lt_obs); + me_lt_obs_na.add(me_lt_obs); } } // end for i @@ -820,7 +856,8 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o // Include in subset: // wgt_na, o_na, cmn_na, csd_na, v_na, r_na, // crps_emp_na, crps_emp_fair_na, crpscl_emp_na, crps_gaus_na, crpscl_gaus_na, - // ign_na, pit_na, var_na, var_oerr_na, var_plus_oerr_na, + // ign_na, pit_na, n_gt_obs_na, me_gt_obs_na, n_lt_obs_na, me_lt_obs_na, + // var_na, var_oerr_na, var_plus_oerr_na, // mn_na, mn_oerr_na, e_na // // Exclude from subset: @@ -841,6 +878,10 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o pd.crpscl_gaus_na.add(crpscl_gaus_na[i]); pd.ign_na.add(ign_na[i]); pd.pit_na.add(pit_na[i]); + pd.n_ge_obs_na.add(n_ge_obs_na[i]); + pd.me_ge_obs_na.add(me_ge_obs_na[i]); + pd.n_lt_obs_na.add(n_lt_obs_na[i]); + pd.me_lt_obs_na.add(me_lt_obs_na[i]); pd.skip_ba.add(false); pd.var_na.add(var_na[i]); pd.var_oerr_na.add(var_oerr_na[i]); @@ -2006,5 +2047,59 @@ double compute_ens_pit(double obs, double m, double s) { return(v); } + +//////////////////////////////////////////////////////////////////////// +// +// Compute the bias ratio terms +// +//////////////////////////////////////////////////////////////////////// + +void compute_bias_ratio_terms(double obs, const NumArray &e_na, + int &n_ge_obs, double &me_ge_obs, + int &n_lt_obs, double &me_lt_obs) { + + // Initialize + n_ge_obs = n_lt_obs = 0; + me_ge_obs = me_lt_obs = 0.0; + + // Loop over ensemble member values + for(int i=0; i= obs) { + n_ge_obs += 1; + me_ge_obs += (e_na[i] - obs); + } + else { + n_lt_obs += 1; + me_lt_obs += (e_na[i] - obs); + } + } + + // Convert sums to means + if(n_ge_obs > 0) me_ge_obs /= n_ge_obs; + else me_ge_obs = bad_data_double; + if(n_lt_obs > 0) me_lt_obs /= n_lt_obs; + else me_lt_obs = bad_data_double; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +double compute_bias_ratio(double me_ge_obs, double me_lt_obs) { + double v; + + // Compute bias ratio + if(is_bad_data(me_ge_obs) || + is_bad_data(me_lt_obs) || is_eq(me_lt_obs, 0.0)) { + v = bad_data_double; + } + else { + v = me_ge_obs / abs(me_lt_obs); + } + + return(v); +} //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_statistics/pair_data_ensemble.h b/src/libcode/vx_statistics/pair_data_ensemble.h index 17501d867a..0dc78e38b5 100644 --- a/src/libcode/vx_statistics/pair_data_ensemble.h +++ b/src/libcode/vx_statistics/pair_data_ensemble.h @@ -88,6 +88,11 @@ class PairDataEnsemble : public PairBase { NumArray ign_na; // Ignorance Score [n_obs] NumArray pit_na; // Probability Integral Transform [n_obs] + NumArray n_ge_obs_na; // Number of ensemble memebers >= obs [n_obs] + NumArray me_ge_obs_na; // Mean error of ensemble members >= obs [n_obs] + NumArray n_lt_obs_na; // Number of ensemble members < obs [n_obs] + NumArray me_lt_obs_na; // Mean error of ensemble members < obs [n_obs] + int n_ens; // Number of ensemble members int n_pair; // Number of valid pairs, n_obs - sum(skip_ba) int ctrl_index; // Index of the control member @@ -123,6 +128,8 @@ class PairDataEnsemble : public PairBase { double me_oerr; // ME for mean of perturbed members double rmse_oerr; // RMSE for mean of perturbed members + double bias_ratio; // Bias ratio + ////////////////////////////////////////////////////////////////// void clear(); @@ -310,6 +317,9 @@ extern double compute_crps_emp(double, const NumArray &); extern double compute_crps_gaus(double, double, double); extern double compute_ens_ign(double, double, double); extern double compute_ens_pit(double, double, double); +extern void compute_bias_ratio_terms(double, const NumArray &, + int &, double &, int &, double &); +extern double compute_bias_ratio(double, double); //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_tc_util/Makefile.am b/src/libcode/vx_tc_util/Makefile.am index 390eb1e5ae..463cb5495a 100644 --- a/src/libcode/vx_tc_util/Makefile.am +++ b/src/libcode/vx_tc_util/Makefile.am @@ -15,6 +15,7 @@ libvx_tc_util_a_SOURCES = \ atcf_line_base.cc atcf_line_base.h atcf_offsets.h \ atcf_track_line.cc atcf_track_line.h \ atcf_prob_line.cc atcf_prob_line.h \ + diag_file.cc diag_file.h \ track_point.cc track_point.h \ track_info.cc track_info.h \ track_pair_info.cc track_pair_info.h \ diff --git a/src/libcode/vx_tc_util/Makefile.in b/src/libcode/vx_tc_util/Makefile.in index b059aaf546..236c0abd6f 100644 --- a/src/libcode/vx_tc_util/Makefile.in +++ b/src/libcode/vx_tc_util/Makefile.in @@ -110,6 +110,7 @@ libvx_tc_util_a_LIBADD = am_libvx_tc_util_a_OBJECTS = libvx_tc_util_a-atcf_line_base.$(OBJEXT) \ libvx_tc_util_a-atcf_track_line.$(OBJEXT) \ libvx_tc_util_a-atcf_prob_line.$(OBJEXT) \ + libvx_tc_util_a-diag_file.$(OBJEXT) \ libvx_tc_util_a-track_point.$(OBJEXT) \ libvx_tc_util_a-track_info.$(OBJEXT) \ libvx_tc_util_a-track_pair_info.$(OBJEXT) \ @@ -144,6 +145,7 @@ am__maybe_remake_depfiles = depfiles am__depfiles_remade = ./$(DEPDIR)/libvx_tc_util_a-atcf_line_base.Po \ ./$(DEPDIR)/libvx_tc_util_a-atcf_prob_line.Po \ ./$(DEPDIR)/libvx_tc_util_a-atcf_track_line.Po \ + ./$(DEPDIR)/libvx_tc_util_a-diag_file.Po \ ./$(DEPDIR)/libvx_tc_util_a-gen_shape_info.Po \ ./$(DEPDIR)/libvx_tc_util_a-genesis_info.Po \ ./$(DEPDIR)/libvx_tc_util_a-pair_data_genesis.Po \ @@ -369,6 +371,7 @@ libvx_tc_util_a_SOURCES = \ atcf_line_base.cc atcf_line_base.h atcf_offsets.h \ atcf_track_line.cc atcf_track_line.h \ atcf_prob_line.cc atcf_prob_line.h \ + diag_file.cc diag_file.h \ track_point.cc track_point.h \ track_info.cc track_info.h \ track_pair_info.cc track_pair_info.h \ @@ -438,6 +441,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_tc_util_a-atcf_line_base.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_tc_util_a-atcf_prob_line.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_tc_util_a-atcf_track_line.Po@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_tc_util_a-diag_file.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_tc_util_a-gen_shape_info.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_tc_util_a-genesis_info.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_tc_util_a-pair_data_genesis.Po@am__quote@ # am--include-marker @@ -516,6 +520,20 @@ libvx_tc_util_a-atcf_prob_line.obj: atcf_prob_line.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_tc_util_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_tc_util_a-atcf_prob_line.obj `if test -f 'atcf_prob_line.cc'; then $(CYGPATH_W) 'atcf_prob_line.cc'; else $(CYGPATH_W) '$(srcdir)/atcf_prob_line.cc'; fi` +libvx_tc_util_a-diag_file.o: diag_file.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_tc_util_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_tc_util_a-diag_file.o -MD -MP -MF $(DEPDIR)/libvx_tc_util_a-diag_file.Tpo -c -o libvx_tc_util_a-diag_file.o `test -f 'diag_file.cc' || echo '$(srcdir)/'`diag_file.cc +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_tc_util_a-diag_file.Tpo $(DEPDIR)/libvx_tc_util_a-diag_file.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='diag_file.cc' object='libvx_tc_util_a-diag_file.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_tc_util_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_tc_util_a-diag_file.o `test -f 'diag_file.cc' || echo '$(srcdir)/'`diag_file.cc + +libvx_tc_util_a-diag_file.obj: diag_file.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_tc_util_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_tc_util_a-diag_file.obj -MD -MP -MF $(DEPDIR)/libvx_tc_util_a-diag_file.Tpo -c -o libvx_tc_util_a-diag_file.obj `if test -f 'diag_file.cc'; then $(CYGPATH_W) 'diag_file.cc'; else $(CYGPATH_W) '$(srcdir)/diag_file.cc'; fi` +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_tc_util_a-diag_file.Tpo $(DEPDIR)/libvx_tc_util_a-diag_file.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='diag_file.cc' object='libvx_tc_util_a-diag_file.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_tc_util_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_tc_util_a-diag_file.obj `if test -f 'diag_file.cc'; then $(CYGPATH_W) 'diag_file.cc'; else $(CYGPATH_W) '$(srcdir)/diag_file.cc'; fi` + libvx_tc_util_a-track_point.o: track_point.cc @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_tc_util_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_tc_util_a-track_point.o -MD -MP -MF $(DEPDIR)/libvx_tc_util_a-track_point.Tpo -c -o libvx_tc_util_a-track_point.o `test -f 'track_point.cc' || echo '$(srcdir)/'`track_point.cc @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_tc_util_a-track_point.Tpo $(DEPDIR)/libvx_tc_util_a-track_point.Po @@ -854,6 +872,7 @@ distclean: distclean-am -rm -f ./$(DEPDIR)/libvx_tc_util_a-atcf_line_base.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-atcf_prob_line.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-atcf_track_line.Po + -rm -f ./$(DEPDIR)/libvx_tc_util_a-diag_file.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-gen_shape_info.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-genesis_info.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-pair_data_genesis.Po @@ -917,6 +936,7 @@ maintainer-clean: maintainer-clean-am -rm -f ./$(DEPDIR)/libvx_tc_util_a-atcf_line_base.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-atcf_prob_line.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-atcf_track_line.Po + -rm -f ./$(DEPDIR)/libvx_tc_util_a-diag_file.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-gen_shape_info.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-genesis_info.Po -rm -f ./$(DEPDIR)/libvx_tc_util_a-pair_data_genesis.Po diff --git a/src/libcode/vx_tc_util/diag_file.cc b/src/libcode/vx_tc_util/diag_file.cc new file mode 100644 index 0000000000..32a37d2f07 --- /dev/null +++ b/src/libcode/vx_tc_util/diag_file.cc @@ -0,0 +1,542 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2022 +// ** 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 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + +using namespace std; + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "vx_util.h" +#include "vx_log.h" + +#include "diag_file.h" + +//////////////////////////////////////////////////////////////////////// + +static const char default_lsdiag_technique[] = "OFCL"; + +static const int lsdiag_wdth[] = { + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5 +}; +static int n_lsdiag_wdth = sizeof(lsdiag_wdth)/sizeof(*lsdiag_wdth); + +static const int tcdiag_fill_value = 9999; +static const int lsdiag_fill_value = 9999; + +static const char lsdiag_skip_str[] = "TIME,DELV,MTPW,IR00,IRM1,IRM3,PC00,PCM1,PCM3,PSLV,IRXX"; + +//////////////////////////////////////////////////////////////////////// +// +// Code for class DiagFile +// +//////////////////////////////////////////////////////////////////////// + +DiagFile::DiagFile() { + init_from_scratch(); +} + +//////////////////////////////////////////////////////////////////////// + +DiagFile::~DiagFile() { + close(); +} + +//////////////////////////////////////////////////////////////////////// + +int DiagFile::lead(int i) const { + + // Check range + if(i < 0 || i >= LeadTime.n()) { + mlog << Error << "\nDiagFile::lead(int) -> " + << "range check error for index value " << i << "\n\n"; + exit(1); + } + + return(LeadTime[i]); +} + +//////////////////////////////////////////////////////////////////////// + +unixtime DiagFile::valid(int i) const { + + // Check range + if(i < 0 || i >= LeadTime.n()) { + mlog << Error << "\nDiagFile::valid(int) -> " + << "range check error for index value " << i << "\n\n"; + exit(1); + } + + return(InitTime == 0 || is_bad_data(LeadTime[i]) ? + 0 : InitTime + LeadTime[i]); +} + +//////////////////////////////////////////////////////////////////////// + +double DiagFile::lat(int i) const { + + // Check range + if(i < 0 || i >= Lat.n()) { + mlog << Error << "\nDiagFile::lat(int) -> " + << "range check error for index value " << i << "\n\n"; + exit(1); + } + + return(Lat[i]); +} + +//////////////////////////////////////////////////////////////////////// + +double DiagFile::lon(int i) const { + + // Check range + if(i < 0 || i >= Lon.n()) { + mlog << Error << "\nDiagFile::lon(int) -> " + << "range check error for index value " << i << "\n\n"; + exit(1); + } + + return(Lon[i]); +} + +//////////////////////////////////////////////////////////////////////// + +bool DiagFile::has_diag(const string &str) const { + return(DiagName.has(str)); +} + +//////////////////////////////////////////////////////////////////////// + +const NumArray & DiagFile::diag_val(const string &str) const { + int i; + + // Find the index of the name + if(!DiagName.has(str, i)) { + mlog << Error << "\nDiagFile::diag_val() -> " + << "requested diagnostic name \"" << str << "\" not found!\n\n"; + exit(1); + } + + return(DiagVal[i]); +} + +//////////////////////////////////////////////////////////////////////// + +DiagFile::DiagFile(const DiagFile &) { + mlog << Error << "\nDiagFile::DiagFile(const DiagFile &) -> " + << "should never be called!\n\n"; + exit(1); +} + +//////////////////////////////////////////////////////////////////////// + +DiagFile & DiagFile::operator=(const DiagFile &) { + mlog << Error << "\nDiagFile::operator=(const DiagFile &) -> " + << "should never be called!\n\n"; + exit(1); + + return(*this); +} + +//////////////////////////////////////////////////////////////////////// + +void DiagFile::init_from_scratch() { + + clear(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void DiagFile::clear() { + + // Initialize values + SourceType = DiagType_None; + StormId.clear(); + Basin.clear(); + Cyclone.clear(); + Technique.clear(); + InitTime = (unixtime) 0; + + NTime = 0; + LeadTime.clear(); + Lat.clear(); + Lon.clear(); + DiagName.clear(); + DiagVal.clear(); + + close(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void DiagFile::set_technique(const StringArray &sa) { + + for(int i=0; i *convert_map) { + + // Read diagnostics baesd on the source type + if(source == TCDiagType) { + read_tcdiag(path, model_names, convert_map); + } + else if(source == LSDiagRTType) { + read_lsdiag_rt(path, model_names, convert_map); + } + else { + mlog << Error << "\nDiagFile::read() -> " + << "diagnostics of type " << diagtype_to_string(source) + << " are not currently supported!\n\n"; + exit(1); + } + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void DiagFile::read_tcdiag(const ConcatString &path, const StringArray &model_names, + const map *convert_map) { + int i; + double v_in, v_out; + NumArray data; + const UserFunc_1Arg *fx_ptr = 0; + + // Initialize + clear(); + + // Store the file type + SourceType = TCDiagType; + + // Store user-specified model names or read it from the file below + if(model_names.n() > 0) set_technique(model_names); + + // Open the file + open(path.c_str()); + + // Parse the header information + DataLine dl; + ConcatString cs; + while(dl.read_line(this)) { + + // Skip empty lines + if(dl.n_items() == 0) continue; + + // First column + cs = dl[0]; + + // Parse header lines + if(cs == "*") { + + // First header line: Technique InitTime (ATCFID YYYMMDDHH) + if(InitTime == 0) { + if(Technique.n() == 0) add_technique(dl[1]); + InitTime = timestring_to_unix(dl[2]); + } + // Second header line: Basin Cyclone Number (BBCC) + else if(StormId.empty()) { + string bbcc = dl[1]; + Basin = bbcc.substr(0, 2); + Cyclone = bbcc.substr(2, 2); + int mon, day, yr, hr, min, sec; + unix_to_mdyhms(InitTime, mon, day, yr, hr, min, sec); + StormId << Basin << Cyclone << yr; + } + } + // Parse time and location info + else if(cs == "NTIME") { + NTime = atoi(dl[1]); + } + else if(cs == "TIME") { + for(i=2; i " + << "the NTIME value (" << NTime + << ") does not match the actual number of times (" + << LeadTime.n() << "), latitudes (" << Lat.n() + << "), or longitudes (" << Lon.n() << "): " + << path << "\n\n"; + exit(1); + } + + // Rewind to beginning in case data rows precede the header + rewind(); + + // Store the diagnostics data + while(dl.read_line(this)) { + + // Skip empty lines + if(dl.n_items() == 0) continue; + + // First column contains the diagnostic name + cs = dl[0]; + + // Skip non-diagnostic data lines + if(cs.startswith("*") || cs.startswith("----") || + cs.startswith("NTIME") || cs.startswith("TIME") || + cs.startswith("NLEV") || cs.startswith("NVAR")) continue; + + // Check for a conversion function based on the diagnostic name or units + if(convert_map) { + if(convert_map->count(cs) > 0) fx_ptr = &convert_map->at(cs); + else if(convert_map->count(dl[1]) > 0) fx_ptr = &convert_map->at(dl[1]); + else fx_ptr = (UserFunc_1Arg *) 0; + } + else { + fx_ptr = (UserFunc_1Arg *) 0; + } + + // Parse the data values + data.erase(); + for(int i=2; i " + << "the number of \"" << cs << "\" diagnostic values (" + << data.n() << ") does not match the expected number (" + << NTime << "): " << path << "\n\n"; + exit(1); + } + // Store the name and values + else { + DiagName.add(cs); + DiagVal.push_back(data); + } + } // end while + + mlog << Debug(4) << "Parsed " << DiagName.n() << " diagnostic values from " + << StormId << " " << write_css(Technique) << " " << unix_to_yyyymmddhh(InitTime) + << " TC diagnostics file: " << path << "\n"; + + // Close the input file + close(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model_names, + const map *convert_map) { + int i, v_int; + double v_dbl; + NumArray data; + const UserFunc_1Arg *fx_ptr = 0; + + // Initialize + clear(); + + // Store the file type + SourceType = LSDiagRTType; + + // Store user-specified model names or the default one + if(model_names.n() > 0) set_technique(model_names); + else add_technique(default_lsdiag_technique); + + // Open the file + open(path.c_str()); + + // Diagnostic names to ignore + StringArray skip_diag_sa; + skip_diag_sa.parse_css(lsdiag_skip_str); + + // Parse the header information from the first line + DataLine dl; + ConcatString cs; + dl.read_line(this); + + // Check for the expected number of items + if(dl.n_items() != 9 || strncasecmp(dl[8], "HEAD", strlen("HEAD") != 0)) { + mlog << Error << "\nDiagFile::read_lsdiag_rt() -> " + << "unexpected header line: " << path << "\n\n"; + exit(1); + } + + // Parse storm information + StormId = dl[7]; + Basin = StormId.string().substr(0, 2); + Cyclone = StormId.string().substr(2, 2); + + // Parse timing info + cs = dl[1]; + int yr = atoi(StormId.string().substr(4, 4).c_str()); + int mon = atoi(cs.string().substr(2, 2).c_str()); + int day = atoi(cs.string().substr(4, 2).c_str()); + int hr = atoi(dl[2]); + InitTime = mdyhms_to_unix(mon, day, yr, hr, 0, 0); + + // Store the location of the beginning of the data + int data_start_location = in->tellg(); + + // Parse time and location info + while(read_fwf_line(dl, lsdiag_wdth, n_lsdiag_wdth)) { + + // Fixed width column 24 has the data name + cs = dl[23]; + + // Strip any whitespace from the fixed-width column + cs.ws_strip(); + + if(cs == "TIME") { + for(i=2; i<23; i++) { + LeadTime.add(atoi(dl[i])*sec_per_hour); + } + NTime = LeadTime.n(); + } + else if(cs == "LAT") { + // Tenths of degree north + for(i=2; i<23; i++) { + Lat.add(atof(dl[i])/10.0); + } + } + else if(cs == "LON") { + for(i=2; i<23; i++) { + // Tenths of degrees west (convert to east) + Lon.add(rescale_lon(-1.0*atof(dl[i])/10.0)); + } + + // Finished parsing the metadata + break; + } + } // end while + + // Rewind to the beginning of the data + in->seekg(data_start_location); + + // Store the diagnostics data + while(read_fwf_line(dl, lsdiag_wdth, n_lsdiag_wdth)) { + + // Skip empty lines + if(dl.n_items() == 0) continue; + + // The 24th column contains the diagnostic name + cs = dl[23]; + + // Strip any whitespace from the fixed-width column + cs.ws_strip(); + + // Check for diagnostic names to skip + if(skip_diag_sa.has(cs)) continue; + + // Quit reading at the LAST line + if(cs == "LAST") break; + + // Check for a conversion function + if(convert_map) { + if(convert_map->count(cs) > 0) fx_ptr = &convert_map->at(cs); + else fx_ptr = (UserFunc_1Arg *) 0; + } + else { + fx_ptr = (UserFunc_1Arg *) 0; + } + + // Parse the data values + data.erase(); + for(i=2; i<23; i++) { + + v_int = atoi(dl[i]); + + // Check for bad data and apply conversions + if(v_int == lsdiag_fill_value) v_dbl = bad_data_double; + else if(fx_ptr) v_dbl = (*fx_ptr)(v_int); + else v_dbl = (double) v_int; + + // Store the value + data.add(v_dbl); + } + + // Check for the expected number of items + if(NTime != data.n()) { + mlog << Error << "\nDiagFile::read_lsdiag_rt() -> " + << "the number of \"" << cs << "\" diagnostic values (" + << data.n() << ") does not match the expected number (" + << NTime << ")!\n\n"; + exit(1); + } + // Store the name and values + else { + DiagName.add(cs); + DiagVal.push_back(data); + } + } // end while + + mlog << Debug(4) << "Parsed " << DiagName.n() << " diagnostic values from " + << StormId << " " << write_css(Technique) << " " << unix_to_yyyymmddhh(InitTime) + << " LS diagnostics file: " << path << "\n"; + + // Close the input file + close(); + + return; +} + +//////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_tc_util/diag_file.h b/src/libcode/vx_tc_util/diag_file.h new file mode 100644 index 0000000000..6db6848870 --- /dev/null +++ b/src/libcode/vx_tc_util/diag_file.h @@ -0,0 +1,147 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2022 +// ** 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 __DIAG_FILE_H__ +#define __DIAG_FILE_H__ + +//////////////////////////////////////////////////////////////////////// + +#include +#include +#include + +#include "vx_config.h" +#include "vx_cal.h" +#include "data_line.h" + +//////////////////////////////////////////////////////////////////////// +// +// TCDIAG files: +// - Add link to sample data +// - Header: +// * ATCF_ID YYYYMMDDHH * +// * BBNN BBNN * +// ATCF_ID is the technique (model) name +// YYYYMMDDHH is the initialization time +// BB is the 2-letter basin name +// CC is the 2-digit cyclone number +// +// Real-time LSDIAG files: +// - https://ftp.nhc.noaa.gov/atcf/lsdiag +// - Header: +// BBCC YYMMDD HH WS LAT LON 9999 BBCCYYYY +// BB is 2-letter basin name +// CC is 2-digit cyclone number +// YYMMDD HH is the initialization time +// WS is the wind speed +// +// Developmental LSDIAG files (not currently supported): +// - https://rammb2.cira.colostate.edu/research/tropical-cyclones/ships +// +//////////////////////////////////////////////////////////////////////// + +class DiagFile : public LineDataFile { + + private: + + DiagFile(const DiagFile &); + DiagFile & operator=(const DiagFile &); + + protected: + + void init_from_scratch(); + + // Diagnostics file type + DiagType SourceType; + + // Storm and model identification + ConcatString StormId; + ConcatString Basin; + ConcatString Cyclone; + StringArray Technique; + unixtime InitTime; + + int NTime; + IntArray LeadTime; + NumArray Lat; + NumArray Lon; + + // Diagnostic names and values + StringArray DiagName; + std::vector DiagVal; + + public: + + DiagFile(); + ~DiagFile(); + + void clear(); + + // + // set stuff + // + + void set_technique(const StringArray &); + void add_technique(const std::string &); + + // + // get stuff + // + + const ConcatString & storm_id() const; + const ConcatString & basin() const; + const ConcatString & cyclone() const; + const StringArray & technique() const; + const ConcatString & initials() const; + unixtime init() const; + + int n_time() const; + int lead(int) const; + unixtime valid(int) const; + double lat(int) const; + double lon(int) const; + + DiagType source() const; + int n_diag() const; + const StringArray & diag_name() const; + bool has_diag(const std::string &) const; + const NumArray & diag_val(const std::string &) const; + + // + // do stuff + // + + void read (const DiagType, + const ConcatString &, const StringArray &, + const std::map *); + void read_tcdiag (const ConcatString &, const StringArray &, + const std::map *); + void read_lsdiag_rt (const ConcatString &, const StringArray &, + const std::map *); + +}; + +//////////////////////////////////////////////////////////////////////// + +inline const ConcatString & DiagFile::storm_id() const { return(StormId); } +inline const ConcatString & DiagFile::basin() const { return(Basin); } +inline const ConcatString & DiagFile::cyclone() const { return(Cyclone); } +inline const StringArray & DiagFile::technique() const { return(Technique); } +inline unixtime DiagFile::init() const { return(InitTime); } +inline int DiagFile::n_time() const { return(NTime); } +inline DiagType DiagFile::source() const { return(SourceType); } +inline int DiagFile::n_diag() const { return(DiagName.n()); } +inline const StringArray & DiagFile::diag_name() const { return(DiagName); } + +//////////////////////////////////////////////////////////////////////// + +#endif /* __DIAG_FILE_H__ */ + +//////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_tc_util/tc_columns.cc b/src/libcode/vx_tc_util/tc_columns.cc index ec73a544d9..4cb90895f4 100644 --- a/src/libcode/vx_tc_util/tc_columns.cc +++ b/src/libcode/vx_tc_util/tc_columns.cc @@ -36,8 +36,7 @@ void open_tc_txt_file(ofstream *&out, const char *file_name) { out->open(file_name); if(!(*out)) { - mlog << Error - << "\nopen_tc_txt_file()-> " + mlog << Error << "\nopen_tc_txt_file()-> " << "can't open the output file \"" << file_name << "\" for writing!\n\n"; exit(1); @@ -108,6 +107,37 @@ void write_tc_mpr_header_row(int hdr_flag, AsciiTable &at, //////////////////////////////////////////////////////////////////////// +void write_tc_diag_header_row(int hdr_flag, int n_diag, AsciiTable &at, + int r, int c) { + int i; + ConcatString s; + char tmp_str[max_str_len]; + + // Write the header column names if requested + if(hdr_flag) { + for(i=0; i i) { - hdr.set_desc((string)p.line(i)->get_item("DESC", false)); + hdr.set_desc((string)p.tcmpr_line(i)->get_item("DESC", false)); } // Write the header columns @@ -171,6 +198,22 @@ void write_tc_mpr_row(TcHdrColumns &hdr, const TrackPairInfo &p, // Increment the row counter i_row++; + + // Write diagnostics line + if(p.adeck()[i].n_diag() > 0) { + + // TCDIAG line type + hdr.set_line_type((string) TCStatLineType_TCDIAG_Str); + + // Write the header columns + write_tc_header_cols(hdr, at, i_row); + + // Write the data columns + write_tc_diag_cols(p, i, at, i_row, n_tc_header_cols); + + // Increment the row counter + i_row++; + } } return; @@ -178,11 +221,11 @@ void write_tc_mpr_row(TcHdrColumns &hdr, const TrackPairInfo &p, //////////////////////////////////////////////////////////////////////// -void write_prob_rirw_row(TcHdrColumns &hdr, const ProbRIRWPairInfo &p, - AsciiTable &at, int &i_row) { +void write_prob_rirw_pair_info(TcHdrColumns &hdr, const ProbRIRWPairInfo &p, + AsciiTable &at, int &i_row) { // PROBRIRW line type - hdr.set_line_type((string)"PROBRIRW"); + hdr.set_line_type((string) TCStatLineType_ProbRIRW_Str); // Timing information hdr.set_init (p.prob_rirw().init()); @@ -243,7 +286,7 @@ void write_tc_mpr_cols(const TrackPairInfo &p, int i, AsciiTable &at, int r, int c) { int j; - // Write tc_mpr columns + // Write TCMPR columns at.set_entry(r, c++, p.n_points()); at.set_entry(r, c++, i+1); at.set_entry(r, c++, cyclonelevel_to_string(p.bdeck()[i].level())); @@ -309,6 +352,35 @@ void write_tc_mpr_cols(const TrackPairInfo &p, int i, //////////////////////////////////////////////////////////////////////// +void write_tc_diag_cols(const TrackPairInfo &p, int i, + AsciiTable &at, int r, int c) { + int j; + + // Write TCDIAG columns + at.set_entry(r, c++, p.n_points()); + at.set_entry(r, c++, i+1); + at.set_entry(r, c++, diagtype_to_string(p.adeck().diag_source())); + at.set_entry(r, c++, p.adeck()[i].n_diag()); + + // Check the number of names and values match + if(p.adeck().n_diag() != p.adeck()[i].n_diag()) { + mlog << Error << "\nwrite_tc_diag_cols()-> " + << "the number of diagnostic names (" << p.adeck().n_diag() + << ") and values (" << p.adeck()[i].n_diag() + << ") do not match!\n\n"; + exit(1); + } + + for(j=0; j 0 ? unix_to_yyyymmdd_hhmmss(MaxValidTime).text() : na_str) << "\n"; out << prefix << "MinWarmCore = \"" << (MinWarmCore > 0 ? unix_to_yyyymmdd_hhmmss(MinWarmCore).text() : na_str) << "\n"; out << prefix << "MaxWarmCore = \"" << (MaxWarmCore > 0 ? unix_to_yyyymmdd_hhmmss(MaxWarmCore).text() : na_str) << "\n"; + out << prefix << "DiagSource = " << diagtype_to_string(DiagSource) << "\n"; + out << prefix << "NDiag = " << DiagName.n() << "\n"; out << prefix << "NPoints = " << NPoints << "\n"; out << prefix << "NAlloc = " << NAlloc << "\n"; - out << prefix << "NTrackLines = " << TrackLines.n_elements() << "\n"; + out << prefix << "NTrackLines = " << TrackLines.n() << "\n"; for(i=0; i " + mlog << Error << "\nvoid TrackInfo::extend(int, bool) -> " << "memory allocation error\n\n"; exit(1); } @@ -309,8 +316,7 @@ void TrackInfo::set_point(int n, const TrackPoint &p) { // Check range if((n < 0) || (n >= NPoints)) { - mlog << Error - << "\nTrackInfo::set_point(int, const TrackPoint &) -> " + mlog << Error << "\nTrackInfo::set_point(int, const TrackPoint &) -> " << "range check error for index value " << n << "\n\n"; exit(1); } @@ -356,8 +362,7 @@ const TrackPoint & TrackInfo::operator[](int n) const { // Check range if((n < 0) || (n >= NPoints)) { - mlog << Error - << "\nTrackInfo::operator[](int) -> " + mlog << Error << "\nTrackInfo::operator[](int) -> " << "range check error for index value " << n << "\n\n"; exit(1); } @@ -391,6 +396,12 @@ int TrackInfo::warm_core_dur() const { //////////////////////////////////////////////////////////////////////// +const char * TrackInfo::diag_name(int i) const { + return(i>=0 && i " + << "the " << StormId << " " << Technique << " " << unix_to_yyyymmddhh(InitTime) + << " lead time " << sec_to_timestring(diag_file.lead(i_time)) + << " track location (" << Point[i_pnt].lat() << ", " + << Point[i_pnt].lon() << ") does not match the diagnostic location (" + << diag_file.lat(i_time) << ", " << diag_file.lon(i_time) << ")\n"; + } + } + + // Store this diagnostic value in the TrackPoint + Point[i_pnt].add_diag_value(diag_val[i_time]); + + } // end for i_time + } // end for i_name + + return(true); +} + +//////////////////////////////////////////////////////////////////////// + +void TrackInfo::add_diag_value(int i_pnt, double val) { + + // Range check + if(i_pnt < 0 || i_pnt >= NPoints) { + mlog << Error << "\nTrackInfo::add_diag_value() -> " + << "range check error for point " << i_pnt << "\n\n"; + exit(1); + } + + Point[i_pnt].add_diag_value(val); + + return; +} + +//////////////////////////////////////////////////////////////////////// + bool TrackInfo::has(const ATCFTrackLine &l) const { return(TrackLines.has(l.get_line())); } @@ -734,8 +817,7 @@ const TrackInfo & TrackInfoArray::operator[](int n) const { // Check range if((n < 0) || (n >= Track.size())) { - mlog << Error - << "\nTrackInfoArray::operator[](int) -> " + mlog << Error << "\nTrackInfoArray::operator[](int) -> " << "range check error for index value " << n << "\n\n"; exit(1); } @@ -758,8 +840,7 @@ void TrackInfoArray::set(int n, const TrackInfo &t) { // Check range if((n < 0) || (n >= Track.size())) { - mlog << Error - << "\nTrackInfoArray::set(int, const TrackInfo &) -> " + mlog << Error << "\nTrackInfoArray::set(int, const TrackInfo &) -> " << "range check error for index value " << n << "\n\n"; exit(1); } @@ -843,6 +924,19 @@ bool TrackInfoArray::erase_storm_id(const ConcatString &s) { return(status); } +//////////////////////////////////////////////////////////////////////// + +int TrackInfoArray::add_diag_data(DiagFile &diag_file, const StringArray &diag_name) { + int n_match = 0; + + // Set the names for each track + for(int i=0; i " + mlog << Error << "\nTrackInfoArray::consensus() -> " << "cannot compute a consensus for zero tracks!\n\n"; exit(1); } @@ -889,8 +982,7 @@ TrackInfo consensus(const TrackInfoArray &tracks, if(tavg.basin() != tracks[i].basin() || tavg.cyclone() != tracks[i].cyclone() || tavg.init() != tracks[i].init()) { - mlog << Error - << "\nTrackInfoArray::consensus() -> " + mlog << Error << "\nTrackInfoArray::consensus() -> " << "the basin, cyclone number, and init time must " << "remain constant.\n\n"; exit(1); @@ -920,7 +1012,7 @@ TrackInfo consensus(const TrackInfoArray &tracks, } // Loop through the lead times and construct a TrackPoint for each - for(i=0, skip=false; i 180.0 && plon[j] < 0.0 ? 360.0 : 0.0); lon_avg += (plon[j] + lon_shift); } @@ -1087,7 +1179,7 @@ bool has_storm_id(const StringArray &storm_id, bool match = false; // Loop over the storm id entries - for(i=0; i= n) return; - - // Compute the allocation size - if(!exact) { - k = n/TrackPairLineAllocInc; - if(n%TrackPairLineAllocInc) k++; - n = k*TrackPairLineAllocInc; - } - - // Allocate a new TCStatLine array of the required length - new_line = new TCStatLine [n]; - - if(!new_line) { - mlog << Error - << "\nvoid TrackPairInfo::extend(int, bool) -> " - << "memory allocation error\n\n"; - exit(1); - } - - // Copy the array contents and delete the old one - if(Line) { - for(j=0; j= NPoints) { - mlog << Error - << "\nTrackPairInfo::set_keep(int, int) -> " + mlog << Error << "\nTrackPairInfo::set_keep(int, int) -> " << "range check error for index value " << i << "\n\n"; exit(1); } @@ -322,6 +273,10 @@ void TrackPairInfo::set_keep(int i, int val) { return; } +//////////////////////////////////////////////////////////////////////// +// +// Add derived TrackPoint pairs such as for consensus tracks +// //////////////////////////////////////////////////////////////////////// void TrackPairInfo::add(const TrackPoint &a, const TrackPoint &b, @@ -354,6 +309,17 @@ void TrackPairInfo::add(const TrackPoint &a, const TrackPoint &b, //////////////////////////////////////////////////////////////////////// void TrackPairInfo::add(const TCStatLine &l) { + + // Check the line type + if(l.type() == TCStatLineType_TCMPR) add_tcmpr_line(l); + else if(l.type() == TCStatLineType_TCDIAG) add_tcdiag_line(l); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void TrackPairInfo::add_tcmpr_line(const TCStatLine &l) { TrackPoint apoint, bpoint; TrackPoint *tp = (TrackPoint *) 0; QuadInfo wind; @@ -363,8 +329,14 @@ void TrackPairInfo::add(const TCStatLine &l) { // Check the line type if(l.type() != TCStatLineType_TCMPR) return; - // Increment the count + // Store the input TCMPR line and TCDIAG placeholder + TCMPRLine.push_back(l); + TCStatLine empty_line; + TCDIAGLine.push_back(empty_line); + + // Increment the point and line count NPoints++; + NLines++; // Initialize the ADECK/BDECK tracks if(NPoints == 1) initialize(l); @@ -457,9 +429,88 @@ void TrackPairInfo::add(const TCStatLine &l) { BDeckPrvInt.add(bad_data_double); Keep.add(1); - // Store the input line - extend(NLines + 1, false); - Line[NLines++] = l; + return; +} + +//////////////////////////////////////////////////////////////////////// + +void TrackPairInfo::add_tcdiag_line(const TCStatLine &l) { + int n_diag, i; + ConcatString cs; + + // Check the line type + if(l.type() != TCStatLineType_TCDIAG) return; + + // Should have already parsed TCMPR + if(NPoints == 0) { + mlog << Error << "\nTrackPairInfo::add_tcdiag_line() -> " + << "each TCDIAG line must be preceded by the TCMPR line " + << "to which it corresponds:\n" + << " " << l << "\n\n"; + exit(1); + } + + // Check for a match + if(ADeck.storm_id() != l.storm_id() || + ADeck.technique() != l.amodel() || + ADeck.init() != l.init() || + ADeck[NPoints-1].lead() != l.lead()) { + mlog << Error << "\nTrackPairInfo::add_tcdiag_line() -> " + << "the TCDIAG data does not match the TCMPR data:\n" + << " TCMPR Line: " << TCMPRLine[NPoints-1] << "\n" + << " TCDIAG Line: " << l << "\n\n"; + exit(1); + } + + // Update the TCDIAGLine for this TrackPoint + TCDIAGLine[NPoints-1] = l; + + // Increment the line count + NLines++; + + // Name of diagnostics read + StringArray diag_name; + + // Diagnostic source type + DiagType t = string_to_diagtype(l.get_item("SOURCE")); + + // Make sure the source type does not change + if(ADeck.diag_source() != DiagType_None && + ADeck.diag_source() != t) { + mlog << Error << "\nTrackPairInfo::add_tcdiag_line() -> " + << "the diagnostic source type has changed (" + << diagtype_to_string(ADeck.diag_source()) << " != " + << diagtype_to_string(t) << ")!\n\n"; + exit(1); + } + + // Store the source type + ADeck.set_diag_source(t); + + // Number of diagnostics + n_diag = atoi(l.get_item("N_DIAG")); + + // Parse all the diagnostics entries + for(i=0; i " + << "the number of TCDIAG diagnostics have changed (" + << ADeck.diag_name().n() << " != " << diag_name.n() + << ")!\n\n"; + exit(1); + } return; } @@ -575,7 +626,7 @@ int TrackPairInfo::check_rirw(const TrackType track_type, if(!exact_adeck && st_adeck.get_type() != thresh_lt && st_adeck.get_type() != thresh_le && st_adeck.get_type() != thresh_gt && st_adeck.get_type() != thresh_ge) { - mlog << Error << "\ncheck_rirw() -> " + mlog << Error << "\nTrackPairInfo::check_rirw() -> " << "for non-exact intensity differences the RI/RW threshold (" << st_adeck.get_str() << ") must be of type <, <=, >, or >=.\n\n"; exit(1); @@ -585,7 +636,7 @@ int TrackPairInfo::check_rirw(const TrackType track_type, if(!exact_bdeck && st_bdeck.get_type() != thresh_lt && st_bdeck.get_type() != thresh_le && st_bdeck.get_type() != thresh_gt && st_bdeck.get_type() != thresh_ge) { - mlog << Error << "\ncheck_rirw() -> " + mlog << Error << "\nTrackPairInfo::check_rirw() -> " << "for non-exact intensity differences the RI/RW threshold (" << st_bdeck.get_str() << ") must be of type <, <=, >, or >=.\n\n"; exit(1); @@ -801,8 +852,9 @@ TrackPairInfo TrackPairInfo::keep_subset() const { if(Keep[i] == 0) continue; // Add from TC-Stat line data - if(NLines == NPoints) { - tpi.add(Line[i]); + if(TCMPRLine.size() > 0) { + tpi.add_tcmpr_line(TCMPRLine[i]); + tpi.add_tcdiag_line(TCDIAGLine[i]); } // Otherwise, add from TC-Pairs track pair else { @@ -963,8 +1015,7 @@ void TrackPairInfoArray::extend(int n, bool exact) { new_info = new TrackPairInfo [n]; if(!new_info) { - mlog << Error - << "\nvoid TrackPairInfoArray::extend(int, bool) -> " + mlog << Error << "\nTrackPairInfoArray::extend(int, bool) -> " << "memory allocation error\n\n"; exit(1); } @@ -991,8 +1042,7 @@ const TrackPairInfo & TrackPairInfoArray::operator[](int n) const { // Check range if((n < 0) || (n >= NPairs)) { - mlog << Error - << "\nTrackPairInfoArray::operator[](int) -> " + mlog << Error << "\nTrackPairInfoArray::operator[](int) -> " << "range check error for index value " << n << "\n\n"; exit(1); } @@ -1012,6 +1062,18 @@ int TrackPairInfoArray::n_points() const { //////////////////////////////////////////////////////////////////////// +int TrackPairInfoArray::max_n_diag() const { + int i, n; + + for(i=0,n=0; i n) n = Pair[i].adeck().n_diag(); + } + + return(n); +} + +//////////////////////////////////////////////////////////////////////// + void TrackPairInfoArray::add(const TrackPairInfo &p) { extend(NPairs + 1, false); diff --git a/src/libcode/vx_tc_util/track_pair_info.h b/src/libcode/vx_tc_util/track_pair_info.h index b4e92739ab..69db58d9e8 100644 --- a/src/libcode/vx_tc_util/track_pair_info.h +++ b/src/libcode/vx_tc_util/track_pair_info.h @@ -14,6 +14,7 @@ //////////////////////////////////////////////////////////////////////// #include +#include #include "track_point.h" #include "track_info.h" @@ -42,10 +43,10 @@ class TrackPairInfo { void init_from_scratch(); void assign(const TrackPairInfo &); - void extend(int, bool exact = true); - // Number of track points + // Number of track points and input lines int NPoints; + int NLines; // ADECK and BDECK tracks TrackInfo ADeck; @@ -63,9 +64,8 @@ class TrackPairInfo { NumArray CrossTrackErr; // TCStatLines used to construct this track - int NLines; - int NAlloc; - TCStatLine * Line; + std::vector TCMPRLine; // sized by NPoint + std::vector TCDIAGLine; // sized by NPoint // Status for whether RI/RW occurred NumArray ADeckRIRW; @@ -121,7 +121,8 @@ class TrackPairInfo { int bdeck_prv_int(int) const; int n_lines() const; - const TCStatLine * line(int i) const; + const TCStatLine * tcmpr_line(int i) const; + const TCStatLine * tcdiag_line(int i) const; int i_init() const; bool keep(int) const; @@ -134,6 +135,8 @@ class TrackPairInfo { void add(const TrackPoint &, const TrackPoint &, double, double, double, double, double, double, double); void add(const TCStatLine &); + void add_tcmpr_line(const TCStatLine &); + void add_tcdiag_line(const TCStatLine &); void add_watch_warn(const ConcatString &, WatchWarnType, unixtime); int check_water_only(); @@ -163,7 +166,8 @@ inline int TrackPairInfo::bdeck_rirw(int i) const { return(n inline int TrackPairInfo::adeck_prv_int(int i) const { return(nint(ADeckPrvInt[i])); } inline int TrackPairInfo::bdeck_prv_int(int i) const { return(nint(BDeckPrvInt[i])); } inline int TrackPairInfo::n_lines() const { return(NLines); } -inline const TCStatLine * TrackPairInfo::line(int i) const { return(&Line[i]); } +inline const TCStatLine * TrackPairInfo::tcmpr_line(int i) const { return(&TCMPRLine[i]); } +inline const TCStatLine * TrackPairInfo::tcdiag_line(int i) const { return(&TCDIAGLine[i]); } inline bool TrackPairInfo::keep(int i) const { return(Keep[i] != 0); } //////////////////////////////////////////////////////////////////////// @@ -206,8 +210,9 @@ class TrackPairInfoArray { // const TrackPairInfo & operator[](int) const; - int n_pairs() const; - int n_points() const; + int n_pairs() const; + int n_points() const; + int max_n_diag() const; // // do stuff diff --git a/src/libcode/vx_tc_util/track_point.cc b/src/libcode/vx_tc_util/track_point.cc index 9937034156..0badcb827c 100644 --- a/src/libcode/vx_tc_util/track_point.cc +++ b/src/libcode/vx_tc_util/track_point.cc @@ -66,8 +66,7 @@ QuadInfo & QuadInfo::operator+=(const QuadInfo &t) { // Check intensity if(is_bad_data(Intensity)) Intensity = t.intensity(); else if(Intensity != t.intensity()) { - mlog << Error - << "\nQuadInfo::operator+=(const QuadInfo &t) -> " + mlog << Error << "\nQuadInfo::operator+=(const QuadInfo &t) -> " << "cannot call += for two different intensity values (" << Intensity << " != " << t.intensity() << ").\n\n"; exit(1); @@ -265,8 +264,7 @@ void QuadInfo::set_quad_vals(QuadrantType ref_quad, break; default: - mlog << Error - << "\nQuadInfo::set_quad_vals() -> " + mlog << Error << "\nQuadInfo::set_quad_vals() -> " << "unexpected quadrant type encountered \"" << quadranttype_to_string(ref_quad) << "\".\n\n"; exit(1); @@ -343,8 +341,7 @@ TrackPoint & TrackPoint::operator+=(const TrackPoint &p) { // Check valid time if(ValidTime == (unixtime) 0) ValidTime = p.valid(); else if(ValidTime != p.valid()) { - mlog << Error - << "\nTrackPoint::operator+=(const TrackPoint &p) -> " + mlog << Error << "\nTrackPoint::operator+=(const TrackPoint &p) -> " << "cannot call += for two different valid times (" << unix_to_yyyymmdd_hhmmss(ValidTime) << " != " << unix_to_yyyymmdd_hhmmss(p.valid()) << ").\n\n"; @@ -354,8 +351,7 @@ TrackPoint & TrackPoint::operator+=(const TrackPoint &p) { // Check lead time if(is_bad_data(LeadTime)) LeadTime = p.lead(); else if(LeadTime != p.lead()) { - mlog << Error - << "\nTrackPoint::operator+=(const TrackPoint &p) -> " + mlog << Error << "\nTrackPoint::operator+=(const TrackPoint &p) -> " << "cannot call += for two different lead times (" << sec_to_hhmmss(LeadTime) << " != " << sec_to_hhmmss(p.lead()) << ").\n\n"; @@ -430,13 +426,15 @@ void TrackPoint::clear() { Depth = NoSystemsDepth; WarmCore = false; WatchWarn = NoWatchWarnType; - + NumMembers = bad_data_int; Spread = bad_data_double; DistMean = bad_data_double; VmaxStdev = bad_data_double; MSLPStdev = bad_data_double; - + + DiagVal.clear(); + // Call clear for each Wind object and then set intensity value for(i=0; i " + mlog << Error << "\nTrackPoint::operator[](int) -> " << "range check error for index value " << n << "\n\n"; exit(1); } @@ -630,6 +630,12 @@ void TrackPoint::set_watch_warn(WatchWarnType ww_type, unixtime ww_ut) { //////////////////////////////////////////////////////////////////////// +double TrackPoint::diag_val(int i) const { + return(i>=0 && i= NWinds)) { - mlog << Error - << "\nQuadInfo::set_wind(int, const QuadInfo) -> " + mlog << Error << "\nQuadInfo::set_wind(int, const QuadInfo) -> " << "range check error (" << n << ").\n\n"; exit(1); } @@ -679,3 +684,13 @@ bool TrackPoint::is_match(const ATCFTrackLine &l) const { } //////////////////////////////////////////////////////////////////////// + +void TrackPoint::add_diag_value(double val) { + + // Store the diagnostic value + DiagVal.add(val); + + return; +} + +//////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_tc_util/track_point.h b/src/libcode/vx_tc_util/track_point.h index 95c7fc3921..27ae223dae 100644 --- a/src/libcode/vx_tc_util/track_point.h +++ b/src/libcode/vx_tc_util/track_point.h @@ -167,13 +167,15 @@ class TrackPoint { // Wind Radii QuadInfo Wind[NWinds]; - // Consensus track variables - int NumMembers; - double Spread; // nautical miles - double DistMean; // nautical miles - double VmaxStdev; // knots - double MSLPStdev; // millibars - + // Consensus track variables + int NumMembers; + double Spread; // nautical miles + double DistMean; // nautical miles + double VmaxStdev; // knots + double MSLPStdev; // millibars + + // Diagnostic values + NumArray DiagVal; public: @@ -213,7 +215,7 @@ class TrackPoint { void set_warm_core(bool); void set_watch_warn(WatchWarnType); void set_watch_warn(WatchWarnType, unixtime); - + void set_num_members(const int); void set_spread(const double); void set_dist_mean(const double); @@ -249,6 +251,9 @@ class TrackPoint { double dist_mean() const; double v_max_stdev() const; double mslp_stdev() const; + + int n_diag() const; + double diag_val(int) const; // // do stuff @@ -257,6 +262,7 @@ class TrackPoint { void set_wind(int, const QuadInfo &); bool set(const ATCFTrackLine &); bool is_match(const ATCFTrackLine &) const; + void add_diag_value(double); }; @@ -286,31 +292,31 @@ inline void TrackPoint::set_dist_mean(const double v) { DistMean = v; } inline void TrackPoint::set_v_max_stdev(const double v) { VmaxStdev = v; } inline void TrackPoint::set_mslp_stdev(const double v) { MSLPStdev = v; } - -inline unixtime TrackPoint::valid() const { return(ValidTime); } +inline unixtime TrackPoint::valid() const { return(ValidTime); } inline int TrackPoint::valid_hour() const { return(unix_to_sec_of_day(ValidTime)); } -inline int TrackPoint::lead() const { return(LeadTime); } -inline double TrackPoint::lat() const { return(Lat); } -inline double TrackPoint::lon() const { return(Lon); } -inline double TrackPoint::v_max() const { return(Vmax); } -inline double TrackPoint::mslp() const { return(MSLP); } -inline CycloneLevel TrackPoint::level() const { return(Level); } -inline double TrackPoint::radp() const { return(RadP); } -inline double TrackPoint::rrp() const { return(RRP); } -inline double TrackPoint::mrd() const { return(MRD); } -inline double TrackPoint::gusts() const { return(Gusts); } -inline double TrackPoint::eye() const { return(Eye); } -inline double TrackPoint::direction() const { return(Direction); } -inline double TrackPoint::speed() const { return(Speed); } -inline SystemsDepth TrackPoint::depth() const { return(Depth); } -inline bool TrackPoint::warm_core() const { return(WarmCore); } -inline WatchWarnType TrackPoint::watch_warn() const { return(WatchWarn); } - -inline int TrackPoint::num_members() const { return(NumMembers); } -inline double TrackPoint::spread() const { return(Spread); } -inline double TrackPoint::dist_mean() const { return(DistMean); } -inline double TrackPoint::v_max_stdev() const { return(VmaxStdev); } -inline double TrackPoint::mslp_stdev() const { return(MSLPStdev); } +inline int TrackPoint::lead() const { return(LeadTime); } +inline double TrackPoint::lat() const { return(Lat); } +inline double TrackPoint::lon() const { return(Lon); } +inline double TrackPoint::v_max() const { return(Vmax); } +inline double TrackPoint::mslp() const { return(MSLP); } +inline CycloneLevel TrackPoint::level() const { return(Level); } +inline double TrackPoint::radp() const { return(RadP); } +inline double TrackPoint::rrp() const { return(RRP); } +inline double TrackPoint::mrd() const { return(MRD); } +inline double TrackPoint::gusts() const { return(Gusts); } +inline double TrackPoint::eye() const { return(Eye); } +inline double TrackPoint::direction() const { return(Direction); } +inline double TrackPoint::speed() const { return(Speed); } +inline SystemsDepth TrackPoint::depth() const { return(Depth); } +inline bool TrackPoint::warm_core() const { return(WarmCore); } +inline WatchWarnType TrackPoint::watch_warn() const { return(WatchWarn); } + +inline int TrackPoint::num_members() const { return(NumMembers); } +inline double TrackPoint::spread() const { return(Spread); } +inline double TrackPoint::dist_mean() const { return(DistMean); } +inline double TrackPoint::v_max_stdev() const { return(VmaxStdev); } +inline double TrackPoint::mslp_stdev() const { return(MSLPStdev); } +inline int TrackPoint::n_diag() const { return(DiagVal.n()); } //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_tc_util/vx_tc_util.h b/src/libcode/vx_tc_util/vx_tc_util.h index f0c1825cf2..902aa3f74c 100644 --- a/src/libcode/vx_tc_util/vx_tc_util.h +++ b/src/libcode/vx_tc_util/vx_tc_util.h @@ -14,6 +14,7 @@ //////////////////////////////////////////////////////////////////////// #include "atcf_track_line.h" +#include "diag_file.h" #include "track_point.h" #include "track_info.h" #include "track_pair_info.h" diff --git a/src/tools/core/stat_analysis/aggr_stat_line.cc b/src/tools/core/stat_analysis/aggr_stat_line.cc index e2061d30c2..d7e64316c6 100644 --- a/src/tools/core/stat_analysis/aggr_stat_line.cc +++ b/src/tools/core/stat_analysis/aggr_stat_line.cc @@ -2632,6 +2632,10 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, m[key].ens_pd.crps_gaus_na.add(cur.crps_gaus); m[key].ens_pd.crpscl_gaus_na.add(cur.crpscl_gaus); m[key].ens_pd.ign_na.add(cur.ign); + m[key].ens_pd.n_ge_obs_na.add(cur.n_ge_obs); + m[key].ens_pd.me_ge_obs_na.add(cur.me_ge_obs); + m[key].ens_pd.n_lt_obs_na.add(cur.n_lt_obs); + m[key].ens_pd.me_lt_obs_na.add(cur.me_lt_obs); m[key].ens_pd.var_na.add(square(cur.spread)); m[key].ens_pd.var_oerr_na.add(square(cur.spread_oerr)); m[key].ens_pd.var_plus_oerr_na.add(square(cur.spread_plus_oerr)); @@ -2668,18 +2672,18 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, it->second.ens_pd.n_pair = it->second.ens_pd.wgt_na.sum(); // Compute ME and RMSE as weighted averages - it->second.ens_pd.me = it->second.me_na.wmean(it->second.ens_pd.wgt_na); - v = it->second.mse_na.wmean(it->second.ens_pd.wgt_na); - it->second.ens_pd.rmse = (is_bad_data(v) ? bad_data_double : sqrt(v)); - it->second.ens_pd.me_oerr = it->second.me_oerr_na.wmean(it->second.ens_pd.wgt_na); - v = it->second.mse_oerr_na.wmean(it->second.ens_pd.wgt_na); - it->second.ens_pd.rmse_oerr = (is_bad_data(v) ? bad_data_double : sqrt(v)); - - crps_emp = it->second.ens_pd.crps_emp_na.wmean(it->second.ens_pd.wgt_na); - crps_emp_fair = it->second.ens_pd.crps_emp_fair_na.wmean(it->second.ens_pd.wgt_na); - crpscl_emp = it->second.ens_pd.crpscl_emp_na.wmean(it->second.ens_pd.wgt_na); - crps_gaus = it->second.ens_pd.crps_gaus_na.wmean(it->second.ens_pd.wgt_na); - crpscl_gaus = it->second.ens_pd.crpscl_gaus_na.wmean(it->second.ens_pd.wgt_na); + it->second.ens_pd.me = it->second.me_na.wmean(it->second.ens_pd.wgt_na); + v = it->second.mse_na.wmean(it->second.ens_pd.wgt_na); + it->second.ens_pd.rmse = (is_bad_data(v) ? bad_data_double : sqrt(v)); + it->second.ens_pd.me_oerr = it->second.me_oerr_na.wmean(it->second.ens_pd.wgt_na); + v = it->second.mse_oerr_na.wmean(it->second.ens_pd.wgt_na); + it->second.ens_pd.rmse_oerr = (is_bad_data(v) ? bad_data_double : sqrt(v)); + + crps_emp = it->second.ens_pd.crps_emp_na.wmean(it->second.ens_pd.wgt_na); + crps_emp_fair = it->second.ens_pd.crps_emp_fair_na.wmean(it->second.ens_pd.wgt_na); + crpscl_emp = it->second.ens_pd.crpscl_emp_na.wmean(it->second.ens_pd.wgt_na); + crps_gaus = it->second.ens_pd.crps_gaus_na.wmean(it->second.ens_pd.wgt_na); + crpscl_gaus = it->second.ens_pd.crpscl_gaus_na.wmean(it->second.ens_pd.wgt_na); // Compute aggregated empirical CRPSS it->second.ens_pd.crpss_emp = @@ -3195,6 +3199,17 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job, m[key].ens_pd.ign_na.add(compute_ens_ign(cur.obs, cur.ens_mean, cur.spread)); m[key].ens_pd.pit_na.add(compute_ens_pit(cur.obs, cur.ens_mean, cur.spread)); + // Store BIAS_RATIO terms + int n_ge_obs, n_lt_obs; + double me_ge_obs, me_lt_obs; + compute_bias_ratio_terms(cur.obs, cur.ens_na, + n_ge_obs, me_ge_obs, + n_lt_obs, me_lt_obs); + m[key].ens_pd.n_ge_obs_na.add(n_ge_obs); + m[key].ens_pd.me_ge_obs_na.add(me_ge_obs); + m[key].ens_pd.n_lt_obs_na.add(n_lt_obs); + m[key].ens_pd.me_lt_obs_na.add(me_lt_obs); + // // Increment the RHIST counts // diff --git a/src/tools/core/stat_analysis/parse_stat_line.cc b/src/tools/core/stat_analysis/parse_stat_line.cc index bf0c7fbc16..d39e4f4d27 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.cc +++ b/src/tools/core/stat_analysis/parse_stat_line.cc @@ -377,6 +377,12 @@ void parse_ecnt_line(STATLine &l, ECNTData &e_data) { e_data.spread_plus_oerr = atof(l.get_item("SPREAD_PLUS_OERR")); + e_data.bias_ratio = atof(l.get_item("BIAS_RATIO")); + e_data.n_ge_obs = atoi(l.get_item("N_GE_OBS")); + e_data.me_ge_obs = atof(l.get_item("ME_GE_OBS")); + e_data.n_lt_obs = atoi(l.get_item("N_LT_OBS")); + e_data.me_lt_obs = atof(l.get_item("ME_LT_OBS")); + return; } diff --git a/src/tools/core/stat_analysis/parse_stat_line.h b/src/tools/core/stat_analysis/parse_stat_line.h index 0ad26c65e9..8cfc6efdc7 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.h +++ b/src/tools/core/stat_analysis/parse_stat_line.h @@ -66,6 +66,9 @@ struct ECNTData { double ign, me, rmse, spread; double me_oerr, rmse_oerr, spread_oerr; double spread_plus_oerr; + double bias_ratio; + int n_ge_obs, n_lt_obs; + double me_ge_obs, me_lt_obs; }; // Ranked Histogram (RHIST) data structure diff --git a/src/tools/other/ascii2nc/airnow_locations.cc b/src/tools/other/ascii2nc/airnow_locations.cc index 10810c8727..17ddcdb666 100644 --- a/src/tools/other/ascii2nc/airnow_locations.cc +++ b/src/tools/other/ascii2nc/airnow_locations.cc @@ -150,6 +150,7 @@ bool AirnowLocations::lookupLatLonElev(const string aqsid, double &lat, double & { string method_name = "AirnowLocations::lookupLatLonElev()"; + int index = -1; vector::const_iterator it = find(monitoringSiteAqsid.begin(), monitoringSiteAqsid.end(), aqsid); if (it == monitoringSiteAqsid.end()) { it = find(monitoringSiteStationId.begin(), monitoringSiteStationId.end(), aqsid); @@ -157,10 +158,19 @@ bool AirnowLocations::lookupLatLonElev(const string aqsid, double &lat, double & it = find(monitoringSiteFullAqsid.begin(), monitoringSiteFullAqsid.end(), aqsid); if (it == monitoringSiteFullAqsid.end()) { return false; + } else { + index = (int)(it - monitoringSiteStationId.begin()); } + } else { + index = (int)(it - monitoringSiteFullAqsid.begin()); } + } else { + index = (int)(it - monitoringSiteAqsid.begin()); + } + if (index < 0) { + return false; } - int index = (int)(it - monitoringSiteAqsid.begin()); + lat = monitoringSiteLat[index]; lon = monitoringSiteLon[index]; elev = monitoringSiteElev[index]; diff --git a/src/tools/other/grid_diag/grid_diag.cc b/src/tools/other/grid_diag/grid_diag.cc index 96e3be2a94..1ffa59d028 100644 --- a/src/tools/other/grid_diag/grid_diag.cc +++ b/src/tools/other/grid_diag/grid_diag.cc @@ -20,6 +20,7 @@ // 003 08/20/21 Halley Gotway Bugfix #1886 for integer overflow. // 004 07/06/22 Howard Soh METplus-Internal #19 Rename main to met_main // 005 10/03/22 Prestopnik MET #2227 Remove using namespace std and netCDF from header files +// 006 10/26/22 Linden MET #2232 Refine the Grid-Diag output variable names when specifying two input data sources // //////////////////////////////////////////////////////////////////////// @@ -264,7 +265,7 @@ void process_series(void) { // Process the 1d histograms for(int i_var=0; i_varmagic_str_attr() << " histogram with " << n_bins << " bins from " << min << " to " << max << ".\n"; + histograms[i_var_str] = vector(); init_pdf(n_bins, histograms[i_var_str]); + + // Keep track of unique output variable names + if(nc_var_sa.has( data_info->magic_str_attr() )) unique_variable_names = false; + nc_var_sa.add(data_info->magic_str_attr()); + } // for i_var } @@ -455,14 +462,14 @@ void setup_joint_histograms(void) { for(int i_var=0; i_varn_bins(); for(int j_var=i_var+1; j_varn_bins(); @@ -484,48 +491,48 @@ void setup_joint_histograms(void) { //////////////////////////////////////////////////////////////////////// void setup_nc_file(void) { - ConcatString cs, i_var_str; + ConcatString cs, i_var_str, j_var_str; - // Create NetCDF file - nc_out = open_ncfile(out_file.c_str(), true); + // Create NetCDF file + nc_out = open_ncfile(out_file.c_str(), true); - if(IS_INVALID_NC_P(nc_out)) { + if(IS_INVALID_NC_P(nc_out)) { mlog << Error << "\nsetup_nc_file() -> " - << "trouble opening output NetCDF file " - << out_file << "\n\n"; - exit(1); - } - - // Add global attributes - write_netcdf_global(nc_out, out_file.c_str(), program_name, - NULL, NULL, conf_info.desc.c_str()); - add_att(nc_out, "mask_grid", (conf_info.mask_grid_name.nonempty() ? - (string)conf_info.mask_grid_name : - na_str)); - add_att(nc_out, "mask_poly", (conf_info.mask_poly_name.nonempty() ? - (string)conf_info.mask_poly_name : - na_str)); - - // Add time range information to the global attributes - add_att(nc_out, "init_beg", (string)unix_to_yyyymmdd_hhmmss(init_beg)); - add_att(nc_out, "init_end", (string)unix_to_yyyymmdd_hhmmss(init_end)); - add_att(nc_out, "valid_beg", (string)unix_to_yyyymmdd_hhmmss(valid_beg)); - add_att(nc_out, "valid_end", (string)unix_to_yyyymmdd_hhmmss(valid_end)); - add_att(nc_out, "lead_beg", (string)sec_to_hhmmss(lead_beg)); - add_att(nc_out, "lead_end", (string)sec_to_hhmmss(lead_end)); - - // Write the grid size, mask size, and series length - write_nc_var_int("grid_size", "number of grid points", grid.nxy()); - write_nc_var_int("mask_size", "number of mask points", conf_info.mask_area.count()); - write_nc_var_int("n_series", "length of series", n_series); - - // Compression level - int deflate_level = compress_level; - if(deflate_level < 0) deflate_level = conf_info.conf.nc_compression(); - - for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - - i_var_str << cs_erase << "VAR" << i_var; + << "trouble opening output NetCDF file " + << out_file << "\n\n"; + exit(1); + } + + // Add global attributes + write_netcdf_global(nc_out, out_file.c_str(), program_name, + NULL, NULL, conf_info.desc.c_str()); + add_att(nc_out, "mask_grid", (conf_info.mask_grid_name.nonempty() ? + (string)conf_info.mask_grid_name : + na_str)); + add_att(nc_out, "mask_poly", (conf_info.mask_poly_name.nonempty() ? + (string)conf_info.mask_poly_name : + na_str)); + + // Add time range information to the global attributes + add_att(nc_out, "init_beg", (string)unix_to_yyyymmdd_hhmmss(init_beg)); + add_att(nc_out, "init_end", (string)unix_to_yyyymmdd_hhmmss(init_end)); + add_att(nc_out, "valid_beg", (string)unix_to_yyyymmdd_hhmmss(valid_beg)); + add_att(nc_out, "valid_end", (string)unix_to_yyyymmdd_hhmmss(valid_end)); + add_att(nc_out, "lead_beg", (string)sec_to_hhmmss(lead_beg)); + add_att(nc_out, "lead_end", (string)sec_to_hhmmss(lead_end)); + + // Write the grid size, mask size, and series length + write_nc_var_int("grid_size", "number of grid points", grid.nxy()); + write_nc_var_int("mask_size", "number of mask points", conf_info.mask_area.count()); + write_nc_var_int("n_series", "length of series", n_series); + + // Compression level + int deflate_level = compress_level; + if(deflate_level < 0) deflate_level = conf_info.conf.nc_compression(); + + for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { + + i_var_str << cs_erase << "VAR" << i_var+1; VarInfo *data_info = conf_info.data_info[i_var]; @@ -534,11 +541,16 @@ void setup_nc_file(void) { var_name.add("_"); var_name.add(data_info->level_attr()); + if(multiple_data_sources && !unique_variable_names) { + var_name.add("_"); + var_name.add(i_var_str); + } + // Define histogram dimensions NcDim var_dim = add_dim(nc_out, var_name, (long) data_info->n_bins()); data_var_dims.push_back(var_dim); - + // Define histogram bins ConcatString var_min_name = var_name; ConcatString var_max_name = var_name; @@ -575,6 +587,8 @@ void setup_nc_file(void) { // Define histograms for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { + i_var_str << cs_erase << "VAR" << i_var+1; + VarInfo *data_info = conf_info.data_info[i_var]; // Set variable NetCDF name @@ -582,6 +596,11 @@ void setup_nc_file(void) { var_name.add("_"); var_name.add(data_info->level_attr()); + if(multiple_data_sources && !unique_variable_names) { + var_name.add("_"); + var_name.add(i_var_str); + } + ConcatString hist_name("hist_"); hist_name.add(var_name); NcDim var_dim = data_var_dims[i_var]; @@ -597,21 +616,36 @@ void setup_nc_file(void) { // Define joint histograms for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { + i_var_str << cs_erase << "VAR" << i_var+1; + VarInfo *data_info = conf_info.data_info[i_var]; for(int j_var=i_var+1; j_varname_attr()); hist_name.add("_"); hist_name.add(data_info->level_attr()); + + if(multiple_data_sources && !unique_variable_names) { + hist_name.add("_"); + hist_name.add(i_var_str); + } + hist_name.add("_"); hist_name.add(joint_info->name_attr()); hist_name.add("_"); hist_name.add(joint_info->level_attr()); + if(multiple_data_sources && !unique_variable_names) { + hist_name.add("_"); + hist_name.add(j_var_str); + } + NcDim var_dim = data_var_dims[i_var]; NcDim joint_dim = data_var_dims[j_var]; vector dims; @@ -662,7 +696,7 @@ void write_histograms(void) { for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - i_var_str << cs_erase << "VAR" << i_var; + i_var_str << cs_erase << "VAR" << i_var+1; VarInfo *data_info = conf_info.data_info[i_var]; NcVar hist_var = hist_vars[i_var]; @@ -690,8 +724,8 @@ void write_joint_histograms(void) { VarInfo *joint_info = conf_info.data_info[j_var]; ij_var_str << cs_erase - << "VAR" << i_var << "_" - << "VAR" << j_var; + << "VAR" << i_var+1 << "_" + << "VAR" << j_var+1; long long *hist = joint_histograms[ij_var_str].data(); @@ -816,6 +850,7 @@ void usage() { void set_data_files(const StringArray & a) { data_files.push_back(a); + if(data_files.size() > 0) multiple_data_sources = true; } //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/other/grid_diag/grid_diag.h b/src/tools/other/grid_diag/grid_diag.h index 18b18d3326..91f6bcc797 100644 --- a/src/tools/other/grid_diag/grid_diag.h +++ b/src/tools/other/grid_diag/grid_diag.h @@ -90,6 +90,12 @@ vector data_var_dims; vector hist_vars; vector joint_hist_vars; +static bool multiple_data_sources = false; +static bool unique_variable_names = true; + +// List of output NetCDF variable names +static StringArray nc_var_sa; + //////////////////////////////////////////////////////////////////////// // // Miscellaneous Variables diff --git a/src/tools/other/ioda2nc/ioda2nc.cc b/src/tools/other/ioda2nc/ioda2nc.cc index 1aaef67ba1..ad0c90353f 100644 --- a/src/tools/other/ioda2nc/ioda2nc.cc +++ b/src/tools/other/ioda2nc/ioda2nc.cc @@ -63,6 +63,15 @@ static const char * DEF_CONFIG_NAME = "MET_BASE/config/IODA2NCConfig_default"; static const char *program_name = "ioda2nc"; static const int REJECT_DEBUG_LEVEL = 9; +static const int string_data_len = 512; + +static const char *metadata_group_name = "MetaData"; +static const char *qc_group_name = "QCFlags"; +static const char *qc_postfix = "PreQC"; +static const char *obs_group_name = "ObsValue"; +static const char *derived_obs_group_name = "DerivedObsValue"; + +enum e_ioda_format { ioda_v1, ioda_v2 }; //////////////////////////////////////////////////////////////////////// @@ -74,7 +83,8 @@ static const int REJECT_DEBUG_LEVEL = 9; static StringArray ioda_files; static StringArray core_dims; -static StringArray core_vars; +static StringArray core_dims_v1; +static StringArray core_meta_vars; // Output NetCDF file name static ConcatString ncfile; @@ -154,13 +164,15 @@ static void set_valid_beg_time(const StringArray &); static void set_valid_end_time(const StringArray &); static void set_verbosity(const StringArray &); -static bool check_core_data(const bool, const bool, StringArray &, StringArray &); +static bool check_core_data(const bool, const bool, StringArray &, StringArray &, e_ioda_format); static bool check_missing_thresh(float value); static ConcatString find_meta_name(StringArray, StringArray); -static bool get_meta_data_float(NcFile *, StringArray &, const char *, float *, const int); -static bool get_meta_data_strings(NcFile *, const ConcatString, char *); -static bool get_obs_data_float(NcFile *, const ConcatString, - NcVar *, float *, int *, const int); +static bool get_meta_data_float(NcFile *, StringArray &, const char *, float *, + const int); +static bool get_meta_data_strings(NcVar &, char *); +static bool get_meta_data_strings(NcVar &, char **); +static bool get_obs_data_float(NcFile *, const ConcatString, NcVar *, + float *, int *, const int, const e_ioda_format); static bool has_postfix(std::string const &, std::string const &); //////////////////////////////////////////////////////////////////////// @@ -214,15 +226,18 @@ void initialize() { nc_point_obs.init_buffer(); core_dims.clear(); - core_dims.add("nvars"); core_dims.add("nlocs"); - //core_dims.add("nstring"); - //core_dims.add("ndatetime"); + + core_dims_v1.clear(); + core_dims_v1.add("nvars"); + core_dims_v1.add("nlocs"); + core_dims_v1.add("nstring"); + //core_dims_v1.add("ndatetime"); - core_vars.clear(); - core_vars.add("datetime"); - core_vars.add("latitude"); - core_vars.add("longitude"); + core_meta_vars.clear(); + core_meta_vars.add("datetime"); + core_meta_vars.add("latitude"); + core_meta_vars.add("longitude"); summary_obs = new SummaryObs(); return; @@ -349,6 +364,9 @@ void process_ioda_file(int i_pb) { int rej_elv, rej_nobs; double x, y; + bool status; + bool is_time_offset = false; + bool is_time_string = false; unixtime file_ut; unixtime adjusted_file_ut; unixtime msg_ut, beg_ut, end_ut; @@ -411,39 +429,76 @@ void process_ioda_file(int i_pb) { file_name << ioda_files[i_pb]; int nrecs = 0; - StringArray var_names, dim_names; + int nstring, nvars; + StringArray dim_names; StringArray metadata_vars; StringArray obs_value_vars; - get_var_names(f_in, &var_names); + bool error_out = true; + e_ioda_format ioda_format = ioda_v2; + int nlocs = get_dim_value(f_in, "nlocs", error_out); // number of locations + + nvars = bad_data_int ; + nstring = string_data_len; + if (! has_nc_group(f_in, obs_group_name)) ioda_format = ioda_v1; + get_dim_names(f_in, &dim_names); - for(idx=0; idx= 8) { + for(idx=0; idx= 6) { - for(idx=0; idx 0 && nmsg < npbmsg) { npbmsg = (nmsg_percent > 0 && nmsg_percent <= 100) - ? (npbmsg * nmsg_percent / 100) : nmsg; + ? (npbmsg * nmsg_percent / 100) : nmsg; } - long lengths[2] = { nlocs, ndatetime }; - long offsets[2] = { 0, 0 }; float *hdr_lat_arr = new float[nlocs]; float *hdr_lon_arr = new float[nlocs]; float *hdr_elv_arr = new float[nlocs]; float *obs_pres_arr = new float[nlocs]; float *obs_hght_arr = new float[nlocs]; + float *hdr_time_arr = new float[nlocs]; char *hdr_vld_block = new char[nlocs*ndatetime]; - char *hdr_msg_types = 0; - char *hdr_station_ids = 0; + char *hdr_msg_types = NULL; + char *hdr_station_ids = NULL; + char **hdr_vld_block2 = NULL; + char **hdr_msg_types2 = NULL; + char **hdr_station_ids2 = NULL; vector v_qc_data; vector v_obs_data; + + if (is_time_string) { + hdr_vld_block2 = (char**) calloc(nlocs, sizeof(char*)); + for (int i=0; i= debug_level_for_performance) { @@ -605,6 +711,7 @@ void process_ioda_file(int i_pb) { for(idx=0; idx 0) { @@ -620,11 +727,25 @@ void process_ioda_file(int i_pb) { } } - char valid_time[ndatetime+1]; - m_strncpy(valid_time, (const char *)(hdr_vld_block + (i_read * ndatetime)), - ndatetime, method_name_s, "valid_time", true); - valid_time[ndatetime] = 0; - msg_ut = yyyymmddThhmmss_to_unix(valid_time); + if (is_time_offset) { + msg_ut = add_to_unixtime(base_ut, sec_per_unit, + hdr_time_arr[i_read], no_leap_year); + } + else if (is_time_string) { + char valid_time[nstring+1]; + + m_strncpy(valid_time, (const char *)hdr_vld_block2[i_read], + nstring, method_name_s, "valid_time", true); + valid_time[nstring] = 0; + msg_ut = yyyymmddThhmmss_to_unix(valid_time); + } + else { + char valid_time[ndatetime+1]; + m_strncpy(valid_time, (const char *)(hdr_vld_block + (i_read * ndatetime)), + ndatetime, method_name_s, "valid_time", true); + valid_time[ndatetime] = 0; + msg_ut = yyyymmddThhmmss_to_unix(valid_time); + } // Check to make sure that the message time hasn't changed // from one IODA message to the next @@ -638,8 +759,7 @@ void process_ioda_file(int i_pb) { // Check if valid_beg_ut and valid_end_ut were set on the // command line. If so, use them. If not, use beg_ds and // end_ds. - if(valid_beg_ut != (unixtime) 0 || - valid_end_ut != (unixtime) 0) { + if(valid_beg_ut != (unixtime) 0 || valid_end_ut != (unixtime) 0) { beg_ut = valid_beg_ut; end_ut = valid_end_ut; } @@ -671,7 +791,13 @@ void process_ioda_file(int i_pb) { } if(has_msg_type) { - m_strncpy(hdr_typ, hdr_msg_types+(i_read*nstring), nstring, method_name_s, "hdr_typ"); + if (NULL != hdr_msg_types2) { + m_strncpy(hdr_typ, hdr_msg_types2[i_read], nstring, method_name_s, "hdr_typ2"); + } + else { + m_strncpy(hdr_typ, hdr_msg_types+(i_read*nstring), nstring, method_name_s, "hdr_typ"); + + } m_rstrip(hdr_typ, nstring); // If the message type is not listed in the configuration @@ -697,7 +823,12 @@ void process_ioda_file(int i_pb) { if(has_station_id) { char tmp_sid[nstring+1]; - m_strncpy(tmp_sid, hdr_station_ids+(i_read*nstring), nstring, method_name_s, "tmp_sid"); + if (NULL != hdr_station_ids2) { + m_strncpy(tmp_sid, hdr_station_ids2[i_read], nstring, method_name_s, "tmp_sid2"); + } + else { + m_strncpy(tmp_sid, hdr_station_ids+(i_read*nstring), nstring, method_name_s, "tmp_sid"); + } m_rstrip(tmp_sid, nstring, false); m_replace_char(tmp_sid, ' ', '_'); hdr_sid = tmp_sid; @@ -881,6 +1012,19 @@ void process_ioda_file(int i_pb) { delete [] obs_hght_arr; if (hdr_msg_types) delete [] hdr_msg_types; if (hdr_station_ids) delete [] hdr_station_ids; + if (NULL != hdr_msg_types2) { + for (int i=0; i " - << "core dimension \"" << core_dims[idx] << "\" is missing.\n\n"; + << "core dimension \"" << t_core_dims[idx] << "\" is missing.\n\n"; is_netcdf_ready = false; } } - if(has_msg_type || has_station_id) { - if(!dim_names.has("nstring")) { - mlog << Error << "\n" << method_name << "-> " - << "core dimension \"nstring\" is missing.\n\n"; - is_netcdf_ready = false; + if (ioda_format == ioda_v1) { + if(has_msg_type || has_station_id) { + if(!dim_names.has("nstring")) { + mlog << Error << "\n" << method_name << "-> " + << "core dimension \"nstring\" is missing.\n\n"; + is_netcdf_ready = false; + } } } - for(int idx=0; idx 0) { + for (int idx2=0; idx2 < alt_names.n(); idx2++) { + if (core_meta_vars[idx] != alt_names[idx2]) { + found = metadata_vars.has(alt_names[idx2]); + if (found) break; + } + } + } + } + if(!found) { mlog << Error << "\n" << method_name << "-> " - << "core variable \"" << core_vars[idx] << "\" is missing.\n\n"; + << "core variable \"" << core_meta_vars[idx] << "\" is missing.\n\n"; is_netcdf_ready = false; } } @@ -1123,10 +1284,9 @@ bool get_meta_data_float(NcFile *f_in, StringArray &metadata_vars, ConcatString metadata_name = find_meta_name( metadata_vars, conf_info.metadata_map[metadata_key]); + if(metadata_name.length() > 0) { - ConcatString ioda_name = metadata_name; - ioda_name.add("@MetaData"); - NcVar meta_var = get_var(f_in, ioda_name.c_str()); + NcVar meta_var = get_var(f_in, metadata_name.c_str(), metadata_group_name); if(IS_VALID_NC(meta_var)) { status = get_nc_data(&meta_var, metadata_buf, nlocs); if(!status) mlog << Debug(3) << method_name @@ -1148,19 +1308,15 @@ bool get_meta_data_float(NcFile *f_in, StringArray &metadata_vars, //////////////////////////////////////////////////////////////////////// -bool get_meta_data_strings(NcFile *f_in, const ConcatString metadata_name, - char *metadata_buf) { +bool get_meta_data_strings(NcVar &var, char *metadata_buf) { bool status = false; static const char *method_name = "get_meta_data_strings() -> "; - ConcatString ioda_name = metadata_name; - ioda_name.add("@MetaData"); - NcVar meta_var = get_var(f_in, ioda_name.c_str()); - if(IS_VALID_NC(meta_var)) { - status = get_nc_data(&meta_var, metadata_buf); + if(IS_VALID_NC(var)) { + status = get_nc_data(&var, metadata_buf); if(!status) { mlog << Error << "\n" << method_name << " -> " - << "trouble getting " << metadata_name << "\n\n"; + << "trouble getting " << GET_NC_NAME(var) << "\n\n"; exit(1); } } @@ -1171,7 +1327,7 @@ bool get_meta_data_strings(NcFile *f_in, const ConcatString metadata_name, bool get_obs_data_float(NcFile *f_in, const ConcatString var_name, NcVar *obs_var, float *obs_buf, int *qc_buf, - const int nlocs) { + const int nlocs, const e_ioda_format ioda_format) { bool status = false; static const char *method_name = "get_obs_data_float() -> "; @@ -1185,21 +1341,37 @@ bool get_obs_data_float(NcFile *f_in, const ConcatString var_name, << "trouble getting " << var_name << "\n\n"; } else mlog << Error << "\n" << method_name - << var_name << "@ObsValue does not exist!\n\n"; + << var_name << " does not exist!\n\n"; if(!status) exit(1); status = false; if(var_name.length() > 0) { - ConcatString ioda_name = var_name; - ioda_name.add("@PreQC"); - NcVar qc_var = get_var(f_in, ioda_name.c_str()); + ConcatString qc_name = var_name; + ConcatString qc_group = qc_postfix; + if (ioda_format == ioda_v2) { + NcGroup nc_grp = get_nc_group(f_in, qc_postfix); + if (IS_INVALID_NC(nc_grp)) qc_group = qc_group_name; + StringArray qc_names = conf_info.obs_to_qc_map[var_name]; + if (0 < qc_names.n()) { + for (int idx=0; idx obs_name_map; map message_type_map; map metadata_map; + map obs_to_qc_map; StringArray surface_message_types; TimeSummaryInfo timeSummaryInfo; diff --git a/src/tools/tc_utils/tc_gen/Makefile.am b/src/tools/tc_utils/tc_gen/Makefile.am index f66a1577ef..529bb89e3f 100644 --- a/src/tools/tc_utils/tc_gen/Makefile.am +++ b/src/tools/tc_utils/tc_gen/Makefile.am @@ -30,7 +30,7 @@ tc_gen_LDADD = -lvx_stat_out \ -lvx_statistics \ -lvx_gis \ -lvx_data2d \ - -lvx_seeps \ + -lvx_seeps \ -lvx_nc_util \ -lvx_regrid \ -lvx_grid \ diff --git a/src/tools/tc_utils/tc_gen/Makefile.in b/src/tools/tc_utils/tc_gen/Makefile.in index 0bda9b9032..35cd8d681e 100644 --- a/src/tools/tc_utils/tc_gen/Makefile.in +++ b/src/tools/tc_utils/tc_gen/Makefile.in @@ -337,10 +337,10 @@ tc_gen_LDADD = -lvx_stat_out \ -lvx_statistics \ -lvx_gis \ -lvx_data2d \ + -lvx_seeps \ -lvx_nc_util \ -lvx_regrid \ -lvx_grid \ - -lvx_seeps \ -lvx_config \ -lvx_pb_util \ -lvx_cal \ diff --git a/src/tools/tc_utils/tc_pairs/tc_pairs.cc b/src/tools/tc_utils/tc_pairs/tc_pairs.cc index 70469e584d..1861e16eb0 100644 --- a/src/tools/tc_utils/tc_pairs/tc_pairs.cc +++ b/src/tools/tc_utils/tc_pairs/tc_pairs.cc @@ -35,6 +35,7 @@ // the gen_vx_mask tool. // 012 07/06/22 Howard Soh METplus-Internal #19 Rename main to met_main // 013 09/28/22 Prestopnik MET #2227 Remove namespace std from header files +// 014 10/06/22 Halley Gotway MET #392 Incorporate diagnostics // //////////////////////////////////////////////////////////////////////// @@ -123,6 +124,8 @@ static void process_prob_files (const StringArray &, static bool is_keeper (const ATCFLineBase *); static void filter_tracks (TrackInfoArray &); static void filter_probs (ProbInfoArray &); +static void process_diags (TrackInfoArray &); + static bool check_masks (const MaskPoly &, const Grid &, const MaskPlane &, double lat, double lon); @@ -147,8 +150,9 @@ static void usage (); static void set_adeck (const StringArray &); static void set_edeck (const StringArray &); static void set_bdeck (const StringArray &); -static void set_atcf_source (const StringArray &, +static void set_atcf_path (const StringArray &, StringArray &, StringArray &); +static void set_diag (const StringArray &); static void set_config (const StringArray &); static void set_out (const StringArray &); @@ -194,23 +198,21 @@ void process_command_line(int argc, char **argv) { cline.set_usage(usage); // Add function calls for the arguments - cline.add(set_adeck, "-adeck", -1); - cline.add(set_edeck, "-edeck", -1); - cline.add(set_bdeck, "-bdeck", -1); - cline.add(set_config, "-config", 1); - cline.add(set_out, "-out", 1); + cline.add(set_adeck, "-adeck", -1); + cline.add(set_edeck, "-edeck", -1); + cline.add(set_bdeck, "-bdeck", -1); + cline.add(set_diag, "-diag", -1); + cline.add(set_config, "-config", 1); + cline.add(set_out, "-out", 1); // Parse the command line cline.parse(); // Check for the minimum number of arguments - if((adeck_source.n() == 0 && - edeck_source.n() == 0) || - bdeck_source.n() == 0 || - config_file.length() == 0) { - mlog << Error - << "\nprocess_command_line(int argc, char **argv) -> " - << "You must specify at least one source of ADECK or EDECK " + if((adeck_path.n() == 0 && edeck_path.n() == 0) || + bdeck_path.n() == 0 || config_file.length() == 0) { + mlog << Error << "\nprocess_command_line(int argc, char **argv) -> " + << "You must specify at least one path of ADECK or EDECK " << "data, BDECK data, and the config file using the " << "\"-adeck\", \"-edeck\", \"-bdeck\", and \"-config\" " << "command line options.\n\n"; @@ -218,29 +220,38 @@ void process_command_line(int argc, char **argv) { } // List the input ADECK track files - for(i=0; i 0) { + if(adeck_path.n() > 0) { process_adecks(bdeck_tracks); } // Process EDECK files - if(edeck_source.n() > 0) { + if(edeck_path.n() > 0) { process_edecks(bdeck_tracks); } @@ -290,7 +301,7 @@ void process_bdecks(TrackInfoArray &bdeck_tracks) { bdeck_tracks.clear(); // Get the list of track files - get_atcf_files(bdeck_source, bdeck_model_suffix, + get_atcf_files(bdeck_path, bdeck_model_suffix, files, files_model_suffix); mlog << Debug(2) @@ -313,7 +324,7 @@ void process_adecks(const TrackInfoArray &bdeck_tracks) { int i, j, n_match; // Get the list of track files - get_atcf_files(adeck_source, adeck_model_suffix, + get_atcf_files(adeck_path, adeck_model_suffix, files, files_model_suffix); mlog << Debug(2) @@ -364,6 +375,9 @@ void process_adecks(const TrackInfoArray &bdeck_tracks) { << " ADECK tracks based on config file settings.\n"; filter_tracks(adeck_tracks); + // Append diagnostic data to the tracks + process_diags(adeck_tracks); + // // Loop through the ADECK tracks and find a matching BDECK track // @@ -430,7 +444,7 @@ void process_edecks(const TrackInfoArray &bdeck_tracks) { int n_match, i, j; // Get the list of ATCF files - get_atcf_files(edeck_source, edeck_model_suffix, + get_atcf_files(edeck_path, edeck_model_suffix, files, files_model_suffix); mlog << Debug(2) @@ -510,17 +524,16 @@ void process_edecks(const TrackInfoArray &bdeck_tracks) { //////////////////////////////////////////////////////////////////////// -void get_atcf_files(const StringArray &source, +void get_atcf_files(const StringArray &path, const StringArray &model_suffix, StringArray &files, StringArray &files_model_suffix) { - StringArray cur_source, cur_files; + StringArray cur_path, cur_files; int i, j; - if(source.n() != model_suffix.n()) { - mlog << Error - << "\nget_atcf_files() -> " - << "the source and suffix arrays must be equal length!\n\n"; + if(path.n() != model_suffix.n()) { + mlog << Error << "\nget_atcf_files() -> " + << "the path and suffix arrays must be equal length!\n\n"; exit(1); } @@ -529,10 +542,10 @@ void get_atcf_files(const StringArray &source, files_model_suffix.clear(); // Build list of files and corresponding model suffix list - for(i=0; i " + mlog << Error << "\nprocess_track_files() -> " << "unable to open file \"" << files[i] << "\"\n\n"; exit(1); } @@ -664,8 +676,7 @@ void process_prob_files(const StringArray &files, // Open the current file if(!f.open(files[i].c_str())) { - mlog << Error - << "\nprocess_prob_files() -> " + mlog << Error << "\nprocess_prob_files() -> " << "unable to open file \"" << files[i] << "\"\n\n"; exit(1); } @@ -995,12 +1006,76 @@ void filter_probs(ProbInfoArray &probs) { //////////////////////////////////////////////////////////////////////// +void process_diags(TrackInfoArray &tracks) { + StringArray cur_path, cur_model_name; + StringArray files, files_model_name, model_names; + DiagFile diag_file; + int i, j, n; + map n_diag_map; + + // Process the diagnostic inputs + for(i=0; i * convert_ptr = 0; + if(conf_info.DiagConvertMap.count(diag_source[i]) > 0) { + convert_ptr = &conf_info.DiagConvertMap.at(diag_source[i]); + } + + // Process the diagnostic files + diag_file.read(diag_source[i], files[j], model_names, convert_ptr); + + // Add diagnostics data to existing tracks + n += tracks.add_diag_data(diag_file, conf_info.DiagName); + + } // end for j + + // Update the diagnostic track counts + if(n_diag_map.count(diag_source[i]) > 0) n_diag_map[diag_source[i]] += n; + else n_diag_map[diag_source[i]] = n; + + } // end for i + + // Print the diagnostic track counts + for(map::iterator it = n_diag_map.begin(); + it != n_diag_map.end(); it++) { + mlog << Debug(3) + << "Added " << diagtype_to_string(it->first) + << " diagnostics to " << it->second << " tracks.\n"; + } + + return; +} + +//////////////////////////////////////////////////////////////////////// + bool check_masks(const MaskPoly &mask_poly, const Grid &mask_grid, const MaskPlane &mask_area, double lat, double lon) { double grid_x, grid_y; // - // Check polyline masking. + // Check polyline masking // if(mask_poly.n_points() > 0) { if(!mask_poly.latlon_is_inside_dege(lat, lon)) { @@ -1009,7 +1084,7 @@ bool check_masks(const MaskPoly &mask_poly, const Grid &mask_grid, } // - // Check grid masking. + // Check grid masking // if(mask_grid.nx() > 0 || mask_grid.ny() > 0) { mask_grid.latlon_to_xy(lat, -1.0*lon, grid_x, grid_y); @@ -1019,7 +1094,7 @@ bool check_masks(const MaskPoly &mask_poly, const Grid &mask_grid, } // - // Check area mask. + // Check area mask // if(mask_area.nx() > 0 || mask_area.ny() > 0) { if(!mask_area.s_is_on(nint(grid_x), nint(grid_y))) { @@ -1451,8 +1526,7 @@ void derive_baseline_model(const ConcatString &model, // Check bounds if(i_start < 0 || i_start >= ti.n_points()) { - mlog << Error - << "\n" << method_name + mlog << Error << "\n" << method_name << "range check error for i_start = " << i_start << "\n\n"; exit(1); } @@ -1568,8 +1642,7 @@ void derive_baseline_model(const ConcatString &model, bl_lat, bl_lon, bl_vmax); } else { - mlog << Error - << "\n" << method_name + mlog << Error << "\n" << method_name << "unsupported baseline model type \"" << model << "\".\n\n"; exit(1); @@ -1811,8 +1884,7 @@ void compute_track_err(const TrackInfo &adeck, const TrackInfo &bdeck, // Check for too many track points if(n_ut > mxp) { - mlog << Error - << "\ncompute_track_err() -> " + mlog << Error << "\ncompute_track_err() -> " << "exceeded the maximum number of allowable track points (" << n_ut << " > " << mxp << ")\n\n"; exit(1); @@ -1955,7 +2027,7 @@ void process_watch_warn(TrackPairInfoArray &p) { //////////////////////////////////////////////////////////////////////// void write_tracks(const TrackPairInfoArray &p) { - int i_row, i; + int i_row, i, n_row, n_col, n_diag; TcHdrColumns tchc; ConcatString out_file; AsciiTable out_at; @@ -1968,8 +2040,18 @@ void write_tracks(const TrackPairInfoArray &p) { // Create the output file open_tc_txt_file(out, out_file.c_str()); + // Number of output rows and columns + n_row = p.n_points() + 1; + n_col = n_tc_mpr_cols; + + // Update rows and columns for diagnosics output + if((n_diag = p.max_n_diag()) > 0) { + n_row += p.n_points(); + n_col = max(n_col, get_n_tc_diag_cols(n_diag)); + } + // Initialize the output AsciiTable - out_at.set_size(p.n_points() + 1, n_tc_header_cols + n_tc_mpr_cols); + out_at.set_size(n_row, n_tc_header_cols + n_col); setup_table(out_at); // Write the header row @@ -1999,7 +2081,7 @@ void write_tracks(const TrackPairInfoArray &p) { tchc.set_storm_name(p[i].bdeck().storm_name()); // Write the current TrackPairInfo object - write_tc_mpr_row(tchc, p[i], out_at, i_row); + write_track_pair_info(tchc, p[i], out_at, i_row); } // Write the AsciiTable contents and clean up @@ -2077,7 +2159,7 @@ void write_prob_rirw(const ProbRIRWPairInfoArray &p) { tchc.set_storm_name(p[i].bdeck()->storm_name()); // Write the current ProbRIRWPairInfo object - write_prob_rirw_row(tchc, p[i], out_at, i_row); + write_prob_rirw_pair_info(tchc, p[i], out_at, i_row); } // Write the AsciiTable contents and clean up @@ -2121,24 +2203,25 @@ void usage() { << ") ***\n\n" << "Usage: " << program_name << "\n" - << "\t-adeck source and/or -edeck source\n" - << "\t-bdeck source\n" + << "\t-adeck path and/or -edeck path\n" + << "\t-bdeck path\n" << "\t-config file\n" + << "\t[-diag source path]\n" << "\t[-out base]\n" << "\t[-log file]\n" << "\t[-v level]\n\n" - << "\twhere\t\"-adeck source\" is used one or more times to " + << "\twhere\t\"-adeck path\" is used one or more times to " << "specify a file or top-level directory containing ATCF " << "model output \"" << atcf_suffix << "\" data to process (required if no -edeck).\n" - << "\t\t\"-edeck source\" is used one or more times to " + << "\t\t\"-edeck path\" is used one or more times to " << "specify a file or top-level directory containing ATCF " << "ensemble model output \"" << atcf_suffix << "\" data to process (required if no -adeck).\n" - << "\t\t\"-bdeck source\" is used one or more times to " + << "\t\t\"-bdeck path\" is used one or more times to " << "specify a file or top-level directory containing ATCF " << "best track \"" << atcf_suffix << "\" data to process (required).\n" @@ -2147,6 +2230,12 @@ void usage() { << "TCPairsConfig file containing the desired configuration " << "settings (required).\n" + << "\t\t\"-diag source path\" is used one or more times to " + << "specify a file or top-level directory containing tropical " + << "cyclone diagnostics \"" << atcf_suffix + << "\" data to process. The supported formats are TCDIAG, " + << "LSDIAG_RT, LSDIAG_DEV (optional).\n" + << "\t\t\"-out base\" overrides the default output file base " << "(" << out_base << ") (optional).\n" @@ -2158,33 +2247,37 @@ void usage() { << "\tNote: The \"-adeck\", \"-edeck\", and \"-bdeck\" options " << "may include \"suffix=string\" to modify the model names " - << "from that source.\n\n"; + << "from that path.\n\n" + + << "\tNote: The \"-diag\" option may include \"model=string\" to " + << "override the model name of the tracks to which those diagnostics " + << "correspond.\n\n"; exit(1); } //////////////////////////////////////////////////////////////////////// -void set_adeck(const StringArray & a) { - set_atcf_source(a, adeck_source, adeck_model_suffix); +void set_adeck(const StringArray &a) { + set_atcf_path(a, adeck_path, adeck_model_suffix); } //////////////////////////////////////////////////////////////////////// -void set_edeck(const StringArray & a) { - set_atcf_source(a, edeck_source, edeck_model_suffix); +void set_edeck(const StringArray &a) { + set_atcf_path(a, edeck_path, edeck_model_suffix); } //////////////////////////////////////////////////////////////////////// -void set_bdeck(const StringArray & a) { - set_atcf_source(a, bdeck_source, bdeck_model_suffix); +void set_bdeck(const StringArray &a) { + set_atcf_path(a, bdeck_path, bdeck_model_suffix); } //////////////////////////////////////////////////////////////////////// -void set_atcf_source(const StringArray & a, - StringArray &source, StringArray &model_suffix) { +void set_atcf_path(const StringArray &a, + StringArray &path, StringArray &model_suffix) { int i; StringArray sa; ConcatString cs, suffix; @@ -2195,8 +2288,7 @@ void set_atcf_source(const StringArray & a, if(cs.startswith("suffix")) { sa = cs.split("="); if(sa.n() != 2) { - mlog << Error - << "\nset_atcf_source() -> " + mlog << Error << "\nset_atcf_path() -> " << "the model suffix must be specified as " << "\"suffix=string\".\n\n"; usage(); @@ -2207,11 +2299,11 @@ void set_atcf_source(const StringArray & a, } } - // Parse the remaining sources + // Parse the remaining paths for(i=0; i " + << "the \"-diag source path\" command line option " + << "must have length >= 2.\n\n"; + exit(1); + } + + // Parse the DiagType source + source = string_to_diagtype(a[0].c_str()); + + // Check for the model sub-argument + for(i=1; i " + << "the model name must be specified as " + << "\"model=string\".\n\n"; + usage(); + } + else { + model_name = sa[1]; + } + } + } + + // Parse the remaining paths + for(i=1; i diag_source; +static StringArray diag_path, diag_model_name; static ConcatString config_file; static TCPairsConfInfo conf_info; diff --git a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc index a378bbbd1d..fa6c9b318c 100644 --- a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc +++ b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc @@ -23,6 +23,11 @@ using namespace std; #include "apply_mask.h" #include "vx_log.h" +//////////////////////////////////////////////////////////////////////// + +void parse_conf_diag_convert_map(Dictionary *, + map > &); + //////////////////////////////////////////////////////////////////////// // // Code for class TCPairsConfInfo @@ -96,6 +101,8 @@ void TCPairsConfInfo::clear() { DLandFile.clear(); WatchWarnFile.clear(); WatchWarnOffset = bad_data_int; + DiagName.clear(); + DiagConvertMap.clear(); BasinMap.clear(); Version.clear(); @@ -305,6 +312,12 @@ void TCPairsConfInfo::process_config() { // Conf: WatchWarnOffset WatchWarnOffset = dict->lookup_int(conf_key_time_offset); + // Conf: DiagName + DiagName = dict->lookup_string_array(conf_key_diag_name); + + // Conf: DiagConvertMap + parse_conf_diag_convert_map(dict, DiagConvertMap); + // Conf: BasinMap BasinMap = parse_conf_key_value_map(dict, conf_key_basin_map); @@ -312,3 +325,73 @@ void TCPairsConfInfo::process_config() { } //////////////////////////////////////////////////////////////////////// +// +// Utility functions +// +//////////////////////////////////////////////////////////////////////// + +void parse_conf_diag_convert_map(Dictionary *dict, + map > &source_map) { + int i, j; + Dictionary *map_dict = (Dictionary *) 0; + map cur_map; + DiagType source; + StringArray sa; + ConcatString key; + UserFunc_1Arg fx; + + const char *method_name = "parse_conf_diag_convert_map() -> "; + + if(!dict) { + mlog << Error << "\n" << method_name + << "empty dictionary!\n\n"; + exit(1); + } + + // Conf: diag_convert_map + map_dict = dict->lookup_array(conf_key_diag_convert_map); + + // Loop through the array entries + for(i=0; in_entries(); i++) { + + // Initialize the current map + cur_map.clear(); + + // Lookup the source, key, and convert function + source = string_to_diagtype( + (*map_dict)[i]->dict_value()->lookup_string(conf_key_source).c_str()); + sa = (*map_dict)[i]->dict_value()->lookup_string_array(conf_key_key); + fx.clear(); + fx.set((*map_dict)[i]->dict_value()->lookup(conf_key_convert)); + + // Check the function + if(!fx.is_set()) { + mlog << Error << "\n" << method_name + << "lookup for \"" << conf_key_convert << "\" failed in the \"" + << conf_key_diag_convert_map << "\" map!\n\n"; + exit(1); + } + + // Add entry to the current map for each string + for(j=0; j(sa[j],fx)); + } + + // Append to the existing source entry + if(source_map.count(source) > 0) { + for(map::iterator it = cur_map.begin(); + it != cur_map.end(); it++) { + source_map.at(source).insert(pair(it->first, it->second)); + } + } + // Add a new source entry + else { + source_map.insert(pair >(source, cur_map)); + } + + } // end for i + + return; +} + +//////////////////////////////////////////////////////////////////////// diff --git a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h index 3fc82f2671..12dc7eea14 100644 --- a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h +++ b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h @@ -111,6 +111,12 @@ class TCPairsConfInfo { // Watch/warnings time offset int WatchWarnOffset; + // Diagnostics to be extracted + StringArray DiagName; + + // Diagnostic conversions + std::map< DiagType, std::map > DiagConvertMap; + // Basin Map std::map BasinMap; diff --git a/src/tools/tc_utils/tc_stat/Makefile.am b/src/tools/tc_utils/tc_stat/Makefile.am index 60827d3d98..248a1d5b99 100644 --- a/src/tools/tc_utils/tc_stat/Makefile.am +++ b/src/tools/tc_utils/tc_stat/Makefile.am @@ -31,7 +31,7 @@ tc_stat_LDADD = -lvx_stat_out \ -lvx_data2d_nccf \ -lvx_statistics \ -lvx_data2d \ - -lvx_seeps \ + -lvx_seeps \ -lvx_nc_util \ -lvx_regrid \ -lvx_grid \ diff --git a/src/tools/tc_utils/tc_stat/Makefile.in b/src/tools/tc_utils/tc_stat/Makefile.in index df8bbe5b6a..620bbce060 100644 --- a/src/tools/tc_utils/tc_stat/Makefile.in +++ b/src/tools/tc_utils/tc_stat/Makefile.in @@ -340,10 +340,10 @@ tc_stat_LDADD = -lvx_stat_out \ -lvx_data2d_nccf \ -lvx_statistics \ -lvx_data2d \ + -lvx_seeps \ -lvx_nc_util \ -lvx_regrid \ -lvx_grid \ - -lvx_seeps \ -lvx_config \ -lvx_gsl_prob \ -lvx_pb_util \ diff --git a/src/tools/tc_utils/tc_stat/out.tcst b/src/tools/tc_utils/tc_stat/out.tcst new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/tools/tc_utils/tc_stat/tc_stat.cc b/src/tools/tc_utils/tc_stat/tc_stat.cc index c536227149..10c97286e5 100644 --- a/src/tools/tc_utils/tc_stat/tc_stat.cc +++ b/src/tools/tc_utils/tc_stat/tc_stat.cc @@ -21,6 +21,7 @@ // the gen_vx_mask tool. // 004 07/06/22 Howard Soh METplus-Internal #19 Rename main to met_main // 005 09/28/22 Prestopnik MET #2227 Remove namespace std from header files +// 006 10/06/22 Halley Gotway MET #392 Incorporate diagnostics // //////////////////////////////////////////////////////////////////////// @@ -167,7 +168,7 @@ void process_jobs() { TCStatJob *cur_job = (TCStatJob *) 0; ConcatString jobstring; int i, n_jobs; - TCLineCounts n; + TCPointCounts n; const char *method_name = "process_jobs() -> "; // Open the output file @@ -210,23 +211,24 @@ void process_jobs() { << cur_job->serialize() << "\n"; // Initialize counts - memset(&n, 0, sizeof(TCLineCounts)); + memset(&n, 0, sizeof(TCPointCounts)); // Do the job cur_job->do_job(tcst_files, n); mlog << Debug(2) << method_name << "Job " << i+1 << " used " << n.NKeep << " out of " - << n.NRead << " lines read.\n"; + << n.NRead << " track points read.\n"; } else mlog << Debug(1) << method_name << "job is missing\n"; mlog << Debug(3) - << "Total lines read = " << n.NRead << "\n" - << "Total lines kept = " << n.NKeep << "\n" + << "Total track points read = " << n.NRead << "\n" + << "Total track points kept = " << n.NKeep << "\n" << "Rejected for track watch/warn = " << n.RejTrackWatchWarn << "\n" << "Rejected for init threshold = " << n.RejInitThresh << "\n" << "Rejected for init string = " << n.RejInitStr << "\n" + << "Rejected for init diag threshold = " << n.RejInitDiagThresh << "\n" << "Rejected for out init mask = " << n.RejOutInitMask << "\n" << "Rejected for water only = " << n.RejWaterOnly << "\n" << "Rejected for rapid inten = " << n.RejRIRW << "\n" @@ -249,6 +251,7 @@ void process_jobs() { << "Rejected for line type = " << n.RejLineType << "\n" << "Rejected for numeric threshold = " << n.RejColumnThresh << "\n" << "Rejected for string matching = " << n.RejColumnStr << "\n" + << "Rejected for diag threshold = " << n.RejDiagThresh << "\n" << "Rejected for match points = " << n.RejMatchPoints << "\n" << "Rejected for event equal = " << n.RejEventEqual << "\n" << "Rejected for out init mask = " << n.RejOutInitMask << "\n" diff --git a/src/tools/tc_utils/tc_stat/tc_stat_conf_info.cc b/src/tools/tc_utils/tc_stat/tc_stat_conf_info.cc index d2edb908df..d211057ef6 100644 --- a/src/tools/tc_utils/tc_stat/tc_stat_conf_info.cc +++ b/src/tools/tc_utils/tc_stat/tc_stat_conf_info.cc @@ -212,6 +212,16 @@ void TCStatConfInfo::process_config() { conf_key_init_str_exc_name, conf_key_init_str_exc_val, Filter.InitStrExcMap); + // Conf: TCStatJob::DiagThreshName, TCStatJob::DiagThreshVal + parse_conf_thresh_map(Conf, + conf_key_diag_thresh_name, conf_key_diag_thresh_val, + Filter.DiagThreshMap); + + // Conf: TCStatJob::InitDiagThreshName, TCStatJob::InitDiagThreshVal + parse_conf_thresh_map(Conf, + conf_key_init_diag_thresh_name, conf_key_init_diag_thresh_val, + Filter.InitDiagThreshMap); + // Conf: TCStatJob::WaterOnly Filter.WaterOnly = Conf.lookup_bool(conf_key_water_only); diff --git a/src/tools/tc_utils/tc_stat/tc_stat_files.cc b/src/tools/tc_utils/tc_stat/tc_stat_files.cc index bd77dad0f4..5f63e0fd57 100644 --- a/src/tools/tc_utils/tc_stat/tc_stat_files.cc +++ b/src/tools/tc_utils/tc_stat/tc_stat_files.cc @@ -137,7 +137,7 @@ bool TCStatFiles::operator>>(TrackPairInfo &pair) { CurFile++; // Check for the last file - if(CurFile == FileList.n_elements()) return(false); + if(CurFile == FileList.n()) return(false); else { // Open the next file for reading @@ -152,7 +152,7 @@ bool TCStatFiles::operator>>(TrackPairInfo &pair) { // List file being read mlog << Debug(3) << "Reading file " << CurFile+1 << " of " - << FileList.n_elements() << ": " << FileList[CurFile] + << FileList.n() << ": " << FileList[CurFile] << "\n"; } // end else @@ -161,16 +161,28 @@ bool TCStatFiles::operator>>(TrackPairInfo &pair) { // Read lines to the end of the track or file while(CurLDF >> line) { - // Skip header and non-TCMPR lines - if(line.is_header() || line.type() != TCStatLineType_TCMPR) continue; + // Skip header and non-TCMPR/TCDIAG lines + if(line.is_header() || + (line.type() != TCStatLineType_TCMPR && + line.type() != TCStatLineType_TCDIAG)) continue; // Add the current point pair.add(line); // Break out of the loop at the end of the track if(atoi(line.get_item("TOTAL")) == - atoi(line.get_item("INDEX"))) break; + atoi(line.get_item("INDEX"))) { + + // Check for a trailing TCDIAG line + if(CurLDF.peek_line(line)) { + if(line.type() == TCStatLineType_TCDIAG) { + pair.add(line); + CurLDF >> line; + } + } + break; + } } // end while return(true); @@ -192,7 +204,7 @@ bool TCStatFiles::operator>>(ProbRIRWPairInfo &pair) { CurFile++; // Check for the last file - if(CurFile == FileList.n_elements()) return(false); + if(CurFile == FileList.n()) return(false); else { // Open the next file for reading @@ -207,7 +219,7 @@ bool TCStatFiles::operator>>(ProbRIRWPairInfo &pair) { // List file being read mlog << Debug(3) << "Reading file " << CurFile+1 << " of " - << FileList.n_elements() << ": " << FileList[CurFile] + << FileList.n() << ": " << FileList[CurFile] << "\n"; } // end else @@ -224,7 +236,7 @@ bool TCStatFiles::operator>>(ProbRIRWPairInfo &pair) { break; - } //end while + } // end while return(status); } @@ -241,7 +253,7 @@ bool TCStatFiles::operator>>(TCStatLine &line) { CurFile++; // Check for the last file - if(CurFile == FileList.n_elements()) return(false); + if(CurFile == FileList.n()) return(false); else { // Open the next file for reading @@ -256,7 +268,7 @@ bool TCStatFiles::operator>>(TCStatLine &line) { // List file being read mlog << Debug(3) << "Reading file " << CurFile+1 << " of " - << FileList.n_elements() << ": " << FileList[CurFile] + << FileList.n() << ": " << FileList[CurFile] << "\n"; } // end else @@ -270,7 +282,7 @@ bool TCStatFiles::operator>>(TCStatLine &line) { break; - } //end while + } // end while return(status); } diff --git a/src/tools/tc_utils/tc_stat/tc_stat_job.cc b/src/tools/tc_utils/tc_stat/tc_stat_job.cc index 79e217ffa4..79f85cb5eb 100644 --- a/src/tools/tc_utils/tc_stat/tc_stat_job.cc +++ b/src/tools/tc_utils/tc_stat/tc_stat_job.cc @@ -224,6 +224,9 @@ void TCStatJob::clear() { InitThreshMap.clear(); InitStrIncMap.clear(); InitStrExcMap.clear(); + DiagThreshMap.clear(); + InitDiagThreshMap.clear(); + PrintDiagWarning.clear(); EventEqualLead.clear(); EventEqualCases.clear(); @@ -307,6 +310,9 @@ void TCStatJob::assign(const TCStatJob & j) { InitThreshMap = j.InitThreshMap; InitStrIncMap = j.InitStrIncMap; InitStrExcMap = j.InitStrExcMap; + DiagThreshMap = j.DiagThreshMap; + InitDiagThreshMap = j.InitDiagThreshMap; + PrintDiagWarning = j.PrintDiagWarning; DumpFile = j.DumpFile; open_dump_file(); @@ -456,6 +462,18 @@ void TCStatJob::dump(ostream & out, int depth) const { str_it->second.dump(out, depth + 1); } + out << prefix << "DiagThreshMap ...\n"; + for(thr_it=DiagThreshMap.begin(); thr_it!= DiagThreshMap.end(); thr_it++) { + out << prefix << thr_it->first << ": \n"; + thr_it->second.dump(out, depth + 1); + } + + out << prefix << "InitDiagThreshMap ...\n"; + for(thr_it=InitDiagThreshMap.begin(); thr_it!= InitDiagThreshMap.end(); thr_it++) { + out << prefix << thr_it->first << ": \n"; + thr_it->second.dump(out, depth + 1); + } + out << prefix << "WaterOnly = " << bool_to_string(WaterOnly) << "\n"; out << prefix << "RIRWTrack = " << tracktype_to_string(RIRWTrack) << "\n"; @@ -507,7 +525,7 @@ void TCStatJob::dump(ostream & out, int depth) const { //////////////////////////////////////////////////////////////////////// bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, - TCLineCounts &n) const { + TCPointCounts &n) { bool keep = true; int i, i_init; double v_dbl; @@ -546,7 +564,7 @@ bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, if(!keep) n.RejTrackWatchWarn += pair.n_points(); // Get the index of the track initialization point - i_init=pair.i_init(); + i_init = pair.i_init(); // Check for bad track initialization point if(is_bad_data(i_init)) { @@ -562,6 +580,10 @@ bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, keep = false; n.RejInitStr += pair.n_points(); } + else if(InitDiagThreshMap.size() > 0) { + keep = false; + n.RejInitDiagThresh += pair.n_points(); + } else if(OutInitMaskName.nonempty()) { keep = false; n.RejOutInitMask += pair.n_points(); @@ -575,7 +597,7 @@ bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, for(thr_it=InitThreshMap.begin(); thr_it!= InitThreshMap.end(); thr_it++) { // Get the numeric init column value - v_dbl = get_column_double(*pair.line(i_init), thr_it->first); + v_dbl = get_column_double(*pair.tcmpr_line(i_init), thr_it->first); // Check the column threshold if(!thr_it->second.check_dbl(v_dbl)) { @@ -592,7 +614,7 @@ bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, for(str_it=InitStrIncMap.begin(); str_it!= InitStrIncMap.end(); str_it++) { // Retrieve the column value - v_str = pair.line(i_init)->get_item(str_it->first.c_str()); + v_str = pair.tcmpr_line(i_init)->get_item(str_it->first.c_str()); // Check the string value if(!str_it->second.has(v_str)) { @@ -609,7 +631,7 @@ bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, for(str_it=InitStrExcMap.begin(); str_it!= InitStrExcMap.end(); str_it++) { // Retrieve the column value - v_str = pair.line(i_init)->get_item(str_it->first.c_str()); + v_str = pair.tcmpr_line(i_init)->get_item(str_it->first.c_str()); // Check the string value if(str_it->second.has(v_str)) { @@ -620,6 +642,27 @@ bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, } } + // Check InitDiagThreshMap + if(keep == true) { + + int i_diag; + + // Loop through the numeric diagnostic thresholds + for(thr_it=InitDiagThreshMap.begin(); thr_it!= InitDiagThreshMap.end(); thr_it++) { + + // Get the numeric diagnostic value + v_dbl = get_diag_double(pair.adeck().diag_name(), + pair.adeck()[i_init], thr_it->first); + + // Check the diagnostic value threshold + if(!thr_it->second.check_dbl(v_dbl)) { + keep = false; + n.RejInitDiagThresh += pair.n_points(); + break; + } + } + } + // Check OutInitMask if(keep == true) { @@ -665,7 +708,11 @@ bool TCStatJob::is_keeper_track(const TrackPairInfo &pair, //////////////////////////////////////////////////////////////////////// bool TCStatJob::is_keeper_line(const TCStatLine &line, - TCLineCounts &n) const { + TCPointCounts &n) const { + + // Does not apply to TCDIAG lines + if(line.type() == TCStatLineType_TCDIAG) return(true); + bool keep = true; double v_dbl, alat, alon, blat, blon; ConcatString v_str; @@ -881,6 +928,62 @@ double TCStatJob::get_column_double(const TCStatLine &line, //////////////////////////////////////////////////////////////////////// +bool TCStatJob::is_keeper_tcdiag(const StringArray &diag_name, + const TrackPoint &point, + TCPointCounts &n) { + bool keep = true; + double v_dbl; + map::const_iterator thr_it; + + // Check DiagThreshMap + + // Loop through the numeric diagnostic thresholds + for(thr_it=DiagThreshMap.begin(); thr_it!= DiagThreshMap.end(); thr_it++) { + + // Get the numeric diagnostic value + v_dbl = get_diag_double(diag_name, point, thr_it->first); + + // Check the diagnostic value threshold + if(!thr_it->second.check_dbl(v_dbl)) { + keep = false; + n.RejDiagThresh++; + break; + } + } + + // Update counts + if(!keep) n.NKeep -= 1; + + return(keep); +} + +//////////////////////////////////////////////////////////////////////// + +double TCStatJob::get_diag_double(const StringArray &diag_name, + const TrackPoint &point, + const ConcatString &diag_cs) { + double v = bad_data_double; + int i_diag; + + // Check whether the diagnostic name exists + if(diag_name.has(diag_cs, i_diag)) { + v = point.diag_val(i_diag); + } + // Print a single warning message if diagnostics are present + // but not the requested one + else if(diag_name.n() > 0 && + !PrintDiagWarning.has(diag_cs)) { + mlog << Warning << "\nTCStatJob::get_diag_double() -> " + << "diagnostic \"" << diag_cs + << "\" not defined for TrackPoint.\n\n"; + PrintDiagWarning.add(diag_cs); + } + + return(v); +} + +//////////////////////////////////////////////////////////////////////// + StringArray TCStatJob::parse_job_command(const char *jobstring) { StringArray a, b; string c; @@ -940,6 +1043,10 @@ StringArray TCStatJob::parse_job_command(const char *jobstring) { a.shift_down(i, 2); } else if(c.compare("-init_str_exc" ) == 0) { parse_string_option(a[i+1].c_str(), a[i+2].c_str(), InitStrExcMap); a.shift_down(i, 2); } + else if(c.compare("-diag_thresh" ) == 0) { parse_thresh_option(a[i+1].c_str(), a[i+2].c_str(), DiagThreshMap); + a.shift_down(i, 2); } + else if(c.compare("-init_diag_thresh" ) == 0) { parse_thresh_option(a[i+1].c_str(), a[i+2].c_str(), InitDiagThreshMap); + a.shift_down(i, 2); } else if(c.compare("-water_only" ) == 0) { WaterOnly = string_to_bool(a[i+1].c_str()); a.shift_down(i, 1); } else if(c.compare("-rirw_track" ) == 0) { RIRWTrack = string_to_tracktype(a[i+1].c_str()); a.shift_down(i, 1); } else if(c.compare("-rirw_time" ) == 0) { RIRWTimeADeck = timestring_to_sec(a[i+1].c_str()); @@ -1102,15 +1209,25 @@ void TCStatJob::dump_pair(const TrackPairInfo &pair, ofstream *out) { TcHdrColumns tchc; AsciiTable out_at; - int i_row, hdr_row; + int i_row, hdr_row, n_row, n_col, n_diag; + + // Number of output rows and columns + n_row = pair.n_points(); + n_col = n_tc_mpr_cols; // Determine if we need to write a header row if((int) out->tellp() == 0) hdr_row = 1; else hdr_row = 0; + // Update rows and columns for diagnosics output + if((n_diag = pair.adeck().n_diag()) > 0) { + n_row += pair.n_points(); + n_col = max(n_col, get_n_tc_diag_cols(n_diag)); + } + // Initialize the output AsciiTable - out_at.set_size(pair.n_points() + hdr_row, - n_tc_header_cols + n_tc_mpr_cols); + out_at.set_size(n_row + hdr_row, + n_tc_header_cols + n_col); // Write the TCMPR header row write_tc_mpr_header_row(1, out_at, 0, 0); @@ -1138,7 +1255,7 @@ void TCStatJob::dump_pair(const TrackPairInfo &pair, ofstream *out) { // Write the TrackPairInfo object i_row = hdr_row; - write_tc_mpr_row(tchc, pair, out_at, i_row); + write_track_pair_info(tchc, pair, out_at, i_row); // Write the AsciiTable to the file *out << out_at; @@ -1253,6 +1370,18 @@ ConcatString TCStatJob::serialize() const { << str_it->second[i] << " "; } } + for(thr_it=DiagThreshMap.begin(); thr_it!= DiagThreshMap.end(); thr_it++) { + for(i=0; isecond.n(); i++) { + s << "-diag_thresh " << thr_it->first << " " + << thr_it->second[i].get_str() << " "; + } + } + for(thr_it=InitDiagThreshMap.begin(); thr_it!= InitDiagThreshMap.end(); thr_it++) { + for(i=0; isecond.n(); i++) { + s << "-init_diag_thresh " << thr_it->first << " " + << thr_it->second[i].get_str() << " "; + } + } if(WaterOnly != default_water_only) s << "-water_only " << bool_to_string(WaterOnly) << " "; if(RIRWTrack != default_rirw_track) { @@ -1309,7 +1438,7 @@ ConcatString TCStatJob::serialize() const { //////////////////////////////////////////////////////////////////////// -void TCStatJob::do_job(const StringArray &file_list, TCLineCounts &n) { +void TCStatJob::do_job(const StringArray &file_list, TCPointCounts &n) { mlog << Error << "\nTCStatJob::do_job() -> " << "the do_job() base class function should never be called!\n\n"; @@ -1326,7 +1455,7 @@ void TCStatJob::do_job(const StringArray &file_list, TCLineCounts &n) { void TCStatJob::event_equalize_tracks() { TrackPairInfo pair; - TCLineCounts n; + TCPointCounts n; ConcatString key, val; StringArray case_list; ConcatString models; @@ -1362,11 +1491,11 @@ void TCStatJob::event_equalize_tracks() { if(skip) continue; // Add event equalization information for each track point - for(i=0; iheader(); // Track line header + key = pair.adeck().technique(); // ADECK model name + val = pair.tcmpr_line(i)->header(); // Track line header // Add a new map entry, if necessary if(case_map.count(key) == 0) { @@ -1414,7 +1543,7 @@ void TCStatJob::event_equalize_tracks() { //////////////////////////////////////////////////////////////////////// void TCStatJob::event_equalize_lines() { - TCLineCounts n; + TCPointCounts n; TCStatLine line; ConcatString key, val; StringArray case_list; @@ -1484,7 +1613,7 @@ void TCStatJob::event_equalize_lines() { //////////////////////////////////////////////////////////////////////// -void TCStatJob::subset_track_pair(TrackPairInfo &pair, TCLineCounts &n) { +void TCStatJob::subset_track_pair(TrackPairInfo &pair, TCPointCounts &n) { int i, n_rej; // Check for no points @@ -1537,13 +1666,19 @@ void TCStatJob::subset_track_pair(TrackPairInfo &pair, TCLineCounts &n) { } // Loop over the track points to check line by line - for(i=0; i " << "the track-based " << TCStatLineType_TCMPR_Str - << " line type cannot be mixed with the " + << " and " << TCStatLineType_TCDIAG_Str + << " line types cannot be mixed with the " << LineType[i] << " line type.\n\n"; exit(1); } @@ -1680,7 +1820,7 @@ void TCStatJobFilter::do_job(const StringArray &file_list, //////////////////////////////////////////////////////////////////////// -void TCStatJobFilter::filter_tracks(TCLineCounts &n) { +void TCStatJobFilter::filter_tracks(TCPointCounts &n) { TrackPairInfo pair; // Apply the event equalization logic to build a list of common cases @@ -1719,7 +1859,7 @@ void TCStatJobFilter::filter_tracks(TCLineCounts &n) { //////////////////////////////////////////////////////////////////////// -void TCStatJobFilter::filter_lines(TCLineCounts &n) { +void TCStatJobFilter::filter_lines(TCPointCounts &n) { TCStatLine line; // Apply the event equalization logic to build a list of common cases @@ -1962,7 +2102,7 @@ ConcatString TCStatJobSummary::serialize() const { //////////////////////////////////////////////////////////////////////// void TCStatJobSummary::do_job(const StringArray &file_list, - TCLineCounts &n) { + TCPointCounts &n) { TrackPairInfo pair; int i; @@ -1979,8 +2119,11 @@ void TCStatJobSummary::do_job(const StringArray &file_list, // Apply logic based on the LineType // - // If not specified, assume TCMPR by adding it to the LineType - if(LineType.n() == 0) LineType.add(TCStatLineType_TCMPR_Str); + // If not specified, assume TCMPR/TCDIAG by adding them to the LineType + if(LineType.n() == 0) { + LineType.add(TCStatLineType_TCMPR_Str); + LineType.add(TCStatLineType_TCDIAG_Str); + } // Add the input file list TCSTFiles.add_files(file_list); @@ -1988,12 +2131,14 @@ void TCStatJobSummary::do_job(const StringArray &file_list, // Process track-based filter job if(LineType.has(TCStatLineType_TCMPR_Str)) { - // TCMPR and non-TCMPR LineTypes cannot be mixed + // TCMPR and non-TCMPR/TCDIAG LineTypes cannot be mixed for(i=0; i " << "the track-based " << TCStatLineType_TCMPR_Str - << " line type cannot be mixed with the " + << " and " << TCStatLineType_TCDIAG_Str + << " line types cannot be mixed with the " << LineType[i] << " line type.\n\n"; exit(1); } @@ -2019,7 +2164,7 @@ void TCStatJobSummary::do_job(const StringArray &file_list, //////////////////////////////////////////////////////////////////////// -void TCStatJobSummary::summarize_tracks(TCLineCounts &n) { +void TCStatJobSummary::summarize_tracks(TCPointCounts &n) { TrackPairInfo pair; // Apply the event equalization logic to build a list of common cases @@ -2061,7 +2206,7 @@ void TCStatJobSummary::summarize_tracks(TCLineCounts &n) { //////////////////////////////////////////////////////////////////////// -void TCStatJobSummary::summarize_lines(TCLineCounts &n) { +void TCStatJobSummary::summarize_lines(TCPointCounts &n) { TCStatLine line; // Apply the event equalization logic to build a list of common cases @@ -2111,8 +2256,8 @@ void TCStatJobSummary::process_pair(TrackPairInfo &pair) { // Initialize the map cur_map.clear(); - // Loop over TCStatLines and construct a summary map - for(i=0; iheader()); - cur_map[key].AModel.add(pair.line(i)->amodel()); - cur_map[key].Init.add(pair.line(i)->init()); - cur_map[key].Lead.add(pair.line(i)->lead()); - cur_map[key].Valid.add(pair.line(i)->valid()); + cur_map[key].Hdr.add(pair.tcmpr_line(i)->header()); + cur_map[key].AModel.add(pair.tcmpr_line(i)->amodel()); + cur_map[key].Init.add(pair.tcmpr_line(i)->init()); + cur_map[key].Lead.add(pair.tcmpr_line(i)->lead()); + cur_map[key].Valid.add(pair.tcmpr_line(i)->valid()); } // end for j } // end for i @@ -2967,7 +3122,7 @@ ConcatString TCStatJobRIRW::serialize() const { //////////////////////////////////////////////////////////////////////// void TCStatJobRIRW::do_job(const StringArray &file_list, - TCLineCounts &n) { + TCPointCounts &n) { TrackPairInfo pair; // Add the input file list @@ -3039,7 +3194,7 @@ void TCStatJobRIRW::process_pair(TrackPairInfo &pair) { for(i=0; iamodel() << sep - << pair.line(i)->bmodel() << sep - << pair.line(i)->storm_id() << sep - << unix_to_yyyymmdd_hhmmss(pair.line(i)->init()) << sep + << pair.tcmpr_line(i)->amodel() << sep + << pair.tcmpr_line(i)->bmodel() << sep + << pair.tcmpr_line(i)->storm_id() << sep + << unix_to_yyyymmdd_hhmmss(pair.tcmpr_line(i)->init()) << sep << (is_bad_data(lead) ? na_str : sec_to_hhmmss(lead).text()) << sep - << unix_to_yyyymmdd_hhmmss(pair.line(i)->valid()) << sep + << unix_to_yyyymmdd_hhmmss(pair.tcmpr_line(i)->valid()) << sep << aprv << sep << acur << sep << adlt << sep << (is_bad_data(a) ? na_str : bool_to_string(a)) << sep << bprv << sep << bcur << sep << bdlt << sep @@ -3619,7 +3774,7 @@ StringArray TCStatJobProbRIRW::parse_job_command(const char *jobstring) { //////////////////////////////////////////////////////////////////////// void TCStatJobProbRIRW::close_dump_file() { - const char *method_name = "TCStatJobProbRIRW::do_job() -> "; + const char *method_name = "TCStatJobProbRIRW::close_dump_file() -> "; // Close the current output dump file stream if(DumpOut) { @@ -3716,7 +3871,7 @@ ConcatString TCStatJobProbRIRW::serialize() const { //////////////////////////////////////////////////////////////////////// void TCStatJobProbRIRW::do_job(const StringArray &file_list, - TCLineCounts &n) { + TCPointCounts &n) { ProbRIRWPairInfo pair; // Check that -probrirw_thresh has been supplied @@ -3967,7 +4122,8 @@ void TCStatJobProbRIRW::do_output(ostream &out) { //////////////////////////////////////////////////////////////////////// -TCLineCounts::TCLineCounts() { +TCPointCounts::TCPointCounts() { + // Read and keep counts NRead = 0; NKeep = 0; @@ -3976,6 +4132,7 @@ TCLineCounts::TCLineCounts() { RejTrackWatchWarn = 0; RejInitThresh = 0; RejInitStr = 0; + RejInitDiagThresh = 0; // Filtering on track attributes RejRIRW = 0; @@ -4000,6 +4157,7 @@ TCLineCounts::TCLineCounts() { RejWaterOnly = 0; RejColumnThresh = 0; RejColumnStr = 0; + RejDiagThresh = 0; RejMatchPoints = 0; RejEventEqual = 0; RejOutInitMask = 0; diff --git a/src/tools/tc_utils/tc_stat/tc_stat_job.h b/src/tools/tc_utils/tc_stat/tc_stat_job.h index c01bcaeb2d..020134b593 100644 --- a/src/tools/tc_utils/tc_stat/tc_stat_job.h +++ b/src/tools/tc_utils/tc_stat/tc_stat_job.h @@ -118,8 +118,8 @@ extern ConcatString tcstatjobtype_to_string(const TCStatJobType); //////////////////////////////////////////////////////////////////////// -// Struct for counts of lines read and rejected -struct TCLineCounts { +// Struct for counts of track points read and rejected +struct TCPointCounts { // Read and keep counts int NRead; @@ -129,6 +129,7 @@ struct TCLineCounts { int RejTrackWatchWarn; int RejInitThresh; int RejInitStr; + int RejInitDiagThresh; // Filtering on track attributes int RejRIRW; @@ -153,13 +154,14 @@ struct TCLineCounts { int RejWaterOnly; int RejColumnThresh; int RejColumnStr; + int RejDiagThresh; int RejMatchPoints; int RejEventEqual; int RejOutInitMask; int RejOutValidMask; int RejLeadReq; - TCLineCounts(); + TCPointCounts(); }; //////////////////////////////////////////////////////////////////////// @@ -205,12 +207,16 @@ class TCStatJob { ////////////////////////////////////////////////////////////////// - bool is_keeper_track(const TrackPairInfo &, TCLineCounts &) const; + bool is_keeper_track(const TrackPairInfo &, TCPointCounts &); - bool is_keeper_line(const TCStatLine &, TCLineCounts &) const; + bool is_keeper_line(const TCStatLine &, TCPointCounts &) const; double get_column_double(const TCStatLine &, const ConcatString &) const; + bool is_keeper_tcdiag(const StringArray &, const TrackPoint &, TCPointCounts &); + + double get_diag_double(const StringArray &, const TrackPoint &, const ConcatString &); + ////////////////////////////////////////////////////////////////// virtual StringArray parse_job_command(const char *); @@ -231,13 +237,13 @@ class TCStatJob { ////////////////////////////////////////////////////////////////// - virtual void do_job(const StringArray &, TCLineCounts &); + virtual void do_job(const StringArray &, TCPointCounts &); void event_equalize_tracks(); void event_equalize_lines(); - void subset_track_pair(TrackPairInfo &, TCLineCounts &); + void subset_track_pair(TrackPairInfo &, TCPointCounts &); ////////////////////////////////////////////////////////////////// @@ -298,6 +304,11 @@ class TCStatJob { std::map InitStrIncMap; std::map InitStrExcMap; + // Numeric diagnostic thresholds + std::map DiagThreshMap; + std::map InitDiagThreshMap; + StringArray PrintDiagWarning; + // Variables to the store the analysis job specification ConcatString DumpFile; // Dump TrackPairInfo used to a file std::ofstream *DumpOut; // Dump output file stream @@ -371,10 +382,10 @@ class TCStatJobFilter : public TCStatJob { void clear(); - void do_job(const StringArray &, TCLineCounts &); // virtual from base class + void do_job(const StringArray &, TCPointCounts &); // virtual from base class - void filter_tracks(TCLineCounts &); - void filter_lines (TCLineCounts &); + void filter_tracks(TCPointCounts &); + void filter_lines (TCPointCounts &); void do_output(std::ostream &); @@ -405,10 +416,10 @@ class TCStatJobSummary : public TCStatJob { ConcatString serialize() const; - void do_job(const StringArray &, TCLineCounts &); // virtual from base class + void do_job(const StringArray &, TCPointCounts &); // virtual from base class - void summarize_tracks(TCLineCounts &); - void summarize_lines (TCLineCounts &); + void summarize_tracks(TCPointCounts &); + void summarize_lines (TCPointCounts &); void process_pair(TrackPairInfo &); void process_line(TCStatLine &); @@ -469,7 +480,7 @@ class TCStatJobRIRW : public TCStatJob { ConcatString serialize() const; - void do_job(const StringArray &, TCLineCounts &); // virtual from base class + void do_job(const StringArray &, TCPointCounts &); // virtual from base class void process_pair(TrackPairInfo &); @@ -519,7 +530,7 @@ class TCStatJobProbRIRW : public TCStatJob { ConcatString serialize() const; - void do_job(const StringArray &, TCLineCounts &); // virtual from base class + void do_job(const StringArray &, TCPointCounts &); // virtual from base class void process_pair(ProbRIRWPairInfo &);