Custom Shift Detection Methods¶

So far TOAD only comes with ASDETECT [Boulton+Lenton2019]. But you can also implement your own shifts detection algorithm. You basically just need a function which transforms the timeseries in each grid cell into a time series of “abruptness”, with floating point values between -1 (negative shift) and 1 (positive shift).

# Prerequisites
import xarray as xr

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["figure.dpi"] = 300
plt.rcParams["figure.figsize"] = (12, 5)

data = (
    xr.open_dataset("test_data/garbe_2020_antarctica.nc")
    .coarsen(x=2, y=2, GMST=2, boundary="trim")
    .reduce(np.mean)
)

data = data.drop_vars("thk_dts")
data

/Users/jakobharteg/miniconda3/envs/toad312/lib/python3.12/site-packages/pyproj/network.py:59: UserWarning: pyproj unable to set PROJ database path.
  _set_context_ca_bundle_path(ca_bundle_path)

<xarray.Dataset> Size: 6MB
Dimensions:  (GMST: 175, y: 95, x: 95)
Coordinates:
  * GMST     (GMST) float64 1kB 0.0501 0.1301 0.2101 ... 13.81 13.89 13.97
  * y        (y) float64 760B -3.016e+06 -2.952e+06 ... 2.936e+06 3e+06
  * x        (x) float64 760B -3.016e+06 -2.952e+06 ... 2.936e+06 3e+06
Data variables:
    thk      (GMST, y, x) float32 6MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes: (12/13)
    CDI:            Climate Data Interface version 1.9.6 (http://mpimet.mpg.d...
    proj4:          +lon_0=0.0 +ellps=WGS84 +datum=WGS84 +lat_ts=-71.0 +proj=...
    CDO:            Climate Data Operators version 1.9.6 (http://mpimet.mpg.d...
    source:         PISM (development v1.0-535-gb3de48787 committed by Julius...
    institution:    PIK / Potsdam Institute for Climate Impact Research
    author:         Julius Garbe (julius.garbe@pik-potsdam.de)
    ...             ...
    title:          garbe_2020_antarctica
    Conventions:    CF-1.9
    projection:     Polar Stereographic South (71S,0E)
    ice_density:    910. kg m-3
    NCO:            netCDF Operators version 4.7.8 (Homepage = http://nco.sf....
    Modifications:  Modified by Jakob Harteg (jakob.harteg@pik-potsdam.de) Se...

Let’s build a simple shifts detection algorithm

def compute_shift_intensity(values_1d: np.ndarray, window_size=5):
    """
    Simple function to compute shift intensity for each timestep [-1, 1]
    where -1 is strong negative shift, +1 is strong positive shift
    """

    data = np.array(values_1d)
    window_size = max(3, len(data) // 15)  # Ensure reasonable window size

    # Initialize results array
    intensity = np.zeros(len(data))

    # Calculate max possible difference for normalization
    max_diff = np.max(data) - np.min(data)
    if max_diff == 0:
        return intensity

    # For each point (except edges), compare before/after windows
    for i in range(window_size, len(data) - window_size):
        window1 = data[i - window_size : i]
        window2 = data[i : i + window_size]

        mean1 = np.mean(window1)
        mean2 = np.mean(window2)

        # Normalize difference to [-1, 1]
        # Positive means increasing, negative means decreasing
        intensity[i] = (mean2 - mean1) / max_diff
        intensity[i] = min(intensity[i], 1)
        intensity[i] = max(intensity[i], -1)

    return intensity

So this function takes a single timeseries with values and timesteps

# get sample timeseries from dataset
ts = data.sel(x=-1032000, y=-1032000, method="nearest").thk

# Calculate shift intensity
intensity = compute_shift_intensity(ts.values)

# Plot original data and shift intensity
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 4), sharex=True)
ax1.plot(ts.GMST, ts.values)
ax1.set_title("Original Time Series")
ax2.plot(ts.GMST, intensity, color="red")
ax2.axhline(y=0, color="k", linestyle="--", alpha=0.3)
ax1.set_ylabel("Ice thickness")
ax2.set_ylabel("Shift Intensity")
ax2.set_xlabel("GMST")
ax2.set_ylim(-1, 1)

(-1.0, 1.0)

../_images/7785730c1047f96db8f11f865e6acc2d7d72de68543bf9d22ab91cf3414aa200.png

Now we can define a new shifts method by abstracting the ShiftsMethod class in toad.shifts.

from toad.shifts import ShiftsMethod


class MY_SHFITS_METHOD(ShiftsMethod):
    # You can pass parameters to the method like this
    def __init__(self, window_size=5):
        self.window_size = window_size

    # This function is called by TOAD on each gridcell independently.
    def fit_predict(
        self,
        values_1d: np.ndarray,
        times_1d: np.ndarray,
    ) -> np.ndarray:
        shifts = compute_shift_intensity(values_1d, window_size=self.window_size)
        return shifts

And then we can pass this method into compute_shifts.

from toad import TOAD

td = TOAD(data, time_dim="GMST")
td.compute_shifts("thk", method=MY_SHFITS_METHOD(window_size=10))

INFO: New shifts variable thk_dts: min/mean/max=-0.981/-0.057/0.937 using 3589 grid cells. Skipped 60.2% grid cells: 0 NaN, 5436 constant.

Let’s inspect the detection time series distribution. He were see we detect mostly negative shifts.

td.plot.shifts_distribution()

(<Figure size 3600x600 with 1 Axes>,
 array([<Axes: ylabel='#thk_dts'>], dtype=object))

../_images/5d658468c2d04746b29a272d1ffb74c55c2eb007b804785165e6605e6c4b4898.png

Let’s compute clusters now and plot the results

from sklearn.cluster import HDBSCAN

td.compute_clusters(
    "thk",
    method=HDBSCAN(min_cluster_size=10),
    overwrite=True,
)

# plot 5 largest clusters
td.plot.overview(
    var="thk",
    cluster_ids=range(5),
    map_style={"projection": "south_pole"},
);

INFO: New cluster variable thk_dts_cluster: Identified 12 clusters in 2,414 pts; Left 7.0% as noise (169 pts).

../_images/67fa8996e30765ebb422d91e87ab561a24f996c14b5f5b303164a1b2a150e5d9.png