diff --git a/data/config/SeriesAnalysisConfig_default b/data/config/SeriesAnalysisConfig_default index 81b2a388f4..0b1c2314f8 100644 --- a/data/config/SeriesAnalysisConfig_default +++ b/data/config/SeriesAnalysisConfig_default @@ -125,6 +125,19 @@ mask = { poly = ""; } +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1 ]; + dy = [ 1 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + // // Number of grid points to be processed concurrently. Set smaller to use less // memory but increase the number of passes through the data. If set <= 0, all @@ -155,6 +168,7 @@ output_stats = { pstd = []; pjc = []; prc = []; + grad = []; } //////////////////////////////////////////////////////////////////////////////// diff --git a/data/table_files/met_header_columns_V12.1.txt b/data/table_files/met_header_columns_V12.1.txt index bbdeaa1e03..b70610e3e8 100644 --- a/data/table_files/met_header_columns_V12.1.txt +++ b/data/table_files/met_header_columns_V12.1.txt @@ -11,7 +11,7 @@ V12.1 : STAT : SEEPS_MPR: VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID V12.1 : 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.1 : 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 V12.1 : STAT : NBRCTS : 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 BASER_NCL BASER_NCU BASER_BCL BASER_BCU FMEAN FMEAN_NCL FMEAN_NCU FMEAN_BCL FMEAN_BCU ACC ACC_NCL ACC_NCU ACC_BCL ACC_BCU FBIAS FBIAS_BCL FBIAS_BCU PODY PODY_NCL PODY_NCU PODY_BCL PODY_BCU PODN PODN_NCL PODN_NCU PODN_BCL PODN_BCU POFD POFD_NCL POFD_NCU POFD_BCL POFD_BCU FAR FAR_NCL FAR_NCU FAR_BCL FAR_BCU CSI CSI_NCL CSI_NCU CSI_BCL CSI_BCU GSS GSS_BCL GSS_BCU HK HK_NCL HK_NCU HK_BCL HK_BCU HSS HSS_BCL HSS_BCU ODDS ODDS_NCL ODDS_NCU ODDS_BCL ODDS_BCU LODDS LODDS_NCL LODDS_NCU LODDS_BCL LODDS_BCU ORSS ORSS_NCL ORSS_NCU ORSS_BCL ORSS_BCU EDS EDS_NCL EDS_NCU EDS_BCL EDS_BCU SEDS SEDS_NCL SEDS_NCU SEDS_BCL SEDS_BCU EDI EDI_NCL EDI_NCU EDI_BCL EDI_BCU SEDI SEDI_NCL SEDI_NCU SEDI_BCL SEDI_BCU BAGSS BAGSS_BCL BAGSS_BCU -V12.1 : STAT : GRAD : 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 FGBAR OGBAR MGBAR EGBAR S1 S1_OG FGOG_RATIO DX DY +V12.1 : STAT : GRAD : 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 FGBAR OGBAR MGBAR EGBAR S1 S1_OG FGOG_RATIO DX DY FGMAG OGMAG MAG_RMSE LAPLACE_RMSE V12.1 : STAT : DMAP : 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 FBIAS BADDELEY HAUSDORFF MED_FO MED_OF MED_MIN MED_MAX MED_MEAN FOM_FO FOM_OF FOM_MIN FOM_MAX FOM_MEAN ZHU_FO ZHU_OF ZHU_MIN ZHU_MAX ZHU_MEAN G GBETA BETA_VALUE V12.1 : STAT : ORANK : 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 OBS PIT RANK N_ENS_VLD (N_ENS) ENS_[0-9]* OBS_QC ENS_MEAN OBS_CLIMO_MEAN SPREAD ENS_MEAN_OERR SPREAD_OERR SPREAD_PLUS_OERR OBS_CLIMO_STDEV FCST_CLIMO_MEAN FCST_CLIMO_STDEV V12.1 : STAT : PCT : 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]* OY_[0-9]* ON_[0-9]* diff --git a/data/table_files/stat_column_description.txt b/data/table_files/stat_column_description.txt index 391f8cebcf..cd5b89b2fe 100644 --- a/data/table_files/stat_column_description.txt +++ b/data/table_files/stat_column_description.txt @@ -272,3 +272,17 @@ PRC_N_THRESH "Number of probability thresholds" PRC_THRESH_I "The i-th probability threshold" PRC_PODY_I "Probability of detecting yes when forecast is between the i-th and (i+1)-th probability thresholds" PRC_POFD_I "Probability of false detection when forecast is between the i-th and (i+1)-th probability thresholds" +GRAD_TOTAL "Total number of matched pairs" +GRAD_FGBAR "Mean of absolute value of forecast gradients" +GRAD_OGBAR "Mean of absolute value of observed gradients" +GRAD_MGBAR "Mean of maximum of absolute values of forecast and observed gradients" +GRAD_EGBAR "Mean of absolute value of forecast minus observed gradients" +GRAD_S1 "S1 score" +GRAD_S1_OG "S1 score with respect to observed gradient" +GRAD_FGOG_RATIO "Ratio of forecast and observed gradients" +GRAD_DX "Gradient size in the X-direction" +GRAD_DY "Gradient size in the Y-direction" +GRAD_FGMAG "Magnitude of the forecast gradient when the X and Y-directions are interpreted as a vector" +GRAD_OGMAG "Magnitude of the observed gradient when the X and Y-directions are interpreted as a vector" +GRAD_MAG_RMSE "Root mean squared difference of the forecast gradient magnitude minus the observed gradient magnitude" +GRAD_LAPLACE_RMSE "Root mean squared difference of the sum of the forecast X and Y-gradients minus the number of the observed X and Y-gradients" diff --git a/docs/Users_Guide/appendixC.rst b/docs/Users_Guide/appendixC.rst index 7d0f944624..39a41773b7 100644 --- a/docs/Users_Guide/appendixC.rst +++ b/docs/Users_Guide/appendixC.rst @@ -739,7 +739,7 @@ These statistics require climatological values for the wind vector components, w Gradient Values --------------- -Called "TOTAL", "FGBAR", "OGBAR", "MGBAR", "EGBAR", "S1", "S1_OG", and "FGOG_RATIO" in GRAD output :numref:`table_GS_format_info_GRAD` +Called "TOTAL", "FGBAR", "OGBAR", "MGBAR", "EGBAR", "S1", "S1_OG", "FGOG_RATIO", "FGMAG", "OGMAG", "MAG_RMSE", and "LAPLACE_RMSE" in GRAD output :numref:`table_GS_format_info_GRAD` These statistics are only computed by the Grid-Stat tool and require vectors. Here :math:`\nabla` is the gradient operator, which in this applications signifies the difference between adjacent grid points in both the grid-x and grid-y directions. TOTAL is the count of grid locations used in the calculations. The remaining measures are defined below: @@ -778,6 +778,39 @@ where the weights are applied at each grid location, with values assigned accord \text{FGOG_RATIO} = \frac{\text{FGBAR}}{\text{OGBAR}} +The following statistics are computed using the magnitude of the vectors formed by the forecast and observation gradients in the grid-x and grid-y directions. + +.. math:: + \text{FGMAG} = \text{Mean}(|| \nabla f ||) = \frac{1}{n} \sum_{i=1}^n \sqrt{\nabla {f_x}_i^2 + \nabla {f_y}_i^2} + + \text{OGMAG} = \text{Mean}(|| \nabla o ||) = \frac{1}{n} \sum_{i=1}^n \sqrt{\nabla {o_x}_i^2 + \nabla {o_y}_i^2} + +.. only:: latex + + .. math:: + + \text{MAG\_RMSE} = \sqrt{ \frac{1}{n} \sum_{i=1}^n {(|| \nabla f_i || - || \nabla o_i ||)}^2 } + +.. only:: html + + .. math:: + + \text{MAG_RMSE} = \sqrt{ \frac{1}{n} \sum_{i=1}^n {(|| \nabla f_i || - || \nabla o_i ||)}^2 } + +Laplace RMSE is very similar to gradient RMSE, but instead of taking the magnitude of the gradient vector at each +point, we compute the divergence of the gradient. + +.. only:: latex + + .. math:: + + \text{LAPLACE\_RMSE} = \sqrt{ \frac{1}{n} \sum_{i=1}^n { ((\nabla {f_x}_i + \nabla {f_y}_i) - (\nabla {o_x}_i + \nabla {o_y}_i))^2 }} + +.. only:: html + + .. math:: + + \text{LAPLACE_RMSE} = \sqrt{ \frac{1}{n} \sum_{i=1}^n { ((\nabla {f_x}_i + \nabla {f_y}_i) - (\nabla {o_x}_i + \nabla {o_y}_i))^2 }} MET Verification Measures for Probabilistic Forecasts ===================================================== diff --git a/docs/Users_Guide/config_options.rst b/docs/Users_Guide/config_options.rst index bdb3f201e4..4f7db166fd 100644 --- a/docs/Users_Guide/config_options.rst +++ b/docs/Users_Guide/config_options.rst @@ -2232,6 +2232,32 @@ This dictionary may include the following entries: prob_cat_thresh = []; } + +.. _gradient: + +gradient +-------- + +:ref:`gradient ` + +The "gradient" entry is a dictionary which specifies the number and size +of gradients to be computed and applies to both Grid-Stat and +Series-Analysis. The "dx" and "dy" entries specify the size of the +gradients in grid units in the X and Y dimensions, respectively. "dx" +and "dy" are arrays of integers (positive or negative) which must have the +same length, and the GRAD output line type will be computed separately for +each entry. When computing gradients, the value at the (x, y) grid point +is replaced by the value at the (x+dx, y+dy) grid point minus the value at +(x, y). This configuration option may be set separately in each "obs.field" +entry. + +.. code-block:: none + + gradient = { + dx = [ 1 ]; + dy = [ 1 ]; + } + output_flag ----------- diff --git a/docs/Users_Guide/grid-stat.rst b/docs/Users_Guide/grid-stat.rst index d7a9a02d41..f69abaadae 100644 --- a/docs/Users_Guide/grid-stat.rst +++ b/docs/Users_Guide/grid-stat.rst @@ -98,6 +98,13 @@ The S1 score has been in historical use for verification of forecasts, particula Differences are computed in both of the horizontal grid directions and is not a true mathematical gradient. Because the S1 score focuses on differences only, any bias in the forecast will not be measured. Further, the score depends on the domain and spacing of the grid, so can only be compared on forecasts with identical grids. +As described in :ref:`Ebert-Uphoff et al., 2024 `, statistics based +on the magnitude of the forecast and observed gradients are also provided. Similiar to +the S1 score, the root-mean-squared error of the magnitude of the gradients and their +divergence quantify the similarity in the texture of the fields, with 0 being a perfect +score. These gradient-based statistics assess the difference in smoothness between the +two fields but not the accuracy of the forecast. + Distance Maps ------------- @@ -263,6 +270,7 @@ __________________________ eclv_points = 0.05; hss_ec_value = NA; rank_corr_flag = TRUE; + gradient = { dx = [ 1 ]; dy = [ 1 ]; } grid_weight_flag = NONE; tmp_dir = "/tmp"; output_prefix = ""; @@ -319,22 +327,6 @@ The available wave numbers start at 0 (the mean across each row of data) and end _____________________ -.. _gradient: - -:ref:`gradient ` - -.. code-block:: none - - gradient = { - dx = [ 1 ]; - dy = [ 1 ]; - } - - -The **gradient** entry is a dictionary which specifies the number and size of gradients to be computed. The **dx** and **dy** entries specify the size of the gradients in grid units in the X and Y dimensions, respectively. **dx** and **dy** are arrays of integers (positive or negative) which must have the same length, and the GRAD output line type will be computed separately for each entry. When computing gradients, the value at the (x, y) grid point is replaced by the value at the (x+dx, y+dy) grid point minus the value at (x, y). This configuration option may be set separately in each **obs.field** entry. - -____________________ - .. _distance_map: :ref:`distance_map ` @@ -832,11 +824,27 @@ The format of the STAT and ASCII output of the Grid-Stat tool are the same as th * - 33 - DX - Gradient size in the X-direction - - Integer + - Integer * - 34 - DY - Gradient size in the Y-direction - - Integer + - Integer + * - 35 + - FGMAG + - Magnitude of the forecast gradient when the X and Y-directions are interpreted as a vector + - Double + * - 36 + - OGMAG + - Magnitude of the observed gradient when the X and Y-directions are intrepreted as a vector + - Double + * - 37 + - MAG_RMSE + - Root mean squared difference of the forecast gradient magnitude minus the observed gradient magnitude + - Double + * - 38 + - LAPLACE_RMSE + - Root mean squared difference of the sum of the forecast X and Y-gradients minus the sum of the observed X and Y-gradients + - Double .. _table_GS_format_info_DMAP: diff --git a/docs/Users_Guide/refs.rst b/docs/Users_Guide/refs.rst index 00e682abad..3ebe0fae7b 100644 --- a/docs/Users_Guide/refs.rst +++ b/docs/Users_Guide/refs.rst @@ -130,6 +130,12 @@ References | a review and proposed framework. *Meteorological Applications*, 15, 51-64. | +.. _Ebert-Uphoff-2024: + +| Ebert-Uphoff, I.,, 2024: An Investigation of Metrics to Evaluate the Sharpness +| in AI-Generated Meteorological Imagery. *Draft version - Jan 26, 2024* +| + .. _Eckel-2012: | Eckel, F. A., M.S. Allen, M. C. Sittel, 2012: Estimation of Ambiguity in @@ -278,7 +284,7 @@ References .. _North-2022: -| North, R.C., M.P. Mittermaier, S.F. Milton, 2022. *Using SEEPS with a* +| North, R.C., M.P. Mittermaier, S.F. Milton, 2022. *Using SEEPS with a* | TRMM-derived Climatology to Assess Global NWP Precipitation Forecast Skill. | *Monthly Weather Review*, 150, 135-155. | https://doi.org/10.1175/MWR-D-20-0347.1 diff --git a/docs/Users_Guide/series-analysis.rst b/docs/Users_Guide/series-analysis.rst index ed9f5578ab..0be4477c90 100644 --- a/docs/Users_Guide/series-analysis.rst +++ b/docs/Users_Guide/series-analysis.rst @@ -116,6 +116,7 @@ ____________________ boot = { interval = PCTILE; rep_prop = 1.0; n_rep = 1000; rng = "mt19937"; seed = ""; } mask = { grid = [ "FULL" ]; poly = []; } + gradient = { dx = [ 1 ]; dy = [ 1 ]; } hss_ec_value = NA; rank_corr_flag = TRUE; tmp_dir = "/tmp"; @@ -139,7 +140,6 @@ ____________________ Ratio of valid matched pairs for the series of values at each grid point required to compute statistics. Set to a lower proportion to allow some missing values. Setting it to 1.0 requires that every data point be valid over the series to compute statistics. - ____________________ .. code-block:: none @@ -157,6 +157,7 @@ ____________________ pstd = []; pjc = []; prc = []; + grad = []; } The output_stats array controls the type of output that the Series-Analysis tool generates. Each flag corresponds to an output line type in the STAT file and is used to specify the comma-separated list of statistics to be computed. Use the column names from the tables listed below to specify the statistics. The output flags correspond to the following types of output line types: @@ -185,4 +186,6 @@ The output_stats array controls the type of output that the Series-Analysis tool 12. PRC for Receiver Operating Characteristic for Probabilistic forecasts (See :numref:`table_PS_format_info_PRC`) -.. note:: When the -input option is used, all partial sum and contingency table count columns are required to aggregate statistics across multiple runs. To facilitate this, the output_stats entries for the CTC, SL1L2, SAL1L2, and PCT line types can be set to "ALL" to indicate that all available columns for those line types should be written. +13. GRAD for Gradient Statistics (See :numref:`table_GS_format_info_GRAD`) + +.. note:: When the -input option is used, all partial sum and contingency table count columns are required to aggregate statistics across multiple runs. To facilitate this, the output_stats entries for the CTC, SL1L2, SAL1L2, PCT, and GRAD line types can be set to "ALL" to indicate that all available columns for those line types should be written. diff --git a/internal/test_unit/config/SeriesAnalysisConfig b/internal/test_unit/config/SeriesAnalysisConfig index 047b454afc..72945b1266 100644 --- a/internal/test_unit/config/SeriesAnalysisConfig +++ b/internal/test_unit/config/SeriesAnalysisConfig @@ -109,6 +109,19 @@ mask = { poly = "${MASK_POLY}"; } +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1, 3 ]; + dy = [ 1, 3 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + // // Number of grid points to be processed concurrently. Set smaller to use less // memory but increase the number of passes through the data. If set <= 0, all @@ -139,6 +152,7 @@ output_stats = { pstd = [ ${PSTD_STATS} ]; pjc = [ ${PJC_STATS} ]; prc = [ ${PRC_STATS} ]; + grad = [ ${GRAD_STATS} ]; } //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/SeriesAnalysisConfig_climo b/internal/test_unit/config/SeriesAnalysisConfig_climo index b32c9e4283..e75afb4ec5 100644 --- a/internal/test_unit/config/SeriesAnalysisConfig_climo +++ b/internal/test_unit/config/SeriesAnalysisConfig_climo @@ -113,6 +113,19 @@ mask = { poly = ""; } +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1 ]; + dy = [ 1 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + // // Number of grid points to be processed concurrently. Set smaller to use less // memory but increase the number of passes through the data. If set <= 0, all @@ -133,16 +146,17 @@ vld_thresh = 0.5; output_stats = { fho = [ "TOTAL", "F_RATE", "H_RATE", "O_RATE" ]; ctc = [ "ALL" ]; - cts = [ ]; - mctc = [ ]; - mcts = [ ]; + cts = []; + mctc = []; + mcts = []; cnt = [ "TOTAL", "RMSE", "ANOM_CORR", "RMSFA", "RMSOA" ]; sl1l2 = [ "ALL" ]; sal1l2 = [ "ALL" ]; - pct = [ ]; - pstd = [ ]; - pjc = [ ]; - prc = [ ]; + pct = []; + pstd = []; + pjc = []; + prc = []; + grad = []; } //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/SeriesAnalysisConfig_climo_prob b/internal/test_unit/config/SeriesAnalysisConfig_climo_prob index ddbee4ec66..aa4e6e7fdc 100644 --- a/internal/test_unit/config/SeriesAnalysisConfig_climo_prob +++ b/internal/test_unit/config/SeriesAnalysisConfig_climo_prob @@ -122,6 +122,19 @@ mask = { poly = ""; } +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1 ]; + dy = [ 1 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + // // Number of grid points to be processed concurrently. Set smaller to use less // memory but increase the number of passes through the data. If set <= 0, all @@ -140,18 +153,19 @@ vld_thresh = 0.5; // Statistical output types // output_stats = { - fho = [ ]; - ctc = [ ]; - cts = [ ]; - mctc = [ ]; - mcts = [ ]; - cnt = [ ]; - sl1l2 = [ ]; - sal1l2 = [ ]; + fho = []; + ctc = []; + cts = []; + mctc = []; + mcts = []; + cnt = []; + sl1l2 = []; + sal1l2 = []; pct = [ "ALL" ]; pstd = [ "TOTAL", "ROC_AUC", "BRIER", "BRIERCL", "BSS", "BSS_SMPL" ]; - pjc = [ ]; - prc = [ ]; + pjc = []; + prc = []; + grad = []; } //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/SeriesAnalysisConfig_conditional b/internal/test_unit/config/SeriesAnalysisConfig_conditional index 813145fe5a..6639bf4a75 100644 --- a/internal/test_unit/config/SeriesAnalysisConfig_conditional +++ b/internal/test_unit/config/SeriesAnalysisConfig_conditional @@ -76,6 +76,19 @@ mask = { poly = ""; } +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1 ]; + dy = [ 1 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + // // Number of grid points to be processed concurrently. Set smaller to use less // memory but increase the number of passes through the data. If set <= 0, all @@ -94,18 +107,19 @@ vld_thresh = 1.0; // Statistical output types // output_stats = { - fho = [ ]; - ctc = [ ]; - cts = [ ]; - mctc = [ ]; - mcts = [ ]; + fho = []; + ctc = []; + cts = []; + mctc = []; + mcts = []; cnt = [ "TOTAL", "ME", "RMSE" ]; sl1l2 = [ "FBAR", "OBAR" ]; - sal1l2 = [ ]; - pct = [ ]; - pstd = [ ]; - pjc = [ ]; - prc = [ ]; + sal1l2 = []; + pct = []; + pstd = []; + pjc = []; + prc = []; + grad = []; } //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/SeriesAnalysisConfig_const_climo b/internal/test_unit/config/SeriesAnalysisConfig_const_climo index 312aaf3189..ffa97b698a 100644 --- a/internal/test_unit/config/SeriesAnalysisConfig_const_climo +++ b/internal/test_unit/config/SeriesAnalysisConfig_const_climo @@ -120,6 +120,19 @@ mask = { poly = ""; } +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1 ]; + dy = [ 1 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + // // Number of grid points to be processed concurrently. Set smaller to use less // memory but increase the number of passes through the data. If set <= 0, all @@ -138,18 +151,19 @@ vld_thresh = 0.5; // Statistical output types // output_stats = { - fho = [ ]; - ctc = [ ]; - cts = [ ]; - mctc = [ ]; - mcts = [ ]; + fho = []; + ctc = []; + cts = []; + mctc = []; + mcts = []; cnt = [ "TOTAL", "RMSE", "ANOM_CORR" ]; - sl1l2 = [ ]; - sal1l2 = [ ]; - pct = [ ]; - pstd = [ ]; - pjc = [ ]; - prc = [ ]; + sl1l2 = []; + sal1l2 = []; + pct = []; + pstd = []; + pjc = []; + prc = []; + grad = []; } //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/SeriesAnalysisConfig_python b/internal/test_unit/config/SeriesAnalysisConfig_python index 200c380fb8..482d84961b 100644 --- a/internal/test_unit/config/SeriesAnalysisConfig_python +++ b/internal/test_unit/config/SeriesAnalysisConfig_python @@ -113,6 +113,19 @@ mask = { poly = ""; } +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1 ]; + dy = [ 1 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + // // Number of grid points to be processed concurrently. Set smaller to use less // memory but increase the number of passes through the data. If set <= 0, all @@ -143,6 +156,7 @@ output_stats = { pstd = []; pjc = []; prc = []; + grad = []; } //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/hdr/met_12_1.hdr b/internal/test_unit/hdr/met_12_1.hdr index f8655a4a47..b9c7f4fa82 100644 --- a/internal/test_unit/hdr/met_12_1.hdr +++ b/internal/test_unit/hdr/met_12_1.hdr @@ -11,7 +11,7 @@ SEEPS_MPR : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_L 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 NBRCTS : 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 BASER_NCL BASER_NCU BASER_BCL BASER_BCU FMEAN FMEAN_NCL FMEAN_NCU FMEAN_BCL FMEAN_BCU ACC ACC_NCL ACC_NCU ACC_BCL ACC_BCU FBIAS FBIAS_BCL FBIAS_BCU PODY PODY_NCL PODY_NCU PODY_BCL PODY_BCU PODN PODN_NCL PODN_NCU PODN_BCL PODN_BCU POFD POFD_NCL POFD_NCU POFD_BCL POFD_BCU FAR FAR_NCL FAR_NCU FAR_BCL FAR_BCU CSI CSI_NCL CSI_NCU CSI_BCL CSI_BCU GSS GSS_BCL GSS_BCU HK HK_NCL HK_NCU HK_BCL HK_BCU HSS HSS_BCL HSS_BCU ODDS ODDS_NCL ODDS_NCU ODDS_BCL ODDS_BCU LODDS LODDS_NCL LODDS_NCU LODDS_BCL LODDS_BCU ORSS ORSS_NCL ORSS_NCU ORSS_BCL ORSS_BCU EDS EDS_NCL EDS_NCU EDS_BCL EDS_BCU SEDS SEDS_NCL SEDS_NCU SEDS_BCL SEDS_BCU EDI EDI_NCL EDI_NCU EDI_BCL EDI_BCU SEDI SEDI_NCL SEDI_NCU SEDI_BCL SEDI_BCU BAGSS BAGSS_BCL BAGSS_BCU -GRAD : 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 FGBAR OGBAR MGBAR EGBAR S1 S1_OG FGOG_RATIO DX DY +GRAD : 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 FGBAR OGBAR MGBAR EGBAR S1 S1_OG FGOG_RATIO DX DY FGMAG OGMAG MAG_RMSE LAPLACE_RMSE DMAP : 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 FBIAS BADDELEY HAUSDORFF MED_FO MED_OF MED_MIN MED_MAX MED_MEAN FOM_FO FOM_OF FOM_MIN FOM_MAX FOM_MEAN ZHU_FO ZHU_OF ZHU_MIN ZHU_MAX ZHU_MEAN G GBETA BETA_VALUE ORANK : 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 OBS PIT RANK N_ENS_VLD N_ENS _VAR_ PCT : 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_ diff --git a/internal/test_unit/xml/unit_series_analysis.xml b/internal/test_unit/xml/unit_series_analysis.xml index 96cda729dd..61075a94eb 100644 --- a/internal/test_unit/xml/unit_series_analysis.xml +++ b/internal/test_unit/xml/unit_series_analysis.xml @@ -40,6 +40,7 @@ PSTD_STATS PJC_STATS PRC_STATS + GRAD_STATS "ALL" \ -fcst &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F006.grib \ @@ -81,6 +82,7 @@ PSTD_STATS PJC_STATS PRC_STATS + GRAD_STATS "ALL" \ -fcst &DATA_DIR_MODEL;/grib1/gfs_hmt/gfs_2012040900_F030.grib \ @@ -131,6 +133,7 @@ PSTD_STATS "TOTAL", "ROC_AUC", "BRIER", "BRIER_NCL", "BRIER_NCU" PJC_STATS "CALIBRATION_1", "REFINEMENT_1" PRC_STATS "PODY_1", "POFD_1" + GRAD_STATS \ -fcst &OUTPUT_DIR;/series_analysis/input_fcst_file_list \ @@ -174,6 +177,7 @@ PSTD_STATS "TOTAL", "ROC_AUC", "BRIER", "BRIER_NCL", "BRIER_NCU" PJC_STATS "CALIBRATION_1", "REFINEMENT_1" PRC_STATS "PODY_1", "POFD_1" + GRAD_STATS \ -fcst &OUTPUT_DIR;/series_analysis/aggregate_fcst_file_list \ @@ -210,6 +214,7 @@ PSTD_STATS PJC_STATS PRC_STATS + GRAD_STATS \ -fcst &DATA_DIR_MODEL;/grib1/arw-fer-gep1/arw-fer-gep1_2012040900_F012.grib \ diff --git a/src/basic/vx_math/is_bad_data.h b/src/basic/vx_math/is_bad_data.h index cb1c4afa8c..680bd5354e 100644 --- a/src/basic/vx_math/is_bad_data.h +++ b/src/basic/vx_math/is_bad_data.h @@ -101,7 +101,7 @@ inline double square(double v) { } inline double square_root(double v) { - return (is_bad_data(v) ? bad_data_double : sqrt(v)); + return (is_bad_data(v) || v < 0 ? bad_data_double : sqrt(v)); } diff --git a/src/basic/vx_util/stat_column_defs.h b/src/basic/vx_util/stat_column_defs.h index 984ef3c40b..7b532719c5 100644 --- a/src/basic/vx_util/stat_column_defs.h +++ b/src/basic/vx_util/stat_column_defs.h @@ -249,7 +249,9 @@ static const char * const grad_columns [] = { "TOTAL", "FGBAR", "OGBAR", "MGBAR", "EGBAR", "S1", "S1_OG", - "FGOG_RATIO", "DX", "DY" + "FGOG_RATIO", "DX", "DY", + "FGMAG", "OGMAG", "MAG_RMSE", + "LAPLACE_RMSE" }; static const char * const dmap_columns [] = { diff --git a/src/libcode/vx_stat_out/stat_columns.cc b/src/libcode/vx_stat_out/stat_columns.cc index c805d447fd..cd2b560c33 100644 --- a/src/libcode/vx_stat_out/stat_columns.cc +++ b/src/libcode/vx_stat_out/stat_columns.cc @@ -3976,7 +3976,9 @@ void write_grad_cols(const GRADInfo &grad_info, // TOTAL, // FGBAR, OGBAR, MGBAR, // EGBAR, S1, S1_OG, - // FGOG_RATIO, DX, DY + // FGOG_RATIO, DX, DY, + // FGMAG, OGMAG, MAG_RMSE, + // LAPLACE_RMSE // at.set_entry(r, c+0, // Total Count grad_info.total); @@ -4008,6 +4010,18 @@ void write_grad_cols(const GRADInfo &grad_info, at.set_entry(r, c+9, // DY grad_info.dy); + at.set_entry(r, c+10, // FGMAG + grad_info.fgmag); + + at.set_entry(r, c+11, // OGMAG + grad_info.ogmag); + + at.set_entry(r, c+12, // MAG_RMSE + grad_info.magnitude_rmse()); + + at.set_entry(r, c+13, // LAPLACE_RMSE + grad_info.laplace_rmse()); + return; } diff --git a/src/libcode/vx_statistics/met_stats.cc b/src/libcode/vx_statistics/met_stats.cc index c1e88d28de..af11437870 100644 --- a/src/libcode/vx_statistics/met_stats.cc +++ b/src/libcode/vx_statistics/met_stats.cc @@ -3639,6 +3639,10 @@ GRADInfo & GRADInfo::operator=(const GRADInfo &c) { GRADInfo & GRADInfo::operator+=(const GRADInfo &c) { GRADInfo g_info; + // Return if nothing to add + if(c.total == 0) return *this; + + // Gradient definition must remain constant if(dx != c.dx || dy != c.dy) { mlog << Error << "\nGRADInfo::operator+=() -> " << "the gradient DX (" << dx << " vs " << c.dx @@ -3656,6 +3660,11 @@ GRADInfo & GRADInfo::operator+=(const GRADInfo &c) { g_info.ogbar = (ogbar*total + c.ogbar*c.total) / g_info.total; g_info.mgbar = (mgbar*total + c.mgbar*c.total) / g_info.total; g_info.egbar = (egbar*total + c.egbar*c.total) / g_info.total; + + g_info.fgmag = (fgmag*total + c.fgmag*c.total) / g_info.total; + g_info.ogmag = (ogmag*total + c.ogmag*c.total) / g_info.total; + g_info.mag_mse = (mag_mse*total + c.mag_mse*c.total) / g_info.total; + g_info.lap_mse = (lap_mse*total + c.lap_mse*c.total) / g_info.total; } assign(g_info); @@ -3676,10 +3685,12 @@ void GRADInfo::init_from_scratch() { void GRADInfo::clear() { - dx = dy = 0; - fgbar = ogbar = 0.0; - mgbar = egbar = 0.0; - total = 0; + total = 0; + dx = dy = 0; + fgbar = ogbar = 0.0; + mgbar = egbar = 0.0; + fgmag = ogmag = 0.0; + mag_mse = lap_mse = 0.0; return; } @@ -3690,6 +3701,8 @@ void GRADInfo::assign(const GRADInfo &c) { clear(); + total = c.total; + // Gradient sizes dx = c.dx; dy = c.dy; @@ -3699,7 +3712,12 @@ void GRADInfo::assign(const GRADInfo &c) { ogbar = c.ogbar; mgbar = c.mgbar; egbar = c.egbar; - total = c.total; + + // Gradient vector partial sums + fgmag = c.fgmag; + ogmag = c.ogmag; + mag_mse = c.mag_mse; + lap_mse = c.lap_mse; return; } @@ -3751,6 +3769,18 @@ double GRADInfo::fgog_ratio() const { //////////////////////////////////////////////////////////////////////// +double GRADInfo::magnitude_rmse() const { + return square_root(mag_mse); +} + +//////////////////////////////////////////////////////////////////////// + +double GRADInfo::laplace_rmse() const { + return square_root(lap_mse); +} + +//////////////////////////////////////////////////////////////////////// + void GRADInfo::set(int grad_dx, int grad_dy, const NumArray &fgx_na, const NumArray &fgy_na, const NumArray &ogx_na, const NumArray &ogy_na, @@ -3801,6 +3831,22 @@ void GRADInfo::set(int grad_dx, int grad_dy, max(fabs(fgy_na[i]), fabs(ogy_na[i]))); egbar += wgt * (fabs(fgx_na[i] - ogx_na[i]) + fabs(fgy_na[i] - ogy_na[i])); + + // Gradient vector magnitude + // Reference: + // Ebert-Uphoff, I.: An Investigation of Metrics to Evaluate the Sharpness in + // AI-Generated Meteorology Imagery + // Draft version, Jan 26, 2024 + + double fmag = square_root(fgx_na[i] * fgx_na[i] + fgy_na[i] * fgy_na[i]); + double omag = square_root(ogx_na[i] * ogx_na[i] + ogy_na[i] * ogy_na[i]); + + // Gradient vector sums + fgmag += wgt * fmag; + ogmag += wgt * omag; + mag_mse += wgt * (fmag - omag)*(fmag - omag); + double diff = (fgx_na[i] + fgy_na[i]) - (ogx_na[i] + ogy_na[i]); + lap_mse += wgt * (diff * diff); total++; } @@ -3813,6 +3859,72 @@ void GRADInfo::set(int grad_dx, int grad_dy, return; } +//////////////////////////////////////////////////////////////////////// + +void GRADInfo::set_stat(const string &stat_name, double v) { + + // Store the statistic by name + if(stat_name == "TOTAL" ) total = nint(v); + else if(stat_name == "FGBAR" ) fgbar = v; + else if(stat_name == "OGBAR" ) ogbar = v; + else if(stat_name == "MGBAR" ) mgbar = v; + else if(stat_name == "EGBAR" ) egbar = v; + else if(stat_name == "DX" ) dx = nint(v); + else if(stat_name == "DY" ) dy = nint(v); + else if(stat_name == "FGMAG" ) fgmag = v; + else if(stat_name == "OGMAG" ) ogmag = v; + else if(stat_name == "MAG_RMSE" ) mag_mse = v*v; + else if(stat_name == "LAPLACE_RMSE") lap_mse = v*v; + else if(stat_name == "S1" || + stat_name == "S1_OG" || + stat_name == "FGOG_RATIO") { + // Ignore derived quantities + } + else { + mlog << Error << "\nGRADInfo::set_stat() -> " + << "unknown gradient statistic name \"" << stat_name + << "\"!\n\n"; + exit(1); + } + + return; +} + +//////////////////////////////////////////////////////////////////////// + +double GRADInfo::get_stat(const string &stat_name) const { + double v = bad_data_double; + + // Find the statistic by name + if(stat_name == "TOTAL" ) v = (double) total; + else if(stat_name == "FGBAR" ) v = fgbar; + else if(stat_name == "OGBAR" ) v = ogbar; + else if(stat_name == "MGBAR" ) v = mgbar; + else if(stat_name == "EGBAR" ) v = egbar; + else if(stat_name == "S1" ) v = s1(); + else if(stat_name == "S1_OG" ) v = s1_og(); + else if(stat_name == "FGOG_RATIO" ) v = fgog_ratio(); + else if(stat_name == "DX" ) v = (double) dx; + else if(stat_name == "DY" ) v = (double) dy; + else if(stat_name == "FGMAG" ) v = fgmag; + else if(stat_name == "OGMAG" ) v = ogmag; + else if(stat_name == "MAG_RMSE" ) v = magnitude_rmse(); + else if(stat_name == "LAPLACE_RMSE") v = laplace_rmse(); + else { + mlog << Error << "\nGRADInfo::get_stat() -> " + << "unknown gradient statistic name \"" << stat_name + << "\"!\n\n"; + exit(1); + } + + // Return bad data for 0 pairs + if(total == 0 && stat_name != "TOTAL") { + v = bad_data_double; + } + + return v; +} + //////////////////////////////////////////////////////////////////////// // // Code for class DMAPInfo diff --git a/src/libcode/vx_statistics/met_stats.h b/src/libcode/vx_statistics/met_stats.h index 1b1d69ec7d..2cb3b343f4 100644 --- a/src/libcode/vx_statistics/met_stats.h +++ b/src/libcode/vx_statistics/met_stats.h @@ -613,22 +613,39 @@ class GRADInfo { int dx; int dy; - // Gradient Partial Sums - int total; - double fgbar, ogbar, mgbar, egbar; + // Number of pairs + int total; + + // Gradient Partial sums + double fgbar; + double ogbar; + double mgbar; + double egbar; + + // Gradient Vector Partial Sums + double fgmag; + double ogmag; + double mag_mse; + double lap_mse; // Gradient Statistics double s1() const; // s1 = 100 * egbar / mgbar double s1_og() const; // s1_og = 100 * egbar / ogbar double fgog_ratio() const; // fgog_ratio = fgbar / ogbar - // Compute sums + double magnitude_rmse() const; + double laplace_rmse() const; + + void clear(); + + // Compute gradient sums void set(int grad_dx, int grad_dy, const NumArray &fgx_na, const NumArray &fgy_na, const NumArray &ogx_na, const NumArray &ogy_na, const NumArray &wgt_na); - void clear(); + void set_stat(const std::string &, double); + double get_stat(const std::string &) const; }; //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/series_analysis/series_analysis.cc b/src/tools/core/series_analysis/series_analysis.cc index e5c3a62fe8..32b7295e2b 100644 --- a/src/tools/core/series_analysis/series_analysis.cc +++ b/src/tools/core/series_analysis/series_analysis.cc @@ -37,6 +37,7 @@ // 016 01/29/24 Halley Gotway MET #2801 Configure time difference warnings. // 017 07/05/24 Halley Gotway MET #2924 Support forecast climatology. // 018 07/26/24 Halley Gotway MET #1371 Aggregate previous output. +// 019 12/09/24 Halley Gotway MET #3030 Add GRAD output. // //////////////////////////////////////////////////////////////////////// @@ -95,6 +96,9 @@ static void do_continuous (int, const PairDataPoint *); static void do_partialsums (int, const PairDataPoint *); static void do_probabilistic (int, const PairDataPoint *); static void do_climo_brier (int, double, int, PCTInfo &); +static void do_gradient (int, int, int, + const PairDataPoint *, + const PairDataPoint *); static int read_aggr_total (int); static void read_aggr_ctc (int, const CTSInfo &, CTSInfo &); @@ -102,6 +106,7 @@ static void read_aggr_mctc (int, const MCTSInfo &, MCTSInfo &); static void read_aggr_sl1l2 (int, const SL1L2Info &, SL1L2Info &); static void read_aggr_sal1l2 (int, const SL1L2Info &, SL1L2Info &); static void read_aggr_pct (int, const PCTInfo &, PCTInfo &); +static void read_aggr_grad (int, const GRADInfo &, GRADInfo &); static void store_stat_categorical(int, STATLineType, const ConcatString &, @@ -118,12 +123,16 @@ static void store_stat_continuous(int, static void store_stat_probabilistic(int, STATLineType, const ConcatString &, const PCTInfo &); +static void store_stat_gradient(int, + STATLineType, const ConcatString &, + const GRADInfo &); static void store_stat_all_ctc (int, const CTSInfo &); static void store_stat_all_mctc (int, const MCTSInfo &); static void store_stat_all_sl1l2 (int, const SL1L2Info &); static void store_stat_all_sal1l2(int, const SL1L2Info &); static void store_stat_all_pct (int, const PCTInfo &); +static void store_stat_all_grad (int, const GRADInfo &); static ConcatString build_nc_var_name_categorical( STATLineType, const ConcatString &, @@ -140,6 +149,9 @@ static ConcatString build_nc_var_name_continuous( static ConcatString build_nc_var_name_probabilistic( STATLineType, const ConcatString &, const PCTInfo &, double); +static ConcatString build_nc_var_name_gradient( + STATLineType, const ConcatString &, + const GRADInfo &); static void setup_nc_file(const VarInfo *, const VarInfo *); static void add_stat_data(const ConcatString &, const ConcatString &, @@ -149,6 +161,7 @@ static void write_stat_data(); static void set_range(const unixtime &, unixtime &, unixtime &); static void set_range(const int &, int &, int &); +static void set_pair_dims(vector &, int, int); static void clean_up(); @@ -182,13 +195,13 @@ int met_main(int argc, char *argv[]) { //////////////////////////////////////////////////////////////////////// -const string get_tool_name() { +string get_tool_name() { return "series_analysis"; } //////////////////////////////////////////////////////////////////////// -void process_command_line(int argc, char **argv) { +static void process_command_line(int argc, char **argv) { int i; CommandLine cline; ConcatString default_config_file; @@ -238,13 +251,13 @@ void process_command_line(int argc, char **argv) { << R"("-obs" or "-both" option.)" << "\n\n"; usage(); } - if(config_file.length() == 0) { + if(config_file.empty()) { mlog << Error << "\nprocess_command_line() -> " << "the configuration file must be set using the " << R"("-config" option.)" << "\n\n"; usage(); } - if(out_file.length() == 0) { + if(out_file.empty()) { mlog << Error << "\nprocess_command_line() -> " << "the output NetCDF file must be set using the " << R"("-out" option.)" << "\n\n"; @@ -388,7 +401,7 @@ void process_command_line(int argc, char **argv) { //////////////////////////////////////////////////////////////////////// -void process_grid(const Grid &fcst_grid, const Grid &obs_grid) { +static void process_grid(const Grid &fcst_grid, const Grid &obs_grid) { // Determine the verification grid grid = parse_vx_grid(conf_info.fcst_info[0]->regrid(), @@ -426,10 +439,10 @@ void process_grid(const Grid &fcst_grid, const Grid &obs_grid) { //////////////////////////////////////////////////////////////////////// -Met2dDataFile *get_mtddf(const StringArray &file_list, - const GrdFileType type) { +static Met2dDataFile *get_mtddf(const StringArray &file_list, + const GrdFileType type) { int i; - Met2dDataFile *mtddf = (Met2dDataFile *) nullptr; + Met2dDataFile *mtddf = nullptr; // Find the first file that actually exists for(i=0; igrid(); // Close the data file - delete mtddf; mtddf = (Met2dDataFile *) nullptr; + delete mtddf; mtddf = nullptr; return found; } //////////////////////////////////////////////////////////////////////// -void open_aggr_file() { +static void open_aggr_file() { mlog << Debug(1) << "Reading aggregate data file: " << aggr_file << "\n"; @@ -791,8 +806,8 @@ void open_aggr_file() { //////////////////////////////////////////////////////////////////////// -DataPlane read_aggr_data_plane(const ConcatString &var_name, - const char *suggestion) { +static DataPlane read_aggr_data_plane(const ConcatString &var_name, + const char *suggestion) { DataPlane aggr_dp; // Setup the data request @@ -833,20 +848,26 @@ DataPlane read_aggr_data_plane(const ConcatString &var_name, //////////////////////////////////////////////////////////////////////// -void process_scores() { - int x; - int y; +static void process_scores() { int i_point = 0; - VarInfo *fcst_info = (VarInfo *) nullptr; - VarInfo *obs_info = (VarInfo *) nullptr; + VarInfo *fcst_info = nullptr; + VarInfo *obs_info = nullptr; DataPlane fcst_dp; DataPlane obs_dp; vector pd_block; - const char *method_name = "process_scores() "; + + // X and Y gradient pairs + bool do_grad = (!conf_info.fcst_info[0]->is_prob() && + conf_info.get_n_grad() > 0 && + conf_info.output_stats[STATLineType::grad].n() > 0); + vector gx_pd_block; + vector gy_pd_block; // Climatology mean and standard deviation - DataPlane fcmn_dp, fcsd_dp; - DataPlane ocmn_dp, ocsd_dp; + DataPlane fcmn_dp; + DataPlane fcsd_dp; + DataPlane ocmn_dp; + DataPlane ocsd_dp; // Open the aggregate file, if needed if(aggr_file.nonempty()) open_aggr_file(); @@ -876,22 +897,33 @@ void process_scores() { // Retrieve the data planes for the current series entry get_series_data(i_series, fcst_info, obs_info, fcst_dp, obs_dp); - // Initialize PairDataPoint vector, if needed - // block_size is defined in get_series_data() + // Set the pair dimensions if(pd_block.empty()) { - pd_block.resize(conf_info.block_size); - for(auto &pd : pd_block) pd.extend(n_series_pair); + set_pair_dims(pd_block, conf_info.block_size, n_series_pair); } + // Set the gradient pair dimensions + if(gx_pd_block.empty() || gy_pd_block.empty()) { + int n_grad_pd = conf_info.block_size * conf_info.get_n_grad(); + set_pair_dims(gx_pd_block, n_grad_pd, n_series_pair); + set_pair_dims(gy_pd_block, n_grad_pd, n_series_pair); + } + // Beginning of each data pass if(i_series == 0) { - // Re-initialize the PairDataPoint objects + // Re-initialize the pairs for(auto &pd : pd_block) { pd.erase(); pd.set_climo_cdf_info_ptr(&conf_info.cdf_info); } + // Re-initialize the gradient pairs + if(do_grad) { + for(auto &pd : gx_pd_block) pd.erase(); + for(auto &pd : gy_pd_block) pd.erase(); + } + // Starting grid point i_point = i_read*conf_info.block_size; @@ -949,10 +981,14 @@ void process_scores() { set_range(obs_dp.lead(), obs_lead_beg, obs_lead_end); // Store matched pairs for each grid point - for(int i=0; if_na.n() + + (aggr_file.empty() ? 0 : read_aggr_total(i_grid)); int n_series = n_series_pair + n_series_aggr; // Check for the required number of matched pairs if(n_valid / (double) n_series < conf_info.vld_data_thresh) { mlog << Debug(4) - << "[" << i+1 << " of " << conf_info.block_size + << "[" << i_block+1 << " of " << conf_info.block_size << "] Skipping point (" << x << ", " << y << ") with " << n_valid << " of " << n_series << " valid matched pairs.\n"; @@ -1001,7 +1095,7 @@ void process_scores() { } else { mlog << Debug(4) - << "[" << i+1 << " of " << conf_info.block_size + << "[" << i_block+1 << " of " << conf_info.block_size << "] Processing point (" << x << ", " << y << ") with " << n_valid << " of " << n_series << " valid matched pairs.\n"; } @@ -1011,27 +1105,27 @@ void process_scores() { (conf_info.output_stats[STATLineType::fho].n() + conf_info.output_stats[STATLineType::ctc].n() + conf_info.output_stats[STATLineType::cts].n()) > 0) { - do_categorical(i_point+i, &pd_block[i]); + do_categorical(i_grid, pd_ptr); } // Compute multi-category contingency table counts and statistics if(!conf_info.fcst_info[0]->is_prob() && (conf_info.output_stats[STATLineType::mctc].n() + conf_info.output_stats[STATLineType::mcts].n()) > 0) { - do_multicategory(i_point+i, &pd_block[i]); + do_multicategory(i_grid, pd_ptr); } // Compute continuous statistics if(!conf_info.fcst_info[0]->is_prob() && conf_info.output_stats[STATLineType::cnt].n() > 0) { - do_continuous(i_point+i, &pd_block[i]); + do_continuous(i_grid, pd_ptr); } // Compute partial sums if(!conf_info.fcst_info[0]->is_prob() && (conf_info.output_stats[STATLineType::sl1l2].n() + conf_info.output_stats[STATLineType::sal1l2].n()) > 0) { - do_partialsums(i_point+i, &pd_block[i]); + do_partialsums(i_grid, pd_ptr); } // Compute probabilistics counts and statistics @@ -1040,11 +1134,27 @@ void process_scores() { conf_info.output_stats[STATLineType::pstd].n() + conf_info.output_stats[STATLineType::pjc].n() + conf_info.output_stats[STATLineType::prc].n()) > 0) { - do_probabilistic(i_point+i, &pd_block[i]); + do_probabilistic(i_grid, pd_ptr); } - } // end for i + // Compute gradient statistics + if(do_grad) { + + // Loop over the gradient options + for(int i_grad=0; i_gradf_na, pd_gy->f_na, + pd_gx->o_na, pd_gy->o_na, pd_gx->wgt_na); + + // Aggregate current gradients with previous statistics + if(aggr_file.nonempty()) { + + // Aggregate gradients + GRADInfo aggr_grad; + read_aggr_grad(n, grad_info, aggr_grad); + grad_info += aggr_grad; + } + + // Add statistic value for each possible GRAD column + for(int j=0; j data(grid.nx()*grid.ny()); // Write output for each stat_data map entry - for(auto &key : stat_data_keys) { + for(const auto &key : stat_data_keys) { - NcVarData *ptr = &stat_data[key]; + const NcVarData *ptr = &stat_data[key]; // Add a new variable to the NetCDF file NcVar nc_var = add_var(nc_out, key, ncFloat, lat_dim, lon_dim, deflate_level); @@ -2287,9 +2525,9 @@ void write_stat_data() { add_att(&nc_var, "_FillValue", bad_data_float); add_att(&nc_var, "name", ptr->name); add_att(&nc_var, "long_name", ptr->long_name); - if(ptr->fcst_thresh.length() > 0) add_att(&nc_var, "fcst_thresh", ptr->fcst_thresh); - if(ptr->obs_thresh.length() > 0) add_att(&nc_var, "obs_thresh", ptr->obs_thresh); - if(!is_bad_data(ptr->alpha)) add_att(&nc_var, "alpha", ptr->alpha); + if(!ptr->fcst_thresh.empty()) add_att(&nc_var, "fcst_thresh", ptr->fcst_thresh); + if(!ptr->obs_thresh.empty()) add_att(&nc_var, "obs_thresh", ptr->obs_thresh); + if(!is_bad_data(ptr->alpha)) add_att(&nc_var, "alpha", ptr->alpha); // Store the data for(int x=0; x &pd_block, + int n_pd, int n_pair) { + + // Resize to the number of objects + pd_block.resize(n_pd); + + // Reserve space for the number of pairs + for(auto &pd : pd_block) pd.extend(n_pair); + + return; +} + + +//////////////////////////////////////////////////////////////////////// + +static void clean_up() { // Close the output NetCDF file if(nc_out) { @@ -2347,7 +2600,7 @@ void clean_up() { << "Output file: " << out_file << "\n"; delete nc_out; - nc_out = (NcFile *) nullptr; + nc_out = nullptr; } // Close the aggregate NetCDF file @@ -2365,7 +2618,7 @@ void clean_up() { //////////////////////////////////////////////////////////////////////// -void usage() { +__attribute__((noreturn)) static void usage() { cout << "\n*** Model Evaluation Tools (MET" << met_version << ") ***\n\n" @@ -2436,58 +2689,59 @@ void usage() { //////////////////////////////////////////////////////////////////////// -void set_fcst_files(const StringArray & a) { +static void set_fcst_files(const StringArray & a) { fcst_files = a; } //////////////////////////////////////////////////////////////////////// -void set_obs_files(const StringArray & a) { +static void set_obs_files(const StringArray & a) { obs_files = a; } //////////////////////////////////////////////////////////////////////// -void set_both_files(const StringArray & a) { +static void set_both_files(const StringArray & a) { set_fcst_files(a); set_obs_files(a); } //////////////////////////////////////////////////////////////////////// -void set_aggr(const StringArray & a) { +static void set_aggr(const StringArray & a) { aggr_file = a[0]; } //////////////////////////////////////////////////////////////////////// -void set_paired(const StringArray & a) { +static void set_paired(const StringArray &) { paired = true; } //////////////////////////////////////////////////////////////////////// -void set_out_file(const StringArray & a) { +static void set_out_file(const StringArray & a) { out_file = a[0]; } //////////////////////////////////////////////////////////////////////// -void set_config_file(const StringArray & a) { +static void set_config_file(const StringArray & a) { config_file = a[0]; } //////////////////////////////////////////////////////////////////////// -void set_compress(const StringArray & a) { +static void set_compress(const StringArray & a) { compress_level = atoi(a[0].c_str()); } //////////////////////////////////////////////////////////////////////// -void parse_long_names() { +static void parse_long_names() { ifstream f_in; - ConcatString line, key; + ConcatString line; + ConcatString key; StringArray sa; ConcatString file_name = replace_path(stat_long_name_file); diff --git a/src/tools/core/series_analysis/series_analysis.h b/src/tools/core/series_analysis/series_analysis.h index 73b2f3d6f6..8cfa318e6f 100644 --- a/src/tools/core/series_analysis/series_analysis.h +++ b/src/tools/core/series_analysis/series_analysis.h @@ -72,8 +72,10 @@ static const char * total_name = "TOTAL"; //////////////////////////////////////////////////////////////////////// // Input files -static StringArray fcst_files, found_fcst_files; -static StringArray obs_files, found_obs_files; +static StringArray fcst_files; +static StringArray found_fcst_files; +static StringArray obs_files; +static StringArray found_obs_files; static GrdFileType ftype = FileType_None; static GrdFileType otype = FileType_None; static ConcatString aggr_file; diff --git a/src/tools/core/series_analysis/series_analysis_conf_info.cc b/src/tools/core/series_analysis/series_analysis_conf_info.cc index fd19bf61bc..b2c73d8b94 100644 --- a/src/tools/core/series_analysis/series_analysis_conf_info.cc +++ b/src/tools/core/series_analysis/series_analysis_conf_info.cc @@ -80,6 +80,8 @@ void SeriesAnalysisConfInfo::clear() { mask_poly_file.clear(); mask_poly_name.clear(); mask_area.clear(); + grad_dx.clear(); + grad_dy.clear(); block_size = bad_data_int; vld_data_thresh = bad_data_double; hss_ec_value = bad_data_double; @@ -131,14 +133,16 @@ void SeriesAnalysisConfInfo::read_config(const char *default_file_name, void SeriesAnalysisConfInfo::process_config(GrdFileType ftype, GrdFileType otype) { - int i, n; + int i; + int n; ConcatString s; StringArray sa; ThreshArray cur_ta; VarInfoFactory info_factory; - Dictionary *fdict = (Dictionary *) nullptr; - Dictionary *odict = (Dictionary *) nullptr; - Dictionary i_fdict, i_odict; + Dictionary *fdict = nullptr; + Dictionary *odict = nullptr; + Dictionary i_fdict; + Dictionary i_odict; BootInfo boot_info; map::iterator it; @@ -188,6 +192,7 @@ void SeriesAnalysisConfInfo::process_config(GrdFileType ftype, bool do_cnt = (output_stats[STATLineType::sl1l2].n() + output_stats[STATLineType::sal1l2].n() + output_stats[STATLineType::cnt].n()) > 0; + bool do_grad = output_stats[STATLineType::grad].n() > 0; // Conf: fcst.field and obs.field fdict = conf.lookup_array(conf_key_fcst_field); @@ -339,6 +344,24 @@ void SeriesAnalysisConfInfo::process_config(GrdFileType ftype, } } // end for i + // Parse gradients + if(do_grad) { + + // Conf: gradient + Dictionary *d = conf.lookup_dictionary(conf_key_gradient); + grad_dx = d->lookup_int_array(conf_key_dx); + grad_dy = d->lookup_int_array(conf_key_dy); + + // Check for the same length + if(grad_dx.n() != grad_dy.n()) { + mlog << Error << "\nSeriesAnalysisConfInfo::process_config() -> " + << "The gradient dx and dy arrays must have " + << "the same length (" << grad_dx.n() << " != " + << grad_dy.n() << ").\n\n"; + exit(1); + } + } + // Conf: block_size block_size = conf.lookup_int(conf_key_block_size); @@ -445,7 +468,8 @@ void SeriesAnalysisConfInfo::process_config(GrdFileType ftype, //////////////////////////////////////////////////////////////////////// void SeriesAnalysisConfInfo::process_masks(const Grid &grid) { - MaskPlane mask_grid, mask_poly; + MaskPlane mask_grid; + MaskPlane mask_poly; ConcatString name; mlog << Debug(2) @@ -461,7 +485,7 @@ void SeriesAnalysisConfInfo::process_masks(const Grid &grid) { mask_poly_file = conf.lookup_string(conf_key_mask_poly); // Parse out the masking grid - if(mask_grid_file.length() > 0) { + if(!mask_grid_file.empty()) { mlog << Debug(3) << "Processing grid mask: " << mask_grid_file << "\n"; parse_grid_mask(mask_grid_file, grid, mask_grid, mask_grid_name); @@ -469,7 +493,7 @@ void SeriesAnalysisConfInfo::process_masks(const Grid &grid) { } // Parse out the masking polyline - if(mask_poly_file.length() > 0) { + if(!mask_poly_file.empty()) { mlog << Debug(3) << "Processing poly mask: " << mask_poly_file << "\n"; parse_poly_mask(mask_poly_file, grid, mask_poly, mask_poly_name); diff --git a/src/tools/core/series_analysis/series_analysis_conf_info.h b/src/tools/core/series_analysis/series_analysis_conf_info.h index 01aff16098..0d2597f12a 100644 --- a/src/tools/core/series_analysis/series_analysis_conf_info.h +++ b/src/tools/core/series_analysis/series_analysis_conf_info.h @@ -33,15 +33,15 @@ class SeriesAnalysisConfInfo { void init_from_scratch(); - // Counts based on the contents of the config file - int n_fcst; // Number of forecast fields - int n_obs; // Number of observation fields - public: // Series-Analysis configuration object MetConfig conf; + // Counts based on the contents of the config file + int n_fcst; // Number of forecast fields + int n_obs; // Number of observation fields + // Store data parsed from the Series-Analysis configuration object ConcatString model; // Model name ConcatString desc; // Description @@ -72,6 +72,9 @@ class SeriesAnalysisConfInfo { ConcatString mask_poly_name; // Name of masking poly area MaskPlane mask_area; + IntArray grad_dx; // Gradient step size in the X direction + IntArray grad_dy; // Gradient step size in the Y direction + int block_size; // Number of grid points to read concurrently double vld_data_thresh; // Minimum valid data ratio for each point double hss_ec_value; // HSS expected correct value @@ -96,13 +99,15 @@ class SeriesAnalysisConfInfo { // Dump out the counts int get_n_fcst() const; int get_n_obs() const; + int get_n_grad() const; }; //////////////////////////////////////////////////////////////////////// -inline int SeriesAnalysisConfInfo::get_n_fcst() const { return n_fcst; } -inline int SeriesAnalysisConfInfo::get_n_obs() const { return n_obs; } inline int SeriesAnalysisConfInfo::get_compression_level() { return conf.nc_compression(); } +inline int SeriesAnalysisConfInfo::get_n_fcst() const { return n_fcst; } +inline int SeriesAnalysisConfInfo::get_n_obs() const { return n_obs; } +inline int SeriesAnalysisConfInfo::get_n_grad() const { return grad_dx.n(); } //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/stat_analysis/parse_stat_line.cc b/src/tools/core/stat_analysis/parse_stat_line.cc index 80a4dd7102..b113b5ba9f 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.cc +++ b/src/tools/core/stat_analysis/parse_stat_line.cc @@ -315,6 +315,18 @@ void parse_grad_line(STATLine &l, GRADInfo &grad_info) { grad_info.egbar = atof(l.get_item("EGBAR")); grad_info.dx = atoi(l.get_item("DX")); grad_info.dy = atoi(l.get_item("DY")); + grad_info.fgmag = atof(l.get_item("FGMAG")); + grad_info.ogmag = atof(l.get_item("OGMAG")); + + double mag_rmse = atof(l.get_item("MAG_RMSE")); + grad_info.mag_mse = (is_bad_data(mag_rmse) ? + bad_data_double : + mag_rmse * mag_rmse); + + double lap_rmse = atof(l.get_item("LAPLACE_RMSE")); + grad_info.lap_mse = (is_bad_data(lap_rmse) ? + bad_data_double : + lap_rmse * lap_rmse); return; }