How to Load GPS Tracks into Python as Trajectories

Problem statement

Loading a GPX file is three lines. Loading it correctly means answering four questions the file does not answer for you:

  1. Which points belong to one journey?
  2. Are the timestamps timezone-aware and monotonic?
  3. What CRS should the coordinates be in before anything is measured?
  4. Which rows are duplicates, and which are genuinely two fixes in one second?

Question one is the one that bites. Downloading 60,000 public trackpoints from the OpenStreetMap API produced 68 <trkseg> elements β€” of which 53 contained roughly ten separate traces each, interleaved. Trusting the file's own structure would have produced 53 tracks that teleport between ten simultaneous journeys.

Quick answer

import numpy as np
import pandas as pd
import geopandas as gpd
import xml.etree.ElementTree as ET

NS = {"g": "http://www.topografix.com/GPX/1/1"}


def read_gpx(path):
    """Every trkpt with its track and segment index preserved."""
    root = ET.parse(path).getroot()
    ns = NS if root.tag.endswith("}gpx") and "1/1" in root.tag else \
        {"g": "http://www.topografix.com/GPX/1/0"}

    rows = []
    for ti, trk in enumerate(root.findall("g:trk", ns)):
        name = trk.findtext("g:name", default="", namespaces=ns)
        for si, seg in enumerate(trk.findall("g:trkseg", ns)):
            for pt in seg.findall("g:trkpt", ns):
                time = pt.findtext("g:time", namespaces=ns)
                rows.append({
                    "track": ti, "segment": si, "name": name,
                    "lat": float(pt.get("lat")), "lon": float(pt.get("lon")),
                    "ele": float(pt.findtext("g:ele", default="nan", namespaces=ns)),
                    "time": time,
                    "hdop": pt.findtext("g:hdop", namespaces=ns),
                })

    df = pd.DataFrame(rows)
    df["time"] = pd.to_datetime(df["time"], utc=True, format="ISO8601")
    df["track_id"] = df["track"].astype(str) + "-" + df["segment"].astype(str)
    return df
Loading a GPX file through parsing, timestamp handling, identity validation, projection and derived columns.
Parsing is the easy part. Identity and time are where files disagree with themselves.

Step-by-step solution

1. Parse without assuming the namespace

GPX exists in versions 1.0 and 1.1 with different namespace URIs. A hard-coded namespace silently returns zero points for the other version β€” findall returns an empty list rather than raising.

Detect the namespace from the root element, or strip namespaces entirely before searching.

2. Make the timestamps timezone-aware immediately

df["time"] = pd.to_datetime(df["time"], utc=True, format="ISO8601")

GPX timestamps are UTC by specification, but exporters write local times, offsets and occasionally nothing at all. utc=True normalises what is there; format="ISO8601" stops pandas inferring a format per row, which is slow and occasionally inconsistent.

Points with no timestamp are not trajectory points. Keep them if you need the geometry, but exclude them from anything temporal.

3. Validate the identity column rather than trusting it

def check_identity(df, id_col="track_id"):
    for key, sub in df.groupby(id_col):
        sub = sub.sort_values("time")
        dt = sub["time"].diff().dt.total_seconds()
        zero = float((dt == 0).mean())
        if zero > 0.1:
            print(f"  {key}: {zero:.1%} of rows share a timestamp β€” "
                  f"about {round(1 / (1 - zero))} traces merged")
  4-0: 89.9% of rows share a timestamp β€” about 10 traces merged
  6-0: 89.9% of rows share a timestamp β€” about 10 traces merged

A large fraction of zero intervals means several traces occupy one element. Filtering to the clean ones left 13 tracks and 6,417 points from a 60,000-point download β€” 10.7%.

4. Project before deriving anything

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat), crs=4326)
gdf = gdf.to_crs(gdf.estimate_utm_crs())
gdf["x"], gdf["y"] = gdf.geometry.x, gdf.geometry.y

estimate_utm_crs() picks the UTM zone containing the data's centroid. For tracks spanning several zones, choose an equidistant projection centred on the data instead.

5. Derive the between-row columns once, in the right order

gdf = gdf.sort_values(["track_id", "time"]).reset_index(drop=True)
g = gdf.groupby("track_id")
gdf["dt"] = g["time"].diff().dt.total_seconds()
gdf["dist"] = np.hypot(g["x"].diff(), g["y"].diff())
gdf["speed"] = gdf["dist"] / gdf["dt"].replace(0, np.nan)

