---
name: trajectory-analysis
description: Analyse and implement GPX/GPS trajectory-processing logic safely and correctly.
---
# Trajectory Analysis Skill
Reusable domain guidance for working with GPX files, geodesic distance calculations, timestamp handling, GPS anomaly detection, and trajectory test design.
---
## 1. GPX Document Structure
A GPX file is an XML document. The canonical hierarchy is:
```
← one or more tracks
← one or more track segments per track
← one or more trackpoints per segment
… ← optional elevation in metres
← optional ISO-8601 timestamp
```
Key rules:
- A **track** (``) groups logically related movement (e.g. a single activity).
- A **track segment** (``) groups *continuously recorded* points. A gap between segments is intentional and meaningful — the recorder was paused or the signal was lost.
- A **trackpoint** (``) is a single position fix with mandatory `lat` and `lon` attributes.
---
## 2. Default XML Namespace Handling
GPX 1.1 files typically declare a default namespace:
```xml
```
When this is present, every element is in that namespace. Queries that ignore it will silently match nothing.
**Always extract the namespace from the root element and use it in all tag lookups.** A robust pattern:
```python
import xml.etree.ElementTree as ET
tree = ET.parse("track.gpx")
root = tree.getroot()
# Extract namespace from the root tag, e.g. "{http://www.topografix.com/GPX/1/1}gpx"
ns = root.tag.split("}")[0].lstrip("{") if "}" in root.tag else ""
ns_map = {"gpx": ns} if ns else {}
def tag(local: str) -> str:
return f"{{{ns}}}{local}" if ns else local
# Usage
for trkpt in root.iter(tag("trkpt")):
lat = float(trkpt.attrib["lat"])
lon = float(trkpt.attrib["lon"])
```
---
## 3. Preserving Segment Boundaries
Segment boundaries encode intentional gaps in recording. **Never flatten all trackpoints from all segments into a single list** before processing distances or speeds — this creates phantom distances across segment boundaries.
Correct approach:
- Process distances and speeds **within** each segment independently.
- When aggregating totals (e.g. total distance), sum per-segment results.
- When reporting anomalies, always include the segment index so results can be mapped back to the source.
```python
for seg_idx, segment in enumerate(track.segments):
for i in range(1, len(segment.points)):
prev = segment.points[i - 1]
curr = segment.points[i]
# compute distance between prev and curr (same segment only)
```
---
## 4. Latitude and Longitude Validation
Always validate coordinate values before feeding them into distance or bearing formulas.
Valid ranges:
- Latitude: **−90.0 ≤ lat ≤ 90.0**
- Longitude: **−180.0 ≤ lon ≤ 180.0**
Additional checks worth applying:
- Reject `NaN` and `±Inf`.
- Treat `(0.0, 0.0)` — the null island — as suspicious unless the track is genuinely in the Gulf of Guinea.
- Reject values that are outside the declared bounding box (``) when one is present in the GPX file.
```python
def is_valid_coord(lat: float, lon: float) -> bool:
return (
-90.0 <= lat <= 90.0
and -180.0 <= lon <= 180.0
and not (math.isnan(lat) or math.isnan(lon))
and not (math.isinf(lat) or math.isinf(lon))
)
```
---
## 5. Distance Calculations: Haversine vs WGS-84 Ellipsoidal Geodesic
Never approximate GPS distance with naive Euclidean degree arithmetic. One degree of latitude is ~111 km, but one degree of longitude varies from ~111 km at the equator to 0 km at the poles. Degree-based Euclidean distance produces wildly incorrect results at higher latitudes and for any east–west movement.
**Two legitimate approaches — choose based on accuracy requirements:**
### Haversine (spherical approximation)
Haversine models the Earth as a perfect sphere. It is simple, dependency-free, and accurate to within ~0.5 % for most practical GPS distances. It is a valid and reasonable choice for anomaly detection, trip summaries, and relative distance comparisons where sub-percent error is acceptable.
```python
import math
EARTH_RADIUS_M = 6_371_000.0 # mean radius in metres
def haversine_m(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
"""Great-circle distance in metres (spherical Earth approximation)."""
phi1, phi2 = math.radians(lat1), math.radians(lat2)
d_phi = math.radians(lat2 - lat1)
d_lam = math.radians(lon2 - lon1)
a = math.sin(d_phi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(d_lam / 2) ** 2
return 2 * EARTH_RADIUS_M * math.asin(math.sqrt(a))
```
### WGS-84 ellipsoidal geodesic
WGS-84 models the Earth as an oblate spheroid. When a project requires accurate geodesic distance — survey-grade measurements, long-haul route totals, or any context where the ~0.5 % spherical error is unacceptable — use a geodesic library rather than Haversine:
```python
# geopy — Vincenty/Karney algorithm on the WGS-84 ellipsoid
from geopy.distance import geodesic
dist_m = geodesic((lat1, lon1), (lat2, lon2)).meters
# pyproj — also ellipsoidal, useful when already using a projection pipeline
from pyproj import Geod
geod = Geod(ellps="WGS84")
_, _, dist_m = geod.inv(lon1, lat1, lon2, lat2)
```
**Guidance:** use Haversine when simplicity and no extra dependencies matter; switch to an ellipsoidal library when the application explicitly requires WGS-84 geodesic accuracy. Do not mix the two approaches within the same distance pipeline.
---
## 6. Timestamp Parsing and UTC-Normalised Comparison
GPX timestamps are ISO-8601 strings, typically `2024-06-01T10:32:45Z` (UTC) or with an explicit offset like `2024-06-01T12:32:45+02:00`.
Rules for valid timestamps:
- Always parse into **timezone-aware** `datetime` objects.
- Normalise to UTC before any arithmetic.
- Never compare a naive datetime to an aware one.
- `timedelta` arithmetic on UTC-normalised datetimes gives correct elapsed seconds regardless of DST or timezone offset in the source file.
```python
from datetime import datetime, timezone
def parse_gpx_time(ts: str) -> datetime:
"""Parse an ISO-8601 GPX timestamp into a UTC-aware datetime."""
dt = datetime.fromisoformat(ts.replace("Z", "+00:00"))
return dt.astimezone(timezone.utc)
def elapsed_seconds(t1: datetime, t2: datetime) -> float:
return (t2 - t1).total_seconds()
```
**Handling timestamp problems — record, don't discard:**
A trackpoint's `