Handling GPS drift in raw trajectory logs

Correcting GPS drift requires a deterministic four-stage pipeline: temporal regularization → kinematic thresholding → spatial smoothing → validation. Raw GNSS receivers routinely produce 2–10 m of positional noise under open-sky conditions, and this compounds into unrealistic zigzag paths, phantom stops, and inflated distance metrics when left uncorrected. The most reliable strategy combines physics-aware constraints (maximum speed, maximum acceleration) with a Savitzky-Golay polynomial smoother, which preserves genuine route geometry while suppressing high-frequency jitter — making the output suitable for downstream mobility analytics, route optimization, and urban flow modeling.

Why this happens

GPS drift is not a single failure mode. It emerges from four concurrent sources: multipath signal reflection in urban canyons and dense foliage, ionospheric and tropospheric propagation delays, satellite geometry degradation expressed as high HDOP or PDOP values, and receiver clock jitter on low-cost hardware. When a device logs at 1 Hz or higher, these micro-errors accumulate as high-frequency spatial noise in the output stream.

The full taxonomy of GNSS error sources, their magnitudes, and their interaction with coordinate reference system mapping is covered in the parent GPS precision and error handling reference. The key production consequence is that uncorrected drift corrupts trip segmentation, dwell-time calculations, and instantaneous speed calculations from discrete GPS points, all of which sit downstream of this cleaning stage.

Core mitigation pipeline

The four stages below follow the order in which each failure mode must be resolved — you cannot smooth before you have a uniform time axis, and you cannot kinematically filter before projecting to a metric space.

GPS Drift Correction Pipeline Four sequential stages: Temporal Regularization, Kinematic Thresholding, Spatial Smoothing, and Validation, connected by arrows. 1. Temporal Regularization 2. Kinematic Thresholding 3. Spatial Smoothing 4. Validation & Metrics
  1. Temporal regularization. Raw logs often arrive at irregular intervals due to OS scheduling, network buffering, or device power-saving states. Resample or linearly interpolate to a fixed cadence (1 s or 2 s is typical) to stabilize the derivative calculations that follow.
  2. Kinematic thresholding. Apply hard constraints based on vehicle physics. Urban passenger vehicles rarely exceed 120 km/h or 4.5 m/s² lateral acceleration. Points violating these limits are flagged as drift artifacts and masked before smoothing.
  3. Spatial smoothing. Project WGS84 coordinates into a local metric CRS (a UTM zone derived from the data centroid), apply a windowed polynomial filter, and reproject back. The metric projection is mandatory — computing distances or velocities in angular degree-space introduces systematic scale distortion that the filter will amplify rather than remove. See CRS transformation best practices for the UTM zone selection logic.
  4. Validation. Confirm the correction did not over-smooth by checking distance delta, turn-angle preservation, and stop-detection accuracy against ground-truth events or telematics logs.

Production-ready Python implementation

The function below applies all four stages. It handles WGS84 input, derives the correct UTM zone from the data centroid, masks kinematic outliers, interpolates across gaps, and applies Savitzky-Golay smoothing. The pyproj transformer uses always_xy=True to guarantee consistent axis order regardless of CRS definition — a common source of silent coordinate swaps. The SciPy Savitzky-Golay implementation provides optimized C-backed filtering that outperforms naive moving averages for trajectory data.

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


