Enhance the Grid-Stat GRAD line type with additional gradient vector-based statistics to measure sharpness

[JohnHalleyGotway]

Describe the Enhancement

As discussed on Nov 18, 2024 by @DanielAdriaansen and @JohnHalleyGotway (see meeting notes), consider adding statistics to quantify "sharpness", especially in AI/ML output, using the methods defined in this paper:
An Investigation of Metrics to Evaluate the Sharpness in AI-Generated Meteorological Imagery

Grid-Stat already writes a gradient (GRAD) output line type that contains the S1 score and it's components. The existing statistics are computed by summing the X- and Y-gradients. The methods described in this paper use the X- and Y-gradients to form a gradient vector, and then compute the magnitude of those vectors in both the forecast and observation fields.

Consider adding 4 new GRAD columns for the following methods:

1. Mean Gradient Magnitude (grad-mag) of the forecast and observation fields is the average of norm of the gradient vectors in each field.
2. Gradient RMSE (grad-rmse) is computed by comparing the forecast and observation gradient vector magnitudes.
3. Laplace RMSE (laplace-rmse) is the RMSE of the divergence of each gradient vector rather than the magnitude.

Also enhance Grid-Stat to include the gradient vector magnitude and divergence fields, if nc_pairs_flag.gradient = TRUE in the Grid-Stat config file.

Also enhance Stat-Analysis to aggregate these new statistics properly.

Prior to adding new columns, update the version number of MET from 12.0 to 12.1.

Recommend NOT adding the following:

1. Image Intensity reports the min/max/mean of the values within each image. The CNT Line Type already includes the mean forecast FBAR and observation OBAR values. While the min/max data values are not actually written, they don't really belong in the gradient line type since they are NOT based on gradients.
2. Root Mean Squared Error is already included in the CNT line type.
3. Total Variation is summed across an image and not normalized by its size. Recommend not including it because its dependence on the number of grid points makes it much less useful.
4. Gradient Total Variation is also summed and dependent on the number of grid points.

Requires further investigation:

1. Structural Similarity Index Measure (ssim) is a similarity measure between two images based on a weighted combination of three simpler comparisons: luminance, contrast, and structure. If this were well-defined, we could consider adding it to the CNT line type.

• See reference on Wikipedia

I'd also recommend that we review the existing implementation of the S1 score to confirm that the Appendix C Equations (and also the WGNE website) actually match the computation of sums starting on this line of code.

      // Gradient sums
      fgbar += wgt * (fabs(fgx_na[i]) + fabs(fgy_na[i]));
      ogbar += wgt * (fabs(ogx_na[i]) + fabs(ogy_na[i]));
      mgbar += wgt * (max(fabs(fgx_na[i]), fabs(ogx_na[i])) +
                      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]));
      total++;

Where the forecast gradient is given by (fgx_na[i], fgy_na[i]) and observed gradient is given by (ogx_na[i], ogy_na[i]).

Time Estimate

2 days?

Sub-Issues

Consider breaking the enhancement down into sub-issues.
None needed.

Relevant Deadlines

NRL charging must be completed on 12/30/24

Funding Source

FY25 Q1 NRL METplus 7730022

Define the Metadata

Assignee

• Select engineer(s) or no engineer required
• Select scientist(s) or no scientist required

Labels

• Review default alert labels
• Select component(s)
• Select priority
• Select requestor(s)

Milestone and Projects

• Select Milestone as a MET-X.Y.Z version, Consider for Next Release, or Backlog of Development Ideas
• For a MET-X.Y.Z version, select the MET-X.Y.Z Development project

Define Related Issue(s)

Consider the impact to the other METplus components.

• METplus, MET, METdataio, METviewer, METexpress, METcalcpy, METplotpy
• Need METdataio issue to update database schema.
• Need METviewer issue to plot these new statistics.

Enhancement Checklist

See the METplus Workflow for details.

• Complete the issue definition above, including the Time Estimate and Funding Source.
• Fork this repository or create a branch of develop.
  Branch name: feature_<Issue Number>_<Description>
• Complete the development and test your changes.
• Add/update log messages for easier debugging.
• Add/update unit tests.
• Add/update documentation.
• Push local changes to GitHub.
• Submit a pull request to merge into develop.
  Pull request: feature <Issue Number> <Description>
• Define the pull request metadata, as permissions allow.
  Select: Reviewer(s) and Development issue
  Select: Milestone as the next official version
  Select: MET-X.Y.Z Development project for development toward the next official release
• Iterate until the reviewer(s) accept and merge your changes.
• Delete your fork or branch.
• Close this issue.

[JohnHalleyGotway]

I think we should validate that the existing contents of the GRAD line type is actually correct.

