CLUSTERING METHODS FOR LAGRANGIAN ANALYSIS#

What do we mean by clustering?#

In ocean Lagrangian analysis, clustering methods are used to group particle or drifter trajectories that share similar characteristics in space, time, or dynamical behaviour. The main question we try to answer is:

How do we classify a spatiotemporal distribution of trajectories?

The choice of method depends on what your particles represent and whether you are interested in spatial aggregation (where particles cluster) or path similarity (how trajectories evolve). Several studies illustrate the diversity of clustering approaches in the field of geophysical fluid dynamics:

From simulated trajectories

Detection of Algulhas rings in the Southern Ocean Wichmann, D., et al. (2021)
Identification of the main pathways of the Labrador Current Jutras, M., Planat, N., Dufour, C. O., & Talbot, L. C. (2024)
Classification of dispersion regimes from clustering of autocorrelation velocities Lagomarsino-Oneto, D., De Leo, A., Stocchino, A., & Cucco, A. (2024).
Theoretical study of inertial particle clustering dynamics in idealised flow Bandopadhyay, T., Tomas-Gil, A., & Villafañe, L. (2025)

From drifter data

Estimation of eddy diffusivities from clustering of surface drifter data I. Koszalka, et al. (2011)
Analysis of the inertial effects of different drifter shapes/size Miron, P. et al. (2020)

Additionally, there are a wide variety of designed clustering techniques (see Rasyid, L. A., and Andayani, S., (2018) for a review). In this tutorial, we will be comparing the results from these algorithms:

K-Means: a spatial partitioning method that minimises within-cluster variance (squared Euclidean distance) to create K clusters
DBSCAN (Sander, J., et al. 1998): a density-based algorithm that groups points that are closely packed (points with many nearby neighbours), and marks as outliers points that lie alone in low-density regions.
OPTICS (Ankerst, M., et al. 1999): extension of DBSCAN algorithm that can handle clusters of varying densities.

I. Selection of the snapshot 🔴#

[ ]:

time_index=... #indicate time index
time= dataset.time.isel(obs=time_index).values[0]
ds=dataset.isel(obs=time_index)

II. Define & apply clustering algorithms 🔴#

💡 Look-up at the webiste of Scikit learn package to understand the physical meaning of the eps, min_samples, and max_eps parameters

[90]:

clustering_dict={
    'K-means': {'algorithm': KMeans(n_clusters=10)},
    'DBSCAN':  {'algorithm': DBSCAN(eps=np.radians(1), min_samples=5,  metric='haversine', algorithm='ball_tree')},
    'OPTICS':  {'algorithm': OPTICS(min_samples= 10, metric='haversine', algorithm='ball_tree', max_eps=np.radians(2))},
}

💡 Go through the cluster function & ensure you understand all the code

[91]:

