Measuring trajectory similarity with Fréchet and DTW

Discrete Fréchet distance and dynamic time warping (DTW) are the two workhorse measures for comparing the shape of two ordered trajectories. Both solve the alignment problem with a small dynamic program over an m×n table, and both must run on projected metric coordinates to return meaningful numbers. This guide gives complete, self-contained implementations of each, the complexity you should expect, the edge cases that break naive code, and a worked comparison that shows exactly when the two measures disagree.

Why this happens: two ways to align two paths

The hard part of trajectory comparison is not measuring distance between two points — it is deciding which points to pair. Two vehicles driving the same road rarely produce samples at the same positions: one logs faster, one starts a little earlier, one pauses at a light. A rigid index-to-index pairing (pᵢ with qᵢ) breaks the moment the point counts differ or the cadence drifts, a failure mode covered in the parent movement similarity and clustering reference.

Fréchet and DTW both solve this by searching over all monotone couplings — pairings that never move backward along either path — and picking the best one. They differ only in how they score a coupling. Fréchet takes the maximum coupled distance (the worst gap on the best alignment); DTW takes the sum. That single choice explains every behavioural difference between them. Because both reduce to Euclidean point distances, the same projection discipline that governs speed and acceleration profiling applies here: measure in metres, never in degrees.

Core comparison pipeline

Trajectory Similarity Computation Steps Four sequential steps: Project both trajectories to a metric CRS, Guard edge cases, Fill the dynamic-programming table, and Read and validate the distance, connected by arrows. 1. Project to metric CRS 2. Guard edge cases 3. Fill DP table 4. Read & validate
  1. Project both trajectories from EPSG:4326 to a shared metric CRS so every point distance is in metres.
  2. Guard edge cases — empty arrays, single-point tracks, and mismatched lengths must be handled before the dynamic program runs.
  3. Fill the DP table — Fréchet propagates a running maximum; DTW propagates a running sum with an optional warping band.
  4. Read and validate — the bottom-right cell is the distance; check it against known invariants.

Production-ready Python implementation

Both functions take already-projected (m, 2) arrays in metres. A small helper projects raw WGS84 input first so the examples are runnable end to end. The Fréchet routine uses an iterative table rather than the classic recursion to avoid Python recursion-depth limits on long tracks.

PYTHON 
from __future__ import annotations

import numpy as np
from pyproj import Transformer


def project_to_metric(track_lonlat: np.ndarray, target_crs_epsg: int = 32633) -> np.ndarray:
    """
    Project an (m, 2) array of (lon, lat) in EPSG:4326 to a metric CRS.

    Distances for Fréchet and DTW are only meaningful in a projected CRS
    (metres). NEVER feed raw EPSG:4326 degrees into the distance functions.
    """
    if track_lonlat.ndim != 2 or track_lonlat.shape[1] != 2:
        raise ValueError(f"expected an (m, 2) array, got shape {track_lonlat.shape}")
    to_metric = Transformer.from_crs("EPSG:4326", f"EPSG:{target_crs_epsg}", always_xy=True)
    x, y = to_metric.transform(track_lonlat[:, 0], track_lonlat[:, 1])  # lon, lat order
    return np.column_stack([x, y])


