Computing bearing and heading change with NumPy

Bearings are angles on a circle, and every bug in turn analysis comes from arithmetic that treats them as numbers on a line. The three rules that prevent all of them: compute the azimuth with arctan2 on projected coordinates, wrap every difference into ±180° at the point of computation, and never take an ordinary mean of a set of headings — use the circular mean, or the answer for a vehicle oscillating around north is south.

Why this happens

arctan2(dy, dx) returns a value in (−π, π], which maps cleanly to a compass bearing after a rotation and a modulo. The trouble starts one line later. Subtracting consecutive bearings gives a value in (−360°, 360°), and roughly half of that range is wrong: a two-degree turn across north reads as −358°.

The same discontinuity breaks averaging. The mean of 359° and 1° is 180° by ordinary arithmetic and 0° in reality. Any rolling mean, any per-segment average heading, any “typical direction on this corridor” computed with np.mean is wrong whenever the data straddles north — which on a north-south corridor is most of it. The consequences for turn detection are set out in directionality and turn analysis.

The failure is silent and the output is a legal bearing, which is why it survives review. Any aggregation over angles needs the vector form.

Core pipeline

Project to a metric CRS so that dx and dy are metres and the azimuth is geometrically meaningful.
Compute the forward azimuth with arctan2(dx, dy) — note the argument order for a compass bearing.
Wrap every difference into ±180° immediately, before it is stored, summed or smoothed.
Aggregate with circular statistics, never with mean or median on the raw degrees.

Production-ready Python implementation

PYTHON

import numpy as np
import pandas as pd


