How to Calculate Speed, Distance and Direction Along a Track #

Problem statement #

Speed, distance and heading are the three derived columns every trajectory analysis needs, and all three are computed between rows rather than within them. That makes them sensitive to four things the raw data does not advertise: the sort order, the group boundaries, the coordinate system and the sampling interval.

The sampling interval alone changes the answer substantially. The same 13 real traces:

native (median 3 s)   37.86 km
every 10 s            33.55 km   -11.4%
every 60 s            29.67 km   -21.6%

And derived speed on the same traces reached 175 m/s — 630 km/h — on a dataset whose median speed was 4.7 km/h.

Quick answer #

import numpy as np

df = df.sort_values(["track_id", "time"]).reset_index(drop=True)
g = df.groupby("track_id")

df["dt"] = g["time"].diff().dt.total_seconds()
df["dx"] = g["x"].diff()
df["dy"] = g["y"].diff()
df["dist"] = np.hypot(df["dx"], df["dy"])
df["speed"] = df["dist"] / df["dt"].replace(0, np.nan)
df["bearing"] = np.degrees(np.arctan2(df["dx"], df["dy"])) % 360

Four details, each preventing a specific wrong answer: sort by identity and time; groupby so track boundaries do not leak; projected coordinates so hypot gives metres; replace(0, np.nan) so a zero interval does not divide.

Note the argument order in arctan2(dx, dy). Compass bearing is measured clockwise from north, which is the opposite convention to the mathematical arctan2(y, x).

Step-by-step solution #

1. Project first #

np.hypot on degrees returns degrees, and one degree of longitude at 53° north is 66 km against 111 km for latitude. Distances computed that way are wrong by a factor that varies with direction.

gdf = gdf.to_crs(gdf.estimate_utm_crs())
gdf["x"], gdf["y"] = gdf.geometry.x, gdf.geometry.y

For tracks spanning more than a UTM zone or two, use a geodesic calculation instead — pyproj.Geod.inv gives distance and forward azimuth on the ellipsoid directly.

2. Group before differencing #

df["x"].diff() on a table containing several tracks computes the difference across every boundary, producing one enormous displacement per track transition.

df.groupby("track_id")["x"].diff() gives NaN at each boundary instead, which is correct — the first point of a track has no predecessor.

3. Handle the zero interval #

Two fixes in one second give dt == 0 and an infinite speed. Replacing with NaN keeps the row and marks the value as unavailable, which is more honest than dropping the row or substituting zero.

4. Prefer the receiver's speed when it exists #

Derived speed inherits position error from both endpoints. With 5 m error and a 1-second interval that is about 7 m/s of noise — 25 km/h, on top of a possibly slower true speed.

Many receivers record a Doppler-derived speed, which is an independent measurement and far less noisy. If your data has it, use it and keep the derived value as a cross-check.

5. Average bearings as vectors, never as numbers #

The mean of 350° and 10° is 180° if you average the numbers, and 0° if you average the directions. Always convert to unit vectors first:

def mean_bearing(degrees):
    radians = np.radians(np.asarray(degrees, dtype=float))
    return float(np.degrees(np.arctan2(np.sin(radians).mean(),
                                       np.cos(radians).mean())) % 360)

6. Report the sampling interval with any distance #

Distance is not a property of the journey alone. State the interval, and resample before comparing datasets.

Code examples #

Example 1 — all three columns, with the traps handled #

import numpy as np
import pandas as pd


def derive_motion(df, id_col="track_id", time_col="time",
                  x="x", y="y", speed_col=None):
    """Distance, speed, bearing and acceleration between consecutive fixes."""
    df = df.sort_values([id_col, time_col]).reset_index(drop=True).copy()
    g = df.groupby(id_col)

    df["dt"] = g[time_col].diff().dt.total_seconds()
    dx = g[x].diff()
    dy = g[y].diff()
    df["dist"] = np.hypot(dx, dy)

    safe_dt = df["dt"].replace(0, np.nan)
    df["speed_derived"] = df["dist"] / safe_dt
    df["speed"] = df[speed_col] if speed_col in df else df["speed_derived"]
    df["bearing"] = np.degrees(np.arctan2(dx, dy)) % 360
    df["accel"] = g["speed"].diff() / safe_dt

    # turn angle: signed change in bearing, wrapped to [-180, 180]
    turn = g["bearing"].diff()
    df["turn"] = ((turn + 180) % 360) - 180

    n_zero = int((df["dt"] == 0).sum())
    if n_zero:
        print(f"  {n_zero:,} zero intervals -> speed left as NaN")
    print(f"  {df['dist'].sum() / 1000:.2f} km, "
          f"median speed {df['speed'].median():.2f} m/s, "
          f"p99 {df['speed'].quantile(0.99):.2f} m/s")
    return df

The turn column is the one people forget to wrap. A heading change from 350° to 10° is a 20° right turn, not a 340° left one.

Example 2 — geodesic distance for long tracks #

import numpy as np
from pyproj import Geod

GEOD = Geod(ellps="WGS84")


def geodesic_motion(df, id_col="track_id", lon="lon", lat="lat"):
    """Ellipsoidal distance and forward azimuth — no projection needed."""
    df = df.sort_values([id_col, "time"]).reset_index(drop=True).copy()
    lon1 = df.groupby(id_col)[lon].shift()
    lat1 = df.groupby(id_col)[lat].shift()

    valid = lon1.notna()
    az, _, dist = GEOD.inv(lon1[valid].values, lat1[valid].values,
                           df[lon][valid].values, df[lat][valid].values)

    df.loc[valid, "dist"] = dist
    df.loc[valid, "bearing"] = az % 360
    df["dt"] = df.groupby(id_col)["time"].diff().dt.total_seconds()
    df["speed"] = df["dist"] / df["dt"].replace(0, np.nan)

    print(f"  geodesic length {df['dist'].sum() / 1000:.3f} km")
    return df

