Directionality & Turn Analysis: Quantifying Heading Changes in Movement Trajectories

Directionality & Turn Analysis is the process of computing forward bearings, angular velocity, and turn severity categories from a time-ordered sequence of GPS coordinates — transforming raw position records into a structured representation of navigational intent. Without this layer, movement pipelines can tell you where an entity was but not how it navigated there: intersection compliance, U-turn frequency, and routing inefficiency remain invisible. This topic sits within Movement Pattern Extraction & Trajectory Analysis and feeds directly into fleet safety scoring, driver behavior modeling, and anomaly detection workflows.

Prerequisites Checklist

Complete these upstream stages before computing any directional metrics:

• Python 3.10+ — pandas 2.x, numpy 1.26+, geopandas 0.14+, scipy 1.12+, pyproj 3.6+
• Input schema: each record must carry track_id (trajectory identifier), ts (timezone-aware datetime), lat and lon (WGS84 decimal degrees). Optional: speed_kmh, heading_sensor (on-board compass reading)
• Upstream stage: GPS precision & error handling must be complete — multipath-flagged points and positional outliers removed — before bearings are computed, because a single erroneous coordinate produces two false turns
• Coordinate reference system: bearing computations use geographic coordinates (EPSG:4326) with spherical trigonometry. Distance and speed calculations must use a projected metric CRS — see coordinate reference system mapping for projection selection
• Sampling regularity: resample trajectories to a uniform interval using time-series synchronization strategies before computing headings; irregular gaps amplify noise into false turn events

Failure-Mode Taxonomy

Deterministic Pipeline Overview

The five stages below must run in order. Skipping or reordering stages amplifies GPS noise into classification errors.

Stage 1 — Segment & Sort: Group records by track_id, sort by ts ascending. Remove duplicate timestamps (keep first) and runs of spatially identical consecutive points. Flag time gaps larger than twice the median interval as segment boundaries; split there rather than compute bearings across a gap.

Stage 2 — Filter Stops: Apply stay-point detection to identify dwell intervals. Mark or drop all points inside identified stay zones before computing headings. This is the single most impactful noise-reduction step: unfiltered stops can contribute 20–40% of all apparent sharp turns.

Stage 3 — Compute Bearings: Apply the spherical forward-azimuth formula (arctan2) to consecutive coordinate pairs. Output is a bearing column in [0, 360) for each point. See the implementation walkthrough for the correct argument ordering.

Stage 4 — Derive Angular Velocity: Compute wrapped heading differences normalized to [−180, 180], then divide by the time delta in seconds. Mask rows where dt ≤ 0. Optionally smooth with a Savitzky-Golay filter (window 5, polynomial 2) before thresholding.

Stage 5 — Classify Turns: Map angular velocity values to categorical turn types using mode-specific threshold tables. See the Calibration section for validated bins by transport mode.

Implementation Walkthrough

The function below handles the complete stages 3–5 in a single vectorized pass. It is schema-validated, defensive against NaN propagation, and handles the 0°/360° wrap-around correctly.

import numpy as np
import pandas as pd
from scipy.signal import savgol_filter


REQUIRED_COLS = {"track_id", "ts", "lat", "lon"}
TURN_THRESHOLDS = {
    "vehicle":    {"minor": 0.5, "standard": 3.0, "sharp": 8.0},
    "pedestrian": {"minor": 1.0, "standard": 5.0, "sharp": 15.0},
    "cycling":    {"minor": 0.8, "standard": 4.0, "sharp": 12.0},
}


