GPS Precision & Error Handling

Raw GPS telemetry is a noisy signal, not ground truth — even under open-sky conditions, consumer and commercial receivers exhibit positional variance of 2–15 metres that must be systematically classified and corrected before any downstream trajectory analysis can be trusted.

Prerequisites Checklist

Before implementing a GPS cleaning pipeline, confirm your environment meets these requirements. This builds directly on the foundations covered in Spatiotemporal Data Foundations & Structures.

• Python ≥ 3.10 with pandas ≥ 2.0 or polars ≥ 0.20 for vectorized telemetry operations
• geopandas ≥ 0.14 for spatial vectorization and shapely ≥ 2.0 for geometry validation
• pyproj ≥ 3.6 for datum transformations — always import with always_xy=True to avoid axis-order bugs
• scipy ≥ 1.11 for Savitzky-Golay filtering and interpolation
• filterpy (optional) for Kalman filter state estimation
• Input schema: each row must have lat (float64), lon (float64), ts (UTC datetime64), hdop (float32), pdop (float32), fix_type (int8), sat_count (int8). Upstream parsing stage must complete before this pipeline runs.
• Projection context: GPS outputs EPSG:4326 geographic coordinates. All distance, velocity, and acceleration calculations require a metric projected CRS — this is flagged explicitly in every code block below.

GPS Error Taxonomy

Positional inaccuracy stems from both deterministic and stochastic sources. Classifying the error type determines which mitigation to apply — over-engineering a filter (applying Kalman smoothing to atmospheric bias) will mask genuine elevation changes, while under-filtering multipath corrupts turn-detection algorithms.

Deterministic Pipeline Overview

The six-stage sequence below is the standard production order. Each stage must preserve data lineage (keep original values in _raw columns), log dropped record counts, and be fully vectorized — no Python-level row iteration.

1. Parse & validate — extract lat, lon, ts, hdop, sat_count, fix_type from NMEA/GPX; enforce UTC, standardize sentinels (0.0, 999.0) to pd.NA.
2. DOP & fix-type filter — reject records where hdop > 4.0 OR sat_count < 6 OR fix_type != 3; write rejection reason to a _qa_flag column.
3. Temporal alignment — enforce monotonic timestamps, detect and remove duplicates, resample to target cadence with linear interpolation, flag interpolated rows.
4. Kinematic thresholding — compute per-row speed and acceleration in a metric CRS; flag records violating mode-specific bounds; interpolate across flagged gaps ≤ 30 s, split trajectory at gaps > 30 s.
5. Spatial smoothing — apply Savitzky-Golay (batch) or Kalman (real-time) filter to projected coordinates to remove residual noise while preserving turn geometry.
6. CRS projection & output — transform to target metric CRS, validate bounding box, write GeoParquet with embedded metadata, cleaning version, and error flags.

Implementation Walkthrough

The function below implements stages 2–5 end-to-end. All distance and velocity calculations are performed in a metric projected CRS (EPSG:32633 as default) — never on raw EPSG:4326 coordinates. Adjust the UTM zone to your data’s geographic extent.

from __future__ import annotations

import logging
from typing import Literal

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

logger = logging.getLogger(__name__)

TransportMode = Literal["pedestrian", "cyclist", "vehicle", "unknown"]

# Mode-specific kinematic thresholds (metric CRS, m/s and m/s²)
KINEMATIC_LIMITS: dict[TransportMode, dict[str, float]] = {
    "pedestrian": {"max_speed": 2.8, "max_accel": 3.0},
    "cyclist":    {"max_speed": 13.0, "max_accel": 4.5},
    "vehicle":    {"max_speed": 50.0, "max_accel": 9.0},
    "unknown":    {"max_speed": 50.0, "max_accel": 9.0},
}

GAP_SPLIT_SECONDS = 30  # split trajectory at gaps longer than this


