Vehicle Speed and Bearing from GTFS-RT

When a GTFS-Realtime feed omits position.speed and position.bearing, derive them from two consecutive fixes for the same vehicle.id: divide the haversine distance between the coordinates by the timestamp delta to get ground speed in metres per second, and feed the same coordinate pair through the initial-bearing formula to get the compass heading in degrees clockwise from north. Keep a per-vehicle record of the previous fix, guard against tiny time deltas that amplify GPS jitter, and cap results at a mode-plausible ceiling.

Root Cause Analysis

The speed and bearing fields are optional in the VehiclePosition message, and a large share of production feeds simply never set them. The cause is upstream of the feed: many AVL units report only a timestamped coordinate, leaving the derived kinematics to the consumer. Because protobuf has no null on the wire, an unset position.speed decodes to 0.0 and an unset position.bearing decodes to 0.0 — a real value meaning “north”, indistinguishable from “unknown” unless the decode layer used HasField guards, as covered in the parent tracking vehicle positions in realtime guide.

Recovering the missing kinematics is pure geometry over time. Two fixes for the same vehicle define a displacement vector; its magnitude divided by the elapsed time is speed, and its orientation is bearing. The only inputs are the two coordinate pairs and the two timestamp values, all of which the feed provides even when it withholds the derived fields. The subtlety is entirely in the guarding: a two-second gap over a three-metre GPS wander implies 5 km/h of motion for a parked bus, so the raw calculation must be wrapped in sanity limits before it is trusted.

Distance at this scale is well served by the haversine great-circle formula, which needs no projection and stays sub-metre-accurate over the short hops between polls. When you are already projecting fixes to a metric coordinate reference system for other spatial work — for example the shape snapping in spatial analysis and route geometry — a planar distance gives an equivalent result and you can reuse the projected coordinates instead.

Production-Ready Python Implementation

The script maintains a per-vehicle_id cache of the previous fix across polls, so it can be called once per snapshot in a live loop. For each vehicle it computes the timestamp delta, the haversine distance, the resulting speed, and the initial bearing, applies jitter and ceiling guards, and returns the snapshot frame with speed and bearing populated only where a valid prior fix existed.

python
import logging
import math
from dataclasses import dataclass

import numpy as np
import pandas as pd

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

EARTH_RADIUS_M = 6_371_000.0
MIN_DT_SECONDS = 3           # reject deltas below this — jitter dominates
MAX_SPEED_MPS = 35.0         # ~126 km/h ceiling; clamp implausible bus/tram speeds
MIN_MOVE_METRES = 5.0        # below this, treat the vehicle as stationary


@dataclass
class Fix:
    lat: float
    lon: float
    timestamp: int


def haversine_metres(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    """Great-circle distance between two WGS84 points, in metres."""
    p1, p2 = math.radians(lat1), math.radians(lat2)
    d_phi = math.radians(lat2 - lat1)
    d_lambda = math.radians(lon2 - lon1)
    a = (
        math.sin(d_phi / 2) ** 2
        + math.cos(p1) * math.cos(p2) * math.sin(d_lambda / 2) ** 2
    )
    return EARTH_RADIUS_M * 2 * math.asin(math.sqrt(a))


def initial_bearing(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    """Initial bearing from point 1 to point 2, degrees clockwise from true north."""
    p1, p2 = math.radians(lat1), math.radians(lat2)
    d_lambda = math.radians(lon2 - lon1)
    x = math.sin(d_lambda) * math.cos(p2)
    y = math.cos(p1) * math.sin(p2) - math.sin(p1) * math.cos(p2) * math.cos(d_lambda)
    return (math.degrees(math.atan2(x, y)) + 360.0) % 360.0


def derive_speed_bearing(
    snapshot: pd.DataFrame, previous: dict[str, Fix] | None = None
) -> tuple[pd.DataFrame, dict[str, Fix]]:
    """
    Fill speed (m/s) and bearing (degrees) from the previous fix per vehicle_id.

    Parameters
    ----------
    snapshot : pd.DataFrame
        Current poll with vehicle_id, latitude, longitude, timestamp columns.
    previous : dict[str, Fix] | None
        Cache of each vehicle's last valid fix, carried across polls.

    Returns
    -------
    (pd.DataFrame, dict[str, Fix])
        The snapshot with derived speed/bearing, and the updated cache.
    """
    previous = dict(previous or {})
    result = snapshot.copy()
    result["speed_derived"] = pd.Series(pd.NA, index=result.index, dtype="float64")
    result["bearing_derived"] = pd.Series(pd.NA, index=result.index, dtype="float64")

    for idx, row in result.iterrows():
        vehicle_id = row["vehicle_id"]
        if pd.isna(row["latitude"]) or pd.isna(row["timestamp"]):
            continue
        current = Fix(float(row["latitude"]), float(row["longitude"]), int(row["timestamp"]))
        prior = previous.get(vehicle_id)

        if prior is not None:
            dt = current.timestamp - prior.timestamp
            if dt >= MIN_DT_SECONDS:
                dist = haversine_metres(prior.lat, prior.lon, current.lat, current.lon)
                if dist >= MIN_MOVE_METRES:
                    speed = min(dist / dt, MAX_SPEED_MPS)
                    result.loc[idx, "speed_derived"] = speed
                    result.loc[idx, "bearing_derived"] = initial_bearing(
                        prior.lat, prior.lon, current.lat, current.lon
                    )
                else:
                    result.loc[idx, "speed_derived"] = 0.0  # confirmed stationary

        previous[vehicle_id] = current  # advance the cache regardless

    filled = int(result["speed_derived"].notna().sum())
    logging.info("Derived speed/bearing for %d/%d vehicles", filled, len(result))
    return result, previous

Step-by-Step Walkthrough

Per-vehicle state across polls. The previous dict keys the last valid fix by vehicle_id. Passing it back into the next call lets a live polling loop compute deltas across snapshots without holding full history in memory. On the very first poll every vehicle is new, so no speed is produced — expected behaviour, not an error.

Delta and jitter guard. dt is the difference of the two per-entity timestamp values, not the feed header. Deltas below MIN_DT_SECONDS are rejected because a short interval turns a metre of GPS noise into a large phantom speed. Likewise, moves below MIN_MOVE_METRES are treated as a stationary vehicle and reported as 0.0 rather than a spurious crawl.

Haversine distance. haversine_metres returns the great-circle distance directly from WGS84 degrees, so no projection step is needed for the speed calculation. Over inter-poll distances the great-circle and planar answers agree to well within GPS error.

Initial-bearing formula. initial_bearing computes the forward azimuth from the prior fix to the current one using atan2, then normalises into [0, 360). The result matches the spec’s definition of position.bearing: degrees clockwise from true north, 0 = north, 90 = east. This is the direction of actual travel between fixes, which is what an operations map wants even when it differs slightly from the vehicle’s instantaneous heading.

Speed ceiling. min(dist / dt, MAX_SPEED_MPS) clamps the result at a mode-plausible maximum, absorbing the occasional GPS teleport where a fix jumps hundreds of metres. Tune the ceiling per mode — a commuter-rail feed warrants a higher cap than a city bus one.

Verification and Output

Validate the derived fields against physical plausibility and against any feed-supplied values that do exist, so you catch unit or formula mistakes before they reach a map.

python
snapshot_one, cache = derive_speed_bearing(poll_one)
snapshot_two, cache = derive_speed_bearing(poll_two, previous=cache)

derived = snapshot_two.dropna(subset=["speed_derived", "bearing_derived"])

# Physical bounds: non-negative, capped speed; bearing on the compass circle
assert (derived["speed_derived"] >= 0).all(), "Negative speed derived"
assert (derived["speed_derived"] <= MAX_SPEED_MPS).all(), "Speed exceeds ceiling"
assert derived["bearing_derived"].between(0.0, 360.0).all(), "Bearing off the circle"

# Where the feed did supply speed, the derived value should be broadly consistent
if "speed" in derived.columns:
    both = derived.dropna(subset=["speed"])
    if not both.empty:
        gap = (both["speed"] - both["speed_derived"]).abs()
        assert gap.median() < 5.0, "Derived speed diverges from feed-reported speed"

print(
    f"Derived {len(derived)} vehicles | "
    f"median {derived['speed_derived'].median() * 3.6:.1f} km/h"
)

Sample output across two consecutive polls 20 seconds apart:

text
INFO: Derived speed/bearing for 0/312 vehicles
INFO: Derived speed/bearing for 297/314 vehicles
Derived 297 vehicles | median 18.4 km/h

The first poll fills nothing because no prior fix exists; the second fills nearly every vehicle. The 17 that stay null are new vehicle_id values that appeared only in the second poll, or vehicles whose delta fell below the jitter threshold.

Gotchas and Edge Cases

  • Stationary GPS wander. A parked vehicle’s coordinate drifts by a few metres between polls. Without MIN_MOVE_METRES, that drift reads as slow motion and its bearing spins randomly. Clamping small moves to zero speed keeps dwelling vehicles visibly stopped.

  • Timestamp gaps from missed polls. If a vehicle drops out for several minutes and returns, the delta spans the whole gap and the straight-line haversine underestimates true path length around curves, deflating speed. Discard pairs whose dt exceeds a staleness limit and treat the return as a fresh track.

  • Bearing at low speed is unreliable. Heading derived from a tiny displacement is dominated by noise. Suppress the bearing_derived value whenever the move is near the MIN_MOVE_METRES floor rather than publishing a jittery arrow.

  • Overnight timestamp continuity. Unlike the schedule times in stop_times.txt, realtime timestamp values are absolute POSIX seconds and roll over midnight cleanly, so no special hours-past-24 handling is needed here — a contrast worth remembering when the same pipeline also parses static schedule times.

Frequently Asked Questions

Why compute speed instead of reading position.speed?

Many GTFS-Realtime feeds never populate position.speed or position.bearing because the agency’s AVL hardware does not report them. When those fields are absent the protobuf decodes them to defaults, so the only way to recover instantaneous speed and heading is to compute them from successive fixes for the same vehicle.id.

Is haversine accurate enough for transit distances?

Yes. Over the tens to hundreds of metres a vehicle moves between polls, the great-circle haversine formula is accurate to well under a metre, which is far tighter than the GPS noise in the fixes themselves. A projected planar distance gives the same answer at this scale.

How do I avoid absurd speeds from GPS jitter?

Reject pairs with a timestamp delta below a few seconds, cap the derived speed at a plausible maximum for the mode, and optionally smooth over a short rolling window. A stationary vehicle whose GPS wanders a few metres will otherwise report a phantom crawl or an impossible sprint.