Movement Pattern Extraction & Trajectory Analysis

Without a disciplined extraction pipeline, raw GPS and telemetry streams are too noisy, irregular, and semantically thin to drive operational decisions — logistics teams misroute vehicles, urban analysts miscount dwell times, and ML models train on artifacts rather than behaviour. This guide covers the data structures, architecture choices, core algorithms, Python tooling, and validation patterns that mobility engineers need to turn coordinate sequences into reliable, queryable movement intelligence.

Prerequisites & Scope

The techniques here assume you are working with ordered spatiotemporal point sequences — each record carrying at minimum a timestamp, latitude, longitude, and a device or entity identifier. Typical sources are GNSS receivers, cellular network positions, AIS transponders, or IoT telematics units.

Python stack assumed: Python 3.10+, movingpandas 0.17+, geopandas 0.14+, pyproj 3.6+, shapely 2.x, scikit-learn 1.4+, pandas 2.x, numpy 1.26+.

Data format assumed: A pandas.DataFrame or GeoParquet file with columns device_id (str), timestamp (datetime64, UTC), lat (float64), lon (float64), and optionally accuracy_m (float32) and altitude_m (float32). All spatial computations must be performed in a metric projected CRS — not in WGS84 (EPSG:4326) — to keep distance and velocity units in metres and metres-per-second.

Before applying any of the algorithms below, the raw stream must pass through the spatiotemporal data foundations layer: GPS precision and error handling to remove jump artifacts, coordinate reference system mapping to project into an appropriate UTM zone, and time-series synchronization strategies to align multi-sensor or multi-device streams before merging.

Core Conceptual Model

The central abstraction in trajectory analytics is the trajectory object: a temporally ordered, spatially continuous line representing a single entity’s path through space over a bounded time interval. In code this is most naturally a movingpandas.Trajectory, which wraps a GeoDataFrame with a DatetimeIndex and enforces temporal ordering.

Three derived concepts sit above the raw trajectory:

• Segments — sub-trajectories split at semantic boundaries (a stay point, a fuel stop, a trip start/end). Segmentation is the gateway to episode-level analytics.
• Stay points — compact spatial regions where movement velocity drops below a threshold for a minimum duration. They represent real-world activities: deliveries, parking, charging sessions, retail visits.
• Kinematic features — scalar time series derived from the trajectory: instantaneous speed, heading (bearing), angular velocity (turn rate), acceleration, and jerk. These are the inputs to behaviour classification models.

The diagram below shows how raw points flow through these three levels of abstraction:

Architecture Decision Map

Choosing the right design patterns early avoids expensive rewrites when data volume or transport-mode diversity grows. The table below covers the four most consequential architectural choices in a trajectory analytics system:

Pipeline Integration

Trajectory analysis does not exist in isolation. It sits in the middle of a larger movement analytics stack. The diagram below maps the upstream and downstream dependencies:

Upstream (Preprocessing): Raw GPS streams enter through a message broker (Kafka, Pub/Sub) and pass through GPS error handling, CRS projection, and sampling rate optimization before arriving as clean (device_id, timestamp, geometry) records.

This layer (Trajectory Analysis): Converts clean point sequences into segmented trajectories, stay points, kinematic time series, and directional features. Outputs land in a columnar feature store (Parquet or Delta Lake) and, for latency-sensitive use cases, in a PostGIS or TimescaleDB read layer.

Downstream (Consumers): Enriched features feed safety scoring models, demand forecasting pipelines, BI dashboards (via materialized views), route optimization APIs, and real-time alerting systems for geofence breaches or anomalous behaviour.

Core Algorithmic Frameworks

Semantic Segmentation and Stay-Point Identification

Raw trajectories are continuous, but movement is episodic. Segmenting trajectories into “move” and “stop” states is the first analytical step and the prerequisite for everything above it. Stay-point detection algorithms use either radius-time thresholds — where a point is flagged as a candidate stop if all subsequent points within a time window T stay within radius R — or density-based spatial clustering via DBSCAN variants. The choice depends on whether you need deterministic, per-trip stops (radius-time) or post-hoc activity-zone mining across a fleet (DBSCAN).

