Use mean RPS for RPSS

notebooks/scripts.py

@@ -235,27 +235,24 @@ def skill_by_year(preds, adapt=False):
         # check inputs
         assert_predictions_2020(obs_p)
         assert_predictions_2020(fct_p)
-    
-    ## RPSS
+        
     # rps_ML
     rps_ML = xs.rps(obs_p, fct_p, category_edges=None, dim=[], input_distributions='p').compute()
     # rps_clim
     rps_clim = xs.rps(obs_p, clim_p, category_edges=None, dim=[], input_distributions='p').compute()
-    
-    # rpss
-    rpss = 1 - (rps_ML / rps_clim)
-    
-    # https://renkulab.io/gitlab/aaron.spring/s2s-ai-challenge-template/-/issues/7
 
-    # penalize
-    penalize = obs_p.where(fct_p!=1, other=-10).mean('category')
-    rpss = rpss.where(penalize!=0, other=-10)
+    ## RPSS
+    # penalize # https://renkulab.io/gitlab/aaron.spring/s2s-ai-challenge-template/-/issues/7
+    expect = obs_p.sum('category')
+    expect = expect.where(expect > 0.98).where(expect < 1.02)  # should be True if not all NaN
+
+    # https://renkulab.io/gitlab/aaron.spring/s2s-ai-challenge-template/-/issues/50
+    rps_ML = rps_ML.where(expect, other=2)  # assign RPS=2 where value was expected but NaN found
 
+    # following Weigel 2007: https://doi.org/10.1175/MWR3280.1
+    rpss = 1 - (rps_ML.groupby('forecast_time.year').mean() / rps_clim.groupby('forecast_time.year').mean())
     # clip
     rpss = rpss.clip(-10, 1)
-
-    # average over all forecasts
-    rpss = rpss.groupby('forecast_time.year').mean()
     
     # weighted area mean
     weights = np.cos(np.deg2rad(np.abs(rpss.latitude)))