Skill scores from traditional dynamical forecast models

Workshops, competitions and events AI Weather Quest
1

As promised, we are sharing the Ranked Probability Skill Scores (RPSSs) for dynamical forecast models. These scores will be updated weekly, throughout the training period, and the latest data is available in the attached file.

JJA25_DYNforecast_evaluation.xls.zip (32.9 KB)

Currently skill scores are provided for forecasts initialised up to week 12 (31st July 2025) of the testing JJA period.

Dynamical forecast data has been downloaded from the WMO lead centre for the sub-seasonal multi-model ensemble database . We are providing RPSSs for six dynamical models. For all models, except WASHINGTON, the reforecast period spans 2006 to 2016. Details regarding the selection of reforecasts and the configuration of each dynamical model can be found here. No bias correction has been performed as we are evaluating forecasted probabilistic values against each model climate.

The attached file is organised by forecast window, forecast variable and skill score type:

· fcwin1 = Week 3 forecast

· fcwin2 = Week 4 forecast

· tas , mslp, pr = Forecast variables

· period = period-aggregated RPSS.

· weekly = single weekly RPSS.

:desktop_computer: You can now evaluate your own forecast and compare it directly to dynamical forecast skill scores. Detailed guidance on forecast self-evaluation can be found in ReadTheDocs Python documentation. :magnifying_glass_tilted_left:

For precipitation, scores are masked where all historical quintiles are equal to zero. See the following forum post for more information.

:exclamation: If you encounter any issues or have questions regarding the self-evaluation process, feel free to reach out - we’re here to help!

2 Likes
6

Dear Joshua_Talib
When calculating probabilistic forecasts using ensemble prediction results the climatological threshold is derived from model hindcast.
Sub-seasonal forecast models - AI Weather Quest - ECMWF Confluence Wiki
in this page, Number of hindcast members per initialisation (full number used to compute climatology). for example EC, 11 (x11x2=242), the first “11” is from Rfc_size, the second “11” is from Hindcast years used, the “2” is from Hindcast date selection.
Models - S2S - ECMWF Confluence Wiki
in this page, HMCR Rfc_size is 1+10, but when calculating climatology only use 10. It is different from other model.

7

Hi @Peng_Lu,

You are correct to raise this error. I had made a mistake.

For HMCR I was only downloading nine perturbed reforecast members. This missed out a reforecast member.

I have rewritten the code so that 11 reforecasts (10 + 1) are downloaded every week. I will change the dynamical forecast scores (above) once all forecast initialisations have been downloaded.

Thanks for noticing the mistake. I’ve also updated Sub-seasonal forecast models - AI Weather Quest - ECMWF Confluence Wiki.

Kind regards,
Josh

8

Thank you for your detailed and clear explanation

9

Update on 8th August.

10

Dear Joshua_Talib

I try to calculate ECMWF s2s dynamical forecast RPSS of pr in the first week (initialised on 20250703), similar as yours, but not include control result. The RPSS is -0.23. I will share my code, can you give me some advice?

11 
from ecmwfapi import ECMWFDataServer

target_month = "07"
target_day = "03"

years = range(2016, 2005, -1)  # From 2016 to 2006
hdate_list = [f"{year}-{target_month}-{target_day}" for year in years]
hdate = "/".join(hdate_list)

server = ECMWFDataServer()

server.retrieve({
    "class": "s2",
    "dataset": "s2s",
    "date": "2025-07-03",
    "expver": "prod",
    "hdate": hdate,
    "levtype": "sfc",
    "model": "glob",
    "number": "1/2/3/4/5/6/7/8/9/10",
    "origin": "ecmf",
    "param": "228228",
    "step": "456/480/504/528/552/576/600",
    "stream": "enfh",
    "time": "00:00:00",
    "type": "pf",
    "target": "pr-0703-h"
})

from ecmwfapi import ECMWFDataServer

target_month = "07"
target_day = "01"

years = range(2016, 2005, -1)  # From 2016 to 2006
hdate_list = [f"{year}-{target_month}-{target_day}" for year in years]
hdate = "/".join(hdate_list)

server = ECMWFDataServer()

server.retrieve({
    "class": "s2",
    "dataset": "s2s",
    "date": "2025-07-01",
    "expver": "prod",
    "hdate": hdate,
    "levtype": "sfc",
    "model": "glob",
    "number": "1/2/3/4/5/6/7/8/9/10",
    "origin": "ecmf",
    "param": "228228",
    "step": "504/528/552/576/600/624/648",
    "stream": "enfh",
    "time": "00:00:00",
    "type": "pf",
    "target": "pr-0701-h"
})

from ecmwfapi import ECMWFDataServer

