Speed Profiling from Raw GPS Coordinates

Raw GPS telemetry does not contain a reliable speed column. Device-reported velocity — when present — is derived from the GPS chipset’s own Doppler estimate, which lags position fixes by one to three seconds and silently resets to zero during satellite reacquisition. Naively trusting it produces odometer errors of 5–15% over a full trip. Computing speed independently from sequential coordinate-timestamp pairs gives you an auditable, reproducible velocity signal you can validate against OBD-II wheel encoder data and tune per vehicle class.

This page covers the full implementation path: from schema validation through vectorized haversine computation, Savitzky-Golay smoothing, and routing engine integration. It is part of the Trajectory Analysis & Map Matching Techniques discipline.

Why Naive Approaches Fail

The two most common mistakes are dividing raw coordinate deltas by elapsed seconds without geodesic correction, and applying a single global smoothing pass across device boundaries.

Planar distance errors accumulate at scale. At 45° latitude, one degree of latitude spans 111 km but one degree of longitude spans only 78 km. A simple Euclidean distance in degree-space will underestimate east-west travel by up to 30% at high latitudes, generating persistent speed underreads for any fleet operating above 40°N or below 40°S.

Cross-device smoothing artefacts occur when a rolling filter or Savitzky-Golay window straddles the last ping of one vehicle and the first ping of another. The artificial speed spike at the join can exceed 500 km/h and will not be caught by physical plausibility checks if the window is small.

Silent dt=0 divisions happen whenever two records share an identical timestamp — common in telematics units that buffer pings and flush them with the same epoch when connectivity resumes. NumPy divides by zero silently, producing inf values that propagate through all downstream aggregations.

Prerequisites

Before proceeding, ensure timestamp synchronization has been applied across all devices. Clock drift between telematics units produces artificial speed spikes at the point where synchronized records are stitched together.

Step-by-Step Workflow

1. Ingestion and Schema Validation

Parse the raw stream, enforce strict dtype contracts, and remove records that will corrupt downstream arithmetic.

import pandas as pd
import numpy as np

def load_and_validate(path: str) -> pd.DataFrame:
    df = pd.read_parquet(path)

    # Enforce required columns
    required = {"device_id", "lat", "lon", "timestamp"}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing columns: {missing}")

    # Type coercion
    df["lat"] = pd.to_numeric(df["lat"], errors="coerce")
    df["lon"] = pd.to_numeric(df["lon"], errors="coerce")
    df["timestamp"] = pd.to_datetime(df["timestamp"], utc=True)

    # Coordinate range validation
    mask = (
        df["lat"].between(-90, 90)
        & df["lon"].between(-180, 180)
        & df["lat"].notna()
        & df["lon"].notna()
        & df["timestamp"].notna()
    )
    dropped = (~mask).sum()
    if dropped:
        print(f"Dropped {dropped} records failing coordinate/timestamp checks")
    df = df[mask].copy()

    # Optional: HDOP gate — drop fixes with poor satellite geometry
    if "hdop" in df.columns:
        df = df[df["hdop"] <= 10.0]

    return df

Expected output shape: same columns, fewer rows. The HDOP gate alone typically removes 0.5–3% of records in urban environments; higher fractions signal a hardware fault worth flagging to the operations team.

2. Temporal Sorting and Monotonicity Verification

Group by device and sort strictly by timestamp. Non-monotonic sequences usually indicate GPS cold-start resets or network reordering — both require explicit handling before distance computation.

def sort_and_verify(df: pd.DataFrame) -> pd.DataFrame:
    df = df.sort_values(["device_id", "timestamp"]).reset_index(drop=True)

    # Flag backward time jumps per device
    df["dt_seconds"] = (
        df.groupby("device_id")["timestamp"]
        .diff()
        .dt.total_seconds()
    )
    backward = df["dt_seconds"] < 0
    if backward.any():
        print(f"Warning: {backward.sum()} backward time jumps detected — "
              "check for GPS cold-start resets")
        # Nullify distances that would span the reset
        df.loc[backward, "dt_seconds"] = np.nan

    return df

3. Vectorized Haversine Distance

The haversine formula computes the great-circle distance between two points on a sphere. It introduces a maximum error of ~0.5% compared to the Vincenty ellipsoidal formula, which is acceptable for velocity at fleet telematics sampling rates. Vectorizing over the whole DataFrame avoids Python-level loops on datasets of tens of millions of rows.

The formula for angular distance c between consecutive coordinates (φ₁, λ₁) and (φ₂, λ₂) in radians:

