Add RPS and RPSS to the output of Ensemble-Stat and Point-Stat.

[dwfncar]

Enhance MET statistics tools for climatological data by doing the following:
   - CRPSS (Continuous Ranked Probability Skill Score) to the RHIST line type.
   - RPSS (Ranked Probability Skill Score) to the RHIST line type.
   - Consider defining a new line type for Ensemble-Stat for continuous ensemble statistics (CRPS, IGN, CRPSS, RPSS). [MET-681] created by johnhg

[dwfncar]

Needed for global ensemble and AF deterministic by end of AOP2016 by jensen

[dwfncar]

AOP2017 VX1 NOAA or NGGPS FY17 by jensen

[dwfncar]

From Tara:

Here's one example from Space Weather that is gives the equation used for RPS and RPSS and some examples of the data. I'm still looking for other examples. by johnhg

[dwfncar]

From CPC - sample data are on dakota:/d3/projects/MET/NCEP_unification/cpc_data

To read the files in Python:

import numpy as np
data = np.fromfile('gefs-legacy-00z_rfcst-cal_tmean_20170508_6-10day.bin', 'f').reshape(19, 181, 360)

where there are 19 records, each a global 1-degree grid with 181 points in the y direction, 360 points in the x direction.

Hi Tara,
Regarding RPSS calculations - we apply it to categorical forecasts. A description of the methodology for how we do that is here: http://www.cpc.ncep.noaa.gov/products/verification/summary/index.php?page=tutorial
Just drop down for RPSS.

Our dynamic verification with categorical RPSS can be performed here:
http://www.vwt.ncep.noaa.gov/
Our static verification page, with spatial maps of RPSS are here
http://www.cpc.ncep.noaa.gov/products/verification/summary/index.php?page=map
'
Since you would also be applying RPSS to various percentiles, a quick general reference for the equation is here: https://journals.ametsoc.org/doi/full/10.1175/MWR3280.1
which just references the RPSS definition all of us use in the Wilks book:

Wilks, D. S., 1995: Statistical Methods in the Atmospheric Sciences. International Geophysics Series, Vol. 59, Academic Press, 467 pp.

I don't believe there are any other special assumptions we would do for RPSS applied to percentiles so I think it is straightforward. Regarding extremes, we do subset the forecasts to be more useful for our own application. Since we typically issue hazards based on a minimum forecast percent of exceeding or being lower than a specific percentile, we subset our model extremes tool as such, e.g. for tmax extremes, we would verify those GEFS calibrated forecasts from our tool where the forecast probability was > 20% of the forecast percentile being greater than or equal to the 85th percentile.

We have some example results here of extremes:
https://docs.google.com/presentation/d/1AaxbOBk2cUo4-Rm9F_V0wwHC4hxc20U2G126XvouSjE/edit#slide=id.g3f12034e57_0_579

If you have any other questions, please let us know!

NOTE: I ADDED Prob Extremes Tool Verification presentation PDF to issue as well. by jensen

[TaraJensen]

Charge 2780541

[lisagoodrich]

Original JIRA attachment location: https://sdg.rap.ucar.edu/jira/browse/MET-681?filter=12903

Prob Extremes Tool Verification - All vars 092018.pdf

RPS_RPSS_Sharpe_SpaceWx.pptx

[TaraJensen]

Data related to description of computing RPS/RPSS above are sitting on dakota: /d3/projects/MET/NCEP_unification/cpc_data/data-for-tara/

Description of the data were included in an email in Aug 2018. See below:

all of 2017 GEFS-reforecast-calibration 2m Temperature forecasts for the 6-10day period into the FTP directory you described above.

The data are in 32-bit floating point binary, and contain 19 records (probability of exceeding the 1st percentile, 2nd percentile, ..., 99th percentile). These are the percentiles if you need them:

[1, 2, 5, 10, 15, 20, 25, 33, 40, 50, 60, 67, 75, 80, 85, 90, 95, 98, 99]

