EnsembleStat: compute gridded CRPSS, CRPS, correlation using two sets of reforecasts (GEFSv12 and GEFSv13) and AORC

[haksulee]

Hi METPlus support:

I am new to METPlus configuration.
My goal is to evaluate GEFSv13 precipitation reforecast data in reference to GEFSv12 reforecasts on a US domain. The observation data used is AORC.
I want to produce the US map plots of gridded ensemble statistics (CRPS, CRPSS, etc) and ensemble mean statistics (correlation, etc) for different lead hours on a 6-hr and 1-day time step.

question 1) how can I configure using two different sets of reforecast data (GEFSv12 and GEFSv13) ?
I see only one line for "FCST_ENSEMBLE_STAT_INPUT_DIR"; Is there an option to define one more forecast data to be used as a reference in computing CRPSS ?

question 2) how can I compute lead time specific statistics from, say, 100 different initializations?
I tried LOOP_BY=INIT or VALID, but looks to me they work for each INIT or VALID separately.

question 3) how can I produce gridded statistics outputs, for example, gridded CRPS, CRPSS, or correlation on a US map ?

question 4) Can EnsembleStat configuration be configured to compute 1-day APCP by summing 6-hrly APCP from 4 different lead times (6/12/18/24Z) from the same initialization time ?
If not, what other METPlus command should I use to compute 1-day APCP from 6-hr APCP?

Could you check my EnsembleStat.conf and user.conf in EnsembleStat.zip (also copied below), and let me know what needs to be modified ?
gefsv13EP5r2.tgz at the google folder contains the input files (GEFSv12 reforecasts, GEFSv13 reforecasts, AORC data): https://drive.google.com/drive/folders/1iVSAWBVb5Lsezf3A3dBZcFdDwx1Pt0EW

I use docker/container metplus:5.1-latest, and the commandline command used is run_metplus.py EnsembleStat.conf user.conf

I checked following websites but didn't find answers to my questions:
#1740
#1725
https://dtcenter.org/metplus-practical-session-guide-feb-2019/metplus-practical-session-guide-feb-2019/session-3-ensemble-pqpf/met-tool-ensemble-stat/ensemble-stat-tool-1
https://metplus.readthedocs.io/en/main_v5.1/Users_Guide/systemconfiguration.html
https://metplus.readthedocs.io/en/latest/generated/met_tool_wrapper/EnsembleStat/EnsembleStat.html#sphx-glr-generated-met-tool-wrapper-ensemblestat-ensemblestat-py

EnsembleStat.conf:
[config]

PROCESS_LIST = EnsembleStat

LOOP_BY = INIT
INIT_TIME_FMT = %Y%m%d%H
#INIT_BEG=2009123112
#INIT_END=2009123112
INIT_BEG=2018010400
INIT_END=2018022200

#INIT_INCREMENT=3600
INIT_INCREMENT=604800

#LEAD_SEQ = 24H
#z LEAD_SEQ = 6H
LEAD_SEQ = begin_end_incr(6,120,6)

#z FCST_ENSEMBLE_STAT_INPUT_DIR = {INPUT_BASE}/met_test/data/sample_fcst
#FCST_ENSEMBLE_STAT_INPUT_DIR = {INPUT_BASE}/sample_fcst
FCST_ENSEMBLE_STAT_INPUT_DIR = {INPUT_BASE}/gefsv12

#FCST_ENSEMBLE_STAT_INPUT_TEMPLATE = {init?fmt=%Y%m%d%H}/arw-???-gep?/d01_{init?fmt=%Y%m%d%H}_0{lead?fmt=%HH}00.grib

FCST_ENSEMBLE_STAT_INPUT_TEMPLATE = {init?fmt=%Y%m%d}/p25/ge???.t00z.pgrb2.0p25.f{lead?fmt=%3H}

#ENSEMBLE_STAT_CTRL_INPUT_DIR = {INPUT_BASE}/met_test/data/sample_fcst
#ENSEMBLE_STAT_CTRL_INPUT_TEMPLATE = {init?fmt=%Y%m%d%H}/arw-fer-gep1/d01_{init?fmt=%Y%m%d%H}_0{lead?fmt=%HH}00.grib

#z OBS_ENSEMBLE_STAT_POINT_INPUT_DIR = {INPUT_BASE}/ascii2nc
#z OBS_ENSEMBLE_STAT_POINT_INPUT_TEMPLATE = precip24_{valid?fmt=%Y%m%d%H}.nc

