Trajectory Object Design Patterns

A trajectory object is a validated, metadata-enriched abstraction over a chronologically ordered GPS point sequence — not a bare LineString — and the design choices you make here propagate to every downstream stage of a movement analytics pipeline.

Prerequisites

Implement these patterns after completing the upstream Spatiotemporal Data Foundations & Structures ingestion stage. You need:

• Python 3.10+, with geopandas >= 0.14, shapely >= 2.0, pandas >= 2.0, and pyarrow >= 14
• Input data as a GeoDataFrame with a DatetimeTZDtype index and a geometry column of Point objects
• A projected CRS applied before object construction — see Coordinate Reference System Mapping for automated projection routing
• Upstream GPS Precision & Error Handling complete: multipath outliers removed, velocity-threshold filter applied, and CRS confirmed as metric

Raw (lat, lon, timestamp) tuples without these preconditions produce trajectory objects that silently report incorrect velocities, nonsensical path lengths, or corrupt geometry the moment they cross a projection boundary.

Failure Mode Taxonomy

Pipeline Overview

The five-stage construction pipeline below converts raw telemetry into a queryable trajectory object. Each stage is a discrete validation checkpoint — abort early rather than propagate corrupt state.

Stage 1 — Temporal sort and deduplication. Sort the raw GeoDataFrame by its DatetimeTZDtype index, drop exact duplicate timestamps, and assert strict monotonic increase before proceeding.

Stage 2 — Velocity threshold filter. Compute successive displacement / time delta for each segment. Discard points whose implied speed exceeds a mode-appropriate ceiling (see calibration table below). Use GPS Precision & Error Handling validation results as the source of that ceiling.

Stage 3 — CRS projection. Re-project from EPSG:4326 to a locally optimal metric CRS. For continental data, UTM zone selection based on trajectory centroid longitude minimises cumulative Mercator distortion. All LineString.length and velocity calculations must occur in this projected space — never in raw degrees.

Stage 4 — Object construction. Wrap the validated, projected GeoDataFrame in an immutable TrajectoryObject instance that enforces temporal ordering, derives geometry, and exposes analytical methods.

Stage 5 — Quality annotation. Attach a per-point quality series (boolean or [0, 1] float) encoding signal confidence. Downstream aggregations can filter on this series without re-parsing the full geometry.

Implementation Walkthrough

Core trajectory object

from __future__ import annotations
from typing import Optional
import pandas as pd
import geopandas as gpd
from shapely.geometry import LineString
from shapely.validation import make_valid


