Movement Similarity & Clustering

Movement similarity and clustering is the practice of quantifying how alike two trajectories are with a distance measure such as discrete Fréchet, dynamic time warping, Hausdorff, or the longest common subsequence, then grouping many trajectories into route families by clustering the resulting pairwise distance matrix.

Prerequisites Checklist

Before building a similarity and clustering workflow, confirm your environment and inputs. This topic sits inside Movement Pattern Extraction & Trajectory Analysis and assumes each trajectory has already been cleaned and split into coherent trips.

• Python ≥ 3.10 with numpy ≥ 1.24 and pandas ≥ 2.0 for vectorized array work.
• scipy ≥ 1.11 for scipy.spatial.distance.directed_hausdorff and scipy.cluster.hierarchy agglomerative linkage.
• scikit-learn ≥ 1.3 for DBSCAN(metric="precomputed") and AgglomerativeClustering(metric="precomputed").
• pyproj ≥ 3.6 for projecting from EPSG:4326 to a metric CRS — import transformers with always_xy=True.
• Optional accelerators: dtaidistance or tslearn provide C-backed DTW with warping-band support if the pure-Python inner loop is too slow.
• Input contract: each trajectory is an ordered array of (x, y) points in a metric projected CRS, one array per trip. Trips are already segmented — a single array must not span two unrelated journeys.
• Projection context: trajectories usually arrive as EPSG:4326. Every similarity computation below assumes projected metres; the transform is flagged in each code block.

Similarity & Clustering Problem Taxonomy

Trajectory similarity looks deceptively simple — sum some point distances — but four distinct problems corrupt naive implementations. Diagnosing which one you face determines the measure, the pre-processing, and the compute budget. Feeding unsegmented or unprojected input into a distance measure produces clusters that reflect artifacts rather than genuine movement patterns.

Deterministic Pipeline Overview

The five stages below are the standard order. Each stage is fully vectorized where possible, records the parameters used, and preserves the mapping from matrix row index back to the source trajectory ID so cluster labels can be joined back to trips.

1. Project traces. Transform every trajectory from EPSG:4326 to a shared metric CRS (a UTM zone covering the study area). All subsequent distances are in metres.
2. Resample & align. Resample each trajectory to uniform arc-length spacing so sampling-rate differences and point-count mismatch do not dominate the distance. Optionally normalise the origin if you want shape-only matching.
3. Build the pairwise distance matrix. Compute a symmetric N×N matrix with the chosen measure. Only the upper triangle is computed and mirrored; the diagonal is zero.
4. Cluster the matrix. Run DBSCAN(metric="precomputed") or agglomerative linkage directly on the distance matrix. No re-derivation of features — the matrix is the input.
5. Extract a representative. For each detected cluster, select the medoid: the trajectory whose summed distance to the other cluster members is minimal. The medoid is a real observed path, unlike a synthetic average.

Implementation Walkthrough

The function below builds a pairwise DTW distance matrix on projected coordinates (a metric CRS) and clusters it with DBSCAN. DTW is chosen because it absorbs sampling-rate differences; swap in Fréchet or Hausdorff by replacing the inner _dtw_distance call. The distance matrix is precomputed, so DBSCAN’s eps is a real distance in metres. For the full derivation of the DTW and Fréchet inner loops, see measuring trajectory similarity with Fréchet and DTW.

from __future__ import annotations

import logging
from dataclasses import dataclass

import numpy as np
from pyproj import Transformer
from sklearn.cluster import DBSCAN

logger = logging.getLogger(__name__)


@dataclass
class ClusterResult:
    labels: np.ndarray          # cluster id per trajectory; -1 = noise
    distance_matrix: np.ndarray  # symmetric N x N, metres
    medoids: dict[int, int]     # cluster id -> index of representative trajectory


def _resample_arclength(track: np.ndarray, n_points: int) -> np.ndarray:
    """Resample an (m, 2) metric-CRS track to n_points of uniform arc-length spacing."""
    if track.shape[0] < 2:
        return np.repeat(track, n_points, axis=0)[:n_points]
    seg = np.sqrt(((np.diff(track, axis=0)) ** 2).sum(axis=1))
    cum = np.concatenate([[0.0], np.cumsum(seg)])
    total = cum[-1]
    if total == 0.0:                      # degenerate: all points identical
        return np.repeat(track[:1], n_points, axis=0)
    targets = np.linspace(0.0, total, n_points)
    x = np.interp(targets, cum, track[:, 0])
    y = np.interp(targets, cum, track[:, 1])
    return np.column_stack([x, y])


def _dtw_distance(a: np.ndarray, b: np.ndarray, band: int) -> float:
    """
    Banded dynamic time warping between two (m, 2) and (n, 2) metric-CRS tracks.
    A Sakoe-Chiba band of half-width `band` caps the inner cost at O(m * band).
    """
    m, n = len(a), len(b)
    inf = np.inf
    prev = np.full(n + 1, inf)
    prev[0] = 0.0
    for i in range(1, m + 1):
        curr = np.full(n + 1, inf)
        j_lo = max(1, i - band)
        j_hi = min(n, i + band)
        for j in range(j_lo, j_hi + 1):
            cost = np.hypot(a[i - 1, 0] - b[j - 1, 0], a[i - 1, 1] - b[j - 1, 1])
            curr[j] = cost + min(prev[j], curr[j - 1], prev[j - 1])
        prev = curr
    return float(prev[n])