Use this when a track crosses UTM zones, or when you want distances that do not depend on a projection choice at all. It is slower than hypot on projected coordinates and exact on the ellipsoid.

Example 3 — robust speed, smoothing before differencing #

import numpy as np
import pandas as pd


def robust_speed(df, id_col="track_id", window_s=9.0):
    """Speed over a longer baseline, which suppresses position noise."""
    out = []
    for key, sub in df.groupby(id_col):
        sub = sub.sort_values("time").set_index("time")

        # displacement over a window, rather than between adjacent fixes
        x0 = sub["x"].rolling(f"{int(window_s)}s", min_periods=2).apply(
            lambda s: s.iloc[0], raw=False)
        y0 = sub["y"].rolling(f"{int(window_s)}s", min_periods=2).apply(
            lambda s: s.iloc[0], raw=False)
        t0 = pd.Series(sub.index, index=sub.index).rolling(
            f"{int(window_s)}s", min_periods=2).apply(lambda s: s.iloc[0])

        elapsed = (sub.index.astype("int64") / 1e9) - t0
        sub["speed_window"] = np.hypot(sub["x"] - x0, sub["y"] - y0) / \
            elapsed.replace(0, np.nan)
        out.append(sub.reset_index())

    result = pd.concat(out, ignore_index=True)
    print(f"  windowed speed p99 {result['speed_window'].quantile(0.99):.2f} m/s "
          f"vs adjacent-fix p99 {result['speed'].quantile(0.99):.2f} m/s")
    return result

Measuring displacement over a nine-second window instead of a three-second one cuts the noise contribution by a factor of three, because the position error is unchanged while the baseline is longer. The cost is that brief accelerations are smoothed away.

Explanation #

Why distance depends on the sampling interval #

Between two fixes, the path is assumed straight, so every deviation shorter than the interval is discarded. Measured length therefore falls monotonically as the interval grows: 37.86 km at 3 seconds, 29.67 km at 60.

There is no interval-free "true" distance. Report the interval alongside the number, and resample datasets to a common interval before comparing them.

Why derived speed is noisier than it looks #

Position error enters twice — once for each endpoint — and is then divided by the interval. Roughly, σ_speed ≈ σ_position × √2 / dt.

Two consequences that surprise people. Shorter intervals give noisier speeds, so a 1 Hz log looks more erratic than a 10-second log of the same journey. And the extreme percentiles of a derived-speed distribution are almost entirely noise: the traces measured here had a median of 1.32 m/s and a maximum of 175 m/s.

Why bearings need circular statistics #

Bearings live on a circle, where 359° and 1° are two degrees apart. Ordinary arithmetic treats them as 358 apart, which breaks means, standard deviations, differences and any binning that crosses north.

The general fix is to convert to unit vectors, do the arithmetic there, and convert back. For differences, wrap into [−180, 180] with ((d + 180) % 360) - 180.

The resultant length of the mean vector is also useful in its own right: near 1 means consistent heading, near 0 means the directions cancel — a vehicle circling a roundabout, or a receiver stationary and jittering.

Why acceleration is the most useful derived column #

Speed alone cannot distinguish a fast vehicle from a displaced point. Acceleration can, because real vehicles are limited to a few metres per second squared while a position error produces hundreds.

It is also the basis for the most effective outlier test: a value that is extreme in both speed and acceleration is an error, while one extreme only in speed may be genuine fast travel.

Edge cases or notes #

• Sort by identity and time, then groupby before diff.
• Project first. Degrees are not metres, and not equally so in each direction.
• dt == 0 gives infinite speed. Replace with NaN.
• arctan2(dx, dy) for compass bearing, not arctan2(dy, dx).
• Wrap turn angles into [−180, 180].
• Average bearings as vectors, never as numbers.
• Prefer the receiver's Doppler speed over derived speed if present.
• Do not derive gradient from GPS altitude; vertical error is roughly twice horizontal.

Internal links #

• Trajectories explained: what makes a track different from points — why these are between-row columns
• GPS error explained: why your track wanders — the noise in derived speed
• Sampling rate and gaps explained — why distance depends on the interval
• How to clean a GPS track: outliers, duplicates and jumps — using acceleration to find errors
• How to load GPS tracks into Python as trajectories — the columns these depend on
• How to measure distance accurately in Python: geodesic vs projected — the projection question in general
• How to resample and interpolate a track to a fixed interval — making distances comparable
• My track has impossible speeds or teleports — when these columns go wrong

FAQ #

How do I calculate speed from GPS points in Python? #

Sort by track and time, group before differencing, compute distance with np.hypot on projected coordinates, and divide by the time difference with zero intervals replaced by NaN.

Why is my calculated speed so noisy? #

Because it divides the error of two positions by a short interval. With 5 m error at 1-second sampling that is about 7 m/s of noise regardless of the true speed.

Should I use the speed column from the device? #

Yes, if it exists. It is usually Doppler-derived, which is an independent measurement and much less noisy than a difference of positions.

How do I calculate a bearing? #

np.degrees(np.arctan2(dx, dy)) % 360 for compass bearing measured clockwise from north — note that dx comes first, unlike the mathematical convention.

Why is the average of my bearings wrong? #

Because bearings are circular. Average 350° and 10° numerically and you get 180°. Convert to unit vectors, average those, and convert back.

Why does my track distance change when I resample? #

Because straight lines between fixes cut corners. Measured length fell 21.6% between 3-second and 60-second sampling on real traces.

Should I use geodesic or projected distance? #

Projected is fine and fast for local tracks. Use pyproj.Geod.inv for tracks spanning several UTM zones or when you want no projection dependence at all.