Fix bug in the computation of the climatological CDF value.

[JohnHalleyGotway]

Describe the Problem

When testing the verification of the climatology-based probabilities for NOAA/CPC with Grid-Stat, John Opatz discovered an apparent bug in MET's computation of the climatological cumulative distribution function values.

While the issue was identified using CPC data, it is also apparent in the existing unit test output:
climatology/grid_stat_PROB_GFS_CLIMO_1.0DEG_150000L_20120409_120000V_pairs.nc

I also found a bug in derive_climo_prob(). When the observation threshold is defined relative to climatology (e.g. >CDP67), the climatological probability value should be constant. However, the constant value should depend on the inequality type, where <CDP{N} has climo probability p=N/100 while >CDP{N} has probability 1.0-p.

This task is to fix both of these bugs.

Expected Behavior

The climo CDF value is computed as the area (between 0 and 1) of the climo distribution to the left of the observation value. However, MET is not computing that correctly.

Environment

Describe your runtime environment:
1. Machine: NB output on kiowa.

To Reproduce

Describe the steps to reproduce the behavior:
1. Inspect the file: climatology/grid_stat_PROB_GFS_CLIMO_1.0DEG_150000L_20120409_120000V_pairs.nc
2. Using ncview, note the values for i=8, j=77:
obs = 13.1514, climo_mean = 3.21593, climo_stdev = 7.9448, climo_cdf = 0.359584
3. Note that the obs > climo_mean but climo_cdf < 0.5. When the obs value greater than the climo mean, the CDF should be greater than 0.5.

So there's a bug in this computation!

Relevant Deadlines

List relevant project deadlines here or state NONE.

Funding Source

2788881

Define the Metadata

Assignee

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

Labels

• Select component(s)
• Select priority
• Select requestor(s)

Projects and Milestone

• Review projects and select relevant Repository and Organization ones or add "alert:NEED PROJECT ASSIGNMENT" label
• Select milestone to relevant bugfix version

Define Related Issue(s)

Consider the impact to the other METplus components.

• METplus, MET, METdatadb, METviewer, METexpress, METcalcpy, METplotpy

Bugfix 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 main_<Version>.
  Branch name: bugfix_<Issue Number>_main_<Version>_<Description>
• Fix the bug 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 main_<Version>.
  Pull request: bugfix <Issue Number> main_<Version> <Description>
• Define the pull request metadata, as permissions allow.
  Select: Reviewer(s), Project(s), Milestone, and Linked issues
• Iterate until the reviewer(s) accept and merge your changes.
• Delete your fork or branch.
• Complete the steps above to fix the bug on the develop branch.
  Branch name: bugfix_<Issue Number>_develop_<Description>
  Pull request: bugfix <Issue Number> develop <Description>
• Close this issue.

[JohnHalleyGotway]

Good news, the bug was easy to find... passing function arguments in the wrong order in Grid-Stat:

