crps_hersbach

hydrostats.ens_metrics.crps_hersbach(obs, fcst_ens, remove_neg=False, remove_zero=False)

Calculate the the continuous ranked probability score (CRPS) as per equation 25-27 in Hersbach et al. (2000).

It is strongly recommended to use the hydrostats.ens_metric.ens_crps() function for the improved performance, instead of using this particular implementation. This is meant more for a proof of concept of the algorithm presented in the literature.

Parameters:
  • obs (1D ndarry) – Array of observations for each start date
  • fcst_ens (2D ndarray) – Array of ensemble forecast of dimension n x M, where n = number of start dates and M = number of ensemble members.
  • remove_neg (bool) – If True, when a negative value is found at the i-th position in the observed OR ensemble array, the i-th value of the observed AND ensemble array are removed before the computation.
  • remove_zero (bool) – If true, when a zero value is found at the i-th position in the observed OR ensemble array, the i-th value of the observed AND ensemble array are removed before the computation.
Returns:

Dictionary contains a number of experimental outputs including:
  • [“crps”] 1D ndarray of crps values per n start dates.
  • [“crpsMean1”] arithmetic mean of crps values.
  • [“crpsMean2”] mean crps using eqn. 28 in Hersbach (2000).
Return type:

dict

Notes

NaN and inf treatment: If any value in obs or fcst_ens is NaN or inf, then the corresponding row in both fcst_ens (for all ensemble members) and in obs will be deleted.

References

  • Hersbach, H. (2000) Decomposition of the Continuous Ranked Porbability Score for Ensemble Prediction Systems, Weather and Forecasting, 15, 559-570.

Examples

>>> import numpy as np
>>> import hydrostats.ens_metrics as em
>>> np.random.seed(3849590438)

Creating an observed 1D array and an ensemble 2D array with all random numbers

>>> ens_array_random = (np.random.rand(15, 52) + 1) * 100
>>> obs_array_random = (np.random.rand(15) + 1) * 100

Creating an observed 1D array and an ensemble 2D array with noise.

>>> noise = np.random.normal(scale=1, size=(15, 52))
>>> x = np.linspace(1, 10, 15)
>>> observed_array = np.sin(x) + 10
>>> ensemble_array_noise = (np.ones((15, 52)).T * observed_array).T + noise

Computing the Hersbach CRPS values between the ensemble mean and the observed data with the random data.

>>> crps_dictionary_rand = em.crps_hersbach(obs_array_random, ens_array_random)
>>> print(crps_dictionary_rand['crps'])
[ 7.73360237  9.59248626 34.46719655 30.10271075  7.451665   16.07882352
 14.59543529  8.55181637 15.4833089   8.32422363 16.55108154 19.20821296
  8.39452279 12.59949378 27.82543302]
>>> crps_dictionary_rand['crpsMean1']
15.797334183277723
>>> crps_dictionary_rand['crpsMean2']
15.797334183277725

Computing the Hersbach CRPS values between the ensemble mean and the observed data with noise in the ensemble data.

>>> crps_dictionary_noise = em.crps_hersbach(obs=observed_array, fcst_ens=ensemble_array_noise)
>>> print(crps_dictionary_noise['crps'])
[0.26921152 0.21388687 0.24927151 0.26047667 0.30234843 0.1996493
 0.2779844  0.29478927 0.275383   0.25682693 0.21485236 0.22824711
 0.2813889  0.21264652 0.18141063]
>>> crps_dictionary_noise['crpsMean1']
0.24789156041214705
>>> crps_dictionary_noise['crpsMean2']
0.24789156041214705