#z OBS_ENSEMBLE_STAT_GRID_INPUT_DIR = {INPUT_BASE}/met_test/data/sample_obs/ST4
#z OBS_ENSEMBLE_STAT_GRID_INPUT_DIR = {INPUT_BASE}/sample_obs/ST4
#z OBS_ENSEMBLE_STAT_GRID_INPUT_TEMPLATE = ST4.{valid?fmt=%Y%m%d%H}.24h

OBS_ENSEMBLE_STAT_GRID_INPUT_DIR = {INPUT_BASE}/aorc_regridded/{valid?fmt=%Y%m}
OBS_ENSEMBLE_STAT_GRID_INPUT_TEMPLATE = AORC-OWP_{valid?fmt=%Y%m%d%H}z.nc4

#ENSEMBLE_STAT_ENS_MEAN_INPUT_DIR =
#ENSEMBLE_STAT_ENS_MEAN_INPUT_TEMPLATE =

ENSEMBLE_STAT_CLIMO_MEAN_INPUT_DIR =
ENSEMBLE_STAT_CLIMO_MEAN_INPUT_TEMPLATE =

ENSEMBLE_STAT_CLIMO_STDEV_INPUT_DIR =
ENSEMBLE_STAT_CLIMO_STDEV_INPUT_TEMPLATE =

ENSEMBLE_STAT_OUTPUT_DIR = {OUTPUT_BASE}/ensemble
ENSEMBLE_STAT_OUTPUT_TEMPLATE = {init?fmt=%Y%m%d%H%M}/ensemble_stat

#z MODEL = WRF
MODEL = GEFS
OBTYPE = MC_PCP

FCST_VAR1_NAME = APCP
#z FCST_VAR1_LEVELS = A24
FCST_VAR1_LEVELS = Z0

FCST_VAR1_OPTIONS = ens_ssvar_bin_size = 0.1; ens_phist_bin_size = 0.05;

#z OBS_VAR1_NAME = {FCST_VAR1_NAME}
#z OBS_VAR1_LEVELS = {FCST_VAR1_LEVELS}

OBS_VAR1_NAME = APCP_surface
#z OBS_VAR1_LEVELS = (0,,)
OBS_VAR1_LEVELS = (,)

OBS_VAR1_OPTIONS = {FCST_VAR1_OPTIONS}

#LOG_ENSEMBLE_STAT_VERBOSITY = 2

#ENSEMBLE_STAT_N_MEMBERS = 6
ENSEMBLE_STAT_N_MEMBERS = 5

ENSEMBLE_STAT_CONFIG_FILE = {PARM_BASE}/met_config/EnsembleStatConfig_wrapped

ENSEMBLE_STAT_DESC = NA

ENSEMBLE_STAT_OBS_WINDOW_BEGIN = -5400
ENSEMBLE_STAT_OBS_WINDOW_END = 5400

ENSEMBLE_STAT_ENS_THRESH = 1.0

ENSEMBLE_STAT_VLD_THRESH = 1.0

ENSEMBLE_STAT_OUTPUT_PREFIX =

#ENSEMBLE_STAT_MET_OBS_ERR_TABLE =

ENSEMBLE_STAT_REGRID_TO_GRID = NONE
ENSEMBLE_STAT_REGRID_METHOD = NEAREST
ENSEMBLE_STAT_REGRID_WIDTH = 1
ENSEMBLE_STAT_REGRID_VLD_THRESH = 0.5
ENSEMBLE_STAT_REGRID_SHAPE = SQUARE
#ENSEMBLE_STAT_REGRID_CONVERT =
#ENSEMBLE_STAT_REGRID_CENSOR_THRESH =
#ENSEMBLE_STAT_REGRID_CENSOR_VAL =

ENSEMBLE_STAT_CENSOR_THRESH =
ENSEMBLE_STAT_CENSOR_VAL =

#ENSEMBLE_STAT_PROB_CAT_THRESH =
#ENSEMBLE_STAT_PROB_PCT_THRESH = ==0.25
#ENSEMBLE_STAT_ECLV_POINTS = 0.05

ENSEMBLE_STAT_MESSAGE_TYPE = ADPSFC
#ENSEMBLE_STAT_OBS_THRESH =
ENSEMBLE_STAT_DUPLICATE_FLAG = NONE
ENSEMBLE_STAT_SKIP_CONST = False

ENSEMBLE_STAT_OBS_ERROR_FLAG = FALSE

ENSEMBLE_STAT_ENS_SSVAR_BIN_SIZE = 1.0
ENSEMBLE_STAT_ENS_PHIST_BIN_SIZE = 0.05