a = sin²(Δφ/2) + cos(φ₁)·cos(φ₂)·sin²(Δλ/2)
c = 2·arctan2(√a, √(1−a))
distance = R·c          (R = 6 371 000 m)

def add_haversine_distance(df: pd.DataFrame) -> pd.DataFrame:
    # Shift within each device group to get the previous coordinate
    prev_lat = df.groupby("device_id")["lat"].shift(1)
    prev_lon = df.groupby("device_id")["lon"].shift(1)

    lat1 = np.radians(prev_lat)
    lat2 = np.radians(df["lat"])
    dlat = lat2 - lat1
    dlon = np.radians(df["lon"]) - np.radians(prev_lon)

    a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))

    df["distance_m"] = 6_371_000 * c  # meters
    return df

The shift operation is applied inside groupby so that the first ping of each device correctly returns NaN rather than measuring the distance from the last ping of the previous device in the sorted DataFrame.

4. Velocity Derivation and Unit Conversion

Divide distance by elapsed seconds. Apply a jitter gate before the division: displacements below 0.5 m are indistinguishable from GPS coordinate noise at rest and should not generate a non-zero speed reading.

def derive_velocity(df: pd.DataFrame) -> pd.DataFrame:
    # Jitter gate: sub-0.5 m displacements are noise, not movement
    df.loc[df["distance_m"] < 0.5, "distance_m"] = np.nan

    # Velocity in m/s, then convert to km/h
    df["speed_kmh"] = (df["distance_m"] / df["dt_seconds"]) * 3.6

    # Mask division-by-zero (duplicate timestamps) and NaN propagation
    df.loc[df["dt_seconds"].le(0) | df["dt_seconds"].isna(), "speed_kmh"] = np.nan

    # Physical plausibility cap: >220 km/h is implausible for road vehicles
    df.loc[df["speed_kmh"] > 220, "speed_kmh"] = np.nan

    return df

The 220 km/h cap should be adjusted for specialist fleets: use 130 km/h for heavy goods vehicles (HGVs), 300+ km/h for rail, and 600+ km/h for airborne assets.

5. Signal Smoothing per Device Group

A Savitzky-Golay filter fits a low-degree polynomial to a sliding window, preserving peaks (sharp braking, hard acceleration) better than a simple moving average. Apply it strictly within each device group to avoid the cross-device boundary artefacts described above.

from scipy.signal import savgol_filter

def smooth_speed(
    df: pd.DataFrame,
    window_length: int = 7,
    poly_order: int = 2
) -> pd.DataFrame:
    """
    window_length: must be odd and > poly_order.
                   Increase for smoother highways; decrease for stop-and-go routes.
    poly_order:    2 preserves acceleration ramps; 3 better fits sharp cornering.
    """
    def _smooth(series: pd.Series) -> pd.Series:
        valid_mask = series.notna()
        valid_vals = series[valid_mask].values
        if len(valid_vals) < window_length:
            return series  # not enough points — return unsmoothed
        smoothed = savgol_filter(
            valid_vals,
            window_length=window_length,
            polyorder=poly_order,
            mode="nearest"  # avoids edge artefacts at trip start/end
        )
        result = series.copy()
        result[valid_mask] = np.clip(smoothed, 0, None)  # no negative speeds
        return result

    df["speed_kmh_smoothed"] = df.groupby("device_id")["speed_kmh"].transform(_smooth)
    return df

Mathematical Model: Haversine vs. Vincenty

For most fleet telematics applications the haversine formula is sufficient, but it is worth knowing where the error comes from. The haversine treats the Earth as a perfect sphere of radius 6 371 000 m. The actual Earth is an oblate spheroid: the equatorial radius is 6 378 137 m and the polar radius is 6 356 752 m (WGS84 ellipsoid). The maximum relative error in haversine distance is approximately:

ε_max ≈ (a − b) / a ≈ 0.33%

where a is the equatorial radius and b the polar radius. Over a 10 km segment this is 33 m — well within GPS horizontal accuracy for consumer-grade chipsets (±5–15 m CEP). For sub-meter accuracy requirements (precision agriculture, surveying), use the Vincenty formula via the geographiclib package.

Numerical stability note: when consecutive points are very close together (< 1 m), a approaches zero and arctan2(√a, √(1−a)) approaches zero as expected. No log-space transformation is needed for haversine; the formula is numerically stable for all valid coordinate pairs on Earth.

Routing Engine Integration Notes

OSRM

OSRM expects coordinates in [longitude, latitude] order (GeoJSON convention), not [latitude, longitude]. Passing coordinates in the wrong order produces a valid HTTP 200 response with a plausible-looking but geographically wrong route. Always verify by checking that the snapped point is within 50 m of your input coordinate.

