Interpolating Missing GPS Points with Kalman Filters

A Kalman filter fills GPS dropouts by replacing naive straight-line assumptions with a recursive Bayesian estimator that models vehicle kinematics and sensor uncertainty. During signal loss only the prediction step runs, advancing the trajectory while inflating positional covariance; when a valid ping returns, the update step fuses prediction and measurement weighted by their relative confidence. The result is a statistically optimal, temporally continuous trace with an explicit uncertainty column that downstream analytics can consume directly.

Why This Happens

Raw telemetry streams in urban mobility pipelines rarely arrive at fixed intervals. Multipath reflections in dense street canyons, tunnel outages, and aggressive device power-saving modes create irregular gaps that break standard resampling routines. These dropouts are the core problem addressed by Gap Filling in Sparse Trajectories, which sits at the foundation of Temporal Aggregation & Window Mapping.

Geometric methods — linear and cubic splines — know nothing about momentum. They connect the last known point to the next observed point as if the vehicle teleported along a smooth curve, producing unrealistic artefacts across even 10-second gaps. A state-space model respects physics: the filter’s motion model encodes that position changes with velocity and that velocity changes slowly, so predictions drift realistically under uncertainty rather than drawing implausible straight lines.

The diagram below shows the predict-update cycle across a 3-ping dropout.

Kalman filter predict-update cycle across a GPS dropout Timeline showing five ticks: ticks 1, 5, and 6 receive GPS measurements and run both predict and update steps; ticks 2, 3, and 4 are GPS dropouts and run only the predict step, with growing uncertainty ellipses. t₁ t₂ t₃ t₄ t₅ t₆ Predict+Update Predict only Predict only Predict only Predict+Update Predict+Update GPS fix received Predicted (no fix) Uncertainty envelope — GPS dropout (t₂ – t₄) —

Core Mitigation Pipeline

  1. Project to a metric CRS — reproject WGS84 coordinates to EPSG:3857 or a local UTM zone so that state velocities carry metre-per-second units rather than meaningless decimal-degree deltas.
  2. Build the constant-velocity state-space model — define the state vector, transition matrix parameterised by elapsed Δt, observation matrix, and noise covariances Q and R.
  3. Run the predict-update loop — for every timestamp: always execute the prediction step; execute the update step only when a non-null measurement is present.
  4. Attach uncertainty metadata and reproject — store the combined positional standard deviation per row, then reproject output back to WGS84 for storage or downstream joins.

Production-Ready Python Implementation

The implementation below uses filterpy and projects through pyproj so all state units are metres and m/s. It handles irregular timestamps, skips missing measurements, returns a dense trajectory aligned to the original index, and raises informative errors on degenerate inputs.

PYTHON 
import numpy as np
import pandas as pd
from filterpy.kalman import KalmanFilter
from pyproj import Transformer

# NOTE: All distance/velocity arithmetic is performed in EPSG:3857 (metres).
# Raw WGS84 degrees are never fed into the filter — doing so makes Q and R
# dimensionally meaningless and silently degrades accuracy.

def interpolate_gps_kalman(
    df: pd.DataFrame,
    lat_col: str = "lat",
    lon_col: str = "lon",
    time_col: str = "timestamp",
    process_noise: float = 0.5,    # Q scale: m²/s⁴ — tune per transport mode
    measurement_noise: float = 25.0,  # R: m² — ~5 m CEP50 for consumer GNSS
) -> pd.DataFrame:
    """
    Fill missing GPS fixes in *df* using a 2D constant-velocity Kalman filter.

    Parameters
    ----------
    df : pd.DataFrame
        Must have columns *lat_col*, *lon_col*, *time_col*.
        Rows where lat/lon are NaN are treated as dropout ticks.
        The DataFrame must be sorted by *time_col* and contain at least
        two non-null fixes for initialisation.
    process_noise : float
        Scale factor for the process noise covariance Q (m²/s⁴).
        Higher values allow the filter to track sharper manoeuvres.
    measurement_noise : float
        GPS measurement variance in m² (R diagonal).
        Derive from receiver CEP50: R ≈ (CEP50)².

    Returns
    -------
    pd.DataFrame
        Copy of *df* with 'lat', 'lon' filled and two new columns:
        'is_interpolated' (bool) and 'kf_uncertainty_m' (float, metres std-dev).
    """
    if df.empty:
        raise ValueError("Input DataFrame is empty.")

    required = {lat_col, lon_col, time_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")

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

    # Identify the first non-null fix for initialisation
    valid_mask = df[lat_col].notna() & df[lon_col].notna()
    if valid_mask.sum() < 2:
        raise ValueError("Need at least 2 non-null GPS fixes to initialise the filter.")

    first_valid = df.loc[valid_mask].iloc[0]

    # --- Project to metric CRS (EPSG:3857) ---
    to_metric = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
    to_wgs84  = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)

    init_x, init_y = to_metric.transform(first_valid[lon_col], first_valid[lat_col])

    # --- Initialise 2D constant-velocity filter ---
    # State: [x, y, vx, vy]  (metres, m/s in EPSG:3857)
    kf = KalmanFilter(dim_x=4, dim_z=2)
    kf.x = np.array([init_x, init_y, 0.0, 0.0], dtype=float)
    kf.P = np.eye(4) * 50.0          # Initial state uncertainty (m² / (m/s)²)
    kf.H = np.array([[1, 0, 0, 0],   # Observe x position
                     [0, 1, 0, 0]])  # Observe y position
    kf.R = np.eye(2) * measurement_noise
    kf.Q = np.eye(4) * process_noise
    # F is set per-tick to handle variable Δt

    out_lats = df[lat_col].values.copy().astype(float)
    out_lons = df[lon_col].values.copy().astype(float)
    out_unc  = np.full(len(df), np.nan)
    is_interp = np.zeros(len(df), dtype=bool)

    timestamps = pd.to_datetime(df[time_col])

    for i in range(len(df)):
        if i == 0:
            # No prediction at the very first tick; just record init state
            out_x, out_y = kf.x[0], kf.x[1]
            out_unc[i] = float(np.sqrt(kf.P[0, 0] + kf.P[1, 1]))
            # If first row has no fix, mark interpolated
            if not valid_mask.iloc[i]:
                is_interp[i] = True
            lon_out, lat_out = to_wgs84.transform(out_x, out_y)
            out_lats[i], out_lons[i] = lat_out, lon_out
            continue

        dt = (timestamps.iloc[i] - timestamps.iloc[i - 1]).total_seconds()
        dt = max(dt, 0.1)  # Guard against duplicate timestamps

        # Time-varying transition matrix
        kf.F = np.array([
            [1, 0, dt,  0],
            [0, 1,  0, dt],
            [0, 0,  1,  0],
            [0, 0,  0,  1],
        ], dtype=float)

        # Prediction step — always runs
        kf.predict()

        # Update step — only when a valid fix exists
        if valid_mask.iloc[i]:
            meas_x, meas_y = to_metric.transform(
                df[lon_col].iloc[i], df[lat_col].iloc[i]
            )
            kf.update(np.array([meas_x, meas_y]))
        else:
            is_interp[i] = True

        out_unc[i] = float(np.sqrt(kf.P[0, 0] + kf.P[1, 1]))
        lon_out, lat_out = to_wgs84.transform(kf.x[0], kf.x[1])
        out_lats[i], out_lons[i] = lat_out, lon_out

    df[lat_col]           = out_lats
    df[lon_col]           = out_lons
    df["is_interpolated"] = is_interp
    df["kf_uncertainty_m"] = out_unc

    return df