- write_nc((string)"CLIMO_CDF", normal_cdf(cmn_dp, csd_dp, obs_dp),
+ write_nc((string)"CLIMO_CDF", normal_cdf(obs_dp, cmn_dp, csd_dp),

Bad news, this won't have any impact on the computation of Brier score and BSS that @j-opatz is investigating. The bug is confined to the contents of the NetCDF matched pairs Grid-Stat output file but those values aren't actually used in any of the statistics computations.

[JohnHalleyGotway]

@j-opatz and @bgbrowntollerud please take a look at the following code:

MET/met/src/libcode/vx_statistics/pair_base.cc

Line 1029 in 76253a9

// For CDP threshold types, the climo_prob is constant.

   // For CDP threshold types, the climo_prob is constant.
   if(othresh.get_ptype() == perc_thresh_climo_dist) {
      climo_prob.add_const(othresh.get_pvalue()/100.0, n_mn);
   }

I think this logic is not sufficient. This code means that thresholds of >=CDP67 and <=CDP67 would both result in a constant climo probability of 0.67. I do think it should be a constant climo probability, but it should depend on the threshold type. Am I correct in thinking the following?

>CDP{N} should be a constant climo probability of 1.0 - N/100 (e.g. >CDP67 has 0.33 probability).
<CDP{N} should be a constant climo probability of N/100 (e.g. <CDP67 has 0.67 probability).
And ==CDP{N} would have a probability of 0 while !=CDP{N} would have a probability of 1.0.

Mathematically, the probability of being exactly equal to a number is 0 while being not-equal is 1.0. But in practice, that isn't very useful. Should I instead error out if the user specifies CDP thresholds which are not >, >=, <, or <=?

Here's the updated logic:

   // For CDP threshold types, the climo probability is constant
   if(othresh.get_ptype() == perc_thresh_climo_dist) {

      // Climo probability varies based on the threshold type
      switch(othresh.get_type()) {

         case thresh_lt:
         case thresh_le:
            prob = othresh.get_pvalue()/100.0;
            break;

         case thresh_eq:
            prob = 0.0;
            break;

         case thresh_ne:
            prob = 1.0;
            break;

         case thresh_gt:
         case thresh_ge:
            prob = 1.0 - othresh.get_pvalue()/100.0;
            break;

         default:
            mlog << Error << "\n\nderive_climo_prob() -> "
                 << "climatological threshold \"" << othresh.get_str()
                 << "\" cannot be converted to a probability!\n\n";
            exit(1);
            break;
      }

      // Add constant climo probability value
      climo_prob.add_const(prob, n_mn);
   }

[JohnHalleyGotway]

All of the units test ran fine but they did produce differences with the regression test output, as expected:
/d1/projects/MET/MET_pull_requests/met-10.0.0_beta3/bugfix_1638/comp_dir.log

file1: MET-bugfix_1638_develop_climo_cdf_into_develop/test_output/climatology/ensemble_stat_NCEP_1.0DEG_20120410_120000V.stat
ERROR: found     6 differences in row type RPS          column RPSS          - max abs: 0.41248 
ERROR: failed tests for RPS: num
--
file1: MET-bugfix_1638_develop_climo_cdf_into_develop/test_output/climatology/grid_stat_PROB_GFS_CLIMO_1.0DEG_150000L_20120409_120000V_pairs.nc
ERROR: NetCDF headers differ:
--
file1: MET-bugfix_1638_develop_climo_cdf_into_develop/test_output/climatology/grid_stat_WMO_CLIMO_1.5DEG_240000L_20120410_000000V_pairs.nc
ERROR: NetCDF headers differ:
--
file1: MET-bugfix_1638_develop_climo_cdf_into_develop/test_output/climatology/point_stat_WMO_CLIMO_1.5DEG_120000L_20120409_120000V_rps.txt
ERROR: found    16 differences in row type RPS          column RPSS          - max abs: 0.37573 
ERROR: failed tests for RPS: num
--
file1: MET-bugfix_1638_develop_climo_cdf_into_develop/test_output/climatology/point_stat_WMO_CLIMO_1.5DEG_120000L_20120409_120000V.stat
ERROR: found    16 differences in row type RPS          column RPSS          - max abs: 0.37573 
ERROR: failed tests for RPS: num

The changes to the RPS output do make sense but the timestamps in the NetCDF files don't really:

COMPARING climatology/grid_stat_PROB_GFS_CLIMO_1.0DEG_150000L_20120409_120000V_pairs.nc
file1: MET-bugfix_1638_develop_climo_cdf_into_develop/test_output/climatology/grid_stat_PROB_GFS_CLIMO_1.0DEG_150000L_20120409_120000V_pairs.nc
file2: /d1/projects/MET/MET_regression/develop/NB20210125/MET-develop/test_output/climatology/grid_stat_PROB_GFS_CLIMO_1.0DEG_150000L_20120409_120000V_pairs.nc
ERROR: NetCDF headers differ:
75,76c75,76
< 		CLIMO_CDF_WIND_Z10_FULL:init_time = "20120409_120000" ;
< 		CLIMO_CDF_WIND_Z10_FULL:init_time_ut = "1333972800" ;
---
> 		CLIMO_CDF_WIND_Z10_FULL:init_time = "19590409_120000" ;
> 		CLIMO_CDF_WIND_Z10_FULL:init_time_ut = "-338644800" ;
145,146c145,146
--
COMPARING climatology/grid_stat_WMO_CLIMO_1.5DEG_240000L_20120410_000000V_pairs.nc
file1: MET-bugfix_1638_develop_climo_cdf_into_develop/test_output/climatology/grid_stat_WMO_CLIMO_1.5DEG_240000L_20120410_000000V_pairs.nc
file2: /d1/projects/MET/MET_regression/develop/NB20210125/MET-develop/test_output/climatology/grid_stat_WMO_CLIMO_1.5DEG_240000L_20120410_000000V_pairs.nc
ERROR: NetCDF headers differ:
989,990c989,990
< 		CLIMO_CDF_TMP_P850_FULL:init_time = "20120410_000000" ;
< 		CLIMO_CDF_TMP_P850_FULL:init_time_ut = "1334016000" ;
---
> 		CLIMO_CDF_TMP_P850_FULL:init_time = "256550410_000000" ;
> 		CLIMO_CDF_TMP_P850_FULL:init_time_ut = "747434995200" ;
1297,1298c1297,1298

Need to debug these timestamps further.

[JohnHalleyGotway]

Found/fixed the time differences in the Grid-Stat NetCDF matched pairs file for climo CDF variables. Instead of initializing the output variables with the observation field, using the climo mean field instead. That results in the output time stamps being consistent for the climo mean, standard deviation, and cdf fields... which is desirable.