def discrete_frechet(p: np.ndarray, q: np.ndarray) -> float:
    """
    Discrete Fréchet distance between two (m, 2) and (n, 2) metric-CRS tracks.

    Uses an iterative dynamic-programming table:
        ca[i, j] = max( ‖p_i - q_j‖ ,
                        min( ca[i-1, j], ca[i-1, j-1], ca[i, j-1] ) )
    The result ca[m-1, n-1] is the minimum leash length over all monotone
    couplings — the classic dog-walking formulation.

    Edge cases
    ----------
    - Empty input on either side -> ValueError (distance undefined).
    - Single point on one side -> reduces to the max nearest-vertex distance.

    Complexity: O(m * n) time, O(m * n) space.
    """
    if p.size == 0 or q.size == 0:
        raise ValueError("discrete_frechet is undefined for an empty trajectory.")
    if p.ndim != 2 or q.ndim != 2 or p.shape[1] != 2 or q.shape[1] != 2:
        raise ValueError("both inputs must be (k, 2) arrays of projected coordinates.")

    m, n = len(p), len(q)
    # Pairwise Euclidean distances in metres, fully vectorized.
    diff = p[:, None, :] - q[None, :, :]          # (m, n, 2)
    dist = np.sqrt((diff ** 2).sum(axis=2))       # (m, n)

    ca = np.empty((m, n), dtype=float)
    ca[0, 0] = dist[0, 0]
    for i in range(1, m):                         # first column: only downward moves
        ca[i, 0] = max(ca[i - 1, 0], dist[i, 0])
    for j in range(1, n):                         # first row: only rightward moves
        ca[0, j] = max(ca[0, j - 1], dist[0, j])
    for i in range(1, m):
        for j in range(1, n):
            ca[i, j] = max(min(ca[i - 1, j], ca[i - 1, j - 1], ca[i, j - 1]), dist[i, j])
    return float(ca[m - 1, n - 1])