class TrajectoryObject:
    """Immutable, metadata-enriched trajectory with temporal indexing.

    Parameters
    ----------
    device_id:     Unique identifier for the originating device/entity.
    points:        GeoDataFrame with a monotonic DatetimeTZDtype index and
                   Point geometry column. Must already be in a METRIC CRS —
                   never pass raw EPSG:4326 here.
    quality_flags: Optional per-point quality series aligned to points.index.
                   True / 1.0 = high confidence; False / 0.0 = suspect.
    """

    def __init__(
        self,
        device_id: str,
        points: gpd.GeoDataFrame,
        quality_flags: Optional[pd.Series] = None,
    ) -> None:
        if points.empty:
            raise ValueError(f"[{device_id}] Empty point GeoDataFrame — cannot construct trajectory.")
        if len(points) < 2:
            raise ValueError(f"[{device_id}] At least 2 points required; got {len(points)}.")
        if "geometry" not in points.columns:
            raise ValueError(f"[{device_id}] GeoDataFrame missing 'geometry' column.")
        if points.crs is None or points.crs.is_geographic:
            raise ValueError(
                f"[{device_id}] CRS is geographic (or unset). Project to a metric CRS before "
                "constructing a TrajectoryObject — distances computed in degrees are meaningless."
            )
        if not isinstance(points.index, pd.DatetimeTZDtype.__mro__[0]) and not hasattr(points.index, "tz"):
            # Accept any DatetimeTZDtype-like index; warn if timezone-naive
            pass

        # Enforce chronological ordering — never mutate caller's frame
        if not points.index.is_monotonic_increasing:
            points = points.sort_index()

        # Drop exact duplicate timestamps (keep first)
        dup_mask = points.index.duplicated(keep="first")
        if dup_mask.any():
            points = points[~dup_mask]

        # Build and validate LineString in projected space
        coords = list(zip(points.geometry.x, points.geometry.y))
        line = LineString(coords)
        valid_geom = make_valid(line)
        if valid_geom.is_empty:
            raise ValueError(f"[{device_id}] Geometry is empty after make_valid — insufficient valid coords.")

        self.device_id: str = device_id
        self._points: gpd.GeoDataFrame = points
        self.geometry: LineString = valid_geom
        self.crs: str = str(points.crs)
        self.quality: pd.Series = (
            quality_flags.reindex(points.index, fill_value=True)
            if quality_flags is not None
            else pd.Series(True, index=points.index, dtype=bool)
        )

    # ------------------------------------------------------------------
    # Temporal properties
    # ------------------------------------------------------------------

    @property
    def start_time(self) -> pd.Timestamp:
        return self._points.index[0]

    @property
    def end_time(self) -> pd.Timestamp:
        return self._points.index[-1]

    def duration(self) -> float:
        """Total elapsed time in seconds."""
        return (self.end_time - self.start_time).total_seconds()

    # ------------------------------------------------------------------
    # Kinematic properties (requires metric CRS — asserted in __init__)
    # ------------------------------------------------------------------

    def path_length(self) -> float:
        """Total path length in the object's projected CRS units (metres for UTM)."""
        return self.geometry.length

    def mean_speed(self) -> float:
        """Mean speed in metres per second (path length / total duration).

        Note: this is a trip-level mean. It includes stationary periods.
        For mode-aware speed profiles, use segment_speeds() and filter by
        quality flag before aggregating.
        """
        d = self.duration()
        return self.path_length() / d if d > 0 else 0.0

    def segment_speeds(self) -> pd.Series:
        """Per-segment speed in m/s between consecutive GPS fixes.

        Returns a Series aligned to the point index (NaN for the first row).
        All distances are computed in the projected CRS — never in WGS84 degrees.
        """
        pts = self._points
        # Vectorised displacement using shifted geometry
        prev_geom = pts.geometry.shift(1)
        distances = pts.geometry.distance(prev_geom)          # metres (metric CRS)
        dt_seconds = pts.index.to_series().diff().dt.total_seconds()
        speeds = distances / dt_seconds                        # m/s
        speeds.iloc[0] = float("nan")                         # no predecessor for first point
        return speeds

    # ------------------------------------------------------------------
    # Quality helpers
    # ------------------------------------------------------------------

    def high_quality_ratio(self) -> float:
        """Fraction of points with quality flag True / >= 0.5."""
        q = self.quality.astype(float)
        return float((q >= 0.5).sum() / len(q)) if len(q) else 0.0

    def __repr__(self) -> str:
        return (
            f"TrajectoryObject(device={self.device_id!r}, "
            f"points={len(self._points)}, "
            f"duration={self.duration():.0f}s, "
            f"crs={self.crs!r})"
        )

Lazy-loading wrapper for large datasets

For trajectories with more than 500,000 points, or when loading a fleet of thousands of trajectories concurrently, materialising a full GeoDataFrame per object exhausts heap memory. The LazyTrajectory wrapper defers loading until an operation actually needs the points:

from __future__ import annotations
from pathlib import Path
import pandas as pd
import geopandas as gpd
import pyarrow.parquet as pq


class LazyTrajectory:
    """Memory-mapped trajectory backed by a Parquet file partition.

    The GeoDataFrame is loaded on first access and cached. For worker-pool
    use, call .evict() after processing to release the reference.
    """

    def __init__(self, parquet_path: Path, device_id: str, target_crs: str = "EPSG:32633") -> None:
        if not parquet_path.exists():
            raise FileNotFoundError(f"Parquet partition not found: {parquet_path}")
        self._path = parquet_path
        self.device_id = device_id
        self._target_crs = target_crs
        self._cache: gpd.GeoDataFrame | None = None

    def _load(self) -> gpd.GeoDataFrame:
        """Load and project the point GeoDataFrame from Parquet."""
        table = pq.read_table(self._path, filters=[("device_id", "=", self.device_id)])
        df = table.to_pandas()
        if df.empty:
            raise ValueError(f"No rows for device_id={self.device_id!r} in {self._path}.")
        gdf = gpd.GeoDataFrame(
            df,
            geometry=gpd.points_from_xy(df["lon"], df["lat"]),
            crs="EPSG:4326",
        )
        gdf = gdf.set_index(pd.to_datetime(df["timestamp"], utc=True)).sort_index()
        return gdf.to_crs(self._target_crs)   # project to metric CRS before any distance op

    @property
    def points(self) -> gpd.GeoDataFrame:
        if self._cache is None:
            self._cache = self._load()
        return self._cache

    def evict(self) -> None:
        """Release cached GeoDataFrame, allowing GC to reclaim memory."""
        self._cache = None

    def to_trajectory_object(self) -> TrajectoryObject:
        return TrajectoryObject(device_id=self.device_id, points=self.points)

Mathematical Grounding

Segment speed and the metric CRS requirement

Speed between two consecutive fixes is:

v_i = d(p_{i-1}, p_i) / (t_i - t_{i-1})