#ENSEMBLE_STAT_CLIMO_MEAN_FILE_NAME =
#ENSEMBLE_STAT_CLIMO_MEAN_FIELD =
#ENSEMBLE_STAT_CLIMO_MEAN_REGRID_METHOD =
#ENSEMBLE_STAT_CLIMO_MEAN_REGRID_WIDTH =
#ENSEMBLE_STAT_CLIMO_MEAN_REGRID_VLD_THRESH =
#ENSEMBLE_STAT_CLIMO_MEAN_REGRID_SHAPE =
#ENSEMBLE_STAT_CLIMO_MEAN_TIME_INTERP_METHOD =
#ENSEMBLE_STAT_CLIMO_MEAN_MATCH_MONTH =
#ENSEMBLE_STAT_CLIMO_MEAN_DAY_INTERVAL = 31
#ENSEMBLE_STAT_CLIMO_MEAN_HOUR_INTERVAL = 6

#ENSEMBLE_STAT_CLIMO_STDEV_FILE_NAME =
#ENSEMBLE_STAT_CLIMO_STDEV_FIELD =
#ENSEMBLE_STAT_CLIMO_STDEV_REGRID_METHOD =
#ENSEMBLE_STAT_CLIMO_STDEV_REGRID_WIDTH =
#ENSEMBLE_STAT_CLIMO_STDEV_REGRID_VLD_THRESH =
#ENSEMBLE_STAT_CLIMO_STDEV_REGRID_SHAPE =
#ENSEMBLE_STAT_CLIMO_STDEV_TIME_INTERP_METHOD =
#ENSEMBLE_STAT_CLIMO_STDEV_MATCH_MONTH =
#ENSEMBLE_STAT_CLIMO_STDEV_DAY_INTERVAL = 31
#ENSEMBLE_STAT_CLIMO_STDEV_HOUR_INTERVAL = 6

ENSEMBLE_STAT_CLIMO_CDF_BINS = 1
ENSEMBLE_STAT_CLIMO_CDF_CENTER_BINS = False
ENSEMBLE_STAT_CLIMO_CDF_WRITE_BINS = True
#ENSEMBLE_STAT_CLIMO_CDF_DIRECT_PROB =

ENSEMBLE_STAT_MASK_GRID = FULL
ENSEMBLE_STAT_MASK_POLY =
MET_BASE/poly/HMT_masks/huc4_1605_poly.nc,
MET_BASE/poly/HMT_masks/huc4_1803_poly.nc,
MET_BASE/poly/HMT_masks/huc4_1804_poly.nc,
MET_BASE/poly/HMT_masks/huc4_1805_poly.nc,
MET_BASE/poly/HMT_masks/huc4_1806_poly.nc

ENSEMBLE_STAT_CI_ALPHA = 0.05

ENSEMBLE_STAT_INTERP_FIELD = BOTH
ENSEMBLE_STAT_INTERP_VLD_THRESH = 1.0
ENSEMBLE_STAT_INTERP_SHAPE = SQUARE
ENSEMBLE_STAT_INTERP_METHOD = NEAREST
ENSEMBLE_STAT_INTERP_WIDTH = 1

ENSEMBLE_STAT_OUTPUT_FLAG_ECNT = BOTH
ENSEMBLE_STAT_OUTPUT_FLAG_RPS = NONE
ENSEMBLE_STAT_OUTPUT_FLAG_RHIST = BOTH
ENSEMBLE_STAT_OUTPUT_FLAG_PHIST = BOTH
ENSEMBLE_STAT_OUTPUT_FLAG_ORANK = BOTH
ENSEMBLE_STAT_OUTPUT_FLAG_SSVAR = BOTH
ENSEMBLE_STAT_OUTPUT_FLAG_RELP = BOTH
#ENSEMBLE_STAT_OUTPUT_FLAG_PCT = BOTH
#ENSEMBLE_STAT_OUTPUT_FLAG_PSTD = BOTH
#ENSEMBLE_STAT_OUTPUT_FLAG_PJC = BOTH
#ENSEMBLE_STAT_OUTPUT_FLAG_PRC = BOTH
#ENSEMBLE_STAT_OUTPUT_FLAG_ECLV = BOTH