To read the files in Python:

import numpy as np data = np.fromfile('gefs-legacy-00z_rfcst-cal_tmean_20170508_6-10day.bin', 'f').reshape(19, 181, 360)

where there are 19 records, each a global 1-degree grid with 181 points in the y direction, 360 points in the x direction.

NOTE: John put together a Python script to read this data and pass it into MET. That script is located at: https://dtcenter.org/community-code/model-evaluation-tools-met/sample-analysis-scripts (the first one)

[TressaFowler]

To compute RPS from this data, it must first be formatted into Mutually Exclusive and Collectively Exhaustive categories, (MECE). The data given to us are cumulative, so the probability in each of the 19 categories must be subtracted from that of the previous category to just get the individual category percentage. Obviously, the final category is 1-(sum of the others).

Next, observations of the event in question must be collected. These are typically all ones and zeros, with a 1 being placed in the bin for the temperature that was observed. I guess we would need to know what the climatology is to place a temperature value in a bin.

These are placed in the squared error equation from Wilks.

[TaraJensen]

From Binbin on 8/13/19 - NOTE: I have requested paths to input and output files. Will update issue when we get those

As we discussed at yesterday's MET meeting, you like to know the detailed computation procedure for RPS and PAC scores. I reviewed code and found the locations of PAC and RPS/RPSS subroutines for you. The entire VSDB verification package has been copied on theia, so you can copy the code and see the detailed computation procedures. The sorc code directory is in /scratch4/NCEPDEV/fv3-cam/save/Binbin.Zhou/work/grid2grid/verf_g2g.v3.0.13/sorc/verf_g2g_grid2grid_grib2.fd

The main program is grid2grid.f, which handles both single model and ensemble verification. For ensemble case, please just focus on "call EFS", (not EFS1, which is for regional ensemble). EFS sub calls another 2 major subs dist ( ) and prob ( ).

These two subs are Yuejian's codes. dist is for distribution scores (e.g. rmse, spread, PAC etc), saved in DIST.f file, and prob is for probabilistic scores (e.g. BS, BSS, RPS, RPSS, reliability, etc) saved in PROB.f file in the same directory.

In DIST.f, search "faccs" which is PAC, you can see how to calculate PAC. As I can see, PAC is similar to conventional AC using average value of 10 climatology bins as reference.

In PROB, search "call RRPS" , which is for RPS computation. Subroutine RRPS code is also in PROB.f file.

Please also note that the this code is for both regional and global ensemble verification. For regional ensemble verif, single event is used, while for global ensemble verif, 10 clim bins are used as 11 multiple events. In RRPS sub, ib is used to specify both cases:

 ib=1 for single event (regional, single event case)), 
 ib.ge.2 (e.g. ib=10,  for multiple events in global ensemble case) 

Please focus on "ib.ge.2" in RRPS for global ensemble case. RRPS sub is pretty short and Yuejian has provided detailed description/comments in the code.

Entire ensemble verification is grid-wide. In other words, loop on all grids and accumulate the results at all grids and out put a single value. The RPS/RPSS computation procedure has following structure in PROB.f:

For specific forecast hour, domain, and field:
Loop over all grid points: At each grid point,
Call RRPS( ) sub
Input: 20 member fcst, 10 clim bins, and 1 analysis Note: RPS is accumulated and averaged over all clim bins in rrps ( ) sub
Output: Bin averaged RPS, for fcst and clim, respectively Accumulate RPS, for fcst and clim, respectively End of grid loop Average RPS over all grids Divided by total number of grid), for fcst and clim, respectively Compute RPSS:
Using fcst RPS and clim RPS. Final output to VSDB file: grid-averaged RPS for fcst, RPS for clim, and RPSS