def cluster(ds, algorithm, algorithm_name):
    """
    Perform clustering on particle trajectories using a specified clustering algorithm.

    Parameters
    ----------
    ds : xarray.Dataset
        dataset containing 'lon', 'lat', and 'traj' (trajectory index) variables.
    algorithm : clustering object
initialised clustering algorithm (e.g., KMeans, DBSCAN, OPTICS) from scikit-learn.
    algorithm_name : str
        name of the algorithm ('KMeans', 'DBSCAN', 'OPTICS') to handle algorithm-specific outputs.

    Returns
    -------
    df_cluster : pd.DataFrame
        DataFrame with original trajectory points and clustering results including:
        - cluster_label: cluster assignment (-1 for noise in DBSCAN/OPTICS)
        - is_noise: boolean flag for noise points
        - reachability: OPTICS reachability distance or zeros for other algorithms
        - ordering: OPTICS ordering of points or zeros for others
        - core_distances: OPTICS core distances or zeros for others
    """

    # convert longitude and latitude to radians and stack into Nx2 array for clustering
    points = np.column_stack([np.radians(ds.lon.values), np.radians(ds.lat.values)])

    # fit the clustering algorithm to the points
    clustering = algorithm.fit(points)

    # get cluster labels assigned by the algorithm
    labels = clustering.labels_

    # ensure noise points are labelled as -1, otherwise keep the cluster label
    unique_labels = np.where(labels == -1, -1, labels)

    # algorithm-specific handling for OPTICS
    if algorithm_name == 'OPTICS':
        # initialize reachability array with NaNs
        reachability_values = np.full(len(ds.lon.values), np.nan)
        # fill reachability in the order determined by OPTICS
        reachability_values[clustering.ordering_] = clustering.reachability_[clustering.ordering_]

        # initialise ordering array with -1
        ordering_values = np.full(len(ds.lon.values), -1)
        # store the position in the OPTICS ordering for each original point
        for order_pos, orig_idx in enumerate(clustering.ordering_):
            ordering_values[orig_idx] = order_pos

        # get core distances from OPTICS
        core_values = clustering.core_distances_

    else:
        # for non-OPTICS algorithms, reachability, ordering, and core distances are not defined
        reachability_values = np.full_like(unique_labels, 0)
        ordering_values = np.full_like(unique_labels, 0)
        core_values = np.full_like(unique_labels, 0)

    # Create a DataFrame with clustering results and original trajectory info
    df_cluster = pd.DataFrame({
        'traj_index': ds.traj.values,
        'lon': ds.lon.values,
        'lat': ds.lat.values,
        'cluster_label': unique_labels,
        'is_noise': (labels == -1),
        'reachability': reachability_values,
        'ordering': ordering_values,
        'core_distances': core_values,
    })

    return df_cluster

[92]:

#execute clustering for all algorithms
for clustering_method in tqdm(clustering_dict):
    algorithm= clustering_dict[clustering_method]['algorithm']
    clustering_dict[clustering_method]['result']= cluster(ds, algorithm, clustering_method)

100%|██████████| 3/3 [00:00<00:00, 80.48it/s]

III. Visualise results 🔴#

💡 Make a plot with the following characteristics:

Noise points appear in grey
Clustered points are colour-coded by cluster label
Includes land mask & coastline

🚀 CHALLENGE: Include the bathymetry of the region in the background!

[ ]:

# for clustering_method in clustering_dict:

🚀 CHALLENGE: Try to calculate the centre of mass of each cluster & add it to the plots above!

💡 What to do now?

Experiment with changing the parameters defining the algorithms & observe the differences
Change the time index & assess the stability of the clusters in time

IV. Statistical analysis of clusters 🔴#

In order to compare the results from different clustering algorithms, we can focus on the following statistical properties:

Number of clusters
Number of points per cluster
Within-cluster concentration [points/\(km^2\)]

[110]:

def calculate_number_clusters(clustering_dict):
    number_clusters={}
    for clustering_method in clustering_dict:
        number_clusters[clustering_method]=len(np.unique(clustering_dict[clustering_method]['result']['cluster_label']))
    return number_clusters

calculate_number_clusters(clustering_dict)

[110]:

{'K-means': 10, 'DBSCAN': 9, 'OPTICS': 6}

🚀 Your turn! Create functions to calculate the number of points per cluster for each algorithm & calculate the concentration of each cluster

Hint: to calculate the area covered by points, use scipy.spatial.ConvexHull

[ ]:

#def cardinality_cluster(clustering_dict):

#def within_cluster_concentration(clustering_dict):

🚀 Choose two options to visualise your results of points per cluster and concentration of each cluster:

Make a box plot with the mean & std property per algorithm
Make a histogram per cluster with the distribution points per cluster/ concentration

Still looking for a challenge ? 🔴#

Try using this other dataset: ‘toy_data_02.zarr’

Hint: make changes to the hyperparameters that define the algorithms
Take a look at other clustering algorithms included in Scikit learn package & try experimenting with other algorithms!