ENSEMBLE_STAT_NC_ORANK_FLAG_LATLON = TRUE
ENSEMBLE_STAT_NC_ORANK_FLAG_MEAN = TRUE
ENSEMBLE_STAT_NC_ORANK_FLAG_RAW = TRUE
ENSEMBLE_STAT_NC_ORANK_FLAG_RANK = TRUE
ENSEMBLE_STAT_NC_ORANK_FLAG_PIT = TRUE
ENSEMBLE_STAT_NC_ORANK_FLAG_VLD_COUNT = TRUE
ENSEMBLE_STAT_NC_ORANK_FLAG_WEIGHT = FALSE

#ENSEMBLE_STAT_OBS_QUALITY_INC =
#ENSEMBLE_STAT_OBS_QUALITY_EXC =

#ENSEMBLE_STAT_ENS_MEMBER_IDS =
#ENSEMBLE_STAT_CONTROL_ID =

#ENSEMBLE_STAT_GRID_WEIGHT_FLAG =

user.conf:
[dir]
INPUT_BASE=/test/EnsembleStat/gefsv13EP5r2/input
OUTPUT_BASE=/test/EnsembleStat/gefsv13EP5r2/output

[JohnHalleyGotway]

@haksulee I see you have a question about how to compute gridded ensemble statistics with METplus.

I want to produce the US map plots of gridded ensemble statistics (CRPS, CRPSS, etc) and ensemble mean statistics (correlation, etc) for different lead hours on a 6-hr and 1-day time step.

First, let me point out a couple of details:

• The Point-Stat, Grid-Stat, and Ensemble-Stat tools compute statistics by aggregating matched pairs over one or more spatial areas at a specific point in time.
• The Series-Analysis tool computes statistics separately for each grid point by aggregating matched pairs at individual grid points over a series of times.

Since you're looking for US map plots of statistics, Series-Analysis is the tool you'll want to run. However, Series-Analysis reads in a time series of gridded forecast data and a time series of gridded analysis data. Note though that it DOES NOT support reading a time series of ensemble inputs. In fact, this dtcenter/MET#770 issue recommends creating an Ensemble-Series-Analysis tool to compute ensemble statistics, like CRPS, but that tool does not exist yet.

So unfortunately, there is no easy/obvious way to create spatial of maps of CRPS and CRPSS with MET. I see two options to consider:

Option 1:
Stick with Ensemble-Stat and compute spatial summary statistics for each model output time. Then run Stat-Analysis to aggregate the ECNT line type (which contains CRPS) across multiple output times. You won't get a spatial map of CRPS but you'll get it for spatial masking region you define.

Option 2:
Run the Gen-Ens-Prod tool to generate ensemble products for the fields you specify, such as the ensemble mean and also derive threshold-based probabilities. Run Series-Analysis to read the Gen-Ens-Prod output and compare to the AORC analysis and compute statistics separately for each grid point. For the ensemble mean, you'd compute continuous statistics (CNT line type) and/or categorical statistics (CTC and CTS line types). For the ensemble probabilities, you'd compute probabilistic statistics (PCT, PSTD, PRC, and PJC line types).

How would you like to proceed? Running Ensemble-Stat/Stat-Analysis? Or running Gen-Ens-Prod/Series-Analysis?

Thanks,
John

[haksulee]

@JohnHalleyGotway
Thanks for the detailed explanation. I would like to proceed with Option 2 that generates gridded statistics.
Could you teach me how to do that?

[JohnHalleyGotway]

@haksulee, I downloaded all the data you provided... and that's a lot! I understand that you're primarily interested in accumulated precipitation, and ideally 24-hour accumulations. So really, you'd like advice setting up a use case to compare precipitation from two different GEFS ensembles to the AORC analysis data.

I see the following data:

1. gefsv12 ensemble forecast data in GRIB2 format.

• 8 different initializations from 20180104 to 2018022, every 7 days.
• Each initialization has data for two different resolutions: p25 and p5 for 1/4 and 1/2 lat/lon grid spacing.
• Each initialization has data for 1 control member (gec00) and 4 forecast members (gep1, gep2, gep3, gep4).
• The p25 members have data every 3 hours from forecast hour 3 to 240.
• The p5 members have data every 6 hours from forecast hour 246 to 384.
• The APCP variable contains 3 or 6 hour accumulations.

2. gefsv13EP5r2 ensemble forecast data in GRIB2 format.

• Same 8 initialization times.
• The resolution sub-directories are named p25 and sp25, containing six-hourly output from forecast hours 6 to 240.
• Each initialization has data for 1 control member (gec00) and 4 forecast members (gep1, gep2, gep3, gep4).

3. aorc_regridded contains hourly data in NetCDF format.

