diff --git a/examples/example_parcels.py b/examples/example_parcels.py index 8b5a8fb..2638ec5 100644 --- a/examples/example_parcels.py +++ b/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() diff --git a/tests/test_skill_score.py b/tests/test_skill_score.py index 7e9c94a..7754dec 100644 --- a/tests/test_skill_score.py +++ b/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 diff --git a/trajan/skill/__init__.py b/trajan/skill/__init__.py index 812cc56..76d9419 100644 --- a/trajan/skill/__init__.py +++ b/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 diff --git a/trajan/traj1d.py b/trajan/traj1d.py index 95c6c9f..795c686 100644 --- a/trajan/traj1d.py +++ b/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: