From 81513602f9b762c74e7842e8e7e1ee23a36f8815 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Wed, 9 Oct 2024 13:21:10 -0500 Subject: [PATCH] Feature #2882 seeps qa (#2987) * Update seeps.h Change variable names to reduce ambiguity for interpretation and aid useability. * Update seeps.cc Pull through variable name changes and renaming of functions to aid legibility and clarity. Introduced some additional debug print statements. * Update grid-stat.rst Add documentation about the location of the gridded climatology files for SEEPS and which environment variable to use. * Replace read_seeps_scores() with get_seeps_climo_grid() * Manually merging Rachel's patch-1 changes. * Getting close to getting these seeps changes to compile. But it's failing in pair_data_point.cc * Per #2882, get branch feature_2882_seeps_qa compiling again. Recommend revisiting the volume of SEEPS-related Debug log messages and reducing them once its fully tested. * Per #2882, need to update the handling of the PPT24_seepsweights_grid.nc file name. Rename as _v12.0.nc for the updated version with the new names so that the existing regressions tests and nightly builds for main_v11.1 and develop continue to work. We can remove the _v12.0 once this feature branch is merged into develop but for the time being, we need both versions to exist. * Per #2882, rename the SEEPS columns from S12, S13, S21, S23, S31, S32 to the more descriptive ODFL, ODFH, OLFD, OLFH, OHFD, OHFL names. * Per #2882, update SEEPS details * Per #2882, store and report the weighted mean fcst and mean obs, just like the SEEPS score itself so that they're handled in a consistent manner. Note however that it's hard-coded to NOT write the weighted means/score, only the unweighted ones. * Per #2882, change SEEPS debug log levels and correct the storage of mean_fcst and mean_obs values. * Per #2882, correct SEEPS column name lookups * Per #2882, call is_bad_data() instead of is_eq(..., -9999.0) to get rid of compiler warning message. * Per #2882, add 2 more variations of the is_eq() function with mixed float and double inputs to satisfy compiler pb2nc compiler warnings. * Per #2882, switch from dynamically allocated arrays to std::vector * Per #2882, enhance Stat-Analysis to write the SEEPS line type to an output .stat file. * Per #2882, update the aggregated seeps computation to use better-initialized vectors. * Per #2882, resolve a few more SonarQube code smells. * Per #2882, now that this PR is ready to merge, remove the v12.0 version number from the gridded SEEPS climo file name ci-skip-all --------- Co-authored-by: mpm-meto <64001904+mpm-meto@users.noreply.github.com> --- data/table_files/met_header_columns_V12.0.txt | 2 +- docs/Users_Guide/grid-stat.rst | 2 + docs/Users_Guide/point-stat.rst | 44 +- internal/test_unit/hdr/met_12_0.hdr | 2 +- internal/test_unit/xml/unit_grid_stat.xml | 2 +- src/basic/vx_util/stat_column_defs.h | 6 +- src/libcode/vx_analysis_util/stat_job.cc | 2 + src/libcode/vx_seeps/seeps.cc | 443 ++++++++++-------- src/libcode/vx_seeps/seeps.h | 109 ++--- src/libcode/vx_stat_out/stat_columns.cc | 71 +-- src/libcode/vx_statistics/compute_stats.cc | 335 ++++++++----- src/libcode/vx_statistics/pair_data_point.cc | 19 +- .../core/stat_analysis/parse_stat_line.cc | 23 +- 13 files changed, 622 insertions(+), 438 deletions(-) diff --git a/data/table_files/met_header_columns_V12.0.txt b/data/table_files/met_header_columns_V12.0.txt index 5a7f27978e..3f4785edf5 100644 --- a/data/table_files/met_header_columns_V12.0.txt +++ b/data/table_files/met_header_columns_V12.0.txt @@ -6,7 +6,7 @@ V12.0 : STAT : ISC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID V12.0 : STAT : MCTC : 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_CAT) F[0-9]*_O[0-9]* EC_VALUE V12.0 : STAT : MCTS : 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_CAT ACC ACC_NCL ACC_NCU ACC_BCL ACC_BCU HK HK_BCL HK_BCU HSS HSS_BCL HSS_BCU GER GER_BCL GER_BCU HSS_EC HSS_EC_BCL HSS_EC_BCU EC_VALUE V12.0 : STAT : MPR : 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 INDEX OBS_SID OBS_LAT OBS_LON OBS_LVL OBS_ELV FCST OBS OBS_QC OBS_CLIMO_MEAN OBS_CLIMO_STDEV OBS_CLIMO_CDF FCST_CLIMO_MEAN FCST_CLIMO_STDEV -V12.0 : STAT : SEEPS : 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 S12 S13 S21 S23 S31 S32 PF1 PF2 PF3 PV1 PV2 PV3 MEAN_FCST MEAN_OBS SEEPS +V12.0 : STAT : SEEPS : 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 ODFL ODFH OLFD OLFH OHFD OHFL PF1 PF2 PF3 PV1 PV2 PV3 MEAN_FCST MEAN_OBS SEEPS V12.0 : STAT : SEEPS_MPR : 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 OBS_SID OBS_LAT OBS_LON FCST OBS OBS_QC FCST_CAT OBS_CAT P1 P2 T1 T2 SEEPS V12.0 : STAT : NBRCNT : 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 FBS FBS_BCL FBS_BCU FSS FSS_BCL FSS_BCU AFSS AFSS_BCL AFSS_BCU UFSS UFSS_BCL UFSS_BCU F_RATE F_RATE_BCL F_RATE_BCU O_RATE O_RATE_BCL O_RATE_BCU V12.0 : STAT : NBRCTC : 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 FY_OY FY_ON FN_OY FN_ON diff --git a/docs/Users_Guide/grid-stat.rst b/docs/Users_Guide/grid-stat.rst index cac539e497..b10b1b3431 100644 --- a/docs/Users_Guide/grid-stat.rst +++ b/docs/Users_Guide/grid-stat.rst @@ -80,6 +80,8 @@ The Stable Equitable Error in Probability Space (SEEPS) was devised for monitori The capability to calculate the SEEPS has also been added to Grid-Stat. This follows the method described in :ref:`North et al, 2022 `, which uses the TRMM 3B42 v7 gridded satellite product for the climatological values and interpolates the forecast and observed products onto this grid for evaluation. A 24-hour TRMM climatology (valid at 00 UTC) constructed from data over the time period 1998-2015 is supplied with the release. Expansion of the capability to other fields will occur as well vetted examples and funding allow. +The gridded climatology required to compute SEEPS is not distributed as part of the code release and can be downloaded from `Zenodo `. The path to the file needs to be specified using MET_SEEPS_GRID_CLIMO_NAME. + Fourier Decomposition --------------------- diff --git a/docs/Users_Guide/point-stat.rst b/docs/Users_Guide/point-stat.rst index d6de2d32b1..73b7ce9a33 100644 --- a/docs/Users_Guide/point-stat.rst +++ b/docs/Users_Guide/point-stat.rst @@ -1549,10 +1549,10 @@ The first set of header columns are common to all of the output files generated - Quality control flag for observation * - 31 - FCST_CAT - - Forecast category to 3 by 3 matrix + - Forecast category (dry, light, or heavy) * - 32 - OBS_CAT - - Observationtegory to 3 by 3 matrix + - Observation category (dry, light, or heavy) * - 33 - P1 - Climo-derived probability value for this station (dry) @@ -1561,10 +1561,10 @@ The first set of header columns are common to all of the output files generated - Climo-derived probability value for this station (dry + light) * - 35 - T1 - - Threshold 1 for p1 + - Threshold 1 for P1 (dry) * - 36 - T2 - - Threshold 2 for p2 + - Threshold 2 for P2 (dry + light) * - 37 - SEEPS - SEEPS (Stable Equitable Error in Probability Space) score @@ -1589,41 +1589,41 @@ The first set of header columns are common to all of the output files generated - TOTAL - Total number of SEEPS matched pairs * - 26 - - S12 - - Counts multiplied by the weights for FCST_CAT 1 and OBS_CAT 2 + - ODFL + - Counts multiplied by the weights for the observation dry, forecast light category * - 27 - - S13 - - Counts multiplied by the weights for FCST_CAT 1 and OBS_CAT 3 + - ODFH + - Counts multiplied by the weights for the observation dry, forecast heavy category * - 28 - - S21 - - Counts multiplied by the weights for FCST_CAT 2 and OBS_CAT 1 + - OLFD + - Counts multiplied by the weights for the observation light, forecast dry category * - 29 - - S23 - - Counts multiplied by the weights for FCST_CAT 2 and OBS_CAT 3 + - OLFH + - Counts multiplied by the weights for the observation light, forecast heavy category * - 30 - - S31 - - Counts multiplied by the weights for FCST_CAT 3 and OBS_CAT 1 + - OHFD + - Counts multiplied by the weights for the observation heavy, forecast dry category * - 31 - - S32 - - Counts multiplied by the weights for FCST_CAT 3 and OBS_CAT 2 + - OHFL + - Counts multiplied by the weights for the observation heavy, forecast light category * - 32 - PF1 - - marginal probabilities of the forecast values (FCST_CAT 1) + - Marginal probabilities of the forecast dry (FCST_CAT 0) * - 33 - PF2 - - marginal probabilities of the forecast values (FCST_CAT 2) + - Marginal probabilities of the forecast light (FCST_CAT 1) * - 34 - PF3 - - marginal probabilities of the forecast values (FCST_CAT 3) + - Marginal probabilities of the forecast heavy (FCST_CAT 2) * - 35 - PV1 - - marginal probabilities of the observed values (OBS_CAT 1) + - Marginal probabilities of the observed dry (OBS_CAT 0) * - 36 - PV2 - - marginal probabilities of the observed values (OBS_CAT 2) + - Marginal probabilities of the observed light (OBS_CAT 1) * - 37 - PV3 - - marginal probabilities of the observed values (OBS_CAT 3) + - Marginal probabilities of the observed heavy (OBS_CAT 2) * - 38 - SEEPS - Averaged SEEPS (Stable Equitable Error in Probability Space) score diff --git a/internal/test_unit/hdr/met_12_0.hdr b/internal/test_unit/hdr/met_12_0.hdr index f3785b7b48..f8655a4a47 100644 --- a/internal/test_unit/hdr/met_12_0.hdr +++ b/internal/test_unit/hdr/met_12_0.hdr @@ -6,7 +6,7 @@ ISC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_L MCTC : 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_CAT _VAR_ MCTS : 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_CAT ACC ACC_NCL ACC_NCU ACC_BCL ACC_BCU HK HK_BCL HK_BCU HSS HSS_BCL HSS_BCU GER GER_BCL GER_BCU HSS_EC HSS_EC_BCL HSS_EC_BCU EC_VALUE MPR : 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 INDEX OBS_SID OBS_LAT OBS_LON OBS_LVL OBS_ELV FCST OBS OBS_QC OBS_CLIMO_MEAN OBS_CLIMO_STDEV OBS_CLIMO_CDF FCST_CLIMO_MEAN FCST_CLIMO_STDEV -SEEPS : 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 S12 S13 S21 S23 S31 S32 PF1 PF2 PF3 PV1 PV2 PV3 MEAN_FCST MEAN_OBS SEEPS +SEEPS : 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 ODFL ODFH OLFD OLFH OHFD OHFL PF1 PF2 PF3 PV1 PV2 PV3 MEAN_FCST MEAN_OBS SEEPS SEEPS_MPR : 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 OBS_SID OBS_LAT OBS_LON FCST OBS OBS_QC FCST_CAT OBS_CAT P1 P2 T1 T2 SEEPS NBRCNT : 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 FBS FBS_BCL FBS_BCU FSS FSS_BCL FSS_BCU AFSS AFSS_BCL AFSS_BCU UFSS UFSS_BCL UFSS_BCU F_RATE F_RATE_BCL F_RATE_BCU O_RATE O_RATE_BCL O_RATE_BCU NBRCTC : 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 FY_OY FY_ON FN_OY FN_ON diff --git a/internal/test_unit/xml/unit_grid_stat.xml b/internal/test_unit/xml/unit_grid_stat.xml index ad070dff60..e0de304f80 100644 --- a/internal/test_unit/xml/unit_grid_stat.xml +++ b/internal/test_unit/xml/unit_grid_stat.xml @@ -296,7 +296,7 @@ SEEPS_FLAG BOTH SEEPS_P1_THRESH NA OUTPUT_PREFIX SEEPS - SEEPS_GRID_CLIMO_NAME&DATA_DIR_CLIMO;/seeps/PPT24_seepsweights_grid.nc + SEEPS_GRID_CLIMO_NAME&DATA_DIR_CLIMO;/seeps/PPT24_seepsweights_grid_v12.0.nc \ &DATA_DIR_MODEL;/seeps/gpm_2021120100_2021120200_trmmgrid.nc \ diff --git a/src/basic/vx_util/stat_column_defs.h b/src/basic/vx_util/stat_column_defs.h index a7b5427143..984ef3c40b 100644 --- a/src/basic/vx_util/stat_column_defs.h +++ b/src/basic/vx_util/stat_column_defs.h @@ -383,9 +383,9 @@ static const char * const seeps_mpr_columns [] = { }; static const char * const seeps_columns [] = { - "TOTAL", "S12", "S13", - "S21", "S23", "S31", - "S32", "PF1", "PF2", + "TOTAL", "ODFL", "ODFH", + "OLFD", "OLFH", "OHFD", + "OHFL", "PF1", "PF2", "PF3", "PV1", "PV2", "PV3", "MEAN_FCST", "MEAN_OBS", "SEEPS" diff --git a/src/libcode/vx_analysis_util/stat_job.cc b/src/libcode/vx_analysis_util/stat_job.cc index 0552bc6d49..c10f392b85 100644 --- a/src/libcode/vx_analysis_util/stat_job.cc +++ b/src/libcode/vx_analysis_util/stat_job.cc @@ -1994,6 +1994,7 @@ void STATAnalysisJob::setup_stat_file(int n_row, int n) { case STATLineType::prc: c = get_n_prc_columns(n); break; case STATLineType::eclv: c = get_n_eclv_columns(n); break; case STATLineType::mpr: c = n_mpr_columns; break; + case STATLineType::seeps: c = n_seeps_columns; break; case STATLineType::nbrctc: c = n_nbrctc_columns; break; case STATLineType::nbrcts: c = n_nbrcts_columns; break; case STATLineType::nbrcnt: c = n_nbrcnt_columns; break; @@ -2065,6 +2066,7 @@ void STATAnalysisJob::setup_stat_file(int n_row, int n) { case STATLineType::prc: write_prc_header_row (1, n, stat_at, 0, 0); break; case STATLineType::eclv: write_eclv_header_row (1, n, stat_at, 0, 0); break; case STATLineType::mpr: write_header_row (mpr_columns, n_mpr_columns, 1, stat_at, 0, 0); break; + case STATLineType::seeps: write_header_row (seeps_columns, n_seeps_columns, 1, stat_at, 0, 0); break; case STATLineType::nbrctc: write_header_row (nbrctc_columns, n_nbrctc_columns, 1, stat_at, 0, 0); break; case STATLineType::nbrcts: write_header_row (nbrcts_columns, n_nbrcts_columns, 1, stat_at, 0, 0); break; case STATLineType::nbrcnt: write_header_row (nbrcnt_columns, n_nbrcnt_columns, 1, stat_at, 0, 0); break; diff --git a/src/libcode/vx_seeps/seeps.cc b/src/libcode/vx_seeps/seeps.cc index 4204f59370..a2d5184ec6 100644 --- a/src/libcode/vx_seeps/seeps.cc +++ b/src/libcode/vx_seeps/seeps.cc @@ -50,7 +50,7 @@ double weighted_average(double, double, double, double); //////////////////////////////////////////////////////////////////////// -SeepsClimo *get_seeps_climo(ConcatString seeps_point_climo_name) { +SeepsClimo *get_seeps_climo(const ConcatString &seeps_point_climo_name) { if (! seeps_climo) seeps_climo = new SeepsClimo(seeps_point_climo_name); return seeps_climo; } @@ -63,7 +63,7 @@ void release_seeps_climo() { //////////////////////////////////////////////////////////////////////// -SeepsClimoGrid *get_seeps_climo_grid(int month, ConcatString seeps_grid_climo_name, int hour) { +SeepsClimoGrid *get_seeps_climo_grid(int month, const ConcatString &seeps_grid_climo_name, int hour) { if (seeps_climo_grid_map_00.count(month) == 0) { seeps_climo_grid_map_00[month] = nullptr; @@ -96,19 +96,25 @@ double weighted_average(double v1, double w1, double v2, double w2) { void SeepsAggScore::clear() { + // Syntax used to define SEEPS obs/fcast categories (s_* or c_*) + // o{d|l|h} : obs in {dry(0)|light(1)|heavy(2)} category + // f{d|l|h} : fcsts in {dry(0)|light(1)|heavy(2)} category + n_obs = 0; - c12 = c13 = c21 = c23 = c31 = c32 = 0; - s12 = s13 = s21 = s23 = s31 = s32 = 0.; - pv1 = pv2 = pv3 = 0.; - pf1 = pf2 = pf3 = 0.; - mean_fcst = mean_obs = bad_data_double; - weighted_score = score = bad_data_double; + c_odfl = c_odfh = c_olfd = c_olfh = c_ohfd = c_ohfl = 0; + s_odfl = s_odfh = s_olfd = s_olfh = s_ohfd = s_ohfl = 0.0; + pv1 = pv2 = pv3 = 0.0; + pf1 = pf2 = pf3 = 0.0; + mean_fcst = mean_fcst_wgt = bad_data_double; + mean_obs = mean_obs_wgt = bad_data_double; + score = score_wgt = bad_data_double; } //////////////////////////////////////////////////////////////////////// SeepsAggScore & SeepsAggScore::operator+=(const SeepsAggScore &c) { + const char *method_name = "SeepsAggScore::operator+=() -> "; // Check for degenerate case if(n_obs == 0 && c.n_obs == 0) return *this; @@ -121,20 +127,23 @@ SeepsAggScore & SeepsAggScore::operator+=(const SeepsAggScore &c) { n_obs += c.n_obs; // Increment counts - c12 += c.c12; - c13 += c.c13; - c21 += c.c21; - c23 += c.c23; - c31 += c.c31; - c32 += c.c32; - + c_odfl += c.c_odfl; + c_odfh += c.c_odfh; + c_olfd += c.c_olfd; + c_olfh += c.c_olfh; + c_ohfd += c.c_ohfd; + c_ohfl += c.c_ohfl; + // Compute weighted averages - s12 = weighted_average(s12, w1, c.s12, w2); - s13 = weighted_average(s13, w1, c.s13, w2); - s21 = weighted_average(s21, w1, c.s21, w2); - s23 = weighted_average(s23, w1, c.s23, w2); - s31 = weighted_average(s31, w1, c.s31, w2); - s32 = weighted_average(s32, w1, c.s32, w2); + s_odfl = weighted_average(s_odfl, w1, c.s_odfl, w2); + s_odfh = weighted_average(s_odfh, w1, c.s_odfh, w2); + s_olfd = weighted_average(s_olfd, w1, c.s_olfd, w2); + s_olfh = weighted_average(s_olfh, w1, c.s_olfh, w2); + s_ohfd = weighted_average(s_ohfd, w1, c.s_ohfd, w2); + s_ohfl = weighted_average(s_ohfl, w1, c.s_ohfl, w2); + mlog << Debug(9) << method_name + << "s_odfl, o_odfh => " + << s_odfl << " " << s_odfh << "\n"; pv1 = weighted_average(pv1, w1, c.pv1, w2); pv2 = weighted_average(pv2, w1, c.pv2, w2); @@ -144,18 +153,21 @@ SeepsAggScore & SeepsAggScore::operator+=(const SeepsAggScore &c) { pf2 = weighted_average(pf2, w1, c.pf2, w2); pf3 = weighted_average(pf3, w1, c.pf3, w2); - mean_fcst = weighted_average(mean_fcst, w1, c.mean_fcst, w2); - mean_obs = weighted_average(mean_obs, w1, c.mean_obs, w2); + mean_fcst = weighted_average(mean_fcst, w1, c.mean_fcst, w2); + mean_fcst_wgt = weighted_average(mean_fcst_wgt, w1, c.mean_fcst_wgt, w2); - score = weighted_average(score, w1, c.score, w2); - weighted_score = weighted_average(weighted_score, w1, c.weighted_score, w2); + mean_obs = weighted_average(mean_obs, w1, c.mean_obs, w2); + mean_obs_wgt = weighted_average(mean_obs_wgt, w1, c.mean_obs_wgt, w2); + + score = weighted_average(score, w1, c.score, w2); + score_wgt = weighted_average(score_wgt, w1, c.score_wgt, w2); return *this; } //////////////////////////////////////////////////////////////////////// -SeepsClimoBase::SeepsClimoBase(ConcatString seeps_climo_name) : climo_file_name{seeps_climo_name} { +SeepsClimoBase::SeepsClimoBase(const ConcatString &seeps_climo_name) : climo_file_name{seeps_climo_name} { clear(); seeps_ready = false; @@ -192,7 +204,8 @@ ConcatString SeepsClimoBase::get_climo_filename() { seeps_ready = file_exists(seeps_filename.c_str()); if (seeps_ready) { - mlog << Debug(7) << method_name << "SEEPS climo name=\"" + mlog << Debug(7) << method_name + << "SEEPS climo name=\"" << seeps_filename.c_str() << "\"\n"; } else { @@ -229,11 +242,11 @@ void SeepsClimoBase::set_p1_thresh(const SingleThresh &p1_thresh) { //////////////////////////////////////////////////////////////////////// -SeepsClimo::SeepsClimo(ConcatString seeps_climo_name) : SeepsClimoBase{seeps_climo_name} { +SeepsClimo::SeepsClimo(const ConcatString &seeps_climo_name) : SeepsClimoBase{seeps_climo_name} { clear(); ConcatString seeps_name = get_climo_filename(); - if (file_exists(seeps_name.c_str())) read_seeps_scores(seeps_name); + if (file_exists(seeps_name.c_str())) read_seeps_climo_grid(seeps_name); } @@ -268,13 +281,17 @@ SeepsClimoRecord *SeepsClimo::create_climo_record( double *p1, double *p2, double *t1, double *t2, double *scores) { int offset; SeepsClimoRecord *record = new SeepsClimoRecord(); + const char *method_name = "SeepsClimo::create_climo_record() -> "; record->sid = sid; record->lat = lat; record->lon = lon; record->elv = elv; if (standalone_debug_seeps && SAMPLE_STATION_ID == sid) { - cout << " sid=" << sid << ", lat=" << lat << ", lon=" << lon << ", elv=" << elv << "\n"; + cout << method_name + << "sid=" << sid << ", lat=" << lat + << ", lon=" << lon << ", elv=" << elv + << "\n"; } for (int idx=0; idxp1[idx] = p1[idx]; @@ -283,7 +300,8 @@ SeepsClimoRecord *SeepsClimo::create_climo_record( record->t2[idx] = t2[idx]; if (standalone_debug_seeps && SAMPLE_STATION_ID == sid) { - cout << str_format("\t%2d: %6.3f %6.3f %6.3f %6.3f ", + cout << method_name + << str_format("\t%2d: %6.3f %6.3f %6.3f %6.3f ", (idx+1), record->p1[idx], record->p2[idx], record->t1[idx], record->t2[idx]); } @@ -291,7 +309,8 @@ SeepsClimoRecord *SeepsClimo::create_climo_record( offset = idx*SEEPS_MATRIX_SIZE + idx_m; record->scores[idx][idx_m] = scores[offset]; if (standalone_debug_seeps && SAMPLE_STATION_ID == sid) { - cout << str_format(" %.3f", record->scores[idx][idx_m]); + cout << method_name + << str_format(" %.3f", record->scores[idx][idx_m]); } } if (standalone_debug_seeps && SAMPLE_STATION_ID == sid) cout << "\n"; @@ -330,14 +349,24 @@ SeepsRecord *SeepsClimo::get_record(int sid, int month, int hour) { record->p2 = climo_record->p2[month-1]; record->t1 = climo_record->t1[month-1]; record->t2 = climo_record->t2[month-1]; + mlog << Debug(9) << method_name + << "record info: sid, lat, lon, month, p1, p2, t1, t2 => " + << record->sid << " " << record->lat << " " + << record->lon << " " << record->month << " " + << record->p1 << " " << record->p2 << " " + << record->t1 << " " << record->t2 << "\n"; for (int idx=0; idxscores[idx] = climo_record->scores[month-1][idx]; + record->scores[idx] = climo_record->scores[month-1][idx]; + mlog << Debug(7) << method_name + << "record info (SEEPS matrix): score => " + << record->scores[idx] << "\n"; } } else if (!is_eq(p1, bad_data_double)) { increase_filtered_count(); - mlog << Debug(7) << method_name << " filtered by threshold p1=" - << climo_record->p1[month-1] <<"\n"; + mlog << Debug(7) << method_name + << "filtered by threshold p1=" + << climo_record->p1[month-1] << "\n"; } } } @@ -349,16 +378,21 @@ SeepsRecord *SeepsClimo::get_record(int sid, int month, int hour) { << "or disable output for SEEPS and SEEPS_MPR.\n\n"; exit(1); } + mlog << Debug(9) << method_name + << "sid = " << sid + << ", month = " << month << ", hour = " << hour + << ", filtered_count = " << get_filtered_count() << "\n"; return record; } //////////////////////////////////////////////////////////////////////// -double SeepsClimo::get_score(int sid, double p_fcst, double p_obs, - int month, int hour) { +double SeepsClimo::get_seeps_category(int sid, double p_fcst, double p_obs, + int month, int hour) { double score = bad_data_double; SeepsRecord *record = get_record(sid, month, hour); + const char *method_name = "SeepsClimo::get_seeps_category() -> "; if (nullptr != record) { // Determine location in contingency table @@ -366,6 +400,10 @@ double SeepsClimo::get_score(int sid, double p_fcst, double p_obs, int jc = (p_fcst>record->t1)+(p_fcst>record->t2); score = record->scores[(jc*3)+ic]; + mlog << Debug(9) << method_name + << "ic, jc, score => " + << ic << " " << jc << " " + << score << "\n"; delete record; } @@ -378,6 +416,7 @@ SeepsScore *SeepsClimo::get_seeps_score(int sid, double p_fcst, double p_obs, int month, int hour) { SeepsScore *score = nullptr; SeepsRecord *record = get_record(sid, month, hour); + const char *method_name = "SeepsClimo::get_seeps_score() -> "; if (nullptr != record) { score = new SeepsScore(); @@ -385,18 +424,27 @@ SeepsScore *SeepsClimo::get_seeps_score(int sid, double p_fcst, double p_obs, score->p2 = record->p2; score->t1 = record->t1; score->t2 = record->t2; - + mlog << Debug(9) << method_name + << "p1, p2, t1, t2 => " + << score->p1 << " " << score->p2 << " " + << score->t1 << " " << score->t2 + << "\n"; + score->obs_cat = (p_obs>record->t1)+(p_obs>record->t2); score->fcst_cat = (p_fcst>record->t1)+(p_fcst>record->t2); score->s_idx = (score->fcst_cat*3)+score->obs_cat; score->score = record->scores[score->s_idx]; + mlog << Debug(9) << method_name + << "obs_cat, fc_cat, s_idx, score => " + << score->obs_cat << " " << score->fcst_cat + << " " << score->s_idx << " " << score->score + << "\n"; delete record; } return score; } - //////////////////////////////////////////////////////////////////////// void SeepsClimo::print_all() { @@ -457,9 +505,9 @@ void SeepsClimo::print_record(SeepsRecord *record, bool with_header) { //////////////////////////////////////////////////////////////////////// -void SeepsClimo::read_seeps_scores(ConcatString filename) { +void SeepsClimo::read_seeps_climo_grid(const ConcatString &filename) { clock_t clock_time = clock(); - const char *method_name = "SeepsClimo::read_records() -> "; + const char *method_name = "SeepsClimo::read_seeps_climo_grid() -> "; try { double p1_00_buf[SEEPS_MONTH]; @@ -478,23 +526,27 @@ void SeepsClimo::read_seeps_scores(ConcatString filename) { // dimensions: month = 12 ; nstn = 5293 ; nmatrix = 9 ; get_dim(nc_file, dim_name_nstn, nstn, true); - mlog << Debug(6) << method_name << "dimensions nstn = " << nstn << "\n"; - if (standalone_debug_seeps) cout << "dimensions nstn = " << nstn << "\n"; - - int *sid_array = new int[nstn]; - double *lat_array = new double[nstn]; - double *lon_array = new double[nstn]; - double *elv_array = new double[nstn]; - double *p1_00_array = new double[nstn*SEEPS_MONTH]; - double *p2_00_array = new double[nstn*SEEPS_MONTH]; - double *t1_00_array = new double[nstn*SEEPS_MONTH]; - double *t2_00_array = new double[nstn*SEEPS_MONTH]; - double *p1_12_array = new double[nstn*SEEPS_MONTH]; - double *p2_12_array = new double[nstn*SEEPS_MONTH]; - double *t1_12_array = new double[nstn*SEEPS_MONTH]; - double *t2_12_array = new double[nstn*SEEPS_MONTH]; - double *matrix_00_array = new double[nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE]; - double *matrix_12_array = new double[nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE]; + mlog << Debug(6) << method_name + << "dimensions nstn = " << nstn << "\n"; + if (standalone_debug_seeps) { + cout << method_name + << "dimensions nstn = " << nstn << "\n"; + } + + vector sid_array(nstn); + vector lat_array(nstn); + vector lon_array(nstn); + vector elv_array(nstn); + vector p1_00_array(nstn*SEEPS_MONTH); + vector p2_00_array(nstn*SEEPS_MONTH); + vector t1_00_array(nstn*SEEPS_MONTH); + vector t2_00_array(nstn*SEEPS_MONTH); + vector p1_12_array(nstn*SEEPS_MONTH); + vector p2_12_array(nstn*SEEPS_MONTH); + vector t1_12_array(nstn*SEEPS_MONTH); + vector t2_12_array(nstn*SEEPS_MONTH); + vector matrix_00_array(nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE); + vector matrix_12_array(nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE); NcVar var_sid = get_nc_var(nc_file, var_name_sid); NcVar var_lat = get_nc_var(nc_file, var_name_lat); @@ -511,72 +563,72 @@ void SeepsClimo::read_seeps_scores(ConcatString filename) { NcVar var_matrix_00 = get_nc_var(nc_file, var_name_matrix_00); NcVar var_matrix_12 = get_nc_var(nc_file, var_name_matrix_12); - if (IS_INVALID_NC(var_sid) || !get_nc_data(&var_sid, sid_array)) { + if (IS_INVALID_NC(var_sid) || !get_nc_data(&var_sid, sid_array.data())) { mlog << Error << "\n" << method_name << "Did not get sid\n\n"; exit(1); } - if (IS_INVALID_NC(var_lat) || !get_nc_data(&var_lat, lat_array)) { + if (IS_INVALID_NC(var_lat) || !get_nc_data(&var_lat, lat_array.data())) { mlog << Error << "\n" << method_name << "Did not get lat\n\n"; exit(1); } - if (IS_INVALID_NC(var_lon) || !get_nc_data(&var_lon, lon_array)) { + if (IS_INVALID_NC(var_lon) || !get_nc_data(&var_lon, lon_array.data())) { mlog << Error << "\n" << method_name << "Did not get lon\n\n"; exit(1); } - if (IS_INVALID_NC(var_elv) || !get_nc_data(&var_elv, elv_array)) { + if (IS_INVALID_NC(var_elv) || !get_nc_data(&var_elv, elv_array.data())) { mlog << Error << "\n" << method_name << "Did not get elv\n\n"; exit(1); } - if (IS_INVALID_NC(var_p1_00) || !get_nc_data(&var_p1_00, p1_00_array)) { + if (IS_INVALID_NC(var_p1_00) || !get_nc_data(&var_p1_00, p1_00_array.data())) { mlog << Error << "\n" << method_name << "Did not get p1_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_p2_00) || !get_nc_data(&var_p2_00, p2_00_array)) { + if (IS_INVALID_NC(var_p2_00) || !get_nc_data(&var_p2_00, p2_00_array.data())) { mlog << Error << "\n" << method_name << "Did not get p2_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_t1_00) || !get_nc_data(&var_t1_00, t1_00_array)) { + if (IS_INVALID_NC(var_t1_00) || !get_nc_data(&var_t1_00, t1_00_array.data())) { mlog << Error << "\n" << method_name << "Did not get t1_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_t2_00) || !get_nc_data(&var_t2_00, t2_00_array)) { + if (IS_INVALID_NC(var_t2_00) || !get_nc_data(&var_t2_00, t2_00_array.data())) { mlog << Error << "\n" << method_name << "Did not get t2_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_p1_12) || !get_nc_data(&var_p1_12, p1_12_array)) { + if (IS_INVALID_NC(var_p1_12) || !get_nc_data(&var_p1_12, p1_12_array.data())) { mlog << Error << "\n" << method_name << "Did not get p1_12\n\n"; exit(1); } - if (IS_INVALID_NC(var_p2_12) || !get_nc_data(&var_p2_12, p2_12_array)) { + if (IS_INVALID_NC(var_p2_12) || !get_nc_data(&var_p2_12, p2_12_array.data())) { mlog << Error << "\n" << method_name << "Did not get p2_12\n\n"; exit(1); } - if (IS_INVALID_NC(var_t1_12) || !get_nc_data(&var_t1_12, t1_12_array)) { + if (IS_INVALID_NC(var_t1_12) || !get_nc_data(&var_t1_12, t1_12_array.data())) { mlog << Error << "\n" << method_name << "Did not get t1_12\n\n"; exit(1); } - if (IS_INVALID_NC(var_t2_12) || !get_nc_data(&var_t2_12, t2_12_array)) { + if (IS_INVALID_NC(var_t2_12) || !get_nc_data(&var_t2_12, t2_12_array.data())) { mlog << Error << "\n" << method_name << "Did not get t2_12\n\n"; exit(1); } - if (IS_INVALID_NC(var_matrix_00) || !get_nc_data(&var_matrix_00, matrix_00_array)) { + if (IS_INVALID_NC(var_matrix_00) || !get_nc_data(&var_matrix_00, matrix_00_array.data())) { mlog << Error << "\n" << method_name << "Did not get matrix_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_matrix_12) || !get_nc_data(&var_matrix_12, matrix_12_array)) { + if (IS_INVALID_NC(var_matrix_12) || !get_nc_data(&var_matrix_12, matrix_12_array.data())) { mlog << Error << "\n" << method_name << "Did not get matrix_12\n\n"; exit(1); @@ -611,26 +663,15 @@ void SeepsClimo::read_seeps_scores(ConcatString filename) { seeps_score_12_map[sid] = rec_12; } - if (sid_array) delete [] sid_array; - if (lat_array) delete [] lat_array; - if (lon_array) delete [] lon_array; - if (elv_array) delete [] elv_array; - if (p1_00_array) delete [] p1_00_array; - if (p2_00_array) delete [] p2_00_array; - if (t1_00_array) delete [] t1_00_array; - if (t2_00_array) delete [] t2_00_array; - if (p1_12_array) delete [] p1_12_array; - if (p2_12_array) delete [] p2_12_array; - if (t1_12_array) delete [] t1_12_array; - if (t2_12_array) delete [] t2_12_array; - if (matrix_00_array) delete [] matrix_00_array; - if (matrix_12_array) delete [] matrix_12_array; - nc_file->close(); float duration = (float)(clock() - clock_time)/CLOCKS_PER_SEC; - mlog << Debug(6) << method_name << "took " << duration << " seconds\n"; - if (standalone_debug_seeps) cout << method_name << "took " << duration << " seconds\n"; + mlog << Debug(6) << method_name + << "took " << duration << " seconds\n"; + if (standalone_debug_seeps) { + cout << method_name + << "took " << duration << " seconds\n"; + } } catch(int i_err) { @@ -645,12 +686,13 @@ void SeepsClimo::read_seeps_scores(ConcatString filename) { //////////////////////////////////////////////////////////////////////// -SeepsClimoGrid::SeepsClimoGrid(int month, int hour, ConcatString seeps_climo_name) +SeepsClimoGrid::SeepsClimoGrid(int month, int hour, const ConcatString &seeps_climo_name) : month{month}, hour{hour}, SeepsClimoBase{seeps_climo_name} { - init_from_scratch(); + clear(); + ConcatString seeps_name = get_climo_filename(); - if (file_exists(seeps_name.c_str())) read_seeps_scores(seeps_name); + if (file_exists(seeps_name.c_str())) read_seeps_climo_grid(seeps_name); } //////////////////////////////////////////////////////////////////////// @@ -662,25 +704,25 @@ SeepsClimoGrid::~SeepsClimoGrid() { //////////////////////////////////////////////////////////////////////// void SeepsClimoGrid::init_from_scratch() { - p1_buf = p2_buf = t1_buf = t2_buf = nullptr; - s12_buf = s13_buf = s21_buf = s23_buf = s31_buf = s32_buf = nullptr; clear(); } //////////////////////////////////////////////////////////////////////// void SeepsClimoGrid::clear() { + SeepsClimoBase::clear(); - if (nullptr != p1_buf) { delete [] p1_buf; p1_buf = nullptr; } - if (nullptr != p2_buf) { delete [] p2_buf; p2_buf = nullptr; } - if (nullptr != t1_buf) { delete [] t1_buf; t1_buf = nullptr; } - if (nullptr != t2_buf) { delete [] t2_buf; t2_buf = nullptr; } - if (nullptr != s12_buf) { delete [] s12_buf; s12_buf = nullptr; } - if (nullptr != s13_buf) { delete [] s13_buf; s13_buf = nullptr; } - if (nullptr != s21_buf) { delete [] s21_buf; s21_buf = nullptr; } - if (nullptr != s23_buf) { delete [] s23_buf; s23_buf = nullptr; } - if (nullptr != s31_buf) { delete [] s31_buf; s31_buf = nullptr; } - if (nullptr != s32_buf) { delete [] s32_buf; s32_buf = nullptr; } + + p1_buf.clear(); + p2_buf.clear(); + t1_buf.clear(); + t2_buf.clear(); + s_odfl_buf.clear(); + s_odfh_buf.clear(); + s_olfd_buf.clear(); + s_olfh_buf.clear(); + s_ohfd_buf.clear(); + s_ohfl_buf.clear(); }; //////////////////////////////////////////////////////////////////////// @@ -689,6 +731,7 @@ SeepsScore *SeepsClimoGrid::get_record(int ix, int iy, double p_fcst, double p_obs) { SeepsScore *seeps_record = nullptr; const char *method_name = "SeepsClimoGrid::get_record() -> "; + if (!is_eq(p_fcst, -9999.0) && !is_eq(p_obs, -9999.0)) { int offset = iy * nx + ix; double p1 = p1_buf[offset]; @@ -697,7 +740,7 @@ SeepsScore *SeepsClimoGrid::get_record(int ix, int iy, // Determine location in contingency table int ic = (p_obs>t1_buf[offset])+(p_obs>t2_buf[offset]); int jc = (p_fcst>t1_buf[offset])+(p_fcst>t2_buf[offset]); - double score = get_score(offset, ic, jc); + double score = get_seeps_score(offset, ic, jc); seeps_record = new SeepsScore(); seeps_record->obs_cat = ic; @@ -711,7 +754,9 @@ SeepsScore *SeepsClimoGrid::get_record(int ix, int iy, } else if (~is_eq(p1, bad_data_double)) { increase_filtered_count(); - mlog << Debug(7) << method_name << " filtered by threshold p1=" << p1_buf[offset] <<"\n"; + mlog << Debug(7) << method_name + << " filtered by threshold p1=" + << p1_buf[offset] << "\n"; } } @@ -720,55 +765,46 @@ SeepsScore *SeepsClimoGrid::get_record(int ix, int iy, //////////////////////////////////////////////////////////////////////// -double SeepsClimoGrid::get_score(int offset, int obs_cat, int fcst_cat) { +double SeepsClimoGrid::get_seeps_score(int offset, int obs_cat, int fcst_cat) { double score = bad_data_double; + const char *method_name = "SeepsClimoGrid::get_seeps_score() -> "; - if (offset >= (nx * ny)) { - mlog << Error << "\nSeepsClimoGrid::get_score() --> offset (" << offset - << " is too big (" << (nx*ny) << ")\n"; + if (offset < 0 || offset >= (nx * ny)) { + mlog << Error << method_name + << "offset (" << offset << ") is too big (" + << (nx*ny) << ")\n"; return score; } + // Place climate score depending on obs and forecast categories if (obs_cat == 0) { - if (fcst_cat == 1) score = s12_buf[offset]; - else if (fcst_cat == 2) score = s13_buf[offset]; + if (fcst_cat == 1) score = s_odfl_buf[offset]; + else if (fcst_cat == 2) score = s_odfh_buf[offset]; else score = 0.; } else if (obs_cat == 1) { - if (fcst_cat == 0) score = s21_buf[offset]; - else if (fcst_cat == 2) score = s23_buf[offset]; + if (fcst_cat == 0) score = s_olfd_buf[offset]; + else if (fcst_cat == 2) score = s_olfh_buf[offset]; else score = 0.; } else { - if (fcst_cat == 0) score = s31_buf[offset]; - else if (fcst_cat == 1) score = s32_buf[offset]; + if (fcst_cat == 0) score = s_ohfd_buf[offset]; + else if (fcst_cat == 1) score = s_ohfl_buf[offset]; else score = 0.; } + mlog << Debug(9) << method_name + << "obs_cat = " << obs_cat + << ", fcst_cat = " << fcst_cat + << ", score = " << score << "\n"; return score; } //////////////////////////////////////////////////////////////////////// -double SeepsClimoGrid::get_score(int ix, int iy, double p_fcst, double p_obs) { - double score = bad_data_double; - - if (!is_eq(p_fcst, -9999.0) && !is_eq(p_obs, -9999.0)) { - int offset = iy * nx + ix; - // Determine location in contingency table - int ic = (p_obs>t1_buf[offset])+(p_obs>t2_buf[offset]); - int jc = (p_fcst>t1_buf[offset])+(p_fcst>t2_buf[offset]); - score = get_score(offset, ic, jc); - } - - return score; -} - -//////////////////////////////////////////////////////////////////////// - -void SeepsClimoGrid::read_seeps_scores(ConcatString filename) { +void SeepsClimoGrid::read_seeps_climo_grid(const ConcatString &filename) { clock_t clock_time = clock(); - const char *method_name = "SeepsClimoGrid::read_seeps_scores() -> "; + const char *method_name = "SeepsClimoGrid::read_seeps_climo_grid() -> "; try { NcFile *nc_file = open_ncfile(filename.c_str()); @@ -777,26 +813,33 @@ void SeepsClimoGrid::read_seeps_scores(ConcatString filename) { if (!has_dim(nc_file, dim_name_lat) || !has_dim(nc_file, dim_name_lon)) { mlog << Error << "\n" << method_name << "\"" << filename << "\" is not valid SEEPS climo file\n\n"; - //exit(1); + exit(1); } get_dim(nc_file, dim_name_lat, ny, true); get_dim(nc_file, dim_name_lon, nx, true); - mlog << Debug(6) << method_name << "dimensions lon = " << nx << " lat = " << ny + mlog << Debug(6) << method_name + << "dimensions lon = " << nx << " lat = " << ny << " month=" << month << "\n"; - if (standalone_debug_seeps) cout << "dimensions lon = " << nx << " lat = " << ny - << " month=" << month << "\n";; - - p1_buf = new double[nx*ny]; - p2_buf = new double[nx*ny]; - t1_buf = new double[nx*ny]; - t2_buf = new double[nx*ny]; - s12_buf = new double[nx*ny]; - s13_buf = new double[nx*ny]; - s21_buf = new double[nx*ny]; - s23_buf = new double[nx*ny]; - s31_buf = new double[nx*ny]; - s32_buf = new double[nx*ny]; + if (standalone_debug_seeps) { + cout << method_name + << "dimensions lon = " << nx << " lat = " << ny + << " month=" << month << "\n"; + } + + // Variables in climo file named as s_odfl, s_odfh etc. These then stored + // into new convention s_odfl, s_odfh etc. + + p1_buf.resize(nx*ny); + p2_buf.resize(nx*ny); + t1_buf.resize(nx*ny); + t2_buf.resize(nx*ny); + s_odfl_buf.resize(nx*ny); + s_odfh_buf.resize(nx*ny); + s_olfd_buf.resize(nx*ny); + s_olfh_buf.resize(nx*ny); + s_ohfd_buf.resize(nx*ny); + s_ohfl_buf.resize(nx*ny); LongArray curs; // = { month-1, 0, 0 }; LongArray dims; // = { 1, ny, nx }; @@ -804,12 +847,12 @@ void SeepsClimoGrid::read_seeps_scores(ConcatString filename) { NcVar var_p2_00 = get_nc_var(nc_file, var_name_p2_00); NcVar var_t1_00 = get_nc_var(nc_file, var_name_t1_00); NcVar var_t2_00 = get_nc_var(nc_file, var_name_t2_00); - NcVar var_s12_00 = get_nc_var(nc_file, var_name_s12_00); - NcVar var_s13_00 = get_nc_var(nc_file, var_name_s13_00); - NcVar var_s21_00 = get_nc_var(nc_file, var_name_s21_00); - NcVar var_s23_00 = get_nc_var(nc_file, var_name_s23_00); - NcVar var_s31_00 = get_nc_var(nc_file, var_name_s31_00); - NcVar var_s32_00 = get_nc_var(nc_file, var_name_s32_00); + NcVar var_odfl_00 = get_nc_var(nc_file, var_name_odfl_00); + NcVar var_odfh_00 = get_nc_var(nc_file, var_name_odfh_00); + NcVar var_olfd_00 = get_nc_var(nc_file, var_name_olfd_00); + NcVar var_olfh_00 = get_nc_var(nc_file, var_name_olfh_00); + NcVar var_ohfd_00 = get_nc_var(nc_file, var_name_ohfd_00); + NcVar var_ohfl_00 = get_nc_var(nc_file, var_name_ohfl_00); curs.add(month-1); curs.add(0); @@ -818,73 +861,86 @@ void SeepsClimoGrid::read_seeps_scores(ConcatString filename) { dims.add(ny); dims.add(nx); - if (IS_INVALID_NC(var_p1_00) || !get_nc_data(&var_p1_00, p1_buf, dims, curs)) { + mlog << Debug(9) << method_name + << "var_odfl_00 = " << &var_odfl_00 << "\n"; + + if (IS_INVALID_NC(var_p1_00) || !get_nc_data(&var_p1_00, p1_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name << "Did not get p1_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_p2_00) || !get_nc_data(&var_p2_00, p2_buf, dims, curs)) { + if (IS_INVALID_NC(var_p2_00) || !get_nc_data(&var_p2_00, p2_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name << "Did not get p2_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_t1_00) || !get_nc_data(&var_t1_00, t1_buf, dims, curs)) { + if (IS_INVALID_NC(var_t1_00) || !get_nc_data(&var_t1_00, t1_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name << "Did not get t1_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_t2_00) || !get_nc_data(&var_t2_00, t2_buf, dims, curs)) { + if (IS_INVALID_NC(var_t2_00) || !get_nc_data(&var_t2_00, t2_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name << "Did not get t2_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_s12_00) || !get_nc_data(&var_s12_00, s12_buf, dims, curs)) { + if (IS_INVALID_NC(var_odfl_00) || !get_nc_data(&var_odfl_00, s_odfl_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name - << "Did not get s12_00\n\n"; + << "Did not get odfl_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_s13_00) || !get_nc_data(&var_s12_00, s13_buf, dims, curs)) { + if (IS_INVALID_NC(var_odfh_00) || !get_nc_data(&var_odfh_00, s_odfh_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name - << "Did not get s13_00\n\n"; + << "Did not get odfh_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_s21_00) || !get_nc_data(&var_s21_00, s21_buf, dims, curs)) { + if (IS_INVALID_NC(var_olfd_00) || !get_nc_data(&var_olfd_00, s_olfd_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name - << "Did not get s21_00\n\n"; + << "Did not get olfd_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_s23_00) || !get_nc_data(&var_s23_00, s23_buf, dims, curs)) { + if (IS_INVALID_NC(var_olfh_00) || !get_nc_data(&var_olfh_00, s_olfh_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name - << "Did not get s23_00\n\n"; + << "Did not get olfh_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_s31_00) || !get_nc_data(&var_s31_00, s31_buf, dims, curs)) { + if (IS_INVALID_NC(var_ohfd_00) || !get_nc_data(&var_ohfd_00, s_ohfd_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name - << "Did not get s31_00\n\n"; + << "Did not get ohfd_00\n\n"; exit(1); } - if (IS_INVALID_NC(var_s32_00) || !get_nc_data(&var_s32_00, s32_buf, dims, curs)) { + if (IS_INVALID_NC(var_ohfl_00) || !get_nc_data(&var_ohfl_00, s_ohfl_buf.data(), dims, curs)) { mlog << Error << "\n" << method_name - << "Did not get s32_00\n\n"; + << "Did not get ohfl_00\n\n"; exit(1); } nc_file->close(); + for(int i = 0; i < ny+3; i++) { + mlog << Debug(9) << method_name + << "s_odfl_buf[" << i << "] = " << s_odfl_buf[i] << "\n"; + } + float duration = (float)(clock() - clock_time)/CLOCKS_PER_SEC; - mlog << Debug(6) << method_name << "took " << duration << " seconds\n"; - if (standalone_debug_seeps) cout << method_name << "took " << duration << " seconds\n"; + mlog << Debug(6) << method_name + << "took " << duration << " seconds\n"; + if (standalone_debug_seeps) { + cout << method_name + << "took " << duration << " seconds\n"; + } } catch(...) { set_seeps_ready(false); mlog << Error << "\n" << method_name << "encountered an error on reading " << filename << ".\n\n"; - exit(-1); + } // end catch block } + //////////////////////////////////////////////////////////////////////// void SeepsClimoGrid::print_all() { @@ -896,36 +952,25 @@ void SeepsClimoGrid::print_all() { cout << method_name << " p2_buf[" << offset << "] = " << p2_buf[offset] << "\n"; cout << method_name << " t1_buf[" << offset << "] = " << t1_buf[offset] << "\n"; cout << method_name << " t2_buf[" << offset << "] = " << t2_buf[offset] << "\n"; - cout << method_name << "s12_buf[" << offset << "] = " << s12_buf[offset] << "\n"; - cout << method_name << "s13_buf[" << offset << "] = " << s13_buf[offset] << "\n"; - cout << method_name << "s21_buf[" << offset << "] = " << s21_buf[offset] << "\n"; - cout << method_name << "s23_buf[" << offset << "] = " << s23_buf[offset] << "\n"; - cout << method_name << "s31_buf[" << offset << "] = " << s31_buf[offset] << "\n"; - cout << method_name << "s32_buf[" << offset << "] = " << s32_buf[offset] << "\n"; - - offset = 400; - cout << method_name << " p1_buf[" << offset << "] = " << p1_buf[offset] << "\n"; - cout << method_name << " p2_buf[" << offset << "] = " << p2_buf[offset] << "\n"; - cout << method_name << " t1_buf[" << offset << "] = " << t1_buf[offset] << "\n"; - cout << method_name << " t2_buf[" << offset << "] = " << t2_buf[offset] << "\n"; - cout << method_name << "s12_buf[" << offset << "] = " << s12_buf[offset] << "\n"; - cout << method_name << "s13_buf[" << offset << "] = " << s13_buf[offset] << "\n"; - cout << method_name << "s21_buf[" << offset << "] = " << s21_buf[offset] << "\n"; - cout << method_name << "s23_buf[" << offset << "] = " << s23_buf[offset] << "\n"; - cout << method_name << "s31_buf[" << offset << "] = " << s31_buf[offset] << "\n"; - cout << method_name << "s32_buf[" << offset << "] = " << s32_buf[offset] << "\n"; - + cout << method_name << "s_odfl_buf[" << offset << "] = " << s_odfl_buf[offset] << "\n"; + cout << method_name << "s_odfh_buf[" << offset << "] = " << s_odfh_buf[offset] << "\n"; + cout << method_name << "s_olfd_buf[" << offset << "] = " << s_olfd_buf[offset] << "\n"; + cout << method_name << "s_olfh_buf[" << offset << "] = " << s_olfh_buf[offset] << "\n"; + cout << method_name << "s_ohfd_buf[" << offset << "] = " << s_ohfd_buf[offset] << "\n"; + cout << method_name << "s_ohfl_buf[" << offset << "] = " << s_ohfl_buf[offset] << "\n"; + offset = (nx*ny) - 1; cout << method_name << " p1_buf[" << offset << "] = " << p1_buf[offset] << "\n"; cout << method_name << " p2_buf[" << offset << "] = " << p2_buf[offset] << "\n"; cout << method_name << " t1_buf[" << offset << "] = " << t1_buf[offset] << "\n"; cout << method_name << " t2_buf[" << offset << "] = " << t2_buf[offset] << "\n"; - cout << method_name << "s12_buf[" << offset << "] = " << s12_buf[offset] << "\n"; - cout << method_name << "s13_buf[" << offset << "] = " << s13_buf[offset] << "\n"; - cout << method_name << "s21_buf[" << offset << "] = " << s21_buf[offset] << "\n"; - cout << method_name << "s23_buf[" << offset << "] = " << s23_buf[offset] << "\n"; - cout << method_name << "s31_buf[" << offset << "] = " << s31_buf[offset] << "\n"; - cout << method_name << "s32_buf[" << offset << "] = " << s32_buf[offset] << "\n"; + cout << method_name << "s_odfl_buf[" << offset << "] = " << s_odfl_buf[offset] << "\n"; + cout << method_name << "s_odfh_buf[" << offset << "] = " << s_odfh_buf[offset] << "\n"; + cout << method_name << "s_olfd_buf[" << offset << "] = " << s_olfd_buf[offset] << "\n"; + cout << method_name << "s_olfh_buf[" << offset << "] = " << s_olfh_buf[offset] << "\n"; + cout << method_name << "s_ohfd_buf[" << offset << "] = " << s_ohfd_buf[offset] << "\n"; + cout << method_name << "s_ohfl_buf[" << offset << "] = " << s_ohfl_buf[offset] << "\n"; + } } diff --git a/src/libcode/vx_seeps/seeps.h b/src/libcode/vx_seeps/seeps.h index 808acefa3e..6f93ddffcb 100644 --- a/src/libcode/vx_seeps/seeps.h +++ b/src/libcode/vx_seeps/seeps.h @@ -39,19 +39,18 @@ constexpr char var_name_t1_12[] = "t1_12"; constexpr char var_name_t2_12[] = "t2_12"; constexpr char var_name_matrix_00[] = "matrix_00"; constexpr char var_name_matrix_12[] = "matrix_12"; -constexpr char var_name_s12_00[] = "s12_00"; -constexpr char var_name_s13_00[] = "s13_00"; -constexpr char var_name_s21_00[] = "s21_00"; -constexpr char var_name_s23_00[] = "s23_00"; -constexpr char var_name_s31_00[] = "s31_00"; -constexpr char var_name_s32_00[] = "s32_00"; -constexpr char var_name_s12_12[] = "s12_12"; -constexpr char var_name_s13_12[] = "s13_12"; -constexpr char var_name_s21_12[] = "s21_12"; -constexpr char var_name_s23_12[] = "s23_12"; -constexpr char var_name_s31_12[] = "s31_12"; -constexpr char var_name_s32_12[] = "s32_12"; - +constexpr char var_name_odfl_00[] = "odfl_00"; +constexpr char var_name_odfh_00[] = "odfh_00"; +constexpr char var_name_olfd_00[] = "olfd_00"; +constexpr char var_name_olfh_00[] = "olfh_00"; +constexpr char var_name_ohfd_00[] = "ohfd_00"; +constexpr char var_name_ohfl_00[] = "ohfl_00"; +constexpr char var_name_odfl_12[] = "odfl_12"; +constexpr char var_name_odfh_12[] = "odfh_12"; +constexpr char var_name_olfd_12[] = "olfd_12"; +constexpr char var_name_olfh_12[] = "olfh_12"; +constexpr char var_name_ohfd_12[] = "ohfd_12"; +constexpr char var_name_ohfl_12[] = "ohfl_12"; constexpr char def_seeps_point_filename[] = "MET_BASE/climo/seeps/PPT24_seepsweights.nc"; constexpr char def_seeps_grid_filename[] = @@ -65,8 +64,8 @@ const double density_radius_rad = density_radius * rad_per_deg; //////////////////////////////////////////////////////////////////////// struct SeepsScore { // For SEEPS_MPR - int obs_cat; // i = obs category 0,1,2 - int fcst_cat; // j = model category 0,1,2 + int obs_cat; // i = obs category 0,1,2 (dry, light, heavy) + int fcst_cat; // j = model category 0,1,2 (dry, light, heavy) int s_idx; // index for 3 by 3 matrix as 1 dimensional (fcst_cat*3)+obs_cat double p1; double p2; @@ -81,19 +80,19 @@ struct SeepsAggScore { // For SEEPS void clear(); SeepsAggScore & operator+=(const SeepsAggScore &); - int n_obs; - int c12; - int c13; - int c21; - int c23; - int c31; - int c32; - double s12; - double s13; - double s21; - double s23; - double s31; - double s32; + int n_obs; + int c_odfl; + int c_odfh; + int c_olfd; + int c_olfh; + int c_ohfd; + int c_ohfl; + double s_odfl; + double s_odfh; + double s_olfd; + double s_olfh; + double s_ohfd; + double s_ohfl; double pv1; // marginal probabilities of the observed values double pv2; double pv3; @@ -101,9 +100,11 @@ struct SeepsAggScore { // For SEEPS double pf2; double pf3; double mean_fcst; + double mean_fcst_wgt; double mean_obs; + double mean_obs_wgt; double score; - double weighted_score; + double score_wgt; }; //////////////////////////////////////////////////////////////////////// @@ -156,12 +157,12 @@ class SeepsClimoBase { virtual void clear(); virtual ConcatString get_env_climo_name() { return "not defined"; }; virtual char *get_def_climo_name() { return nullptr; }; - virtual void read_seeps_scores(ConcatString filename) {}; + virtual void read_seeps_climo_grid(const ConcatString &filename) {}; void set_seeps_ready(bool _seeps_ready) { seeps_ready = _seeps_ready; }; public: - SeepsClimoBase(ConcatString seeps_climo_name); + SeepsClimoBase(const ConcatString &seeps_climo_name); virtual ~SeepsClimoBase(); void set_p1_thresh(const SingleThresh &p1_thresh); int get_filtered_count() const; @@ -182,22 +183,25 @@ class SeepsClimo : public SeepsClimoBase { double *p1, double *p2, double *t1, double *t2, double *scores); void print_record(SeepsClimoRecord *record, bool with_header=false); - void read_records(ConcatString filename); + void read_records(const ConcatString &filename); protected: void clear() override; ConcatString get_env_climo_name() override { return MET_ENV_SEEPS_POINT_CLIMO_NAME; }; char *get_def_climo_name() override { return (char *)def_seeps_point_filename; }; - void read_seeps_scores(ConcatString filename) override; + void read_seeps_climo_grid(const ConcatString &filename) override; public: - SeepsClimo(ConcatString seeps_climo_name); + SeepsClimo(const ConcatString &seeps_climo_name); ~SeepsClimo(); SeepsRecord *get_record(int sid, int month, int hour); - double get_score(int sid, double p_fcst, double p_obs, int month, int hour); - SeepsScore *get_seeps_score(int sid, double p_fcst, double p_obs, int month, int hour); + double get_seeps_category(int sid, double p_fcst, double p_obs, + int month, int hour); + SeepsScore *get_seeps_score(int sid, double p_fcst, double p_obs, + int month, int hour); + void print_all(); void print_record(SeepsRecord *record, bool with_header=false); @@ -205,8 +209,6 @@ class SeepsClimo : public SeepsClimoBase { // // - SeepsRecord get_seeps_record(int sid) const; - }; //////////////////////////////////////////////////////////////////////// @@ -219,16 +221,16 @@ class SeepsClimoGrid : public SeepsClimoBase { int hour; // not implemented int nx; int ny; - double *p1_buf; - double *p2_buf; - double *t1_buf; - double *t2_buf; - double *s12_buf; - double *s13_buf; - double *s21_buf; - double *s23_buf; - double *s31_buf; - double *s32_buf; + std::vector p1_buf; + std::vector p2_buf; + std::vector t1_buf; + std::vector t2_buf; + std::vector s_odfl_buf; + std::vector s_odfh_buf; + std::vector s_olfd_buf; + std::vector s_olfh_buf; + std::vector s_ohfd_buf; + std::vector s_ohfl_buf; void init_from_scratch(); @@ -236,16 +238,15 @@ class SeepsClimoGrid : public SeepsClimoBase { void clear() override; ConcatString get_env_climo_name() override { return MET_ENV_SEEPS_GRID_CLIMO_NAME; }; char *get_def_climo_name() override { return (char *)def_seeps_grid_filename; }; - void read_seeps_scores(ConcatString filename) override; + void read_seeps_climo_grid(const ConcatString &filename) override; public: - SeepsClimoGrid(int month, int hour, ConcatString seeps_climo_name); + SeepsClimoGrid(int month, int hour, const ConcatString &seeps_climo_name); ~SeepsClimoGrid(); SeepsScore *get_record(int ix, int iy, double p_fcst, double p_obs); - double get_score(int offset, int obs_cat, int fcst_cat); - double get_score(int ix, int iy, double p_fcst, double p_obs); + double get_seeps_score(int offset, int obs_cat, int fcst_cat); void print_all(); // @@ -261,8 +262,8 @@ inline int SeepsClimoBase::get_filtered_count() const { return filtered_count; } //////////////////////////////////////////////////////////////////////// -extern SeepsClimo *get_seeps_climo(ConcatString seeps_point_climo_name); -extern SeepsClimoGrid *get_seeps_climo_grid(int month, ConcatString seeps_grid_climo_name, int hour=0); +extern SeepsClimo *get_seeps_climo(const ConcatString &seeps_point_climo_name); +extern SeepsClimoGrid *get_seeps_climo_grid(int month, const ConcatString &seeps_grid_climo_name, int hour=0); extern void release_seeps_climo(); extern void release_seeps_climo_grid(); diff --git a/src/libcode/vx_stat_out/stat_columns.cc b/src/libcode/vx_stat_out/stat_columns.cc index f4ccb6354b..d971fe2cf7 100644 --- a/src/libcode/vx_stat_out/stat_columns.cc +++ b/src/libcode/vx_stat_out/stat_columns.cc @@ -4162,35 +4162,40 @@ void write_seeps_cols(const SeepsAggScore *seeps, // // Stable Equitable Error in Probability Space (SEEPS) // Dump out the SEEPS line: - // TOTAL S12, S13, - // S21, S23, S31, - // S32, PF1, PF2, - // PF3, PV1, PV2, - // PV3, MEAN_FCST, MEAN_OBS, - // SEEPS + // TOTAL, + // ODFL, ODFH, OLFD, + // OLFH, OHFD, OHFL, + // PF1, PF2, PF3, + // PV1, PV2, PV3, + // MEAN_FCST, MEAN_OBS, SEEPS // at.set_entry(r, c+0, seeps->n_obs); // Total Number of Pairs - at.set_entry(r, c+1, seeps->s12); // s12 - at.set_entry(r, c+2, seeps->s13); // s13 - at.set_entry(r, c+3, seeps->s21); // s21 - at.set_entry(r, c+4, seeps->s23); // s23 - at.set_entry(r, c+5, seeps->s31); // s31 - at.set_entry(r, c+6, seeps->s32); // s32 - - at.set_entry(r, c+7, seeps->pf1); // pf1 - at.set_entry(r, c+8, seeps->pf2); // pf2 - at.set_entry(r, c+9, seeps->pf3); // pf3 - - at.set_entry(r, c+10, seeps->pv1); // pv1 - at.set_entry(r, c+11, seeps->pv2); // pv2 - at.set_entry(r, c+12, seeps->pv3); // pv3 - - at.set_entry(r, c+13, seeps->mean_fcst); // mean_fcst - at.set_entry(r, c+14, seeps->mean_obs); // mean_obs - - at.set_entry(r, c+15, (use_weighted_seeps ? seeps->weighted_score : seeps->score)); // SEEPS score/weighted score + at.set_entry(r, c+1, seeps->s_odfl); // ODFL + at.set_entry(r, c+2, seeps->s_odfh); // ODFH + at.set_entry(r, c+3, seeps->s_olfd); // OLFD + at.set_entry(r, c+4, seeps->s_olfh); // OLFH + at.set_entry(r, c+5, seeps->s_ohfd); // OHFD + at.set_entry(r, c+6, seeps->s_ohfl); // OHFL + + at.set_entry(r, c+7, seeps->pf1); // PF1 + at.set_entry(r, c+8, seeps->pf2); // PF2 + at.set_entry(r, c+9, seeps->pf3); // PF3 + + at.set_entry(r, c+10, seeps->pv1); // PV1 + at.set_entry(r, c+11, seeps->pv2); // PV2 + at.set_entry(r, c+12, seeps->pv3); // PV3 + + at.set_entry(r, c+13, (use_weighted_seeps ? + seeps->mean_fcst_wgt : + seeps->mean_fcst)); // MEAN_FCST + at.set_entry(r, c+14, (use_weighted_seeps ? + seeps->mean_obs_wgt : + seeps->mean_obs)); // MEAN_OBS + at.set_entry(r, c+15, (use_weighted_seeps ? + seeps->score_wgt: + seeps->score)); // SEEPS return; } @@ -4219,18 +4224,18 @@ void write_seeps_mpr_cols(const PairDataPoint *pd_ptr, int i, at.set_entry(r, c+4, pd_ptr->o_na[i]); // Observation Value - at.set_entry(r, c+5, (string)pd_ptr->o_qc_sa[i]); // Observation Quality Control + at.set_entry(r, c+5, (string)pd_ptr->o_qc_sa[i]); // Observation Quality Control - at.set_entry(r, c+6, pd_ptr->seeps_mpr[i]->fcst_cat); // model category - at.set_entry(r, c+7, pd_ptr->seeps_mpr[i]->obs_cat); // observation category + at.set_entry(r, c+6, pd_ptr->seeps_mpr[i]->fcst_cat); // FCST_CAT + at.set_entry(r, c+7, pd_ptr->seeps_mpr[i]->obs_cat); // OBS_CAT - at.set_entry(r, c+8, pd_ptr->seeps_mpr[i]->p1); // p1 - at.set_entry(r, c+9, pd_ptr->seeps_mpr[i]->p2); // p2 + at.set_entry(r, c+8, pd_ptr->seeps_mpr[i]->p1); // P1 + at.set_entry(r, c+9, pd_ptr->seeps_mpr[i]->p2); // P2 - at.set_entry(r, c+10, pd_ptr->seeps_mpr[i]->t1); // t1 - at.set_entry(r, c+11, pd_ptr->seeps_mpr[i]->t2); // t2 + at.set_entry(r, c+10, pd_ptr->seeps_mpr[i]->t1); // T1 + at.set_entry(r, c+11, pd_ptr->seeps_mpr[i]->t2); // T2 - at.set_entry(r, c+12, pd_ptr->seeps_mpr[i]->score); // SEEPS score + at.set_entry(r, c+12, pd_ptr->seeps_mpr[i]->score); // SEEPS } diff --git a/src/libcode/vx_statistics/compute_stats.cc b/src/libcode/vx_statistics/compute_stats.cc index 4dc90e4aeb..094587b732 100644 --- a/src/libcode/vx_statistics/compute_stats.cc +++ b/src/libcode/vx_statistics/compute_stats.cc @@ -1427,13 +1427,21 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps_agg) } SeepsScore *seeps_mpr = nullptr; - int count, count_diagonal; - int c12, c13, c21, c23, c31, c32; - double score_sum, obs_sum, fcst_sum; + int count = 0; + int count_diagonal = 0; + int c_odfl = 0; + int c_odfh = 0; + int c_olfd = 0; + int c_olfh = 0; + int c_ohfd = 0; + int c_ohfl = 0; + double score_sum = 0.0; + double obs_sum_wgt = 0.0; + double fcst_sum_wgt = 0.0; + double obs_sum = 0.0; + double fcst_sum = 0.0; vector seeps_mprs; - score_sum = obs_sum = fcst_sum = 0.0; - count = count_diagonal = c12 = c13 = c21 = c23 = c31 = c32 = 0; for(int i=0; in_obs; i++) { if (i >= pd->seeps_mpr.size()) break; seeps_mpr = pd->seeps_mpr[i]; @@ -1444,49 +1452,64 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps_agg) fcst_sum += pd->f_na[i]; // Forecast Value score_sum += seeps_mpr->score; if (seeps_mpr->fcst_cat == 0) { - if (seeps_mpr->obs_cat == 1) c12++; - else if(seeps_mpr->obs_cat == 2) c13++; + if (seeps_mpr->obs_cat == 1) c_olfd++; + else if(seeps_mpr->obs_cat == 2) c_ohfd++; else count_diagonal++; } else if (seeps_mpr->fcst_cat == 1) { - if (seeps_mpr->obs_cat == 0) c21++; - else if(seeps_mpr->obs_cat == 2) c23++; + if (seeps_mpr->obs_cat == 0) c_odfl++; + else if(seeps_mpr->obs_cat == 2) c_ohfl++; else count_diagonal++; } else if (seeps_mpr->fcst_cat == 2) { - if (seeps_mpr->obs_cat == 0) c31++; - else if (seeps_mpr->obs_cat == 1) c32++; + if (seeps_mpr->obs_cat == 0) c_odfh++; + else if (seeps_mpr->obs_cat == 1) c_olfh++; else count_diagonal++; } seeps_mprs.push_back(seeps_mpr); } if (count > 0) { - vector density_vector; - double pvf[SEEPS_MATRIX_SIZE]; - double weighted_score, weight_sum, weight[count]; + mlog << Debug(9) << method_name + << "Categories c_odfl, c_odfh, c_olfd, c_olfh, c_ohfd, c_ohfl => " + << c_odfl << " " << c_odfh << " " << c_olfd << " " + << c_olfh << " " << c_ohfd << " " << c_ohfl << "\n"; + + // Unweighted means seeps_agg->n_obs = count; seeps_agg->mean_fcst = fcst_sum / count; seeps_agg->mean_obs = obs_sum / count; seeps_agg->score = score_sum / count; - weighted_score = 0.; - for (int i=0; i " + << seeps_agg->mean_fcst << " " + << seeps_agg->mean_obs << " " + << seeps_agg->score << "\n"; + + double score_sum_wgt = 0.0; + vector pvf(SEEPS_MATRIX_SIZE, 0.0); + vector svf(SEEPS_MATRIX_SIZE, 0.0); + vector density_vector; compute_seeps_density_vector(pd, seeps_agg, density_vector); int density_cnt = density_vector.size(); if(density_cnt > count) density_cnt = count; //IDL: w = 1/d - weight_sum = 0.; - for (int i=0; i weight(count, 0.0); for (int i=0; i " + << i << " " << density_vector[i] << " " + << weight[i] << " " << weight_sum << "\n"; } } - if (!is_eq(weight_sum, 0)) { + if (!is_eq(weight_sum, 0.0)) { //IDL: w = w/sum(w) for (int i=0; iscore * weight[i]; + mlog << Debug(9) << method_name + << "i, seeps_mpr, weight(i), s_idx => " + << i << " " << seeps_mpr->score + << " " << weight[i] << " " << seeps_mpr->s_idx << "\n"; + score_sum_wgt += seeps_mpr->score * weight[i]; + obs_sum_wgt += pd->o_na[i] * weight[i]; + fcst_sum_wgt += pd->f_na[i] * weight[i]; + mlog << Debug(9) << method_name + << "score_sum_wgt (seeps_mpr*weight) => " + << score_sum_wgt << "\n"; //IDL: svf(cat{i)) = svf(cat{i)) + c(4+cat(i) * w{i) //IDL: pvf(cat{i)) = pvf(cat{i)) + w{i) pvf[seeps_mpr->s_idx] += weight[i]; + svf[seeps_mpr->s_idx] += seeps_mpr->score * weight[i]; } else { mlog << Debug(1) << method_name @@ -1510,45 +1543,49 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps_agg) } density_vector.clear(); - seeps_mprs.clear(); - // The weight for s12 to s32 should come from climo file, but not available yet + // The weight for odfl to ohfl come from climo file seeps_agg->pv1 = pvf[0] + pvf[3] + pvf[6]; // sum by column for obs seeps_agg->pv2 = pvf[1] + pvf[4] + pvf[7]; // sum by column for obs seeps_agg->pv3 = pvf[2] + pvf[5] + pvf[8]; // sum by column for obs seeps_agg->pf1 = pvf[0] + pvf[1] + pvf[2]; // sum by row for forecast seeps_agg->pf2 = pvf[3] + pvf[4] + pvf[5]; // sum by row for forecast seeps_agg->pf3 = pvf[6] + pvf[7] + pvf[8]; // sum by row for forecast - seeps_agg->s12 = c12 * seeps_agg->pf1 * seeps_agg->pv2; - seeps_agg->s13 = c13 * seeps_agg->pf1 * seeps_agg->pv3; - seeps_agg->s21 = c21 * seeps_agg->pf2 * seeps_agg->pv1; - seeps_agg->s23 = c23 * seeps_agg->pf2 * seeps_agg->pv3; - seeps_agg->s31 = c31 * seeps_agg->pf3 * seeps_agg->pv1; - seeps_agg->s32 = c32 * seeps_agg->pf3 * seeps_agg->pv2; - seeps_agg->weighted_score = weighted_score; + seeps_agg->s_odfl = (is_eq(svf[3], 0.0) ? 0.0 : svf[3]); + seeps_agg->s_odfh = (is_eq(svf[6], 0.0) ? 0.0 : svf[6]); + seeps_agg->s_olfd = (is_eq(svf[1], 0.0) ? 0.0 : svf[1]); + seeps_agg->s_olfh = (is_eq(svf[7], 0.0) ? 0.0 : svf[7]); + seeps_agg->s_ohfd = (is_eq(svf[2], 0.0) ? 0.0 : svf[2]); + seeps_agg->s_ohfl = (is_eq(svf[5], 0.0) ? 0.0 : svf[5]); + seeps_agg->mean_fcst_wgt = fcst_sum_wgt; + seeps_agg->mean_obs_wgt = obs_sum_wgt; + seeps_agg->score_wgt = score_sum_wgt; mlog << Debug(7) << method_name - << "SEEPS score=" << seeps_agg->score << " weighted_score=" << weighted_score + << "SEEPS score=" << seeps_agg->score << " score_wgt=" << seeps_agg->score_wgt << " pv1=" << seeps_agg->pv1 << " pv2=" << seeps_agg->pv2 << " pv3=" << seeps_agg->pv3 - << " pf1=" << seeps_agg->pf1 << " pf2=" << seeps_agg->pf2 << " pf3=" << seeps_agg->pf3 << "\n"; + << " pf1=" << seeps_agg->pf1 << " pf2=" << seeps_agg->pf2 << " pf3=" << seeps_agg->pf3 + << "\n"; } else { mlog << Debug(5) << method_name << "no SEEPS_MPR available\n"; } - seeps_agg->c12 = c12; - seeps_agg->c13 = c13; - seeps_agg->c21 = c21; - seeps_agg->c23 = c23; - seeps_agg->c31 = c31; - seeps_agg->c32 = c32; - - if (count != (c12+c13+c21+c23+c31+c32+count_diagonal)){ + + seeps_agg->c_odfl = c_odfl; + seeps_agg->c_odfh = c_odfh; + seeps_agg->c_olfd = c_olfd; + seeps_agg->c_olfh = c_olfh; + seeps_agg->c_ohfd = c_ohfd; + seeps_agg->c_ohfl = c_ohfl; + + if (count != (c_odfl+c_odfh+c_olfd+c_olfh+c_ohfd+c_ohfl+count_diagonal)){ mlog << Debug(6) << method_name - << "INFO check count: all=" << count << " s12=" << c12<< " s13=" << c13 - << " s21=" << c21 << " s23=" << c23 - << " s31=" << c31 << " s32=" << c32 << "\n"; + << "INFO check count: all=" << count + << " s_odfl=" << c_odfl << " s_odfh=" << c_odfh + << " s_olfd=" << c_olfd << " s_olfh=" << c_olfh + << " s_ohfd=" << c_ohfd << " s_ohfl=" << c_ohfl << "\n"; } return; @@ -1561,31 +1598,44 @@ void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &ob DataPlane &seeps_dp_ocat, SeepsAggScore *seeps_agg, int month, int hour, const SingleThresh &seeps_p1_thresh, const ConcatString &seeps_climo_name) { - int fcst_cat, obs_cat; - int seeps_count, count_diagonal, nan_count, bad_count; int nx = fcst_dp.nx(); int ny = fcst_dp.ny(); int dp_size = (nx * ny); - int pvf_cnt[SEEPS_MATRIX_SIZE]; - double pvf[SEEPS_MATRIX_SIZE]; - int c12, c13, c21, c23, c31, c32; - double obs_sum, fcst_sum; - double seeps_score, seeps_score_sum, seeps_score_partial_sum; static const char *method_name = "compute_aggregated_seeps_grid() -> "; + int fcst_cat = bad_data_int; + int obs_cat = bad_data_int; + int seeps_count = 0; + int count_diagonal = 0; + int nan_count = 0; + int bad_count = 0; + int c_odfl = 0; + int c_odfh = 0; + int c_olfd = 0; + int c_olfh = 0; + int c_ohfd = 0; + int c_ohfl = 0; + double seeps_score = 0.0; + double seeps_score_sum = 0.0; + double seeps_score_partial_sum = 0.0; + double obs_sum = 0.0; + double fcst_sum = 0.0; + seeps_dp.set_size(nx, ny); seeps_dp_fcat.set_size(nx, ny); seeps_dp_ocat.set_size(nx, ny); - obs_sum = fcst_sum = seeps_score_sum = 0.; - seeps_count = count_diagonal = nan_count = bad_count = 0; - c12 = c13 = c21 = c23 = c31 = c32 = 0; + seeps_agg->clear(); + mlog << Debug(9) << method_name + << "month is " << month << "\n"; + SeepsClimoGrid *seeps_climo = get_seeps_climo_grid(month, seeps_climo_name); seeps_climo->set_p1_thresh(seeps_p1_thresh); - for (int i=0; i pvf_cnt(SEEPS_MATRIX_SIZE, 0); + vector pvf(SEEPS_MATRIX_SIZE, 0.0); + vector svf(SEEPS_MATRIX_SIZE, 0.0); + for (int ix=0; ixget_record(ix, iy, fcst_value, obs_value); if (seeps_mpr != nullptr) { fcst_cat = seeps_mpr->fcst_cat; obs_cat = seeps_mpr->obs_cat; if (fcst_cat == 0) { - if (obs_cat == 1) c12++; - else if(obs_cat == 2) c13++; + if (obs_cat == 1) c_olfd++; + else if(obs_cat == 2) c_odfh++; else count_diagonal++; } else if (fcst_cat == 1) { - if (obs_cat == 0) c21++; - else if(obs_cat == 2) c23++; + if (obs_cat == 0) c_odfl++; + else if(obs_cat == 2) c_ohfd++; else count_diagonal++; } else if (fcst_cat == 2) { - if (obs_cat == 0) c31++; - else if (obs_cat == 1) c32++; + if (obs_cat == 0) c_odfh++; + else if (obs_cat == 1) c_olfh++; else count_diagonal++; } seeps_score = seeps_mpr->score; + mlog << Debug(9) << method_name + << "ix, iy, obs_cat, fcst_cat, seeps_score:" + << ix << " " << iy << " " << obs_cat << " " << fcst_cat + << " " << seeps_score << "\n"; + if (std::isnan(seeps_score)) { nan_count++; seeps_score = bad_data_double; @@ -1629,6 +1689,10 @@ void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &ob //IDL: pvf(cat{i)) = pvf(cat{i)) + w{i) //pvf[seeps_mpr->s_idx] += weight; pvf_cnt[seeps_mpr->s_idx] += 1; + mlog << Debug(9) << method_name + << "obs_sum, fcst_sum, seeps_score_partial_sum, category: " + << obs_sum << " " << fcst_sum << " " + << seeps_score_partial_sum << " " << seeps_mpr->s_idx << "\n"; } if(seeps_mpr) { delete seeps_mpr; seeps_mpr = nullptr; } @@ -1640,43 +1704,52 @@ void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &ob } seeps_score_sum += seeps_score_partial_sum; } + mlog << Debug(9) << method_name + << "dp_size, nan_count, bad_count: " + << dp_size << " " << nan_count << " " << bad_count << "\n"; int cell_count = dp_size - nan_count - bad_count; if (cell_count > 0) { - seeps_agg->weighted_score = seeps_score_sum/cell_count; + seeps_agg->score_wgt = seeps_score_sum/cell_count; for (int i=0; in_obs = seeps_count; - seeps_agg->c12 = c12; - seeps_agg->c13 = c13; - seeps_agg->c21 = c21; - seeps_agg->c23 = c23; - seeps_agg->c31 = c31; - seeps_agg->c32 = c32; - if (seeps_count > 0) { - seeps_agg->mean_fcst = fcst_sum / seeps_count; - seeps_agg->mean_obs = obs_sum / seeps_count; + seeps_agg->c_odfl = c_odfl; + seeps_agg->c_odfh = c_odfh; + seeps_agg->c_olfd = c_olfd; + seeps_agg->c_olfh = c_olfh; + seeps_agg->c_ohfd = c_ohfd; + seeps_agg->c_ohfl = c_ohfl; + if (seeps_count > 0) { seeps_agg->pv1 = pvf[0] + pvf[3] + pvf[6]; // sum by column for obs seeps_agg->pv2 = pvf[1] + pvf[4] + pvf[7]; // sum by column for obs seeps_agg->pv3 = pvf[2] + pvf[5] + pvf[8]; // sum by column for obs seeps_agg->pf1 = pvf[0] + pvf[1] + pvf[2]; // sum by row for forecast seeps_agg->pf2 = pvf[3] + pvf[4] + pvf[5]; // sum by row for forecast seeps_agg->pf3 = pvf[6] + pvf[7] + pvf[8]; // sum by row for forecast - seeps_agg->s12 = c12 * seeps_agg->pf1 * seeps_agg->pv2; - seeps_agg->s13 = c13 * seeps_agg->pf1 * seeps_agg->pv3; - seeps_agg->s21 = c21 * seeps_agg->pf2 * seeps_agg->pv1; - seeps_agg->s23 = c23 * seeps_agg->pf2 * seeps_agg->pv3; - seeps_agg->s31 = c31 * seeps_agg->pf3 * seeps_agg->pv1; - seeps_agg->s32 = c32 * seeps_agg->pf3 * seeps_agg->pv2; - seeps_agg->score = seeps_score_sum / seeps_count; + + seeps_agg->s_odfl = (is_eq(svf[3], 0.0) ? 0.0 : svf[3]); + seeps_agg->s_odfh = (is_eq(svf[6], 0.0) ? 0.0 : svf[6]); + seeps_agg->s_olfd = (is_eq(svf[1], 0.0) ? 0.0 : svf[1]); + seeps_agg->s_olfh = (is_eq(svf[7], 0.0) ? 0.0 : svf[7]); + seeps_agg->s_ohfd = (is_eq(svf[2], 0.0) ? 0.0 : svf[2]); + seeps_agg->s_ohfl = (is_eq(svf[5], 0.0) ? 0.0 : svf[5]); + + seeps_agg->mean_fcst = fcst_sum / seeps_count; + seeps_agg->mean_obs = obs_sum / seeps_count; + seeps_agg->score = seeps_score_sum / seeps_count; } + mlog << Debug(6) << method_name - << "SEEPS score=" << seeps_agg->score << " weighted_score=" << seeps_agg->weighted_score + << "SEEPS score=" << seeps_agg->score + << " score_wgt=" << seeps_agg->score_wgt << " pv1=" << seeps_agg->pv1 << " pv2=" << seeps_agg->pv2 << " pv3=" << seeps_agg->pv3 - << " pf1=" << seeps_agg->pf1 << " pf2=" << seeps_agg->pf2 << " pf3=" << seeps_agg->pf3 << "\n"; + << " pf1=" << seeps_agg->pf1 << " pf2=" << seeps_agg->pf2 << " pf3=" << seeps_agg->pf3 + << "\n"; + if(mlog.verbosity_level() >= detailed_debug_level) { char buffer[100]; ConcatString log_message; @@ -1688,8 +1761,9 @@ void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &ob snprintf ( buffer, 100, " nan: %d ", nan_count); log_message.add(buffer); } - mlog << Debug(7) << method_name << "pvf = " << log_message - << " weight=" << (1. / cell_count) << " (1/" << cell_count << ")" << "\n"; + mlog << Debug(7) << method_name + << "pvf = " << log_message << " weight=" << (1. / cell_count) + << " (1/" << cell_count << ")" << "\n"; } } @@ -1725,7 +1799,8 @@ void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &ob // ; PV-WAVE prints: 2.00000 4.00000 //////////////////////////////////////////////////////////////////////// -void compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps, vector &density_vector) { +void compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps, + vector &density_vector) { int seeps_idx; SeepsScore *seeps_mpr; int seeps_cnt = seeps->n_obs; @@ -1751,7 +1826,14 @@ void compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps, for(int i=0; in_obs; i++) { if (i >= pd->seeps_mpr.size()) break; seeps_mpr = pd->seeps_mpr[i]; + mlog << Debug(9) << method_name + << "seeps_idx, seeps_mpr => " + << seeps_idx << " " << seeps_mpr << "\n"; + if (!seeps_mpr || is_eq(seeps_mpr->score, bad_data_double)) continue; + mlog << Debug(4) << method_name + << "lat, long => " << pd->lat_na[i] << " " + << fmod((pd->lon_na[i] + 360.), 360.) << "\n"; rlat[seeps_idx] = pd->lat_na[i] * rad_per_deg; // radian of lat rlon[seeps_idx] = fmod((pd->lon_na[i] + 360.), 360.) * rad_per_deg; // radian of long @@ -1759,15 +1841,22 @@ void compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps, clon[seeps_idx] = cos(rlon[seeps_idx]); slat[seeps_idx] = sin(rlat[seeps_idx]); slon[seeps_idx] = sin(rlon[seeps_idx]); + mlog << Debug(9) << method_name + << "clat, clon, slat, slon => " + << clat[seeps_idx] << " " << clon[seeps_idx] + << " " << slat[seeps_idx] << " " << slon[seeps_idx] << "\n"; seeps_idx++; } - // prooducs n by n matrix by multipling transposed vector - int v_count; - double density, mask1, mask2, mask3, temp; + // produces n by n matrix by multiplying transposed vector + // Initialize - v_count = 0; + int v_count = 0; + mlog << Debug(9) << method_name + << "seeps_idx, seeps_cnt => " + << seeps_idx << " " << seeps_cnt << "\n"; + if (seeps_idx < seeps_cnt) seeps_cnt = seeps_idx; density_vector.reserve(seeps_cnt); for(int j=0; j " + << j << " " << i << " " << clat_m[i][j] << " " << slat_m[i][j] << " " + << clon_m[i][j] << " " << slon_m[i][j] << "\n"; + //IDL: r=(clat#transpose(clat))*(slon#transpose(slon)) + (clon#transpose(clon))*(slat#transpose(slat)) - density = clat_m[i][j] * slon_m[i][j] + clon_m[i][j] * slat_m[i][j]; + double density = clat_m[i][j] * slon_m[i][j] + clon_m[i][j] * slat_m[i][j]; + double density2 = (clat_m[i][j] * (slon_m[i][j] + clon_m[i][j])) + slat_m[i][j]; + mlog << Debug(9) << method_name + << "density, density2 => " << density << " " << density2 << "\n"; + //IDL: r * ((r lt 1.) and (r gt -1.)) + (r ge 1.) - (r le -1.) - mask1 = (density < 1.0 && density > -1.0) ? 1. : 0.; - mask2 = (density >= 1.0 ) ? 1. : 0.; - mask3 = (density <= -1.0) ? 1. : 0.; + double mask1 = (density < 1.0 && density > -1.0) ? 1. : 0.; + double mask2 = (density >= 1.0 ) ? 1. : 0.; + double mask3 = (density <= -1.0) ? 1. : 0.; + double mask5 = (density2 < 1.0 && density > -1.0) ? 1. : 0.; + double mask6 = (density2 >= 1.0 ) ? 1. : 0.; + double mask7 = (density2 <= -1.0) ? 1. : 0.; density = density * mask1 + mask2 - mask3; + density2 = density2 * mask5 + mask6 - mask7; + mlog << Debug(9) << method_name + << "density, density2 => " << density << " " << density2 << "\n"; + //IDL: r = acos(r) density = acos(density); //IDL: if r0 gt 0.0 then r = exp(-(r/r0)^2) * (r le 4. * r0) else r = (r*0.)+1. - if (density_radius_rad <= 0.) density = 1.0; - else { - mask3 = (density <= 4.0) ? 1. : 0.; - temp = density / density_radius_rad; - density = exp(-(temp * temp)) * mask3 * density_radius_rad; + if (density_radius_rad > 0.) { + if (density < 4.0 * density_radius_rad) { + mask3 = (density <= 4.0) ? 1. : 0.; + double temp = density / density_radius_rad; + density = exp(-(temp * temp)) * mask3 * density_radius_rad; + } + else { + density = 0.; + } + if (density2 < 4.0 * density_radius_rad) { + density2 = exp(-(pow(density2 / density_radius_rad,2))); + } + else { + density2 = 0.; + } + mlog << Debug(4) << method_name + << "final density, density2 => " + << density << " " << density2 << "\n"; + } else { + density = 1.; } - density_vector[j] += density; + mlog << Debug(4) << method_name + << "For Info - Feeding density2 (not density) back as vector " + << "as density all zeros in final output" << "\n"; + density_vector[j] += density2; } if (!is_eq(density_vector[j], 0.)) v_count++; } @@ -1803,8 +1926,8 @@ void compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps, << "no non-zero values at density_vector\n"; } if (seeps_cnt > 0) { - mlog << Debug(10) << method_name - << " non zero count=" << v_count + mlog << Debug(9) << method_name + << " non zero count=" << v_count << " seeps_cnt=" << seeps_cnt << " density_vector[0]=" << density_vector[0] << " density_vector[" << (seeps_cnt-1) << "]=" << density_vector[seeps_cnt-1] << "\n"; } diff --git a/src/libcode/vx_statistics/pair_data_point.cc b/src/libcode/vx_statistics/pair_data_point.cc index 0e7f6203e7..a22194a297 100644 --- a/src/libcode/vx_statistics/pair_data_point.cc +++ b/src/libcode/vx_statistics/pair_data_point.cc @@ -223,8 +223,10 @@ void PairDataPoint::set_seeps_score(SeepsScore *seeps, int index) { } } } - else mlog << Warning << "\nPairDataPoint::set_seeps_score(" - << index << ") is out of range " << seeps_count << "\n\n"; + else { + mlog << Warning << "\nPairDataPoint::set_seeps_score(" + << index << ") is out of range " << seeps_count << "\n\n"; + } } } @@ -317,14 +319,17 @@ SeepsScore *PairDataPoint::compute_seeps(const char *sid, double f, int month, day, year, hour, minute, second; int sid_no = atoi(sid); - if (sid_no && nullptr != seeps_climo) { + if (sid_no && seeps_climo) { unix_to_mdyhms(ut, month, day, year, hour, minute, second); seeps = seeps_climo->get_seeps_score(sid_no, f, o, month, hour); - if (mlog.verbosity_level() >= seeps_debug_level - && seeps && !is_eq(seeps->score, bad_data_double) - && !is_eq(seeps->score, 0) && seeps_record_count < 10) { + if (mlog.verbosity_level() >= seeps_debug_level && + seeps && + !is_bad_data(seeps->score) && + !is_eq(seeps->score, 0) && + seeps_record_count < 10) { mlog << Debug(seeps_debug_level) - << "PairDataPoint::compute_seeps() score = " << seeps->score << "\n"; + << "PairDataPoint::compute_seeps() score = " + << seeps->score << "\n"; seeps_record_count++; } } diff --git a/src/tools/core/stat_analysis/parse_stat_line.cc b/src/tools/core/stat_analysis/parse_stat_line.cc index 02283617fe..23cf791f55 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.cc +++ b/src/tools/core/stat_analysis/parse_stat_line.cc @@ -600,12 +600,12 @@ void parse_seeps_line(STATLine &l, SeepsAggScore &agg_score) { agg_score.n_obs = atoi(l.get_item("TOTAL")); - agg_score.s12 = atoi(l.get_item("S12")); - agg_score.s13 = atoi(l.get_item("S13")); - agg_score.s21 = atof(l.get_item("S21")); - agg_score.s23 = atof(l.get_item("S23")); - agg_score.s31 = atof(l.get_item("S31")); - agg_score.s32 = atof(l.get_item("S32")); + agg_score.s_odfl = atoi(l.get_item("ODFL")); + agg_score.s_odfh = atoi(l.get_item("ODFH")); + agg_score.s_olfd = atof(l.get_item("OLFD")); + agg_score.s_olfh = atof(l.get_item("OLFH")); + agg_score.s_ohfd = atof(l.get_item("OHFD")); + agg_score.s_ohfl = atof(l.get_item("OHFL")); agg_score.pf1 = atof(l.get_item("PF1")); agg_score.pf2 = atof(l.get_item("PF2")); @@ -615,11 +615,12 @@ void parse_seeps_line(STATLine &l, SeepsAggScore &agg_score) { agg_score.pv2 = atof(l.get_item("PV2")); agg_score.pv3 = atof(l.get_item("PV3")); - agg_score.mean_fcst = atof(l.get_item("MEAN_FCST")); - agg_score.mean_obs = atof(l.get_item("MEAN_OBS")); - - agg_score.score = atof(l.get_item("SEEPS")); - agg_score.weighted_score = agg_score.score; + agg_score.mean_fcst = atof(l.get_item("MEAN_FCST")); + agg_score.mean_fcst_wgt = agg_score.mean_fcst; + agg_score.mean_obs = atof(l.get_item("MEAN_OBS")); + agg_score.mean_obs_wgt = agg_score.mean_fcst; + agg_score.score = atof(l.get_item("SEEPS")); + agg_score.score_wgt = agg_score.score; return; }