Cumulative skillscores along trajectory are now calculated when providing cumulative=skill to skill method (e.g. Liu-Weisberg). example_parcels updated to illustrate

Author: knutfrode

examples/example_parcels.py

@@ -53,3 +53,15 @@
 plt.colorbar(mappable=mappable, orientation='horizontal', label='Skillscore per trajectory')
 plt.legend()
 plt.show()
+
+#%%
+# Illustration of cumulative skillscore:
+# Plotting another random trajectory ("modeled"), colored by the cumulative skillscore
+# where the first trajectory is regarded as truth
+ds_model = ds.isel(trajectory=4)
+skillscore = ds_model.traj.skill(expected=ds_true, method='liu-weissberg', cumulative=True)
+ds_true.traj.plot(color='k', linewidth=3, label='"True" trajectory', land='mask')
+mappable = ds_model.traj.plot(linewidth=3, color=skillscore, label='"Modeled" trajectory')
+plt.colorbar(mappable=mappable, orientation='horizontal', label='Cumulative skillscore')
+plt.legend()
+plt.show()

tests/test_skill_score.py

@@ -71,6 +71,78 @@ def test_skillscores():
     skill_lw = ta.skill.liu_weissberg(lon_obs, lat_obs, lon_model, lat_model)
     np.testing.assert_almost_equal(skill_lw, 0.99099, 5)
 
+def test_skillscores_cumulative():
+    lon_obs = np.array([0, 1, 2, 3, 4, 5], dtype=float)
+    lat_obs = np.array([0, 0, 0, 0, 0, 0], dtype=float)
+    lon_model = lon_obs.copy()
+    km2deg = 111
+    lon_model[-1] = lon_obs[-1] + 1.8/km2deg
+    lat_model = np.array([0, 1.2/km2deg, -3.4/km2deg, 6.3/km2deg, 4.2/km2deg, 0])
+
+    skill_cum = ta.skill.liu_weissberg(lon_obs, lat_obs, lon_model, lat_model, cumulative=True)
+
+    # Output length must match input length
+    assert len(skill_cum) == len(lon_obs)
+    # First value is NaN (no arc length at t=0); second value is always valid
+    assert np.isnan(skill_cum[0])
+    assert not np.isnan(skill_cum[1])
+    # All other values are in [0, 1]
+    assert np.all(skill_cum[1:] >= 0) and np.all(skill_cum[1:] <= 1)
+    # Last value equals the non-cumulative score
+    skill_scalar = ta.skill.liu_weissberg(lon_obs, lat_obs, lon_model, lat_model)
+    np.testing.assert_almost_equal(skill_cum[-1], skill_scalar, 10)
+    # Score decreases when divergence accelerates (d_k/L_k exceeds running average)
+    lon_div = np.array([0, 1, 2, 3, 4, 5], dtype=float)
+    lat_div = np.zeros(6)
+    lat_model_div = np.array([0, 0.5/111, 2.0/111, 5.0/111, 10.0/111, 18.0/111])  # super-linear divergence
+    skill_div = ta.skill.liu_weissberg(lon_div, lat_div, lon_div, lat_div + lat_model_div, cumulative=True)
+    assert skill_div[-1] < skill_div[-2], "Score should decrease for accelerating divergence"
+
+
+def test_skillscores_cumulative_2d():
+    lon_obs = np.array([[0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5]], dtype=float).T
+    lat_obs = np.zeros_like(lon_obs)
+    lon_model = lon_obs.copy()
+    km2deg = 111
+    lon_model[:, -1] = lon_obs[:, -1] + 1.8/km2deg
+    lat_model = np.array([[0, 1.2/km2deg, -3.4/km2deg, 6.3/km2deg, 4.2/km2deg, 0],
+                          [0, 1.2/km2deg, -3.4/km2deg, 6.3/km2deg, 4.2/km2deg, 0]]).T
+
+    skill_cum = ta.skill.liu_weissberg(lon_obs, lat_obs, lon_model, lat_model, cumulative=True)
+
+    assert skill_cum.shape == (6, 2)
+    assert np.all(np.isnan(skill_cum[0, :]))
+    assert not np.any(np.isnan(skill_cum[1, :]))
+    assert np.all(skill_cum[1:, :] >= 0) and np.all(skill_cum[1:, :] <= 1)
+    # Last value matches non-cumulative score for each trajectory
+    skill_scalar = ta.skill.liu_weissberg(lon_obs, lat_obs, lon_model, lat_model)
+    np.testing.assert_array_almost_equal(skill_cum[-1, :], skill_scalar)
+
+
+def test_skillscores_cumulative_xarray(barents):
+    barents = barents.traj.gridtime('1h')
+    b0 = barents.isel(trajectory=0).dropna('time')
+    b1 = barents.isel(trajectory=1).sel(time=slice(b0.time[0], b0.time[-1]))
+    b1 = b1.traj.gridtime(b0.time)
+
+    skill_cum = b0.traj.skill(b1, cumulative=True)
+
+    # Cumulative result has exactly the same dims and coords as self's internal dataset
+    assert skill_cum.dims == b0.traj.ds.lon.dims
+    assert skill_cum.sizes == {d: b0.traj.ds.sizes[d] for d in b0.traj.ds.lon.dims}
+    np.testing.assert_array_equal(skill_cum.coords['time'].values, b0.coords['time'].values)
+    assert np.all(np.isnan(skill_cum.isel(time=0).values))
+    assert not np.any(np.isnan(skill_cum.isel(time=1).values))
+    # Cumulative array values are in [0, 1] beyond t=0
+    assert np.all(skill_cum.isel(time=slice(1, None)).values >= 0)
+    assert np.all(skill_cum.isel(time=slice(1, None)).values <= 1)
+    # Last value equals the non-cumulative score
+    skill_scalar = b0.traj.skill(b1)
+    np.testing.assert_array_almost_equal(
+        skill_cum.isel(time=-1).values.squeeze(),
+        skill_scalar.values.squeeze(), 5)
+
+
 def test_skillscores_2d():
     lon_obs = np.array([[0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5]]).T
     lat_obs = np.array([[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]]).T

