How to Calculate Speed, Distance and Direction Along a Track
On this page
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. Measured on 46 real OpenStreetMap traces over Brighton and Hove, 55,316 points logged at a median interval of one second (Example 2):
| Keeping one fix | Total length | Change |
|---|---|---|
| every fix | 163.17 km | — |
| every 10 s | 158.95 km | −2.6% |
| every 60 s | 143.27 km | −12.2% |
And derived speed on the same traces reached 48 m/s — 172 km/h, above the UK's 113 km/h motorway limit — on a dataset whose median speed was 4.4 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"])) % 360Four 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).
What you get — from The data used on this page
60,000 points downloaded: 0 duplicates within a track, 2 anywhere
46 tracks, 55,316 points
Project, group, guard zero intervals, average bearings as vectors
1. Project first
np.hypot on degrees returns degrees, and one degree of longitude at 53° north is 67 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.yFor 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.
All three columns, and geodesic distance for long tracks
The data used on this page
Public GPS traces from the OpenStreetMap trackpoints API for a box over central Brighton and Hove, −0.20 to −0.10° E, 50.81 to 50.86° N: the first 12 pages of 5,000 points. The setup loads them the way How to load GPS tracks in Python does — one track per <trkseg>, exact duplicates dropped, elements with more than 10% repeated timestamps set aside, tracks under ten points dropped — and projects them to UTM zone 30N, adding x, y and the interval dt in seconds. The API serves the newest uploads first, so the setup reads the pages as saved on 9 October 2026. Nothing on this page prints a location or a username.
import glob
import io
import os
import xml.etree.ElementTree as ET
import zipfile
import geopandas as gpd
import numpy as np
import pandas as pd
import requests
# The OSM trackpoints API for this box, pages 0-11, as saved on 9 October 2026. The API serves the
# newest uploads first, so a live download changes whenever someone uploads a trace here. Live:
# requests.get("https://api.openstreetmap.org/api/0.6/trackpoints",
# params={"bbox": "-0.20,50.81,-0.10,50.86", "page": page})
if not os.path.exists("traces/page_11.gpx"):
reply = requests.get("https://data.spatialworkflow.io/osm-gps-traces-brighton-2026-10-09-v2.zip",
timeout=120)
reply.raise_for_status()
with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
archive.extractall("traces")
NS = {"g": "http://www.topografix.com/GPX/1/0"}
rows = []
for n, path in enumerate(sorted(glob.glob("traces/*.gpx"))):
for t, trk in enumerate(ET.parse(path).getroot().findall("g:trk", NS)):
for s, seg in enumerate(trk.findall("g:trkseg", NS)):
rows += [(f"p{n}:{t}-{s}", float(p.get("lat")), float(p.get("lon")),
p.findtext("g:time", namespaces=NS)) for p in seg.findall("g:trkpt", NS)]
raw = pd.DataFrame(rows, columns=["track_id", "lat", "lon", "time"]).dropna()
raw["time"] = pd.to_datetime(raw["time"], utc=True, format="ISO8601")
within = raw.duplicated(["track_id", "time", "lat", "lon"]).sum()
anywhere = raw.duplicated(["time", "lat", "lon"]).sum()
print(f"{len(raw):,} points downloaded: {within:,} duplicates within a track, {anywhere:,} anywhere")
raw = raw.drop_duplicates(["time", "lat", "lon"]).sort_values(["track_id", "time"], kind="stable")
zero = raw.groupby("track_id")["time"].apply(lambda s: (s.diff().dt.total_seconds() == 0).mean())
size = raw.groupby("track_id").size()
keep = zero.index[(zero <= 0.1) & (size >= 10)]
df = gpd.GeoDataFrame(raw[raw.track_id.isin(keep)], crs=4326,
geometry=gpd.points_from_xy(raw.lon, raw.lat)[raw.track_id.isin(keep).to_numpy()])
df = df.to_crs(32630).reset_index(drop=True)
df["x"], df["y"] = df.geometry.x, df.geometry.y
df["dt"] = df.groupby("track_id")["time"].diff().dt.total_seconds()
print(f"{df.track_id.nunique()} tracks, {len(df):,} points")60,000 points downloaded: 0 duplicates within a track, 2 anywhere
46 tracks, 55,316 points
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 dfmotion = derive_motion(df.drop(columns="dt"))
fastest = motion["speed"].max()
print(f" median {motion['speed'].median() * 3.6:.1f} km/h, "
f"fastest step {fastest:.0f} m/s ({fastest * 3.6:.0f} km/h)") 1 zero intervals -> speed left as NaN
163.17 km, median speed 1.23 m/s, p99 12.03 m/s
median 4.4 km/h, fastest step 48 m/s (172 km/h)
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 — what the sampling interval does to distance
def length_at(df, every_s=None, id_col="track_id"):
"""Track length keeping only the first fix in each every_s bucket."""
if every_s:
elapsed = df.groupby(id_col)["time"].transform(lambda s: s - s.iloc[0])
df = df.groupby([df[id_col], elapsed.dt.total_seconds() // every_s]).head(1)
g = df.groupby(id_col)
return np.hypot(g["x"].diff(), g["y"].diff()).sum() / 1000
native = length_at(df)
print(f"native (median {df['dt'].median():.0f} s) {native:8.2f} km")
for every in (10, 60):
km = length_at(df, every)
print(f"every {every:3d} s {km:8.2f} km {km / native - 1:+.1%}")native (median 1 s) 163.17 km
every 10 s 158.95 km -2.6%
every 60 s 143.27 km -12.2%
Example 3 — 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 dfgeo = geodesic_motion(df.drop(columns="dt"))
print(f" projected (UTM 30N) length {motion['dist'].sum() / 1000:.3f} km") geodesic length 163.157 km
projected (UTM 30N) length 163.173 km
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. For a city-sized box inside one zone the choice hardly matters: the two totals above differ by 16 m in 163 km.
Example 4 — 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
secs = pd.Series((sub.index - sub.index[0]).total_seconds(), index=sub.index)
first = lambda col: col.rolling(f"{int(window_s)}s", min_periods=2).apply(
lambda a: a[0], raw=True)
x0, y0, t0 = first(sub["x"]), first(sub["y"]), first(secs)
elapsed = secs - 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 resultwindowed = robust_speed(motion) windowed speed p99 11.82 m/s vs adjacent-fix p99 12.03 m/s
Measuring displacement over a nine-second window instead of between fixes a second apart divides the noise contribution by nine, because the position error is unchanged while the baseline is longer. The cost is that brief accelerations are smoothed away. On these traces the 99th percentile barely moves — 12.03 to 11.82 m/s — because the top 1% is real driving rather than noise; the window changes the individual values, not the tail.
Noisy speed, circular bearings and acceleration
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 as the interval grows: 163.17 km at the native one second, 143.27 km keeping one fix a minute.
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 values of a derived-speed distribution are where the noise collects: the traces measured here had a median of 1.23 m/s, a 99th percentile of 12.03 m/s — 43 km/h, ordinary driving — and a maximum of 48 m/s, which no road in the box allows.
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.
arctan2 order, turn angles and Doppler speed
- Sort by identity and time, then
groupbybeforediff. - Project first. Degrees are not metres, and not equally so in each direction.
dt == 0gives infinite speed. Replace with NaN.arctan2(dx, dy)for compass bearing, notarctan2(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. On 46 real traces logged every second, keeping one fix a minute shortened the total by 12.2%.
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.