Validation Block

After running interpolate_gps_kalman, confirm the output before passing it downstream:

PYTHON 
import numpy as np

def validate_kalman_output(result: pd.DataFrame) -> None:
    # No NaN coordinates should remain
    assert result["lat"].notna().all(), "NaN latitudes remain after interpolation"
    assert result["lon"].notna().all(), "NaN longitudes remain after interpolation"

    # Coordinates must stay in plausible WGS84 range
    assert result["lat"].between(-90, 90).all(),  "Latitude out of range"
    assert result["lon"].between(-180, 180).all(), "Longitude out of range"

    # Uncertainty column must be non-negative and finite
    assert (result["kf_uncertainty_m"] >= 0).all(), "Negative uncertainty values"
    assert np.isfinite(result["kf_uncertainty_m"]).all(), "Non-finite uncertainty"

    # Log interpolation rate for monitoring
    interp_rate = result["is_interpolated"].mean()
    if interp_rate > 0.30:
        import warnings
        warnings.warn(
            f"High interpolation rate: {interp_rate:.1%} of fixes are synthetic. "
            "Verify gap classification before aggregating.",
            stacklevel=2,
        )

    print(
        f"Validation passed — {len(result)} rows, "
        f"{result['is_interpolated'].sum()} interpolated "
        f"({interp_rate:.1%}), "
        f"max uncertainty {result['kf_uncertainty_m'].max():.1f} m"
    )
  • Expected shape: same row count as input, two new columns added.
  • is_interpolated fraction above 30 % signals that gap classification upstream needs review before these rows feed into rolling statistics for mobility metrics.
  • kf_uncertainty_m should collapse after each real fix; a monotonically growing value with no collapses means no valid measurements were processed.
  • Spot-check the projection round-trip: project a known fix to EPSG:3857, run it through the filter for one tick with no dropout, and confirm the output lat/lon matches the input to within 0.5 m.

Common Mistakes & Gotchas

  • Running the filter on raw WGS84 degrees. State velocity units become degrees/second, making Q and R dimensionally incoherent. Always project to a metric CRS before initialising the state vector — this is flagged explicitly in every code comment above.
  • Forgetting to update F per tick. A fixed transition matrix hardcoded with a constant dt silently corrupts trajectories with irregular sampling. The loop above recomputes F at every step from the actual elapsed seconds.
  • Treating interpolated rows the same as observed rows in aggregations. Downstream dynamic time-binning strategies and ETA models should weight or exclude rows where is_interpolated=True. The uncertainty column makes this weighting continuous rather than binary.
  • Using a single global Q across all device types. A pedestrian’s Q (slow acceleration budget, ~0.05 m²/s⁴) destroys accuracy when applied to a courier van (0.5–2.0 m²/s⁴). Segment by transport mode and tune separately; see speed and acceleration profiling for mode-specific velocity priors.
  • Initialising with a NaN fix. If the very first row is a dropout the initial state is undefined. The implementation above finds the first valid fix and initialises from there; anything before that fix is dead-reckoned backward only if you explicitly reverse the time axis — the simpler default is to mark those early rows as interpolated and flag them.
  • Neglecting numerical stability on long dropouts. After many predict-only steps P grows large. If H P Hᵀ + R approaches singularity (rare with diagonal R but possible with poorly conditioned Q), the matrix inversion inside filterpy raises a LinAlgError. Add a maximum-gap cap: classify gaps longer than, say, 5 minutes as unrecoverable in the upstream gap-filling pipeline and reset the filter state after each break.

Back to Gap Filling in Sparse Trajectories