Stay points map directly to semantic locations: delivery addresses, fuel stops, parking zones, or charging sessions. Once identified they can be reverse-geocoded against POI databases to attach human-readable labels, turning a coordinate into “Depot A — 14 min dwell”.

import movingpandas as mpd
import geopandas as gpd
from shapely.geometry import Point
import pandas as pd

def detect_stay_points(
    traj: mpd.Trajectory,
    max_diameter_m: float = 100.0,
    min_duration_s: float = 300.0,
) -> gpd.GeoDataFrame:
    """
    Return stay points for a single trajectory using MovingPandas
    stop-detection. Requires the trajectory CRS to be metric (not WGS84).

    Parameters
    ----------
    traj            : movingpandas.Trajectory — must be in a metric CRS
    max_diameter_m  : spatial threshold — points within this diameter qualify
    min_duration_s  : minimum stationary time in seconds

    Returns
    -------
    GeoDataFrame with columns [geometry, start_time, end_time, duration_s]
    """
    if traj.crs.is_geographic:
        raise ValueError(
            f"Trajectory CRS {traj.crs} is geographic. "
            "Project to a metric CRS (e.g. EPSG:32633) before calling this function."
        )
    if len(traj.df) < 2:
        return gpd.GeoDataFrame(columns=["geometry", "start_time", "end_time", "duration_s"])

    detector = mpd.TrajectoryStopDetector(traj)
    stops = detector.get_stop_points(
        min_duration=pd.Timedelta(seconds=min_duration_s),
        max_diameter=max_diameter_m,
    )
    stops["duration_s"] = (
        stops["end_time"] - stops["start_time"]
    ).dt.total_seconds()
    return stops[["geometry", "start_time", "end_time", "duration_s"]]

Kinematic Feature Engineering

Beyond coordinates, movement intelligence requires derived kinematic metrics. Speed and acceleration profiling computes instantaneous speed, heading, acceleration, and jerk via vectorized finite differences over the trajectory’s metric-projected coordinates. These scalars become the feature vectors for behaviour classification: aggressive-braking detection, highway-cruising band identification, and idle-time distribution analysis. In fleet management they feed safety scoring models and predictive maintenance triggers.

import numpy as np
import geopandas as gpd
import pandas as pd
from pyproj import Transformer

def compute_kinematics(
    gdf: gpd.GeoDataFrame,
    crs_metric: str = "EPSG:32633",
) -> gpd.GeoDataFrame:
    """
    Add speed_mps, heading_deg, accel_mps2 columns to a point GeoDataFrame.
    Input must be sorted by timestamp ascending.

    Parameters
    ----------
    gdf        : GeoDataFrame with geometry (Point) and DatetimeIndex
    crs_metric : a projected metric CRS — NEVER pass EPSG:4326 here

    Returns
    -------
    GeoDataFrame with three new columns; first row has NaN for diff-based cols.
    """
    if gdf.crs.is_geographic:
        gdf = gdf.to_crs(crs_metric)

    if not isinstance(gdf.index, pd.DatetimeIndex):
        raise TypeError("GeoDataFrame must have a DatetimeIndex.")
    if len(gdf) < 2:
        gdf["speed_mps"] = np.nan
        gdf["heading_deg"] = np.nan
        gdf["accel_mps2"] = np.nan
        return gdf

    x = gdf.geometry.x.to_numpy()
    y = gdf.geometry.y.to_numpy()
    dt = np.diff(gdf.index.astype("int64")) / 1e9  # nanoseconds → seconds

    dx = np.diff(x)
    dy = np.diff(y)
    dist = np.hypot(dx, dy)

    speed = np.empty(len(gdf))
    speed[0] = np.nan
    speed[1:] = dist / np.where(dt == 0, np.nan, dt)

    heading = np.empty(len(gdf))
    heading[0] = np.nan
    heading[1:] = np.degrees(np.arctan2(dx, dy)) % 360.0

    accel = np.empty(len(gdf))
    accel[0] = np.nan
    accel[1] = np.nan
    # dt between consecutive speed samples (midpoints)
    dt2 = (dt[:-1] + dt[1:]) / 2.0
    accel[2:] = np.diff(speed[1:]) / np.where(dt2 == 0, np.nan, dt2)

    gdf = gdf.copy()
    gdf["speed_mps"] = speed
    gdf["heading_deg"] = heading
    gdf["accel_mps2"] = accel
    return gdf