server = ECMWFDataServer()

server.retrieve({
    "class": "s2",
    "dataset": "s2s",
    "date": "2025-07-03",
    "expver": "prod",
    "levtype": "sfc",
    "model": "glob",
    "number": "1/2/3/4/5/6/7/8/9/10/11/12/13/14/15/16/17/18/19/20/21/22/23/24/25/26/27/28/29/30/31/32/33/34/35/36/37/38/39/40/41/42/43/44/45/46/47/48/49/50/51/52/53/54/55/56/57/58/59/60/61/62/63/64/65/66/67/68/69/70/71/72/73/74/75/76/77/78/79/80/81/82/83/84/85/86/87/88/89/90/91/92/93/94/95/96/97/98/99/100",
    "origin": "ecmf",
    "param": "228228",
    "step": "456/480/504/528/552/576/600",
    "stream": "enfo",
    "time": "00:00:00",
    "type": "pf",
    "target": "pr-0703-f"
})

import xarray as xr
import numpy as np
ds = xr.open_dataset('pr-0703-h', engine='cfgrib')
#tp[number:10, time:20, step:7, latitude:121, longitude:240]
pr_h_1 = ds['tp'].values
pr_h_2 = np.sum(pr_h_1,axis=2)
pr_h_3 = pr_h_2.reshape(-1,121,240)

ds = xr.open_dataset('pr-0701-h', engine='cfgrib')
#tp[number:10, time:20, step:7, latitude:121, longitude:240]
pr_h_4 = ds['tp'].values
pr_h_5 = np.sum(pr_h_4,axis=2)
pr_h_6 = pr_h_5.reshape(-1,121,240)

pr_h = np.concatenate((pr_h_3, pr_h_6), axis=0)

ds = xr.open_dataset('pr-0703-f', engine='cfgrib')
pr_f_1 = ds['tp'].values  
pr_f = np.sum(pr_f_1,axis=1)

from AI_WQ_package import retrieve_evaluation_data
obs = retrieve_evaluation_data.retrieve_weekly_obs('20250721','pr','xxx')
quintile_clim = retrieve_evaluation_data.retrieve_20yr_quintile_clim('20250721','pr','xxx')

ds = xr.open_dataset('pr_obs_WEEKLYSUM_20250721.nc')
obs_0721 = ds['pr'].values
obs_721 = ds['pr']

ds = xr.open_dataset('pr_20yrCLIM_WEEKLYSUM_quintiles_20250721.nc')
q_0721 = ds['pr'][0,:,:,:].values
q_721 = ds['pr'][0,:,:,:]

# quantile based hintcast
quantiles = [0.2, 0.4, 0.6, 0.8]
q_h = np.quantile(pr_h, q=quantiles, axis=0) 

def calculate_quantile_probability(var_q, var_f):
    bins = np.empty((6, *var_q.shape[1:]))
    bins[0] = -np.inf
    bins[1:-1] = var_q
    bins[-1] = np.inf
    
    # var_f: (n, 1, lat, lon) vs bins: (1, 6, lat, lon)
    in_bin = (var_f[:, np.newaxis] >= bins[np.newaxis, :-1]) & \
             (var_f[:, np.newaxis] < bins[np.newaxis, 1:])
    
    counts = np.sum(in_bin, axis=0)
    
    probs = counts / var_f.shape[0]    
    return counts, probs

counts, probs = calculate_quantile_probability(q_h, pr_f)

from AI_WQ_package import forecast_evaluation
obs_pbs = forecast_evaluation.conditional_obs_probs(obs_721,q_721)
from AI_WQ_package import retrieve_evaluation_data
land_sea_mask = retrieve_evaluation_data.retrieve_land_sea_mask('xxx')
fc_pbs = obs_pbs.copy()
fc_pbs.values = probs
global_RPSS = forecast_evaluation.work_out_RPSS(fc_pbs,obs_pbs,'pr',land_sea_mask)
12

Hi Peng,

Thanks for the code. Good effort.

First quick reply: forecast precipitation is stored as accumulations from forecast lead time so you need to subtract end of day 25 from start of day 19 (i.e. 600 - 456). At the moment, you are summing across all daily time steps for forecast & reforecast. This is different to what is done with the observations.

Second check: have you updated the python package? I made some slight modifications in light of zero precip regions: Updates to Python Package - AI Weather Quest - ECMWF Confluence Wiki

Kind regards,

Josh

1 Like
13

Thank you for your reply.

I change the code and updated the python package, and now the RPSS is -0.0038. I think it is similar with your result and it is reasonable. Thank you very much

14

Much closer. I guess you’re missing the control forecast and reforecast.

Good work. It’s not an easy thing to code.

Josh