where d(·, ·) is Euclidean distance in a projected metric space. In a UTM zone, 1 unit = 1 metre, making the formula numerically correct. In EPSG:4326, the same formula produces degrees-per-second — a unit with no physical meaning. Shapely’s geometry.distance() is always Cartesian; it has no awareness of CRS and will silently compute the wrong value if you pass unprojected coordinates.

Temporal indexing and monotonicity

Trajectory analytics depend on the index being strictly monotonically increasing. Any downstream resample(), rolling(), or shift() operation on a DatetimeTZDtype index assumes this invariant. Violations produce:

• Negative dt values → negative speeds
• Out-of-order resampling bins → duplicated or missing time windows
• merge_asof alignment failures when joining to external signal sources (see Time-Series Synchronization Strategies)

Calibration and Parameter Tuning

Use the transport-mode table below to set velocity threshold ceilings at stage 2 of the pipeline. Thresholds should be enforced before object construction; points exceeding the ceiling become candidates for removal or quality-flag downgrade.

For Sampling Rate Optimization scenarios where points are thinned before object construction, velocity thresholds must account for the increased segment length — a 30-second gap between points will produce a legitimately high apparent speed even for a slow vehicle. Adjust ceilings proportionally or compute displacement over the full gap rather than per-second speed.

Integration and Compatibility

TrajectoryObject is the entry point for all downstream analytical stages:

• Speed and acceleration profiling — segment_speeds() feeds directly into speed-acceleration profiling pipelines; filter by quality before computing acceleration to avoid artifact spikes from low-confidence fixes.
• Stay-point detection — pass _points (the underlying GeoDataFrame) to DBSCAN-based stay detection; the temporal index aligns naturally with dwell-duration thresholds.
• Directionality analysis — derive bearing from consecutive projected coordinates using math.atan2(dy, dx) on the metric point pairs; see directionality and turn analysis for U-turn and heading-shift detection built on top of this object.
• GeoParquet export — use gdf.to_parquet(path, geometry_encoding="WKB") from the _points frame; this produces Parquet-native spatial files that DuckDB’s spatial extension reads without conversion. Partition by device_id hash and date-truncated timestamp for efficient range scans.
• movingpandas compatibility — TrajectoryObject._points is a valid input to mpd.Trajectory(gdf, traj_id=device_id); the metric CRS and monotonic index already satisfy movingpandas preconditions.

Storage Selection

For datasets exceeding 100 million points, combine GeoHash or H3 spatial partitioning with Z-order (Morton) curve indexing on (spatial_partition, date_bucket) to reduce I/O on range scans by 60–80% compared to unpartitioned Parquet.

FAQ

Should a trajectory store a LineString or a GeoDataFrame of points? Both, in layers. Keep the raw point GeoDataFrame for temporal resampling and per-point attribute access; derive the LineString from it for spatial operations like path length and spatial join. Storing only the LineString loses timestamp granularity and makes velocity profiling impossible without re-parsing.

How do I handle gaps or signal loss inside a trajectory? Do not interpolate blindly. Insert explicit NaN-coordinate rows at the gap boundary so downstream consumers can detect the discontinuity. Segment the trajectory at any gap exceeding your transport-mode dwell threshold before computing kinematic metrics — a 10-minute gap in an urban vehicle trace is a parking event, not a tunnel.

Why does mean_speed() differ from averaging segment_speeds()? mean_speed() divides total path length by total elapsed time, weighting by distance. Averaging segment_speeds() weights by number of segments. A trajectory with many short slow segments and one long fast segment will show a higher average when computed from segments. Use mean_speed() for trip-level summaries; use segment_speeds().mean() cautiously and only after filtering stationary periods.

When should I switch to LazyTrajectory? When a single trajectory exceeds roughly 500,000 points, or when you must process more than 1,000 trajectories concurrently in a worker pool. In-memory GeoDataFrame objects at that scale exhaust heap on standard cloud workers. The lazy wrapper loads on demand and evicts after processing.

Can I compute distances in EPSG:4326 for a rough estimate? No. Shapely’s distance() is always Cartesian — it has no CRS awareness. Calling geometry.distance() on EPSG:4326 coordinates gives degree-space distances that are physically meaningless and cannot be converted to metres without knowing the latitude. Always project to a metric CRS first.

Related

• How to Structure Trajectory Data in GeoPandas — extended examples on spatial joins, temporal resampling, and attribute alignment
• GPS Precision & Error Handling — upstream validation this page depends on
• Coordinate Reference System Mapping — automated projection routing and CRS validation workflows
• Sampling Rate Optimization — downsampling strategies that feed into trajectory object construction
• Time-Series Synchronization Strategies — aligning multi-sensor timestamps before building trajectory objects

Back to Spatiotemporal Data Foundations & Structures