Running in seneca:/d1/personal/johnhg/MET/MET_development/MET-feature_3024_GRAD/test/run_gs.sh to compute GRAD output for (dx, dy) = (1, 1) in out/grid_stat_120000L_20050807_120000V_grad.txt, the GRAD line type for TMP/Z2 reports an S1 score of 76.36778:

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.0 FCST  NA   120000    20050807_120000 20050807_120000 000000   20050807_120000 20050807_120000 TMP      K          Z2       TMP     K         Z2      ANALYS FULL    NEAREST     1           NA          NA         NA         NA    GRAD       2401 1.12786 1.1541 1.64431 1.25573 76.36778 108.80549    0.97726  1  1 0.99423 1.01573  1.12915      2.08693

The following R commands report an S1 score of 4.048109:

R
library(SpatialVx)
library(ncdf4)
nc = nc_open("/d1/personal/johnhg/MET/MET_development/MET-feature_3024_GRAD/test/out/grid_stat_120000L_20050807_120000V_pairs.nc", write=F)
fcst = ncvar_get(nc, nc$var[["FCST_TMP_Z2_FULL"]])
obs = ncvar_get(nc, nc$var[["OBS_TMP_Z2_FULL"]])
S1(obs, xhat=fcst)
[1] 4.048109

However, the default gradient function in SpatialVx does a large amount of smoothing:

KernelGradFUN <- function(x,ktype="LoG", nx=10, ny=12, sigma=1) return(kernel2dsmooth(x,kernel.type=ktype, nx=nx, ny=ny, sigma=sigma))

Switching to using a simple difference (as MET uses) produces a much more similar result:

S1(obs, xhat=fcst, gradFUN="diff")
[1] 75.88391

While 75.88391 and MET's 76.36778 are not identical, they are rather close. It seems that the diff function in R only operates in a single dimension rather than both.
So there are likely differences in the algorithms.

[JohnHalleyGotway]

Need to confirm that this is the correct way to compute the laplace_mse for each point that we can take the square root of the MSE prior to writing LAPLACE_RMSE to the output:

      double diff = (fgx_na[i] + fgy_na[i]) - (ogx_na[i] - ogy_na[i]);
      lap_mse += wgt * (diff * diff);

[JohnHalleyGotway]

• Already switched from V12.0 to V12.1 via PR Feature #3024 V12.1 (PR 1 of 2) #3026.
• Updated MET to add 4 new columns to GRAD: FGMAG, OGMAG, MAG_RMSE, LAPLACE_RMSE
• Added the columns to the 2 header tables in data/table_files and internal/test_unit/hdr
• Updated Stat-Analysis to parse the data and aggregate the values.
• Updated the GRAD table in the User's Guide
• Need to confirm that the existing GRAD columns are actually correct. They are rather close, as described in this issue comment.
• Need to also updated Appendix C with the Equations. Recommend asking @bgbrowntollerud to review the result.
• Look for and address TODO comments in the code.

[JohnHalleyGotway]

@JohnHalleyGotway, @KathrynNewman, and @DanielAdriaansen met with @bgbrowntollerud on Dec 4, 2024 to discuss adding these new statistics to the GRAD line type. Barb confirmed that this is a reasonable thing to do and recommended that we proceed. She emphasized that the gradients really just measure the "texture" of the data and say nothing about the accuracy of the forecast.

We decided NOT to include SSIM in this work, but documented that as a potential future feature request via #3031. We also note that these gradient statistics could be of use in the Series-Analysis tool, as described in MET #3030.

[JohnHalleyGotway]

@KathrynNewman and @bgbrowntollerud please advise on the following:

• Compare this excerpt from An Investigation of Metrics to Evaluate the Sharpness in AI-Generated Meteorological Imagery:

Laplace RMSE is very similar to gradient RMSE, but instead of taking the magnitude of the gradient vector at each pixel, we compute the divergence of the gradient at each pixel, which is a way of quantifying the local shape of the gradient vector field.

To this line of code:

      double diff = (fgx_na[i] + fgy_na[i]) - (ogx_na[i] + ogy_na[i]);
      lap_mse += wgt * (diff * diff);

Where fgx, fgy, ogx, and ogy are the forecast/observation gradients at a point in the X- and Y-dimensions.
Is this the correct way to increment the Laplacian Mean-Squared-Error for each grid point?
Is it just the sum of the forecast wind components minus the sum of the observed wind components?
Or is the Laplacian more complex than that?

• Please also review the formatting of the equations I added to this section of Appendix C.

The ones for FGMAG, OGMAG, MAG_RMSE, and LAPLACE_RMSE are new. Is this notation correct?

[bgbrowntollerud]

John, sorry for the delay. Unfortunately I am not an expert by any means on this topic. What you suggest sounds logical to me, but I can't confirm for sure. Do we have any experts on this topic in RAL (perhaps David-John?) who could help us? I'll take another look at the paper to see if I can surmise anything.