def clean_gps_telemetry(
    df: pd.DataFrame,
    mode: TransportMode = "vehicle",
    target_crs_epsg: int = 32633,   # UTM zone 33N — adjust to your region
    hdop_max: float = 4.0,
    sat_min: int = 6,
    resample_hz: float = 1.0,
    sg_window: int = 11,
    sg_poly: int = 3,
) -> pd.DataFrame:
    """
    Apply a deterministic GPS cleaning pipeline to a single-device trajectory.

    Parameters
    ----------
    df : DataFrame with columns lat (float64), lon (float64), ts (datetime64[ns, UTC]),
         hdop (float32), sat_count (int8), fix_type (int8).
    mode : Transport mode — controls kinematic threshold selection.
    target_crs_epsg : EPSG code of the metric projected CRS for distance/velocity computation.
                      NEVER use 4326 here — distances in degrees are not metric.
    hdop_max : Maximum acceptable HDOP; records above this are filtered.
    sat_min : Minimum satellite count; records below this are filtered.
    resample_hz : Target resampling frequency in Hz (1.0 = 1-second cadence).
    sg_window : Savitzky-Golay window length (must be odd, > sg_poly).
    sg_poly : Savitzky-Golay polynomial order.

    Returns
    -------
    Cleaned DataFrame with additional columns:
      x_m, y_m        — projected coordinates in target_crs_epsg (metres)
      speed_ms        — instantaneous speed (m/s), computed in metric CRS
      accel_ms2       — instantaneous acceleration (m/s²)
      _qa_flag        — pipe-separated string of applied QA flags
      _interpolated   — True for rows created by resampling interpolation
      segment_id      — integer ID; incremented at trajectory gaps > GAP_SPLIT_SECONDS
    """
    required = {"lat", "lon", "ts", "hdop", "sat_count", "fix_type"}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Input DataFrame missing columns: {missing}")

    if df.empty:
        logger.warning("clean_gps_telemetry received an empty DataFrame; returning empty.")
        return df.copy()

    out = df.copy()

    # ── Sentinel normalisation ────────────────────────────────────────────────
    # Replace GPS-protocol sentinel values with pd.NA so they propagate cleanly
    for col in ("hdop", "pdop"):
        if col in out.columns:
            out[col] = out[col].replace({0.0: pd.NA, 999.0: pd.NA})
    out["sat_count"] = out["sat_count"].replace({0: pd.NA})

    # ── Stage 2: DOP & fix-type filter ───────────────────────────────────────
    out["_qa_flag"] = ""
    dop_fail = (
        (out["hdop"].isna() | (out["hdop"] > hdop_max))
        | (out["sat_count"].isna() | (out["sat_count"] < sat_min))
        | (out["fix_type"] != 3)
    )
    out.loc[dop_fail, "_qa_flag"] += "dop_fail|"
    n_dop = dop_fail.sum()
    logger.info("DOP filter: dropped %d / %d records (%.1f%%)", n_dop, len(out), 100 * n_dop / len(out))
    out = out.loc[~dop_fail].copy()

    if len(out) < 4:
        logger.warning("Fewer than 4 records after DOP filter; returning empty cleaned frame.")
        return out

    # ── Stage 3: Temporal alignment ──────────────────────────────────────────
    out["ts"] = pd.to_datetime(out["ts"], utc=True)
    out = out.sort_values("ts").drop_duplicates(subset="ts")
    out = out.set_index("ts")

    interval_s = int(1.0 / resample_hz)
    full_index = pd.date_range(out.index.min(), out.index.max(), freq=f"{interval_s}s", tz="UTC")
    out = out.reindex(full_index)
    out["_interpolated"] = out["lat"].isna()

    # Linear interpolation for numeric columns only
    num_cols = out.select_dtypes(include="number").columns
    out[num_cols] = out[num_cols].interpolate(method="time", limit=30, limit_direction="forward")
    out = out.reset_index().rename(columns={"index": "ts"})
    out["_qa_flag"] = out["_qa_flag"].fillna("")

    # ── Stage 4: Kinematic thresholding ──────────────────────────────────────
    # All speed/acceleration computed in metric CRS — NEVER on WGS84 degrees
    transformer = Transformer.from_crs("EPSG:4326", f"EPSG:{target_crs_epsg}", always_xy=True)
    out["x_m"], out["y_m"] = transformer.transform(out["lon"].values, out["lat"].values)

    dt_s = out["ts"].diff().dt.total_seconds().fillna(1.0).clip(lower=0.001)
    dx = out["x_m"].diff().fillna(0.0)
    dy = out["y_m"].diff().fillna(0.0)
    out["speed_ms"] = np.sqrt(dx**2 + dy**2) / dt_s
    out["accel_ms2"] = out["speed_ms"].diff().fillna(0.0) / dt_s

    limits = KINEMATIC_LIMITS[mode]
    kin_fail = (out["speed_ms"] > limits["max_speed"]) | (out["accel_ms2"].abs() > limits["max_accel"])
    out.loc[kin_fail, "_qa_flag"] += "kin_fail|"
    n_kin = kin_fail.sum()
    logger.info("Kinematic filter: flagged %d records for mode '%s'", n_kin, mode)

    # Interpolate across short kinematic gaps; split trajectory at long gaps
    out.loc[kin_fail, ["lat", "lon", "x_m", "y_m"]] = np.nan
    out[["lat", "lon", "x_m", "y_m"]] = out[["lat", "lon", "x_m", "y_m"]].interpolate(
        method="linear", limit=int(GAP_SPLIT_SECONDS * resample_hz)
    )

    gap_seconds = out["ts"].diff().dt.total_seconds().fillna(0)
    out["segment_id"] = (gap_seconds > GAP_SPLIT_SECONDS).cumsum().astype(int)

    # ── Stage 5: Spatial smoothing (Savitzky-Golay) ──────────────────────────
    # Smooth projected x_m / y_m; convert back to lat/lon for output
    # SG preserves peak values better than a moving average
    if len(out) >= sg_window:
        out["x_m"] = savgol_filter(out["x_m"].fillna(method="ffill"), sg_window, sg_poly)
        out["y_m"] = savgol_filter(out["y_m"].fillna(method="ffill"), sg_window, sg_poly)
        inv_transformer = Transformer.from_crs(f"EPSG:{target_crs_epsg}", "EPSG:4326", always_xy=True)
        out["lon"], out["lat"] = inv_transformer.transform(out["x_m"].values, out["y_m"].values)
    else:
        logger.warning("Trajectory too short for Savitzky-Golay (need ≥ %d points); skipping smoothing.", sg_window)

    # Recompute speed after smoothing
    dt_s = out["ts"].diff().dt.total_seconds().fillna(1.0).clip(lower=0.001)
    dx = out["x_m"].diff().fillna(0.0)
    dy = out["y_m"].diff().fillna(0.0)
    out["speed_ms"] = np.sqrt(dx**2 + dy**2) / dt_s

    return out