trajan/skill/__init__.py

@@ -22,13 +22,23 @@ def distance_along_trajectory(lon, lat):
 
     return distance
 
-def liu_weissberg(lon_obs, lat_obs, lon_model, lat_model, tolerance_threshold=1):
+def liu_weissberg(lon_obs, lat_obs, lon_model, lat_model, tolerance_threshold=1, cumulative=False):
     '''
-    Calculate skill score from normalized cumulative seperation distance. Liu and Weisberg 2011.
+    Calculate skill score from normalized cumulative separation distance. Liu and Weisberg 2011.
+
+    Parameters:
+        cumulative: If True, return an array of skill scores at each time step
+                    (shape same as input, time first axis), with NaN at t=0.
+                    At step k the score is max(0, 1 - s_k/n₀) where
+                    s_k = (Σ_{i=0}^{k} d_i) / (Σ_{i=1}^{k} L_i), i.e. the
+                    Liu-Weisberg formula evaluated on the data up to step k.
+                    The last value therefore equals the non-cumulative score.
+                    The score decreases at step k when d_k/L_k exceeds the
+                    running average s_{k-1} (accelerating divergence).
 
     Returns:
 
-        Skill score between 0. and 1.
+        Skill score between 0. and 1., or array thereof if cumulative=True.
     '''
 
     lon_obs = np.array(lon_obs)
@@ -38,12 +48,30 @@ def liu_weissberg(lon_obs, lat_obs, lon_model, lat_model, tolerance_threshold=1)
 
     d = distance_between_trajectories(lon_obs, lat_obs, lon_model, lat_model)
     l = distance_along_trajectory(lon_obs, lat_obs)
-    s = np.nansum(d, axis=0) / np.nansum(np.nancumsum(l, axis=0), axis=0)
 
-    if tolerance_threshold==0:
-        skillscore = 0
+    if cumulative:
+        # Running Liu-Weisberg: at step k, evaluate the formula on data [0..k].
+        # Numerator: Σ d_i, Denominator: Σ L_i (sum of cumulative arc lengths).
+        # At k=0 the denominator is 0 → NaN; at k=N-1 the value equals the
+        # non-cumulative score.  Score decreases when d_k/L_k > running average
+        # (accelerating divergence).
+        cum_d = np.nancumsum(d, axis=0)            # shape (N, ...)
+        cum_L = np.nancumsum(np.nancumsum(l, axis=0), axis=0)  # shape (N-1, ...)
+        # s[k] = Σ_{0..k+1} d / Σ_{1..k+1} L  →  shape (N-1, ...)
+        s = cum_d[1:] / cum_L
+        if tolerance_threshold == 0:
+            skillscore = np.where(np.isfinite(s), 0.0, np.nan)
+        else:
+            skillscore = np.maximum(0, 1 - s / tolerance_threshold)
+        # Prepend NaN for t=0 (denominator undefined, no arc length yet)
+        nan_row = np.full((1,) + skillscore.shape[1:], np.nan)
+        skillscore = np.concatenate([nan_row, skillscore], axis=0)
     else:
-        skillscore = np.maximum(0, 1 - s/tolerance_threshold)
+        s = np.nansum(d, axis=0) / np.nansum(np.nancumsum(l, axis=0), axis=0)
+        if tolerance_threshold == 0:
+            skillscore = 0
+        else:
+            skillscore = np.maximum(0, 1 - s / tolerance_threshold)
 
     return skillscore
 

trajan/traj1d.py

@@ -254,6 +254,20 @@ def skill(self, expected, method='liu-weissberg', **kwargs) -> xr.DataArray:
         s = skill_method(expected.traj.tlon, expected.traj.tlat, ds.traj.tlon,
                          ds.traj.tlat, **kwargs)
 
+        if kwargs.get('cumulative', False):
+            # Build DataArray with explicit dims from the (broadcasted) ds,
+            # then drop any size-1 dims that are not in self, and reorder to
+            # match self.ds.lon dimension order.
+            result = xr.DataArray(s, dims=ds.lon.dims, coords=ds.lon.coords,
+                                  name='Skillscore', attrs={'method': method})
+            self_dims = set(self.ds.lon.dims)
+            for dim in list(result.dims):
+                if dim not in self_dims and result.sizes[dim] == 1:
+                    result = result.squeeze(dim=dim, drop=True)
+            ordered = [d for d in self.ds.lon.dims if d in result.dims]
+            extra = [d for d in result.dims if d not in ordered]
+            return result.transpose(*(ordered + extra))
+
         newcoords = dict(ds.lon.sizes)
         newcoords.pop('time')
         for dim in newcoords: