Calculating Instantaneous Speed from Discrete GPS Points

To compute instantaneous speed from discrete GPS fixes, calculate the geodesic distance between consecutive coordinate pairs, divide by the elapsed seconds between their UTC timestamps, and apply noise masking to suppress multipath spikes. Because a GNSS receiver samples a moving object at discrete moments rather than continuously, the result is a finite-difference approximation — the mean speed across the sampling interval — not a true limit-based derivative. In practice this distinction matters only when analysing rapid events such as emergency braking; for logistics, fleet scoring, and urban mobility analytics the approximation is sufficient once the pipeline handles irregular sampling, timezone hazards, and GPS drift correctly.

Why this calculation is error-prone

Speed is a first derivative of position. Any error in the input positions is amplified when you divide a small noisy distance by a small time interval, so the quality of the output is tightly coupled to the upstream GPS precision and error handling stage. Consumer-grade GNSS receivers in open sky achieve 2–5 m horizontal accuracy at 1 Hz. A vehicle travelling at 30 m/s covers 30 m per sample; positional jitter of ±3 m therefore injects roughly ±10% variance into the raw speed estimate before any smoothing is applied.

Three root causes produce the worst artefacts:

  • Multipath reflection in urban canyons shifts a point 5–30 m sideways for one epoch, then snaps back. This generates a phantom high-speed excursion that looks like a genuine acceleration spike in the derivative.
  • Clock jitter and timezone errors corrupt the denominator. A timestamp in local time accidentally treated as UTC produces a time delta that is off by hours, giving absurd speed values without raising a Python exception.
  • Zero or negative time deltas from device buffering or log merging cause divide-by-zero or negative speed values that propagate silently through pandas arithmetic as inf or NaN.

These failure modes sit within the broader speed and acceleration profiling workflow, which also covers Savitzky-Golay smoothing, entity-level groupby patterns, and CRS projection for acceleration calculations.

Core pipeline steps

  1. Parse and normalise timestamps to UTC. Coerce all timestamp columns to timezone-aware datetime64[ns, UTC] using pd.to_datetime(..., utc=True) so that diff() produces monotonic, daylight-saving-safe deltas.
  2. Sort by entity and time. Multi-entity logs must be sorted within each entity_id group before shifting; out-of-order rows produce negative deltas and phantom reverse movement.
  3. Compute haversine distance between consecutive pairs. Vectorise across the entire frame using shift(1) alignment — no Python-level loop. The haversine formula gives metre-accurate distances for inter-sample gaps under 50 km without requiring a projected CRS.
  4. Divide distance by elapsed seconds, masking invalid intervals. Guard against dt <= 0, NaN on the first row of each entity, and gaps that exceed the expected sampling cadence.
Instantaneous Speed Pipeline Four sequential stages — UTC Normalisation, Haversine Distance, Speed Division, and Noise Filtering — connected left to right by arrows. 1. UTC Normalisation 2. Haversine Distance (m) 3. Speed Division (m/s) 4. Noise Filter & Validate

Production-ready Python implementation

The function below is fully vectorised — no iterrows, no Python-level loops. It handles unsorted inputs, multi-entity frames, missing coordinate values, duplicate timestamps, and trajectory gaps in one pass.

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


def haversine_vectorized(
    lat1: np.ndarray,
    lon1: np.ndarray,
    lat2: np.ndarray,
    lon2: np.ndarray,
) -> np.ndarray:
    """
    Return great-circle distance in metres between two arrays of WGS84 points.

    Uses the spherical haversine formula (error < 0.3% for distances under 50 km).
    Inputs are NumPy arrays of equal length; NaN propagates through cleanly.
    """
    R = 6_371_000.0  # IUGG mean radius in metres
    phi1 = np.radians(lat1)
    phi2 = np.radians(lat2)
    dphi = np.radians(lat2 - lat1)
    dlambda = np.radians(lon2 - lon1)

    a = (
        np.sin(dphi / 2.0) ** 2
        + np.cos(phi1) * np.cos(phi2) * np.sin(dlambda / 2.0) ** 2
    )
    return 2.0 * R * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a))