Sort by identity and time, group before diffing, and guard the zero interval. Each of those three prevents a specific, silent, wrong answer.

A trkseg containing ten interleaved traces shown as repeated timestamps, against a clean single trace with monotonic three-second intervals.
Ninety percent of intervals equal to zero is the fingerprint of about ten merged traces.

Code examples

Example 1 β€” a loader that reports what it found

import numpy as np
import pandas as pd
import geopandas as gpd


def load_tracks(df, id_col="track_id", min_points=10, crs=None):
    """Validated, projected trajectories plus a per-track summary."""
    df = df.dropna(subset=["time", "lat", "lon"]).copy()
    df["time"] = pd.to_datetime(df["time"], utc=True)

    before = len(df)
    df = df.drop_duplicates(subset=[id_col, "time", "lat", "lon"])
    if len(df) < before:
        print(f"  dropped {before - len(df):,} exact duplicate rows "
              f"({1 - len(df) / before:.1%})")

    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat),
                           crs=4326)
    gdf = gdf.to_crs(crs or gdf.estimate_utm_crs())
    gdf["x"], gdf["y"] = gdf.geometry.x, gdf.geometry.y
    gdf = gdf.sort_values([id_col, "time"]).reset_index(drop=True)

    g = gdf.groupby(id_col)
    gdf["dt"] = g["time"].diff().dt.total_seconds()
    gdf["dist"] = np.hypot(g["x"].diff(), g["y"].diff())
    gdf["speed"] = gdf["dist"] / gdf["dt"].replace(0, np.nan)

    summary = g.agg(points=("time", "size"), start=("time", "min"),
                    end=("time", "max"))
    summary["duration_s"] = (summary["end"] - summary["start"]).dt.total_seconds()
    summary["length_km"] = g["dist"].sum() / 1000
    summary["median_dt_s"] = g["dt"].median()
    summary["zero_dt_frac"] = g["dt"].apply(lambda s: float((s == 0).mean()))

    suspect = summary[summary["zero_dt_frac"] > 0.1]
    if len(suspect):
        print(f"  {len(suspect)} of {len(summary)} tracks look interleaved")
    keep = summary[(summary["points"] >= min_points) &
                   (summary["zero_dt_frac"] <= 0.1)]
    print(f"  keeping {len(keep)} tracks, {int(keep['points'].sum()):,} points, "
          f"{keep['length_km'].sum():.2f} km")
    return gdf[gdf[id_col].isin(keep.index)], summary
  dropped 1,388 exact duplicate rows (2.3%)
  53 of 68 tracks look interleaved
  keeping 13 tracks, 6,417 points, 37.86 km

Example 2 β€” reading other formats onto the same shape

import pandas as pd
import geopandas as gpd


def load_any(path, id_col=None, time_col=None, lon=None, lat=None):
    """CSV, GeoPackage, Parquet or GPX into one trajectory table."""
    suffix = path.rsplit(".", 1)[-1].lower()

    if suffix in ("csv", "txt"):
        df = pd.read_csv(path)
    elif suffix == "gpx":
        return read_gpx(path)
    elif suffix in ("parquet", "pq"):
        df = pd.read_parquet(path)
    else:
        gdf = gpd.read_file(path)
        gdf["lon"], gdf["lat"] = gdf.geometry.x, gdf.geometry.y
        df = pd.DataFrame(gdf.drop(columns="geometry"))

    rename = {}
    if time_col: rename[time_col] = "time"
    if lon: rename[lon] = "lon"
    if lat: rename[lat] = "lat"
    if id_col: rename[id_col] = "track_id"
    df = df.rename(columns=rename)

    missing = {"time", "lon", "lat"} - set(df.columns)
    if missing:
        raise ValueError(f"missing columns {sorted(missing)}; "
                         f"available: {sorted(df.columns)[:12]}")
    if "track_id" not in df:
        print("  no identity column β€” treating the whole file as one track")
        df["track_id"] = "all"
    return df

Normalising every source to track_id, time, lon, lat early means the rest of the pipeline has one shape to handle. The alternative β€” special cases threaded through every function β€” is how trajectory code becomes unmaintainable.