def correct_gps_drift(
    df: pd.DataFrame,
    time_col: str = "timestamp",
    lat_col: str = "lat",
    lon_col: str = "lon",
    max_speed_kmh: float = 120.0,
    max_accel_ms2: float = 4.5,
    window_length: int = 11,
    polyorder: int = 3,
) -> pd.DataFrame:
    """
    Correct GPS drift via kinematic masking + Savitzky-Golay smoothing.

    Parameters
    ----------
    df : pd.DataFrame
        Must contain time_col, lat_col, lon_col columns.
    max_speed_kmh : float
        Hard speed ceiling for the vehicle type (km/h).
    max_accel_ms2 : float
        Hard lateral acceleration ceiling (m/s²).
    window_length : int
        Savitzky-Golay window (must be odd; auto-corrected if even).
    polyorder : int
        Polynomial order for the SG filter (2–3 is typical).

    Returns
    -------
    pd.DataFrame
        Input frame with added '{lat_col}_corrected' and
        '{lon_col}_corrected' columns. Original columns are untouched.

    Raises
    ------
    ValueError
        If required columns are missing, the frame is empty, or there
        are fewer rows than window_length after resampling.
    """
    required = {time_col, lat_col, lon_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if df.empty:
        raise ValueError("Input DataFrame is empty.")

    df = df.copy().sort_values(time_col).reset_index(drop=True)

    # ── Stage 1: Temporal regularization ──────────────────────────────
    # Resample to a uniform 1 s cadence. OS scheduling and power-saving
    # modes cause irregular arrival intervals in raw device logs.
    df[time_col] = pd.to_datetime(df[time_col])
    df = (
        df.set_index(time_col)
        .resample("1s")
        .interpolate("linear")
        .reset_index()
    )
    df = df.dropna(subset=[lat_col, lon_col]).reset_index(drop=True)

    if len(df) < window_length:
        raise ValueError(
            f"Only {len(df)} points after resampling; "
            f"need at least window_length={window_length}."
        )

    # ── Stage 2a: Project to local metric CRS (UTM) ───────────────────
    # NEVER compute distances in EPSG:4326. A degree of longitude near
    # the equator is ~111 km but shrinks toward the poles. Derive the
    # UTM zone from the data centroid to stay within the zone's
    # low-distortion corridor (<0.04% scale error).
    centroid_lat = df[lat_col].mean()
    centroid_lon = df[lon_col].mean()
    utm_zone = int((centroid_lon + 180) // 6) + 1
    utm_epsg = (32600 if centroid_lat >= 0 else 32700) + utm_zone

    to_local = Transformer.from_crs(
        "EPSG:4326", f"EPSG:{utm_epsg}", always_xy=True
    )
    to_wgs = Transformer.from_crs(
        f"EPSG:{utm_epsg}", "EPSG:4326", always_xy=True
    )

    # always_xy=True: lon before lat — avoids silent axis-order swaps.
    x, y = to_local.transform(df[lon_col].values, df[lat_col].values)

    # ── Stage 2b: Kinematic thresholding ──────────────────────────────
    # np.gradient uses central differences; first/last points get
    # forward/backward differences, so no NaN sentinel is introduced.
    dt = 1.0  # seconds (fixed cadence from Stage 1)
    vx = np.gradient(x, dt)
    vy = np.gradient(y, dt)
    speed_ms = np.sqrt(vx**2 + vy**2)

    ax = np.gradient(vx, dt)
    ay = np.gradient(vy, dt)
    accel_ms2 = np.sqrt(ax**2 + ay**2)

    speed_mask = speed_ms > (max_speed_kmh / 3.6)
    accel_mask = accel_ms2 > max_accel_ms2
    outlier_mask = speed_mask | accel_mask

    # Replace outlier positions with NaN, then interpolate linearly.
    # Do NOT drop rows — the uniform time axis must be preserved for SG.
    x_clean = np.where(outlier_mask, np.nan, x).astype(float)
    y_clean = np.where(outlier_mask, np.nan, y).astype(float)

    x_interp = pd.Series(x_clean).interpolate(method="linear").values
    y_interp = pd.Series(y_clean).interpolate(method="linear").values

    # ── Stage 3: Spatial smoothing (Savitzky-Golay) ───────────────────
    # SG fits a local polynomial; it suppresses jitter while preserving
    # peak geometry better than a simple moving average.
    # window_length must be odd — auto-correct if caller passed even.
    if window_length % 2 == 0:
        window_length += 1

    x_smooth = savgol_filter(x_interp, window_length, polyorder, mode="interp")
    y_smooth = savgol_filter(y_interp, window_length, polyorder, mode="interp")

    # ── Stage 5: Reproject to WGS84 ───────────────────────────────────
    lon_out, lat_out = to_wgs.transform(x_smooth, y_smooth)

    df[f"{lat_col}_corrected"] = lat_out
    df[f"{lon_col}_corrected"] = lon_out
    return df

Validation block

Run these checks immediately after calling correct_gps_drift. Each check targets a distinct failure mode — over-smoothing, projection errors, and phantom stop injection — so a passing set gives reasonable confidence before the output enters downstream analysis.

PYTHON 
import math

def validate_correction(
    raw: pd.DataFrame,
    corrected: pd.DataFrame,
    lat_col: str = "lat",
    lon_col: str = "lon",
) -> None:
    """
    Sanity-check a corrected trajectory DataFrame.
    Raises AssertionError on failure; prints a summary on success.
    """
    # 1. Shape contract: correction must not drop rows
    assert len(corrected) >= len(raw), (
        f"Row count shrank: {len(raw)} → {len(corrected)}"
    )

    # 2. Corrected columns exist and contain no NaN
    for col in (f"{lat_col}_corrected", f"{lon_col}_corrected"):
        assert col in corrected.columns, f"Missing output column: {col}"
        assert corrected[col].notna().all(), f"NaN values in {col}"

    # 3. Coordinate bounds (basic sanity — catches axis-order swaps)
    assert corrected[f"{lat_col}_corrected"].between(-90, 90).all(), (
        "Corrected latitudes out of range — possible axis-order swap"
    )
    assert corrected[f"{lon_col}_corrected"].between(-180, 180).all(), (
        "Corrected longitudes out of range"
    )

    # 4. Distance delta: correction should reduce trip length.
    #    >15% reduction → successful jitter removal.
    #    >40% reduction → possible over-smoothing; review window_length.
    def haversine_km(lat1, lon1, lat2, lon2):
        R = 6371.0
        phi1, phi2 = math.radians(lat1), math.radians(lat2)
        dphi = math.radians(lat2 - lat1)
        dlam = math.radians(lon2 - lon1)
        a = math.sin(dphi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(dlam / 2) ** 2
        return R * 2 * math.asin(math.sqrt(a))

    raw_dist = sum(
        haversine_km(
            raw.iloc[i][lat_col], raw.iloc[i][lon_col],
            raw.iloc[i + 1][lat_col], raw.iloc[i + 1][lon_col],
        )
        for i in range(len(raw) - 1)
    )
    corr_dist = sum(
        haversine_km(
            corrected.iloc[i][f"{lat_col}_corrected"],
            corrected.iloc[i][f"{lon_col}_corrected"],
            corrected.iloc[i + 1][f"{lat_col}_corrected"],
            corrected.iloc[i + 1][f"{lon_col}_corrected"],
        )
        for i in range(len(corrected) - 1)
    )

    pct_reduction = (raw_dist - corr_dist) / raw_dist * 100
    assert pct_reduction < 40, (
        f"Over-smoothing likely: trip length reduced by {pct_reduction:.1f}%"
    )
    print(
        f"Validation passed. Distance: {raw_dist:.2f} km → "
        f"{corr_dist:.2f} km ({pct_reduction:.1f}% reduction, "
        f"{len(corrected)} points)."
    )

Key assertions:

  • Shape contract: the correction must not silently drop rows; the uniform time axis is required by downstream sampling rate optimization steps.
  • Coordinate bounds: out-of-range values almost always indicate a lat/lon axis-order swap introduced during the pyproj round-trip.
  • Distance delta: a 15–30% reduction indicates successful jitter removal; above 40% suggests over-smoothing — tighten window_length or lower polyorder.

Common mistakes and gotchas

  • Computing gradients in WGS84 degree-space. Dividing a degree-of-latitude difference by seconds gives a value in degrees/second, not m/s. The speed threshold comparison then fires on completely normal movement near the equator and fails silently near the poles. Always project to a metric CRS before computing any derivative.

  • Using iterrows() to compute per-point distances. Pure Python loops over a million-row trajectory log take minutes. np.gradient on projected coordinate arrays runs in milliseconds. Use vectorized NumPy operations throughout.

  • Ignoring NaN sentinels after masking. After replacing outlier coordinates with np.nan, linear interpolation may still leave leading or trailing NaN values if the first or last points were masked. Call .ffill().bfill() after .interpolate() to guarantee a clean array before passing it to savgol_filter.

  • Choosing an even window_length. scipy.signal.savgol_filter raises a ValueError for even windows. The function above auto-corrects this, but if you call savgol_filter directly, ensure the window is odd before the call — silent wrapping in a wrapper function can mask this bug.

  • Smoothing across trajectory segment boundaries. A long gap in the input log (e.g. a 5-minute tunnel blackout) followed by a new GPS fix looks like a sudden teleport. Resampling fills this gap with linearly interpolated phantom positions, and the smoother then propagates them into the valid data on both sides. Split the trajectory at gaps exceeding a threshold (typically 30 s) and correct each segment independently. This integrates naturally with downsampling high-frequency GPS tracks workflows that also operate on segmented tracks.

  • Using the same thresholds for all vehicle types. The default 120 km/h / 4.5 m/s² parameters suit passenger vehicles. For motorcycles raise the acceleration ceiling to 7.0 m/s²; for heavy goods vehicles lower the speed ceiling to 90 km/h; for pedestrian or cycling data drop both significantly (15 km/h / 1.5 m/s²). Mismatched thresholds mask legitimate motion patterns rather than noise.

Back to GPS Precision & Error Handling