Hope this description is clear enough for you to understand the RPS/RPSS computation procedure. For even more detailed questions, I'll have refer to Yuejian who is the original author of the code. Since the structure code is pretty simple, it is not difficult to cover this Fortran code to C/C++ version.

Yuejian, Yan

If you have any suggestions or comments, please let us know.

Binbin

[TaraJensen]

Email from Binbin on 8/21/19: I created a sample file directory on theia for our GEFS verification package located in directory:

/scratch4/NCEPDEV/fv3-cam/save/Binbin.Zhou/work/grid2grid/verf_g2g.v3.0.13/sample_data

In this directory, 20 grib2 GEFS member files are those headed as "fcst", containing only one 24hr forecast (run at 12Z 20190818).

The analysis data file is obsv.grib.gefs2p5.ensbl which is GDAS analysis validated at 12Z 20190819.

The climatology data are those files ended with JUL, AUG and SEP. They are binary, 101 bin data from which 11 bin data can be derived (averaged) For detailed description please look at readme.txt file in the same directory.

The output VSDB file is GEFS2P5_2019081912.vsdb. The RPS and RPSS can be seen in "RPS" record line. The first 6 scores in each RPS record line are RPSf, RPSc, RPSS, CRPSf, CRPSc, CRPSS. "f" refers to forecast, "c" to climatology. So if you like to see the values of RPS and RPSS, just look at the first 3 values on RPS line records.

Our verification is over 6 sub-regions, identified with G2/XX, G2 means grid#2 (2.5x2.5degree), XX is NH, SH, TR, NA, AS, EU, etc (refers to N. and S. hemisphere, Tropical, N. America, Asia and Europe) . The la and lon for these regions are defined in control
file g2g.ctl.ensemble.

If I miss something or you have further questions, please let me know.

Binbin

[JohnHalleyGotway]

This came up during a telecon with AF on 10/4/2019. Bob Craig pointed out that they'd like to compute RPS during the application of HiRA where the categories for wind are defined by the Beaufort scale.

[JohnHalleyGotway]

MET on 10/22/2019 to discuss this with Barb, Eric, Tara, and John. We discussed the attached write-up that Barb put together.

Adding_RPS_to_MET_BGB.pdf

RPS is needed by: (1) WPC to process 19 GRIB records with "nested" probabilities. (2) Global Ensembles group in EMC when verifying an ensemble. (3) Bob Craig at the AF when applying HiRA in Point-Stat.

It is not immediately obvious where/how to add this functionality in MET. It doesn't fit in any other existing line types and will probably require a new one... perhaps two: one for RPS and RPSS and one for the individual RPS matched pairs (RPSMPR).

One option would be ...

• Enhance Ensemble-Stat to compute nested probabilities and write NetCDF output using an additional dimension for the changing thresholds.
  • This enables Ensemble-Stat to compute nested probabilities for relative frequencies, Neighborhood Ensemble Probabilities, and Neighborhood Max Ensemble Probabilities.
  • Would also need to add the same derivation logic to HiRA methodology in Point-Stat.
• Enhance Grid-Stat and Point-Stat by adding a new "is_rps_data = TRUE" (or something like that) flag to read these nested probabilities... either from the NetCDF output of Ensemble-Stat or from multiple GRIB records. So when is_rps_data = TRUE, we'd call data_plane_array() instead of data_plane().
• Further enhance Grid-Stat and Point-Stat to compute RPS/RPSS and write 2 new line types for RPS and RPSMPR (in Point-Stat).

[bgbrowntollerud]

I made a few changes and updates to the RPS document following our meeting yesterday - including a bit of re-organization and adding info about computing RPSS and "overall" accumulated RPS values. Please use this new version, and let me know if you have any questions! Notes on RPS Computation Updated 23Oct19.pdf

[JohnHalleyGotway]

Added new ERPS line type to the output of Ensemble-Stat and Point-Stat. Updated Stat-Analysis to parse and aggregate this new line type. Merged changes into develop and am rerunning the unit tests now.