def dtw_distance(p: np.ndarray, q: np.ndarray, band: int | None = None) -> float:
    """
    Dynamic time warping distance between two (m, 2) and (n, 2) metric-CRS tracks.

    Minimises the SUM of coupled Euclidean distances over all monotone
    warping paths:
        D[i, j] = ‖p_i - q_j‖ + min( D[i-1, j], D[i, j-1], D[i-1, j-1] )

    Parameters
    ----------
    band : optional Sakoe-Chiba half-width. When set, only cells with
           |i - j| <= band are considered, bounding cost at O(m * band)
           and preventing pathological warps. Pass None for the full table.

    Edge cases
    ----------
    - Empty input on either side -> ValueError.
    - Single point on one side -> sum of distances to that point.

    Complexity: O(m * n) unbanded; O(m * band) banded. O(n) rolling space.
    """
    if p.size == 0 or q.size == 0:
        raise ValueError("dtw_distance is undefined for an empty trajectory.")
    if p.ndim != 2 or q.ndim != 2 or p.shape[1] != 2 or q.shape[1] != 2:
        raise ValueError("both inputs must be (k, 2) arrays of projected coordinates.")

    m, n = len(p), len(q)
    if band is None:
        band = max(m, n)                          # effectively unbanded

    prev = np.full(n + 1, np.inf)
    prev[0] = 0.0
    for i in range(1, m + 1):
        curr = np.full(n + 1, np.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(p[i - 1, 0] - q[j - 1, 0], p[i - 1, 1] - q[j - 1, 1])
            curr[j] = cost + min(prev[j], curr[j - 1], prev[j - 1])
        prev = curr
    return float(prev[n])

Two design choices are worth calling out. First, discrete_frechet precomputes the full pairwise distance matrix with NumPy broadcasting, so only the O(m·n) table fill runs in Python — acceptable for tracks resampled to tens or low hundreds of points. Second, dtw_distance keeps only the previous row (prev) rather than the full table, cutting memory to O(n); this is safe because the distance value, unlike the alignment path, does not require backtracking.

Worked comparison

The snippet below builds three short synthetic tracks: a baseline, a slower re-sampled copy of the same corridor, and a genuinely different path with one large excursion. It shows how Fréchet and DTW rank them.

PYTHON 
import numpy as np

# Metric-CRS coordinates (metres) — pretend these came from project_to_metric.
base = np.array([[0, 0], [10, 0], [20, 0], [30, 0], [40, 0]], dtype=float)

# Same corridor, sampled at different density (slower vehicle, more points).
slower = np.array([[0, 0], [5, 0], [12, 1], [20, 0], [28, -1], [35, 0], [40, 0]], dtype=float)

# Different path: identical except one 15 m lateral excursion mid-route.
excursion = np.array([[0, 0], [10, 0], [20, 15], [30, 0], [40, 0]], dtype=float)

print("Fréchet base vs slower   :", round(discrete_frechet(base, slower), 2))
print("DTW     base vs slower   :", round(dtw_distance(base, slower, band=3), 2))
print("Fréchet base vs excursion:", round(discrete_frechet(base, excursion), 2))
print("DTW     base vs excursion:", round(dtw_distance(base, excursion, band=3), 2))

The pattern to expect: for base vs slower, DTW reports a small value because it warps the extra points onto the baseline, while Fréchet stays modest because the worst gap is only the ~1 m wobble. For base vs excursion, Fréchet jumps to roughly 15 (dominated by the single worst-aligned pair at the excursion), whereas DTW rises more gently because the lone excursion contributes one large term to a sum of otherwise-zero terms. This is the core intuition: Fréchet punishes the worst moment, DTW averages over the whole path. When you need route deduplication that must not tolerate a single large deviation, Fréchet is the stricter gate.

Validation block

Run these checks whenever you implement or modify either measure. Each targets a property the algorithm must satisfy, so a passing set gives confidence before the measure feeds a distance matrix.

PYTHON 
import numpy as np


def validate_similarity_measures() -> None:
    """
    Assert the mathematical invariants of both distance functions.
    Raises AssertionError on any violation; prints a summary on success.
    """
    rng = np.random.default_rng(42)
    a = np.column_stack([np.arange(20), rng.normal(0, 1, 20)]).astype(float)
    b = np.column_stack([np.arange(20), rng.normal(0, 1, 20)]).astype(float)

    # 1. Identity: distance from a track to itself is zero.
    assert discrete_frechet(a, a) == 0.0, "Fréchet self-distance must be 0"
    assert dtw_distance(a, a) == 0.0, "DTW self-distance must be 0"

    # 2. Non-negativity.
    assert discrete_frechet(a, b) >= 0.0
    assert dtw_distance(a, b) >= 0.0

    # 3. Fréchet symmetry: d(a, b) == d(b, a).
    assert np.isclose(discrete_frechet(a, b), discrete_frechet(b, a)), "Fréchet must be symmetric"

    # 4. Banded DTW is an upper bound on unbanded DTW (fewer paths allowed).
    full = dtw_distance(a, b, band=None)
    banded = dtw_distance(a, b, band=2)
    assert banded >= full - 1e-9, "Banded DTW cannot be smaller than unbanded"

    # 5. Edge cases raise, not silently return.
    for fn in (discrete_frechet, dtw_distance):
        try:
            fn(np.empty((0, 2)), a)
        except ValueError:
            pass
        else:
            raise AssertionError(f"{fn.__name__} must reject empty input")

    print("All similarity-measure invariants passed.")


validate_similarity_measures()

Common mistakes and gotchas

  • Computing in degrees instead of metres. Passing raw EPSG:4326 coordinates makes the returned “distance” a mixture of longitude and latitude degree units that scale differently with latitude. The number looks plausible but is not comparable across locations, so any clustering threshold silently breaks. Always call project_to_metric (or an equivalent) first — the same discipline required across coordinate reference system mapping.

  • Leaving the DTW warping window unbounded. Without a Sakoe-Chiba band, DTW can align the start of one track to the end of another, producing a small distance for two genuinely different paths and inflating runtime to full O(m·n). Bound the band to a small fraction of the track length once both tracks are resampled to comparable point counts.

  • Assuming DTW is symmetric and metric. DTW does not satisfy the triangle inequality and, with asymmetric step weights, is not even guaranteed symmetric. Do not build a metric-tree index on raw DTW, and do not assume d(a, b) == d(b, a) unless your step pattern is symmetric. Discrete Fréchet is a true metric and is safe for both.

  • Skipping the single-point and empty-track cases. A trajectory trimmed down to one surviving GPS fix, or an empty array after aggressive filtering, will crash a naive table fill or return a misleading zero. Guard these explicitly, as the implementations above do, before segmentation output ever reaches the measure.

  • Comparing unsegmented logs. Running either measure on a full multi-trip day log compares apples to oranges — the alignment tries to warp a morning commute onto an evening errand. Split into single trips first with trajectory segmentation so each comparison is between coherent journeys.

Back to Movement Similarity & Clustering