Geometric and Mathematical Grounding

DOP geometry. The Dilution of Precision family (HDOP, VDOP, PDOP) quantifies how satellite geometry amplifies range-measurement errors. For a simplified 2D case with two satellites at angles θ₁ and θ₂ from the receiver, the horizontal position uncertainty σ_h relates to the pseudorange noise σ_p by:

σ_h = HDOP × σ_p

HDOP is computed from the diagonal of the inverse of the satellite-receiver geometry matrix H (the (HᵀH)⁻¹ normal equations). When satellites cluster near the horizon or all appear in the same azimuthal quadrant, det(H) shrinks and HDOP grows, amplifying any pseudorange error. An HDOP of 1.0 means satellite geometry is ideal; above 4.0, positional error routinely exceeds 15 m on consumer receivers.

Savitzky-Golay vs. moving average. Both are linear FIR filters, but SG fits a least-squares polynomial of degree p across the window rather than weighting coefficients equally. This preserves the amplitude and shape of peaks (velocity spikes from genuine sharp turns) while attenuating high-frequency noise. For trajectory smoothing, a window of 11 points at 1 Hz with polynomial order 3 is a practical starting point that attenuates noise above ~0.1 Hz without smearing turn geometry.

Kalman filter as an alternative. For real-time GPS cleaning, a constant-velocity Kalman model uses state vector [x, y, vx, vy] and a process noise matrix Q tuned to the expected platform dynamics. The Kalman gain K = P·Hᵀ·(H·P·Hᵀ + R)⁻¹ balances prediction trust (P, process covariance) against measurement trust (R, receiver noise covariance). Compared to SG, Kalman is causal and updates incrementally — preferable for streaming pipelines. For batch analytics, SG generally produces cleaner output at lower computational cost.

Calibration & Parameter Tuning

Calibrate max_speed and max_accel against a ground-truth subset of your actual fleet or participant data, not theoretical maximums. Use the 99.9th percentile of observed values as the ceiling to avoid masking genuine edge events while still filtering multipath jumps.

For multi-modal datasets (rideshare, micro-mobility, transit), run separate cleaning passes per mode inferred from speed-acceleration profiling rather than using a single “unknown” configuration, which will either under-filter pedestrian data or over-filter vehicle data.

Integration & Compatibility

The output of this pipeline feeds directly into several downstream stages:

• Coordinate Reference System Mapping: The projected x_m/y_m columns, once validated, feed into spatial join operations that require a consistent CRS. Use the same target_crs_epsg throughout your pipeline to avoid re-projection drift.
• Sampling Rate Optimization: After cleaning, resample_hz has already standardized cadence. Downsampling passes downstream can apply Douglas-Peucker or step-frequency decimation without first fighting irregular timestamps.
• Trajectory Object Design Patterns: Wrap the cleaned DataFrame in a Trajectory object immediately after this pipeline. The segment_id column maps directly to sub-trajectory splitting logic.
• Stay-Point Detection: Stay-point algorithms consume cleaned, projected coordinates. Applying DBSCAN to raw or lightly-cleaned GPS data produces spurious clusters from multipath noise rather than genuine activity locations.
• Kalman/HMM map-matching: The speed_ms and segment_id outputs provide the velocity priors and trajectory boundaries that Hidden Markov Model map-matchers require for emission probability computation.
• Handling GPS drift in raw trajectory logs: For deeper treatment of the specific drift patterns this pipeline handles, see the dedicated deep-dive page.