def cluster_trajectories(
    tracks_wgs84: list[np.ndarray],
    target_crs_epsg: int = 32633,   # UTM zone 33N — adjust to your region
    n_resample: int = 64,
    dtw_band: int = 8,
    eps_m: float = 150.0,
    min_samples: int = 3,
) -> ClusterResult:
    """
    Cluster trajectories by DTW shape similarity.

    Parameters
    ----------
    tracks_wgs84 : list of (m, 2) arrays of (lon, lat) in EPSG:4326, one per trip.
    target_crs_epsg : metric projected CRS for all distance work. NEVER 4326.
    n_resample : uniform arc-length point count per track.
    dtw_band : Sakoe-Chiba half-width bounding the warp; caps inner cost.
    eps_m : DBSCAN neighbourhood radius in METRES (matrix is precomputed).
    min_samples : DBSCAN core-point threshold.

    Returns
    -------
    ClusterResult with labels, the distance matrix, and per-cluster medoid indices.

    Raises
    ------
    ValueError if the input list is empty or contains a non-2D array.
    """
    if not tracks_wgs84:
        raise ValueError("tracks_wgs84 is empty; nothing to cluster.")
    for k, t in enumerate(tracks_wgs84):
        if t.ndim != 2 or t.shape[1] != 2:
            raise ValueError(f"track {k} must be an (m, 2) array, got shape {t.shape}")

    # ── Stage 1: project every trace to a metric CRS ──────────────────────────
    # DTW sums point-to-point distances; those are only metric in projected space.
    to_metric = Transformer.from_crs("EPSG:4326", f"EPSG:{target_crs_epsg}", always_xy=True)
    projected: list[np.ndarray] = []
    for t in tracks_wgs84:
        x, y = to_metric.transform(t[:, 0], t[:, 1])   # always_xy: lon, lat order
        projected.append(np.column_stack([x, y]))

    # ── Stage 2: resample to uniform arc-length spacing ───────────────────────
    resampled = [_resample_arclength(t, n_resample) for t in projected]

    # ── Stage 3: pairwise DTW distance matrix (upper triangle only) ───────────
    n = len(resampled)
    dm = np.zeros((n, n), dtype=float)
    for i in range(n):
        for j in range(i + 1, n):
            d = _dtw_distance(resampled[i], resampled[j], band=dtw_band)
            dm[i, j] = d
            dm[j, i] = d               # enforce symmetry explicitly
    logger.info("Built %d x %d DTW distance matrix", n, n)

    # ── Stage 4: cluster the precomputed matrix ───────────────────────────────
    labels = DBSCAN(eps=eps_m, min_samples=min_samples, metric="precomputed").fit_predict(dm)

    # ── Stage 5: medoid per cluster (min summed distance within cluster) ──────
    medoids: dict[int, int] = {}
    for cid in sorted(set(labels)):
        if cid == -1:
            continue                   # skip the noise label
        members = np.where(labels == cid)[0]
        sub = dm[np.ix_(members, members)]
        medoids[int(cid)] = int(members[sub.sum(axis=1).argmin()])

    n_clusters = len({c for c in labels if c != -1})
    n_noise = int((labels == -1).sum())
    logger.info("DBSCAN: %d clusters, %d noise trajectories", n_clusters, n_noise)
    return ClusterResult(labels=labels, distance_matrix=dm, medoids=medoids)

The medoid is deliberately preferred over a coordinate-wise average. Averaging two trajectories that diverge around a roundabout produces a physically impossible path that cuts through the centre island; the medoid is always a genuine observed trip, which makes it safe to render on a map or feed back into route optimisation.

Geometric and Mathematical Grounding

The four common measures differ in how they pair points between two ordered point sequences P = (p₁, …, pₘ) and Q = (q₁, …, qₙ). Understanding the pairing rule tells you exactly when each is appropriate.

Discrete Fréchet distance. Often described with the dog-and-leash analogy: a person walks the vertices of P while a dog walks the vertices of Q, both monotonically forward, never backtracking. The Fréchet distance is the minimum leash length that permits both to reach the end. Formally it minimises, over all monotone couplings C, the maximum coupled distance:

d_F(P, Q) = min over C of max over (i, j) in C of ‖pᵢ − qⱼ‖

Because it takes a maximum, Fréchet is dominated by the single worst-aligned pair, which makes it sensitive to a lone excursion but faithful to overall path shape and order. Use it when continuity and traversal order matter, such as verifying that two driving routes truly follow the same corridor.

Dynamic time warping (DTW). DTW also allows a monotone alignment but minimises the sum of coupled distances instead of the maximum:

d_DTW(P, Q) = min over warping paths W of Σ over (i, j) in W of ‖pᵢ − qⱼ‖