def calculate_speed(
    df: pd.DataFrame,
    lat_col: str = "lat",
    lon_col: str = "lon",
    time_col: str = "timestamp",
    entity_col: str | None = None,
    max_speed_ms: float = 83.3,       # ~300 km/h hard ceiling
    min_dist_m: float = 0.5,          # below this, treat as stationary
    sg_window: int = 11,              # Savitzky-Golay window (must be odd)
    sg_poly: int = 3,
) -> pd.DataFrame:
    """
    Compute smoothed speed (m/s) from discrete GPS fixes.

    Parameters
    ----------
    df : pd.DataFrame
        Must contain lat_col, lon_col, time_col. If entity_col is provided,
        sorting and derivative calculations are scoped to each entity.
    max_speed_ms : float
        Hard ceiling in m/s. Rows exceeding this are treated as GPS artefacts
        and replaced with NaN before smoothing.
    min_dist_m : float
        Minimum inter-sample distance below which speed is set to 0.0
        (suppresses stationary GPS jitter).
    sg_window : int
        Savitzky-Golay window length. Must be odd and >= sg_poly + 2.
        Auto-incremented by 1 if an even value is passed.
    sg_poly : int
        SG polynomial order. 3 is standard for vehicle kinematics.

    Returns
    -------
    pd.DataFrame
        Original frame with added columns: 'dist_m', 'dt_s',
        'speed_raw_mps', 'speed_mps' (smoothed).

    Raises
    ------
    ValueError
        If required columns are missing or the frame is empty.
    """
    required = {lat_col, lon_col, time_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()

    # ── Step 1: UTC timestamp normalisation ───────────────────────────
    # utc=True coerces naive timestamps to UTC and converts tz-aware ones.
    # This eliminates daylight-saving-time discontinuities in diff().
    df[time_col] = pd.to_datetime(df[time_col], utc=True)

    # Sort within each entity (or globally if no entity column).
    sort_keys = ([entity_col, time_col] if entity_col else [time_col])
    df = df.sort_values(sort_keys).reset_index(drop=True)

    # ── Step 2: Align consecutive pairs with shift(1) ─────────────────
    # Group-aware shift so entity boundaries never bleed into each other.
    if entity_col:
        grp = df.groupby(entity_col, sort=False)
        lat_prev = grp[lat_col].shift(1)
        lon_prev = grp[lon_col].shift(1)
        dt_s = grp[time_col].diff().dt.total_seconds()
    else:
        lat_prev = df[lat_col].shift(1)
        lon_prev = df[lon_col].shift(1)
        dt_s = df[time_col].diff().dt.total_seconds()

    # ── Step 3: Vectorised haversine distance ─────────────────────────
    # IMPORTANT: pass positional arrays, not Series, so NaN propagates
    # correctly through arctan2 without raising a casting warning.
    dist_m = haversine_vectorized(
        lat_prev.to_numpy(),
        lon_prev.to_numpy(),
        df[lat_col].to_numpy(),
        df[lon_col].to_numpy(),
    )

    # Sub-threshold movement → treat as stationary (clamp to 0.0)
    dist_m = np.where(dist_m < min_dist_m, 0.0, dist_m)

    # ── Step 4: Speed division with validity mask ─────────────────────
    # Mask: dt > 0, not NaN, and not the first row of each entity.
    valid = (dt_s > 0) & dt_s.notna()
    speed_raw = np.full(len(df), np.nan)
    speed_raw[valid] = dist_m[valid] / dt_s[valid]

    # Hard-ceiling outlier masking (GPS multipath / buffered flush artefacts)
    speed_raw = np.where(speed_raw > max_speed_ms, np.nan, speed_raw)

    df["dist_m"] = dist_m
    df["dt_s"] = dt_s.to_numpy()
    df["speed_raw_mps"] = speed_raw

    # ── Step 5: Savitzky-Golay smoothing per entity ───────────────────
    # SG preserves genuine acceleration peaks better than a simple EMA.
    # window_length must be odd; auto-correct even caller values.
    wl = sg_window if sg_window % 2 == 1 else sg_window + 1

    def _sg_smooth(series: pd.Series) -> np.ndarray:
        vals = series.fillna(0.0).to_numpy(dtype=float)
        effective_wl = min(wl, len(vals) if len(vals) % 2 == 1 else len(vals) - 1)
        if effective_wl < sg_poly + 2 or effective_wl < 3:
            return vals  # segment too short to filter
        return savgol_filter(vals, window_length=effective_wl, polyorder=sg_poly,
                             mode="nearest")

    if entity_col:
        df["speed_mps"] = (
            df.groupby(entity_col, sort=False)["speed_raw_mps"]
            .transform(_sg_smooth)
        )
    else:
        df["speed_mps"] = _sg_smooth(df["speed_raw_mps"])

    return df

Validation block

Run these checks immediately after calling calculate_speed. Each assertion targets a distinct failure mode: unit confusion, projection errors, and smoothing over-reach.

PYTHON 
def validate_speed_output(df: pd.DataFrame, max_speed_ms: float = 83.3) -> None:
    """
    Sanity-check the output of calculate_speed().
    Raises AssertionError with a descriptive message on failure.
    """
    assert "speed_mps" in df.columns, "Missing 'speed_mps' column — did calculate_speed() run?"

    non_nan = df["speed_mps"].dropna()
    assert len(non_nan) > 0, "All speed values are NaN — check timestamp parsing and entity sorting."

    assert (non_nan >= 0).all(), (
        f"Negative speed values found: min={non_nan.min():.3f} m/s. "
        "Likely cause: unsorted timestamps or swapped lat/lon columns."
    )
    assert (non_nan <= max_speed_ms).all(), (
        f"Speed exceeds ceiling {max_speed_ms} m/s (max found: {non_nan.max():.2f}). "
        "Increase max_speed_ms or investigate GPS multipath artefacts."
    )

    nan_frac = df["speed_mps"].isna().mean()
    assert nan_frac < 0.10, (
        f"{nan_frac:.1%} of speed values are NaN — unusually high. "
        "Check for duplicate timestamps, excessive gaps, or missing entities."
    )

    # Coordinate bounds — catches axis-order swaps introduced upstream
    assert df["lat"].between(-90, 90).all(), (
        "Latitude values out of range — possible lon/lat column swap."
    )
    assert df["lon"].between(-180, 180).all(), (
        "Longitude values out of range."
    )

    print(
        f"Speed validation passed. "
        f"Points: {len(df)}, NaN rate: {nan_frac:.1%}, "
        f"Speed range: {non_nan.min():.2f}–{non_nan.max():.2f} m/s."
    )

Key assertions:

  • Non-negative speeds. Negative values almost always indicate unsorted timestamps or a swapped lat/lon column that produces the coordinate-pair ordering expected by shift(1).
  • Ceiling check. Any speed above the max_speed_ms threshold that survived the masking step indicates uncorrected GPS drift in the input.
  • NaN fraction. Isolated NaN on the first row of each entity is expected and normal. A fraction above 10% points to systematic timestamp problems, duplicate coordinate pairs, or gaps in the telemetry stream that should be handled at the time-series synchronization stage.

Common mistakes and gotchas

  • Computing haversine in degree-space and calling it “metric”. The haversine formula itself returns metres correctly — the mistake is applying Euclidean distance directly to latitude and longitude columns (e.g. np.sqrt((lat2-lat1)**2 + (lon2-lon1)**2)) and treating the result as metres. One degree of longitude near the equator is ~111 km; near 60° N it is ~55 km. The error is systematic and transport-mode-dependent.

  • Using iterrows() or apply(lambda row: ...) for the distance loop. Iterating row-by-row over a million-point GPS log is 100–1000× slower than the vectorised shift(1) + NumPy approach. Even a 10 Hz fleet log with 100 vehicles generating 24 h of data reaches ~86 million rows; iterrows makes this infeasible.

  • Ignoring entity boundaries when shifting. If df contains multiple vehicles and you call df["lat"].shift(1) without a groupby, the last fix of vehicle A becomes the “previous point” for the first fix of vehicle B. The resulting distance is the straight-line gap between two unrelated locations, producing a phantom multi-kilometre excursion in the derivative.

  • Treating naive timestamps as UTC. pandas silently accepts naive timestamps in diff(). If the device logs in local time (e.g. CET = UTC+1) and you forget the utc=True flag, the deltas are still numerically correct — but a log spanning a daylight-saving boundary will contain a one-hour gap or a negative delta that invalidates the entire night segment.

  • Choosing window_length without checking series length. scipy.signal.savgol_filter raises a ValueError if window_length exceeds the series length. Short trajectory segments — such as a vehicle with only three GPS fixes in a log — silently fail if the guard on effective_wl is omitted. The function above guards this explicitly; replicate the guard if you call savgol_filter directly.

  • Skipping the stationary-drift floor. A stopped vehicle still accumulates 0.5–2 m of positional jitter per second from satellite geometry fluctuations. Without the min_dist_m floor, the pipeline reports 0.5–2 m/s for a parked truck. Downstream stay-point detection algorithms and driver behaviour scoring both depend on accurate zero-speed identification.

Back to Speed & Acceleration Profiling