def bearings_and_changes(
    df: pd.DataFrame,
    x_col: str = "x",
    y_col: str = "y",
    entity_col: str = "entity_id",
    min_step_m: float = 1.0,
) -> pd.DataFrame:
    """
    Vectorised forward azimuth and wrapped heading change per fix.

    Parameters
    ----------
    df : pd.DataFrame
        PROJECTED metric coordinates. Computing an azimuth from degrees gives
        an angle in a distorted space that is wrong by up to tens of degrees
        at high latitude.
    min_step_m : float
        Steps shorter than this produce a meaningless bearing — the direction
        of a 20 cm step is pure noise — so they inherit the previous heading.

    Returns
    -------
    pd.DataFrame
        Input plus 'bearing_deg' in [0, 360) and 'dheading_deg' in (-180, 180].

    Raises
    ------
    ValueError
        On missing columns or an empty frame.
    """
    required = {x_col, y_col, entity_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if df.empty:
        raise ValueError("Input DataFrame is empty.")

    out = df.copy()
    dx = out.groupby(entity_col)[x_col].diff().to_numpy()
    dy = out.groupby(entity_col)[y_col].diff().to_numpy()
    step = np.hypot(dx, dy)

    # arctan2(dx, dy) — NOT (dy, dx). Compass bearings measure clockwise
    # from north, so the easting is the "opposite" and the northing the
    # "adjacent". Getting this backwards mirrors every angle about 45°.
    bearing = np.degrees(np.arctan2(dx, dy)) % 360.0

    # A step below the noise floor has no meaningful direction. Carrying the
    # previous bearing forward is better than emitting a random one, because
    # a random bearing feeds a spurious heading change into the accumulator.
    bearing = np.where(step >= min_step_m, bearing, np.nan)
    out["bearing_deg"] = pd.Series(bearing, index=out.index).groupby(
        out[entity_col]
    ).ffill()

    # Wrap at the point of computation. This is the single most important
    # line on the page: an un-wrapped difference stored anywhere contaminates
    # every cumulative statistic built on top of it.
    raw = out.groupby(entity_col)["bearing_deg"].diff().to_numpy()
    out["dheading_deg"] = (raw + 180.0) % 360.0 - 180.0
    return out


def circular_mean_deg(angles_deg: np.ndarray) -> float:
    """
    Mean direction of a set of headings, in [0, 360).

    Returns NaN when the resultant length is near zero — a set of headings
    pointing in all directions has no meaningful mean, and returning an
    arbitrary one is worse than admitting that.
    """
    a = np.radians(np.asarray(angles_deg, dtype=float))
    a = a[np.isfinite(a)]
    if a.size == 0:
        return float("nan")
    s, c = np.sin(a).mean(), np.cos(a).mean()
    if np.hypot(s, c) < 1e-6:
        return float("nan")
    return float(np.degrees(np.arctan2(s, c)) % 360.0)


def circular_std_deg(angles_deg: np.ndarray) -> float:
    """Angular spread in degrees; 0 means perfectly aligned, ~81 means uniform."""
    a = np.radians(np.asarray(angles_deg, dtype=float))
    a = a[np.isfinite(a)]
    if a.size < 2:
        return float("nan")
    r = np.hypot(np.sin(a).mean(), np.cos(a).mean())
    return float(np.degrees(np.sqrt(-2.0 * np.log(max(r, 1e-12)))))

Validation block

PYTHON

def validate_bearings(out: pd.DataFrame) -> None:
    """Four assertions that catch every version of the wrap bug."""
    b = out["bearing_deg"].dropna()
    d = out["dheading_deg"].dropna()

    # 1. Range. A bearing outside [0, 360) means the modulo was skipped.
    assert b.between(0, 360, inclusive="left").all(), "bearing outside [0, 360)"

    # 2. Wrapped range. Anything outside ±180 is an unwrapped difference.
    assert d.between(-180, 180).all(), (
        f"heading change outside ±180 (max {d.abs().max():.0f}°) — the wrap "
        "was applied after storage rather than at computation"
    )

    # 3. A synthetic north-crossing case must give +2, not -358.
    probe = pd.DataFrame({
        "entity_id": ["p", "p"],
        "x": [0.0, 0.35], "y": [0.0, 10.0],
    })
    got = bearings_and_changes(probe)["dheading_deg"].iloc[-1]
    assert abs(got) < 90, f"north-crossing probe returned {got:.0f}°"

    # 4. The circular mean of two headings straddling north is near 0, not 180.
    assert circular_mean_deg(np.array([359.0, 1.0])) < 5.0
    print(f"OK — {len(b)} bearings, max |Δ| {d.abs().max():.1f}°")

The staircase is diagnostic: steps of almost exactly 360° in a cumulative angle mean the wrap is missing, and they appear only on routes that happen to cross north.

The curve is why the step-length gate exists rather than being an optimisation: below a metre or so, the derived bearing carries no information at all and feeds pure noise into every cumulative angle.

Common mistakes and gotchas

arctan2(dy, dx) for a compass bearing. That is the mathematical convention, measured anticlockwise from east. Compass bearings are clockwise from north, so the arguments swap.
Computing azimuths in degrees of latitude and longitude. The angle comes out in a distorted space and is wrong by up to tens of degrees at high latitude. Project first.
Wrapping late. A difference stored unwrapped and wrapped at read time is fine; a difference stored unwrapped and summed before wrapping is unrecoverable.
np.mean on headings. Use the circular mean. The same applies to rolling means, groupby().mean(), and any interpolation across a gap.
Bearings from sub-metre steps. A stationary vehicle’s jitter produces uniformly distributed bearings and a huge cumulative change. Gate on step length, as the code does.
Assuming the device’s own heading field is safe. Many receivers report the last non-zero heading while stationary, which looks like a vehicle facing one way for ten minutes and then snapping around. Derive it or validate it.

FAQ

Should I use the device heading or derive it from positions?

Derive it, and use the device heading as a cross-check. Receiver heading is usually derived from Doppler and is excellent while moving fast, unreliable below walking pace, and often frozen at the last valid value while stationary. A derived bearing gated on step length behaves consistently across the whole speed range.

How do I average headings over a rolling window?

Convert to unit vectors, take the rolling mean of the sine and cosine components separately, and convert back with arctan2. Pandas rolling works on the components; there is no correct way to do it on the degrees directly.

What about the resultant length?

It is the most useful by-product. The magnitude of the mean unit vector is between 0 and 1 and measures how consistent the headings were — near 1 on a straight motorway, near 0 in a car park. It is a better “is this vehicle travelling purposefully” feature than any speed threshold.

Directionality & Turn Analysis — the parent reference and the manoeuvre taxonomy these angles feed.
Smoothing Noisy Heading Signals for Turn Detection — filtering the series this page produces.
Detecting U-turns and Directional Shifts in Fleet Data — the detector built on the accumulated change.
Calculating Instantaneous Speed from Discrete GPS Points — the sibling derivation from the same projected coordinates.
Best Practices for CRS Transformations in Movement Data — the projection this page assumes.

Back to Directionality & Turn Analysis