• Hourly output from 20180101 to 2018022823.
• Presumably, the APCP_surface(lat, lon) variable contains 1-hour accumulations of precip, even though that's not explicitly stated.

Looking in your existing EnsembleStat.conf file, I see that you're looping through the 8 init times, and for each, running Ensemble-Stat for each lead time from 6 to 120, every 6 hours. But since you're asking about 24-hour precip, I assume you'd rather be evaluating 24-hour precip.

Based on your input data and questions, I assume you'd like to do the following:

• Run Series-Analysis once for each 24-hour lead time (so 24, 48, 72, 96, 120) to compare gefsv12 and gefsv13EP5r2 to aorc.
• Assume that we'd only do this for the p25 data.
• For each run, pass as input ensemble products (i.e. probabilities and ensemble mean) for 24-hour precip and AORC analysis for 24-hour precip.
• So for each ensemble member, we'll need to run PCP-Combine to sum up the 6-hour precip into a 24-hour bucket. And then run Gen-Ens-Prod to generate the products.
• And we'll need to run PCP-Combine to sum up the 1-hour analysis precip into a 24-hour bucket.

So you'll get Series-Analysis output for 24, 48, 72, 96, and 120 hour lead times for the gefsv12 and gefsv13EP5r2 models. And you could compare them side-by-side to see the relatively performance. I understand that'd you'd prefer to treat gefsv12 as a "climatology" or "reference" dataset, but would recommend that you start by processing them separately.

Does this sound like what you'd like to do? Or are there any details to be clarified?

I just want to make sure I understand the goal clearly before diving into config file settings.

[haksulee]

@JohnHalleyGotway

Thank you for downloading and checking the data and suggestions:

The following may clarify some details:

1. "for each ensemble member, we'll need to run PCP-Combine to sum up the 6-hour precip into a 24-hour bucket"
   -> let me talk to my colleagues if this is ok.
   I am not sure if each ensemble member valid at different lead hours (e.g., the same initialization time but different valid time, for example, ensemble member 1 at f006 and ensemble member 1 at f012 from the same initialization time) are dynamically connected.

2. aorc_regridded: yes, APCP_surface(lat, lon) variable contains 1-hour accumulations of precip

3. gefsv12 ensemble forecast data : p25 data (ignore p5) will be used and 6 hour accumulations (ignore 3 hour accumulations) will be used.

4. gefsv13EP5r2 ensemble forecast data : sp25 (ignore p25 ) will be used because sp25 has APCP

5. Yes, I want to evaluate 24-hour precip for each of the first 5 days.
   I also want to evaluate a 5-day aggregation of precipitation for the first 5 days

6. Yes, I want to treat gefsv12 as a "reference" dataset

7. I would like to compute following statistics, if possible:
   Fractional bias (bias) = sum(GEFS_i - AORC_i) / sum(AORC_i), i=1,...,n; n is the number of pairs of (GEFS_i, AORC_i),
   Mean absolute error, MAE
   Mean square error, MSE
   Pearson correlation
   Spearman’s rank correlation
   MAE skill score
   MSE skill score
   The Brier Skill Score, BSS
   The ROC score

[haksulee]

@JohnHalleyGotway

Regarding "for each ensemble member, we'll need to run PCP-Combine to sum up the 6-hour precip into a 24-hour bucket", yes, please do so.

Thank you for your support.

[JohnHalleyGotway]

Met with @georgemccabe on Dec 10, 2024 to brainstorm.

• Loop over the 8 initialization times.
• For each, process lead times 24, 48, ..., 120.
• Recommend setting up a use case conf file to control the runtime logic, but have model-specific information in 2 separate conf files.
• Process list:
  1. PcpCombine(aorc) ... add 24 1-hour accumulations
  2. PcpCombine(gefs) ... add 4 6-hour accumulations with a custom loop list to loop over 5 members
  3. GenEnsProd ... to compute gefsv12 ensemble mean and probabilities
  4. SeriesAnalysis(by_lead) ... to compute stats for EACH lead time
  5. SeriesAnalysis(all_leads) ... to compute stats for ALL lead times

Steps for running:

cd /d1/projects/METplus/discussions/2809/METplus
ln -sf /nrit/ral/METplus-5.1.0/ush/run_metplus.py .
conda activate /home/met_test/.conda/envs/metplus_base.v5.1

... more details to come...

[haksulee]

@JohnHalleyGotway
Thanks for the update.
Please let me know if I have to do something.

[haksulee]

@JohnHalleyGotway

Could you continue supporting and updating this ticket ?
Please, let me know if you have any comments.

Thanks
Haksu