Directionality and Geometric Pattern Mining

Route geometry encodes navigational constraints and behavioural preferences. Directionality and turn analysis extracts turn angles, cumulative curvature, and intersection compliance by computing bearing deltas between consecutive trajectory segments. U-turn detection, illegal-manoeuvre flagging, and route deviation alerts all derive from this layer. When combined with road-network topology (OpenStreetMap via OSMnx, or HERE Maps), geometric features can be map-matched to graph edges, enabling network-constrained routing analysis and traffic flow modelling.

Anomaly Detection and Behavioural Drift

Movement patterns are rarely static across time. Seasonal changes, infrastructure disruptions, or fleet policy shifts manifest as statistical drift in trajectory distributions. Distributional distance metrics (Wasserstein distance, KL divergence) applied to sliding windows of speed profiles or stay-point durations can flag when the current behavioural baseline diverges from a historical reference. For real-time alerting, isolation forests or rule-based kinematic thresholds isolate individual-trip anomalies — geofence breaches, route hijacking, sensor degradation — without requiring labelled training data.

Engineering Pitfalls and Production Gotchas

These are the five failure modes most commonly encountered in production trajectory systems, along with their diagnostic signals:

1. Computing distances in geographic coordinates (EPSG:4326) Applying shapely.distance() or NumPy Euclidean distance to raw (lon, lat) tuples gives results in decimal degrees, not metres. Speed values will be off by a factor of roughly 111,000 at the equator (and more at higher latitudes). Diagnostic: check if speed_mps values exceed 200 m/s for ordinary vehicle trips. Fix: always call .to_crs("EPSG:32633") or the appropriate UTM zone before computing distances.

2. Aliasing through short stops at low sampling rates A 5-minute delivery stop sampled at 60-second intervals may generate only 5 points, all within the GPS accuracy circle, making it indistinguishable from stationary noise. The sampling rate optimization layer should be calibrated to the shortest behavioural episode you need to detect — not just the shortest movement segment.

3. Clock drift and timestamp collisions in merged streams When combining data from multiple sensors or network handoffs, the same device_id can appear with two records at identical timestamps after rounding. This breaks temporal ordering and produces zero dt values that cause division-by-zero in speed calculations. Diagnostic: assert not gdf.index.duplicated().any(). Fix: apply the time-series synchronization deduplication step before segmenting.

4. Trajectory object fragmentation from over-aggressive splitting Splitting trajectories at every data gap — regardless of gap duration — creates thousands of micro-segments that overwhelm the feature store and inflate stay-point counts. Use a minimum gap threshold (e.g. 5 minutes for vehicle data) and only split when the implied speed during the gap is implausibly high. MovingPandas’s split_by_observation_gap() accepts a min_gap parameter for this.

5. Projection zone boundary crossing in long-haul routes A fleet operating across multiple UTM zones will produce discontinuous metric coordinates when a single EPSG code is used throughout. Routes crossing zone boundaries exhibit a sudden coordinate jump that appears as an implausible velocity spike. Fix: use a continental equal-area projection (e.g. EPSG:3035 for Europe, EPSG:5070 for the contiguous US) for multi-zone datasets, or split trajectories at zone boundaries and process each zone independently.

Python Tooling Landscape

For the stay-point and segmentation workflows, movingpandas is the most complete starting point. For trajectory object design patterns and schema conventions — including how to build a TrajectoryCollection from raw DataFrame input — see the dedicated page on structuring trajectory data in GeoPandas.

For compute-intensive fleet analytics (10M+ daily points), replacing pandas groupby loops with dask.dataframe and parallelising the kinematic computation across device_id partitions typically reduces runtime by 10–50× on an 8-core machine, with no logic changes beyond wrapping the compute graph.

Validation and Testing Patterns

Trajectory analytics pipelines are stateful and geometrically sensitive. Standard unit tests catch logic errors, but spatiotemporal-specific sanity checks catch the more insidious projection and ordering bugs:

Velocity bounds assertion: After computing speed_mps, assert that no value exceeds the physical maximum for the transport mode (e.g. 70 m/s for road vehicles, 0.5 m/s for pedestrians in a dwell zone). Values above threshold indicate either a projection error or an unfiltered GPS jump artifact.

MAX_ROAD_SPEED_MPS = 70.0  # ~252 km/h — well above legal limits

def assert_speed_bounds(gdf: gpd.GeoDataFrame) -> None:
    """Raise ValueError if any non-NaN speed value exceeds physical maximum."""
    speeds = gdf["speed_mps"].dropna()
    if (speeds > MAX_ROAD_SPEED_MPS).any():
        n_violations = (speeds > MAX_ROAD_SPEED_MPS).sum()
        raise ValueError(
            f"{n_violations} speed values exceed {MAX_ROAD_SPEED_MPS} m/s. "
            "Check CRS projection and GPS outlier filtering."
        )

Temporal monotonicity check: Trajectories must have strictly increasing timestamps within each device_id partition. A violated ordering breaks all finite-difference calculations silently.

def assert_temporal_order(gdf: gpd.GeoDataFrame) -> None:
    """Raise if timestamps are not strictly increasing."""
    if not gdf.index.is_monotonic_increasing:
        raise ValueError(
            "Timestamps are not monotonically increasing. "
            "Sort by index before computing kinematic features."
        )

Inverse-transform round-trip test: After projecting to a metric CRS and back to WGS84, coordinates should round-trip within sub-metre accuracy. A failing round-trip indicates a malformed pyproj pipeline definition.

from pyproj import Transformer

def assert_crs_round_trip(lat: float, lon: float, epsg_metric: int = 32633) -> None:
    """Verify that WGS84 → metric → WGS84 round-trip error is below 1 mm."""
    fwd = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg_metric}", always_xy=True)
    inv = Transformer.from_crs(f"EPSG:{epsg_metric}", "EPSG:4326", always_xy=True)
    x, y = fwd.transform(lon, lat)
    lon2, lat2 = inv.transform(x, y)
    assert abs(lat - lat2) < 1e-8 and abs(lon - lon2) < 1e-8, (
        f"Round-trip error too large: Δlat={abs(lat-lat2):.2e}, Δlon={abs(lon-lon2):.2e}"
    )

Stay-point coverage sanity check: For a typical urban delivery fleet, stay points should account for 30–60% of total trip time. A value below 10% suggests the stop detector threshold is too permissive; above 80% suggests the trajectory was not pre-filtered for parking/overnight idle.

Frequently Asked Questions

Q: When should I use MovingPandas’s built-in speed computation versus writing my own? MovingPandas’s add_speed() method handles CRS detection and computes speed between consecutive points in m/s when the trajectory is in a metric CRS. Write your own only when you need vectorised batch processing across millions of trajectories simultaneously — in that case, NumPy-based finite differences on pre-extracted coordinate arrays are 5–20× faster than the MovingPandas row-level implementation.

Q: How many points does a trajectory need for reliable kinematic analysis? Acceleration and jerk require at least 3 and 4 consecutive points respectively (due to second and third-order differences). For meaningful statistical summaries of a speed profile (mean, 95th percentile), a minimum of 20 points is a practical floor. Trajectories shorter than this should be flagged and excluded from model training.

Q: Can I run DBSCAN for stay-point detection directly on WGS84 coordinates? Technically yes, using the Haversine metric in scikit-learn’s DBSCAN (metric='haversine', eps in radians). But this adds complexity, runs slower than Euclidean DBSCAN, and means your epsilon parameter has no intuitive metre-equivalent — making calibration harder. Project to a metric CRS and use Euclidean DBSCAN; the results are equivalent and the parameter semantics are transparent.

Related

• Stay-Point Detection Algorithms — radius-time and DBSCAN methods for identifying where an entity stopped, and for how long
• Speed & Acceleration Profiling — computing instantaneous speed, jerk, and kinematic signatures from discrete GPS points
• Directionality & Turn Analysis — bearing deltas, U-turn detection, and geometric pattern mining
• GPS Precision & Error Handling — filtering jump artifacts and multipath errors before trajectory analysis
• Coordinate Reference System Mapping — projecting to metric CRS and managing zone boundaries

Back to Spatiotemporal Data Foundations & Structures