def compute_turn_metrics(
    df: pd.DataFrame,
    mode: str = "vehicle",
    smooth: bool = True,
    sg_window: int = 5,
    sg_poly: int = 2,
) -> pd.DataFrame:
    """
    Compute bearing, heading difference, angular velocity, and turn class.

    Parameters
    ----------
    df : pd.DataFrame
        Must contain 'track_id', 'ts' (datetime64[ns, tz]), 'lat', 'lon'.
        Must already be sorted chronologically per track and stops filtered.
    mode : str
        Transport mode key — one of 'vehicle', 'pedestrian', 'cycling'.
    smooth : bool
        Apply Savitzky-Golay smoothing to angular velocity before classifying.
    sg_window, sg_poly : int
        Savitzky-Golay window length (odd) and polynomial order.

    Returns
    -------
    pd.DataFrame with added columns: bearing, heading_diff, dt, angular_velocity,
        angular_velocity_smooth (if smooth=True), turn_class.

    Raises
    ------
    ValueError  if required columns are missing or mode is unknown.
    """
    missing = REQUIRED_COLS - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if mode not in TURN_THRESHOLDS:
        raise ValueError(f"Unknown mode '{mode}'. Choose from {list(TURN_THRESHOLDS)}")
    if df.empty or len(df) < 2:
        # Return with NaN columns — caller decides whether to skip or raise
        for col in ["bearing", "heading_diff", "dt", "angular_velocity", "turn_class"]:
            df = df.copy()
            df[col] = np.nan
        return df

    thresholds = TURN_THRESHOLDS[mode]
    df = df.copy()

    # Convert to radians for spherical formula
    lat1 = np.radians(df["lat"].shift(1))
    lon1 = np.radians(df["lon"].shift(1))
    lat2 = np.radians(df["lat"])
    lon2 = np.radians(df["lon"])
    d_lon = lon2 - lon1

    # Spherical forward-azimuth: arctan2(y, x)
    # y = sin(Δlon) * cos(lat2)
    # x = cos(lat1)*sin(lat2) − sin(lat1)*cos(lat2)*cos(Δlon)
    # IMPORTANT: argument order is arctan2(y, x) — swapping gives bearings offset by 90°
    y = np.sin(d_lon) * np.cos(lat2)
    x = np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(d_lon)
    df["bearing"] = np.degrees(np.arctan2(y, x)) % 360.0

    # Wrap heading difference to [−180, 180] to avoid 360° spike artefacts
    raw_diff = df["bearing"] - df["bearing"].shift(1)
    df["heading_diff"] = (raw_diff + 180.0) % 360.0 - 180.0

    # Time delta in seconds; mask zero/negative intervals (duplicates or out-of-order)
    df["dt"] = df["ts"].diff().dt.total_seconds()
    valid_dt = df["dt"] > 0

    df["angular_velocity"] = np.nan
    df.loc[valid_dt, "angular_velocity"] = (
        df.loc[valid_dt, "heading_diff"].abs() / df.loc[valid_dt, "dt"]
    )

    # Optional Savitzky-Golay smoothing on angular velocity
    av_col = "angular_velocity"
    if smooth and df["angular_velocity"].notna().sum() >= sg_window:
        filled = df["angular_velocity"].fillna(0.0).to_numpy()
        smoothed = savgol_filter(filled, window_length=sg_window, polyorder=sg_poly)
        df["angular_velocity_smooth"] = smoothed
        df.loc[~valid_dt, "angular_velocity_smooth"] = np.nan
        av_col = "angular_velocity_smooth"

    # Turn severity classification
    av = df[av_col]
    conditions = [
        av.between(thresholds["minor"], thresholds["standard"], inclusive="left"),
        av.between(thresholds["standard"], thresholds["sharp"], inclusive="left"),
        av >= thresholds["sharp"],
    ]
    choices = ["minor", "standard", "sharp"]
    df["turn_class"] = np.select(conditions, choices, default="none")
    # Rows without a valid dt remain unclassified
    df.loc[~valid_dt, "turn_class"] = np.nan

    return df

Key reliability points:

• arctan2(y, x) argument order is critical. The formula above assigns y = sin(Δlon)·cos(lat2) and x = cos(lat1)·sin(lat2) − sin(lat1)·cos(lat2)·cos(Δlon). Swapping produces bearings offset by exactly 90°, which passes all unit tests unless you explicitly validate against a known bearing.
• Wrap normalization prevents 360° spikes. The (diff + 180) % 360 − 180 formulation correctly handles transitions across 0°/360° — without it, a vehicle heading from 5° to 355° registers a −350° change instead of the correct −10°.
• Never compute bearing across a trajectory segment boundary. If you split at large time gaps (Stage 1), ensure shift(1) does not reach across the split. Group-level application (.groupby("track_id").apply(...)) enforces this automatically.
• Smoothing angular velocity, not raw bearings. Smoothing the bearing column introduces wrap-around discontinuities at 0°/360°. The Savitzky-Golay filter belongs on the angular velocity (or heading difference) signal, which is already bounded and continuous.

Geometric & Mathematical Grounding

The forward azimuth from point A to point B on a sphere is:

θ = atan2(sin(Δλ)·cos(φ₂),  cos(φ₁)·sin(φ₂) − sin(φ₁)·cos(φ₂)·cos(Δλ))

where φ is latitude in radians and λ is longitude in radians. The result is in (−π, π]; adding 360 mod 360 maps it to [0°, 360°).

Angular velocity is the instantaneous rate of heading change:

ω = |Δθ_wrapped| / Δt   (degrees per second)