The warping path may stay on the same index of one sequence while advancing the other, which is exactly what absorbs sampling-rate and speed differences. Because it sums, DTW is more forgiving of a single outlier than Fréchet but can drift if the warp is left unbounded — hence the Sakoe-Chiba band that restricts |i − j| ≤ band.

Hausdorff distance. The directed Hausdorff distance from P to Q is the greatest distance from any point of P to its nearest point in Q; the symmetric Hausdorff takes the larger of the two directed values:

d_H(P, Q) = max( sup over p in P of inf over q in Q of ‖p − q‖ , sup over q in Q of inf over p in P of ‖p − q‖ )

Hausdorff is a set measure — it ignores ordering entirely. That makes it ideal for comparing spatial coverage (do two trajectories occupy the same streets?) but useless when direction matters, because a route and its exact reverse have Hausdorff distance zero.

Longest common subsequence (LCSS). LCSS counts matched points that fall within a spatial tolerance ε, allowing unmatched points to be skipped on both sides. It converts to a distance via 1 − LCSS/min(m, n). Skipping makes LCSS the most robust to noise and to short spurious detours, at the cost of a tolerance parameter that must be tuned. It is the natural choice when trajectories contain GPS gaps or dropouts.

In practice DTW is the pragmatic default for corridor and commute work, Fréchet for strict route-conformance checks, Hausdorff for coverage comparisons, and LCSS for noisy or gappy data.

Calibration & Parameter Tuning

The right measure and threshold depend entirely on the question being asked. Because the distance matrix is precomputed in metres, DBSCAN’s eps is a literal metric radius you can reason about physically.

Calibrate eps empirically rather than guessing. Sort the k-distance graph (each trajectory’s distance to its min_samples-th nearest neighbour in the matrix) and place eps at the visible knee. Recheck the value whenever the geographic scale changes: a threshold tuned on inner-city delivery trips will over-merge long inter-city hauls.

For multi-modal fleets, cluster within a transport mode rather than across modes. A cyclist corridor and a vehicle corridor may overlap geometrically but represent different behaviour; segment the input by mode inferred from speed and acceleration profiling before building the matrix.

Integration & Compatibility

This topic consumes segmented, cleaned trajectories and produces cluster labels plus representative paths that feed several downstream stages:

• Trajectory Segmentation: Similarity is only meaningful on coherent single-trip segments. Feed the segment boundaries from segmentation in as the unit of comparison — never compare a full multi-day log as one array.
• Stay-Point Detection Algorithms: Stops recovered by stay-point detection let you trim dwell periods before clustering, so the distance reflects the moving corridor rather than time parked at a depot.
• Trajectory Object Design Patterns: Wrap each input trajectory in a consistent object so the matrix row index maps cleanly back to a trip ID and its assigned cluster label can be attached as an attribute.
• Measuring Trajectory Similarity with Fréchet and DTW: The dedicated deep-dive gives the complete dynamic-programming implementations of the two ordered measures used above.
• Directionality & Turn Analysis: Direction-aware measures (Fréchet, DTW) pair naturally with turn features when you need to separate outbound and inbound legs of the same corridor.

FAQ

Which trajectory similarity measure should I use — Fréchet, DTW, or Hausdorff?

Use discrete Fréchet when the ordering and continuity of the path matter, such as comparing two driving routes point by point. Use DTW when two trajectories follow the same corridor at different speeds or sampling rates, because it stretches the time axis to align them. Use Hausdorff when you only care about geometric coverage and not traversal order. Hausdorff ignores direction, so it will call a route and its reverse identical.

Why is my distance matrix so slow to compute?

A pairwise matrix over N trajectories needs N×(N−1)/2 comparisons, and each DTW or Fréchet comparison is itself O(m·n) in the point counts. The combined cost is quadratic in trajectories and quadratic in points. Resample to a bounded point count, add a Sakoe-Chiba band to DTW, prune distant pairs with a centroid pre-filter, and parallelise the upper triangle across cores.

Should I cluster with DBSCAN or agglomerative clustering on a distance matrix?

Both accept a precomputed matrix. Use DBSCAN(metric="precomputed") when you expect noise and an unknown number of route families, because it labels outliers as -1. Use agglomerative clustering with average or complete linkage when you want a dendrogram, a fixed cluster count, or a deterministic hierarchy you can cut at multiple thresholds.

Do I need to project coordinates before measuring trajectory similarity?

Yes. Every measure reduces to summing or maximising point-to-point distances, and those are only metric in a projected CRS. Computing DTW or Fréchet on EPSG:4326 degrees mixes latitude and longitude units that scale differently with latitude, so a threshold tuned in one city breaks in another. Project to a common UTM zone first.

How do I pick the DBSCAN eps threshold for trajectories?

Because the matrix is precomputed in metres, eps is a real distance in metres. Plot the sorted k-distance graph and set eps near the knee, then anchor to your use case: route deduplication wants tens of metres, commute clustering a few hundred, and anomaly detection a wide radius so only unusual paths become noise.

Related

• Measuring Trajectory Similarity with Fréchet and DTW
• Trajectory Segmentation
• Stay-Point Detection Algorithms
• Trajectory Object Design Patterns
• Speed & Acceleration Profiling

Back to Movement Pattern Extraction & Trajectory Analysis