Validation & Continuous Monitoring

After running the cleaning pipeline, verify output quality with these checks before committing data to your analytics store:

def validate_cleaned_trajectory(df: pd.DataFrame, target_crs_epsg: int = 32633) -> dict:
    """
    Post-pipeline sanity checks. Returns a dict of metric names → values.
    Raises AssertionError if critical invariants are violated.
    """
    assert not df.empty, "Cleaned trajectory is empty"
    assert df["ts"].is_monotonic_increasing, "Timestamps not monotonic after cleaning"
    assert df["lat"].between(-90, 90).all(), "Latitude out of range"
    assert df["lon"].between(-180, 180).all(), "Longitude out of range"

    # Speed should be physically plausible — flag but don't raise for outliers
    p999_speed = df["speed_ms"].quantile(0.999)
    drop_rate = df["_qa_flag"].str.contains("dop_fail|kin_fail", na=False).mean()

    return {
        "n_records": len(df),
        "n_segments": df["segment_id"].nunique(),
        "p999_speed_ms": round(p999_speed, 2),
        "interpolated_pct": round(100 * df["_interpolated"].mean(), 1),
        "qa_flagged_pct": round(100 * drop_rate, 1),
        "crs_epsg": target_crs_epsg,
    }

Establish baseline metrics per device cohort and trigger automated alerts when:

• qa_flagged_pct exceeds 3σ of the historical device mean (indicates firmware update, mounting change, or environment shift)
• p999_speed_ms exceeds the mode’s max_speed after cleaning (filter threshold may need recalibration)
• n_segments per trip increases sharply (urban deployments moving into canyon areas)

For rigorous accuracy benchmarking, collect RTK-GPS ground-truth tracks on the same routes and compute Hausdorff distance between cleaned trajectories and reference paths. A well-tuned pipeline should achieve Hausdorff distances below 8 m in open-sky environments and below 20 m in dense urban cores.

FAQ

What HDOP threshold should I use for mobility data?

Filter records where hdop > 4.0 as a general starting point. Multi-constellation receivers (GPS + GLONASS + Galileo) can tolerate HDOP up to 5.0 without noticeable accuracy loss. Pedestrian analysis in urban canyons may require tighter thresholds (hdop < 2.5). Always combine HDOP filtering with satellite count (≥ 6) and 3D fix type validation — each check catches different failure modes.

Should I drop GPS jumps or interpolate across them?

Interpolate for gaps under 30 seconds; split the trajectory for longer outages. Dropping a GPS jump leaves a gap that creates phantom velocity spikes when the next valid record arrives. Instead, flag the offending record as corrected, interpolate across the gap, and preserve the original values in _raw columns for auditability.

Why do my velocity calculations spike near tunnel exits?

GPS receivers lose satellite lock inside tunnels and re-acquire on exit. The coordinate jump between the last pre-tunnel point and the first post-tunnel point is measured over the full tunnel traversal time, producing phantom velocities that can exceed 200 m/s for a car. Detect tunnel-exit events via a sudden HDOP spike followed by a rapid drop; flag the first 2–3 post-acquisition points as unreliable and apply kinematic filtering before any velocity derivation.

Can I compute distances and speeds directly on EPSG:4326 coordinates?

No. EPSG:4326 is an angular coordinate system; treating degree differences as metric distances introduces errors that grow with latitude and trajectory length. Always project to a metric CRS before computing distances, speeds, or accelerations. The haversine formula gives correct great-circle distances but is expensive at scale and cannot produce metric-accurate acceleration derivatives.

What is the right smoothing filter for real-time vs batch GPS cleaning?

For batch processing, Savitzky-Golay filtering preserves peak velocity and turn geometry better than a simple moving average. For real-time streaming, a Kalman filter is preferred because it is causal and can model the receiver’s noise covariance. A fixed-gain alpha-beta tracker is a lightweight alternative when full Kalman state estimation is too expensive for your latency budget.

Related

• Handling GPS Drift in Raw Trajectory Logs
• Coordinate Reference System Mapping
• Sampling Rate Optimization
• Trajectory Object Design Patterns
• Speed & Acceleration Profiling

Back to Spatiotemporal Data Foundations & Structures