Example 3 β€” reconstructing tracks when the file has no identity

import numpy as np


def infer_tracks(df, max_gap_s=600, max_speed_ms=56.0):
    """Build track identity from continuity when the file provides none."""
    df = df.sort_values("time").reset_index(drop=True).copy()
    dt = df["time"].diff().dt.total_seconds()
    dist = np.hypot(df["x"].diff(), df["y"].diff())
    speed = dist / dt.replace(0, np.nan)

    new_track = (dt > max_gap_s) | (speed > max_speed_ms) | dt.isna()
    df["track_id"] = new_track.cumsum()

    sizes = df.groupby("track_id").size()
    print(f"  inferred {len(sizes)} tracks; "
          f"{int((sizes < 10).sum())} have fewer than 10 points")
    return df

This works when a single device logged continuously. It cannot separate simultaneous traces from different devices, which is what the OpenStreetMap trace API produces β€” for that, no continuity rule helps, because the points genuinely interleave in time.

Explanation

Why the file's own structure is not trustworthy

GPX has <trk>, <trkseg> and <trkpt>, and the specification says a <trkseg> is a contiguous run of points. Producers interpret that loosely:

  • Some emit one segment per file regardless of gaps.
  • Some emit a new segment on every GPS loss, so one journey becomes twenty.
  • Anonymising services merge points from many uploads into one segment.

The validation in step 3 costs four lines and catches all three.

Why duplicate rows are common and worth removing

A 60,000-point download contained 1,388 exact (lat, lon, time) duplicates β€” 2.3%. They come from overlapping API pages, from devices writing a point twice, and from files concatenated more than once.

They matter because a duplicate produces a zero interval and therefore an infinite derived speed, and because they double-weight one location in any density or clustering analysis.

Why timezone-naive timestamps break silently

Two tracks recorded either side of a daylight-saving change, stored as naive local times, sort into the wrong order for an hour. diff() produces negative intervals, speeds go negative, and nothing raises.

Parsing with utc=True at load time makes this impossible. If the source really is local time with no offset, attach the zone explicitly with tz_localize and let it raise on the ambiguous hour rather than guessing.

Why to keep points rather than lines

Converting to a LineString per track is convenient and lossy: a line has no time, so speed, stops and duration cannot be recovered from it.

Keep points as the primary representation. Build lines for display, and keep them as a derived product that can be rebuilt after any filtering step.

Five things a trajectory loader should report, ending with which tracks were kept and why.
Of a 60,000-point download, 6,417 points in 13 tracks were usable trajectories.

Edge cases or notes

  • Detect the GPX namespace. Version 1.0 and 1.1 differ, and a wrong namespace returns zero points silently.
  • Parse timestamps with utc=True and format="ISO8601".
  • Validate the identity column by checking the fraction of zero intervals.
  • Drop exact duplicates on (id, time, lat, lon) β€” 2.3% of a real download.
  • Project before measuring. estimate_utm_crs() is a good default for local data.
  • groupby before diff, always.
  • Keep hdop, sat and fix where the file has them; they beat inferred quality.
  • GPS elevation is poor. Sample a DEM instead for anything about gradient.

FAQ

How do I read a GPX file in Python?

Parse it with xml.etree.ElementTree, detecting the namespace from the root element, and collect trkpt elements with their track and segment indices. Then parse timestamps with utc=True and project.

Can I trust each trkseg to be one track?

No. Of 68 segments from the OpenStreetMap trace API, 53 contained about ten interleaved traces each. Validate by checking what fraction of intervals are zero.

Why do I get infinite speeds after loading?

Duplicate rows or two fixes in the same second give a zero interval. Drop exact duplicates and replace zero intervals with NaN before dividing.

What CRS should GPS tracks be in?

WGS 84 for storage, a local projected CRS for any measurement. estimate_utm_crs() picks a sensible zone for local data.

How do I handle points with no timestamp?

Keep them for geometry if you need it, but exclude them from anything temporal. A point without a time is not a trajectory point.

Should I convert tracks to LineStrings on load?

No. A LineString discards the time dimension. Keep points as the primary representation and build lines for display.

What if my file has no track identifier?

Infer one from continuity β€” a new track wherever the time gap or implied speed is impossible. That works for a single device logging over time, but cannot separate simultaneous traces.