Speed profiles feed into the OSRM match endpoint via the radiuses and timestamps parameters. Providing timestamps (Unix seconds, integer) allows OSRM to compute expected travel times per edge and reject GPS fixes that imply impossible speeds for the road class.

import requests

def osrm_match(coords: list[tuple[float, float]], timestamps: list[int]) -> dict:
    """
    coords: list of (lat, lon) tuples
    timestamps: parallel list of Unix timestamps (integer seconds)
    """
    # OSRM requires lon,lat order
    coord_str = ";".join(f"{lon},{lat}" for lat, lon in coords)
    ts_str = ";".join(str(t) for t in timestamps)
    url = (
        f"http://router.project-osrm.org/match/v1/driving/{coord_str}"
        f"?timestamps={ts_str}&overview=full&geometries=geojson"
    )
    resp = requests.get(url, timeout=10)
    resp.raise_for_status()
    return resp.json()

Valhalla

Valhalla’s trace_attributes endpoint accepts a shape_match parameter. Set it to "map_snap" for sparse traces (< 1 Hz) and "edge_walk" for dense traces (≥ 1 Hz). Supplying a speed field per trace point causes Valhalla to weight the emission probability by the ratio of observed speed to the posted speed limit for the candidate edge — this significantly improves match quality on parallel roads with different speed limits.

GraphHopper

The GraphHopper Map Matching API accepts a GPX file with <trkpt> elements. Include <time> elements for each point; without them GraphHopper assumes constant speed and produces worse snapping results. The gps_accuracy parameter (default 40 m) controls the search radius around each GPS fix; reduce it to 15 m in dense urban networks to prevent snapping to parallel streets.

Operational Troubleshooting

Deployment Checklist

All input coordinates confirmed in WGS84 (EPSG:4326) — not EPSG:3857 or any local CRS
Timestamps converted to UTC-aware datetime64[ns, UTC] before sorting
Haversine shift scoped inside groupby("device_id") — not a bare DataFrame shift
dt_seconds <= 0 rows masked before velocity division
Jitter gate (< 0.5 m) applied before speed calculation
Physical plausibility cap set to the correct value for each vehicle class in the fleet
Savitzky-Golay window applied per device group with a small-group fallback
speed_kmh_smoothed compared against OBD-II readings on at least one validation trip
Gap-interpolated rows flagged for downstream consumers
Pipeline benchmarked against largest expected daily batch (target: < 60 s per 10 M rows on 8 cores)

Validation and Quality Assurance

Before deploying to production, run these automated checks:

Physical plausibility tests: assert df["speed_kmh_smoothed"].max() <= vehicle_class_limit per fleet segment. Log percentile distributions (p50, p95, p99) daily — a shift in p99 often signals a new device model with different firmware behavior before it causes operational issues.

Cross-sensor consistency: Compare GPS-derived speed_kmh_smoothed against OBD-II wheel speed for a representative sample of trips. The Pearson correlation should exceed 0.97 for well-tuned pipelines; values below 0.90 indicate a systematic issue with HDOP thresholds or tyre data.

Stop detection agreement: Periods where speed_kmh_smoothed < 2.0 for more than 60 consecutive seconds should align with independently detected stop events. Discrepancies larger than 30 seconds indicate either the speed pipeline or the stop detector needs retuning.

Integration with map matching: After running a Hidden Markov Model map matching pass, compare matched-edge posted speed limits against speed_kmh_smoothed. Segments where GPS-derived speed exceeds the posted limit by more than 30% on a residential road are usually multipath artefacts, not genuine speeding — they merit a second-pass smoothing with a wider window.

For heading consistency, verify that sharp speed changes (> 20 km/h per second) coincide with plausible heading changes. A large speed change with no heading change indicates a GPS coordinate jump rather than a genuine braking event.

Related

• Calculating Instantaneous vs. Average Speed from GPS Traces — extends this pipeline with time-windowed aggregation and segment-level averages
• Directionality & Heading Synchronization — adds bearing vectors to the velocity output for routing engine integration
• Hidden Markov Model Map Matching in Python — road-constrains smoothed speed vectors for legally compliant trip reports
• Kalman Filtering for GPS Noise Reduction — alternative noise suppression for pipelines requiring velocity continuity across signal gaps
• Outlier Removal in Raw Telematics Streams — upstream coordinate cleaning that reduces jitter before the haversine step

Parent topic: Trajectory Analysis & Map Matching Techniques