Correction to crps_climo calculation, CRPSS

j-opatz commented 4 years ago

Describe the Enhancement

During the use of Ensemble Stat, a user observed that the CRPSS from MET did not match the CRPSS from VSDB. Further investigation showed that the CRPS matched between the two calculations. As the climatology CRPS is the added factor for calculating CRPSS, this is the most likely source of the error.

Investigation has revealed that MET's handling of climo data for ensemble verification differs from how it's done at EMC. Ensemble-Stat's use of climo data was made consistent with Point-Stat and Grid-Stat, where individual obs are subset into climo CDF bins, stats are computed within each bin, and the mean of those stats are reported. However, that is NOT how NOAA/EMC uses climo data for ensemble verification.

For met-10.0.0, we should remove the climo CDF subsetting logic. Instead, use the climo mean and standard deviation, assuming a normal distribution to sample from the climo distribution. Process those climo values as an ensemble. By default, NOAA/EMC uses a climo ensemble of size 10, but this should be configurable.

Since the climo CRPS and CRPS skill score are now dependent on both the climo mean and standard deviation, we need to add the climo standard deviation to the ORANK line type. That's necessary for Stat-Analysis to be aggregate ORANK lines together when computing the ECNT output line type.

Time Estimate

Not established; dependent on ID of problem

Sub-Issues

[ ] Discover where in the crps_climo calculation the error is made

Relevant Deadlines

List relevant project deadlines here or state NONE.

Funding Source

2791541

Define the Metadata

Assignee

[x] Select engineer(s) or no engineer required: John HG
[x] Select scientist(s) or no scientist required: John O

Labels

[x] Select component(s)
[x] Select priority
[x] Select requestor(s)

Projects and Milestone

[x] Review projects and select relevant Repository and Organization ones
[x] Select milestone

Define Related Issue(s)

Consider the impact to the other METplus components.

[x] METplus, MET, METdb, METviewer, METexpress, METcalcpy, METplotpy
[x] Add METviewer and/or METdatadb issue for the ORANK line type: renaming CLIMO to CLIMO_MEAN and adding CLIMO_STDEV
[x] dtcenter/METviewer#266
[x] dtcenter/METdatadb#49

Enhancement Checklist

See the METplus Workflow for details.

[x] Complete the issue definition above.
[x] Fork this repository or create a branch of develop. Branch name: feature_<Issue Number>_<Description>
[x] Complete the development and test your changes.
[x] Add/update unit tests.
[x] Add/update documentation.
[x] Push local changes to GitHub.
[x] Submit a pull request to merge into develop. Pull request: feature <Issue Number> <Description>
[x] Iterate until the reviewer(s) accept and merge your changes.
[x] Delete your fork or branch.
[x] Close this issue.

JohnHalleyGotway commented 4 years ago

Please find Binbin's sample test data attached. crpss_met_vsdb.tar.gz

JohnHalleyGotway commented 4 years ago

Met with EMC folks on 9/8/2020 and Binbin requested that CRPS climo to be added to the ECNT line type.

JohnHalleyGotway commented 3 years ago

Currently the CRPS climo is computed using only the climo mean. To make the logic consistent with NCEP (i.e. Binbin's example) we'll need both the climo mean and standard deviation. Use those to extract an ensemble of equally-likely climo values.

One implication is that we'll need to add CLIMO_STDEV to the ORANK line type. That's need to be able to aggregate ORANK line types into ECNT and recompute the Gaussian CRPSS.

JohnHalleyGotway commented 3 years ago

I realize that HiRA in Point-Stat also writes the ECNT line type. Here's a selection of columns taken from:

cat MET_test_output/climatology/point_stat_WMO_CLIMO_1.5DEG_120000L_20120409_120000V_ecnt.txt | awk '{print $10, $12, $18, $19, $24, $27, $28, $37, $39, $40}'

FCST_VAR FCST_LEV INTERP_MTHD INTERP_PNTS LINE_TYPE CRPS CRPSS CRPSCL CRPSCL_EMP CRPSS_EMP
TMP P850 NBRHD_SQUARE 4 ECNT 0.80461 0.61819 2.10735 2.05663 0.57563
TMP P850 NBRHD_SQUARE 9 ECNT 0.85146 0.59596 2.10735 2.05663 0.5682
TMP P850 NBRHD_SQUARE 16 ECNT 0.96009 0.54441 2.10735 2.05663 0.51426
TMP P850 NBRHD_SQUARE 25 ECNT 1.04051 0.50625 2.10735 2.05663 0.47649

Note that the CRPSCL and CRPSCL_EMP columns remain constant. That is because we're using the climo_mean/climo_stdev at the obs location to compute them. So while the neighborhood size of the forecast ensemble members change, the reference climatology remains fixed across all neighborhood sizes.

This is probably just fine, but is interesting to note.

JohnHalleyGotway commented 3 years ago

I compared the result of the updated code to the sample data from Binbin as follows.

(1) Wrote to_orank.py to convert CRPSS.data to 14760 ORANK lines as if they were written by Ensemble-Stat. to_orank.py.txt That produces this result: CRPSS_orank.txt

(2) Ran that through Stat-Analysis:

stat_analysis -lookin CRPSS_orank.txt -job aggregate_stat -line_type ORANK -out_line_type ECNT -v 4 -out_bin_size 0.1 -out orank_to_ecnt.txt
COL_NAME: TOTAL N_ENS    CRPS   CRPSS     IGN      ME   RMSE  SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR  CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP
    ECNT: 14760    20 4.94721 0.49744 3.64504 1.19287 9.2985 10.0277      NA        NA          NA               NA 9.84398  5.10959   10.24482   0.50125

orank_to_ecnt.txt (3) For easier comparison, I tweaked CRPSS.f to use a constant weight of 1.0. I also wrote CRPSS as the mean of the CPRSS for each point and the CRPSS of the avg CRPS scores:

 CRPSF =   5.10960197      CRPSC =   9.77567101    
 Avg of CRPSS =  0.234634370      CRPSS from Avg =  0.477314472

Results:

EMC is using the empirical Hersbach CRPS method.
MET ensemble CRPS_EMP of 5.10959 is ALMOST IDENTICAL to the VSDB ensemble CRPSF of 5.10960197.
MET climo CRPSCL_EMP of 10.24482 is SIMILAR TO to the VSDB climo CRPSC of 9.77567101. MET computes this using the mean and stdev of the 9 climo vals while VSDB uses the 9 climo vals directly. And that will definitely cause a difference.
MET CRPSS_EMP of 0.50125 is SIMILAR TO VSDB CRPSS from the average of 0.477314472.

Notes:

I assume the CRPSS.f code has a bug. Instead of computing the CRPSS for each point and then averaging them, it should really be averaging the CRPSF and CRPSC across all those points and computing CRPSS once using the average scores. But I need Binbin to confirm.
For some reason this simple Stat-Analysis job takes almost 2.5 minutes to run. I should figure out what's slowing it down so much and fix it!

JohnHalleyGotway commented 3 years ago

From NOAA/EMC on 3/16/21: We just completed the retrospective verification of GEFS with MET_10.0.0_beta4 over past 4 weeks, and obtained the averaged CRPSS for different forecast hours. The comparison shows that CRPSS-MET and CRPSS-VSDB are very close. See attached plot for Z500 over N. Hemisphere. This is very nice. So the new CRPSS code in MET_10.0.beta4 is successful. CRPSS-GEFS_Z500-Compare

dtcenter / MET