where Δθ_wrapped is the heading difference normalized to [−180, 180] and Δt is the elapsed time in seconds. High angular velocity at low linear speed is a hallmark of tight parking manoeuvres or pedestrian pivots; the same angular velocity at highway speed indicates a lane-change or hazard avoidance event. Separating these cases requires pairing turn metrics with speed & acceleration profiling.

For distances greater than ~500 km, or for datasets crossing the ±180° meridian, replace the spherical formula with pyproj.Geod(ellps="WGS84").inv(lon1, lat1, lon2, lat2) which returns the geodesic forward azimuth on the WGS84 ellipsoid.

Calibration & Parameter Tuning

Turn severity thresholds are not universal. The table below lists validated starting points by transport mode; recalibrate using the 90th-percentile angular velocity distribution from your own dataset before deploying to production.

Sampling-rate corrections: thresholds derived at 1 Hz do not transfer directly to 5 Hz data. At higher rates, multiply threshold values by the ratio of new_rate / reference_rate (e.g., 1 Hz thresholds × 5 for 5 Hz data) as a first approximation, then validate against labelled manoeuvre samples.

Segment-length minimum: classification is unreliable on segments shorter than three valid bearing intervals. Drop or flag turn events from short segments (fewer than 5 points after stop filtering) and exclude them from aggregate statistics.

Integration & Compatibility

Upstream dependencies:

• GPS precision & error handling — multipath removal and HDOP filtering must precede bearing computation
• Stay-point detection — dwell masks are required inputs to Stage 2
• Sampling rate optimization — uniform resampling stabilizes angular velocity; see downsampling high-frequency GPS tracks for the recommended approach

Downstream consumers:

• Speed & acceleration profiling — combine angular_velocity with speed_ms to build composite manoeuvre risk scores; a sharp turn at 60 km/h is categorically different from the same heading change at 5 km/h
• Detecting U-turns specifically — see Detecting U-turns and directional shifts in fleet data for threshold calibration and false-positive suppression strategies tuned to commercial vehicle dynamics
• Hidden Markov Model map matching — turn class sequences are strong evidence features for HMM state transitions between road segments
• Driver behavior scoring — aggregate turn_class counts per trip form the directional component of composite safety indices

Output schema from compute_turn_metrics: the returned DataFrame adds bearing (°, [0, 360)), heading_diff (°, [−180, 180]), dt (seconds), angular_velocity (°/s), optionally angular_velocity_smooth (°/s), and turn_class (none | minor | standard | sharp | NaN). This schema is consumed directly by downstream aggregation and scoring stages.

FAQ

Why do I see 360° heading spikes even after filtering stops? These are wrap-around artefacts from computing raw differences across the 0°/360° boundary. Normalize heading differences to [−180, 180] using (diff + 180) % 360 − 180 before any further processing. The raw bearing_b − bearing_a approach is the most common implementation mistake in this domain.

Can I compute bearings directly from EPSG:4326 coordinates? Yes — the spherical arctan2 formula operates correctly on geographic lat/lon in radians. However, distance-based metrics (speed, displacement) must use a projected CRS or a haversine/geodesic formula. Never use planar Euclidean distance on raw WGS84 data — see coordinate reference system mapping for projection guidance.

What sampling rate gives the cleanest turn classification? 1–2 Hz is optimal for vehicle-level manoeuvre classification. Above 5 Hz, steering micro-corrections dominate and inflate turn counts by 3–8×. Below 0.2 Hz (typical wildlife GPS collars), sharp turns become aliased into apparent straight segments. Match your threshold table to your actual fix rate.

How do I handle trajectories that cross the ±180° meridian? Use pyproj.Geod(ellps="WGS84").inv(lon1, lat1, lon2, lat2) for geodesic forward azimuth. The standard spherical arctan2 formula produces incorrect bearings for trans-meridian segments because the longitude difference wraps sign unexpectedly.

Should I smooth bearings before or after classification? Smooth angular velocity — not raw bearings — using a Savitzky-Golay filter before thresholding. Smoothing raw bearings introduces discontinuities at the 0°/360° boundary; smoothing the velocity signal (which is already bounded and continuous) avoids this entirely. Use sg_window=5, sg_poly=2 as a starting point and increase the window if noise remains dominant.

Related:

• Detecting U-turns and directional shifts in fleet data
• Speed & Acceleration Profiling
• Stay-Point Detection Algorithms
• GPS Precision & Error Handling
• Sampling Rate Optimization

Back to Movement Pattern Extraction & Trajectory Analysis