Downsampling High-Frequency GPS Tracks Without Losing Path Integrity

Compress high-frequency GPS telemetry without destroying path integrity by pairing the Ramer–Douglas–Peucker (RDP) geometric simplification algorithm with a kinematic stop-point retention pass. Project coordinates to a local metric CRS, run RDP to drop geometrically redundant points, then force-retain any point where the time delta or instantaneous velocity signals a stop or sharp maneuver. The result is a 60–95% smaller dataset that preserves road curvature, dwell locations, and directional shifts—the features that matter most downstream.

Why Naive Decimation Destroys the Analysis

High-frequency telemetry (1–10 Hz) introduces massive geometric redundancy: a vehicle at constant speed on a straight road generates 10 identical coordinates per second. Sampling rate optimization addresses this by treating redundancy as compressible without treating topology as compressible. The two are not the same.

Uniform decimation (df.iloc[::10]) or fixed-interval time-slicing does not distinguish between a redundant straight-line point and an inflection point that defines an intersection. It will discard both at equal probability, systematically biasing cumulative distance, turn-angle distributions, and speed profiles. Those biases propagate to every downstream stage: spatial joins become geometrically inaccurate, stop-detection algorithms miss dwell events, and routing algorithms encounter simplified paths that skip real road geometry.

The root cause sits in the broader context of spatiotemporal data foundations and structures: raw telemetry is only valid when its temporal resolution matches the analytical task and its spatial representation preserves the topological features required by downstream models. Stripping either dimension carelessly breaks the data contract.

The Two-Pass Compression Pipeline

The following four steps define a deterministic pipeline that survives production edge cases.

  1. Project to metric CRS. Convert WGS84 coordinates to a UTM zone so all distance thresholds operate in meters, not degrees.
  2. Apply RDP spatial simplification. Use Shapely’s simplify() with a calibrated meter tolerance to remove points that are geometrically redundant.
  3. Compute temporal and velocity masks. Flag force-retain points where the inter-point time delta exceeds a stop threshold or instantaneous velocity falls below a stationary limit.
  4. Merge index sets and return. Union the RDP-retained indices with the temporal-force indices, deduplicate, and slice the original DataFrame.
Two-pass GPS downsampling pipeline Four sequential stages: Project to UTM, RDP Spatial Simplification, Temporal and Velocity Masking, Merge and Return. Arrows connect each stage left to right. Project to UTM CRS RDP Spatial Simplification Temporal & Velocity Mask Merge & Return

Production-Ready Python Implementation

The function below is complete and runnable. It requires pandas, numpy, geopandas, and shapely. All distance arithmetic happens in projected metric coordinates; no tolerance is ever applied to raw EPSG:4326 degrees.

PYTHON 
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import LineString


def downsample_gps_track(
    df: pd.DataFrame,
    tolerance_m: float = 10.0,
    min_stop_seconds: float = 30.0,
    velocity_threshold_ms: float = 0.5,
) -> pd.DataFrame:
    """
    Downsample a high-frequency GPS track while preserving path integrity.

    Uses two-pass compression: RDP spatial simplification in a projected
    metric CRS, followed by temporal and kinematic stop-point retention.

    Parameters
    ----------
    df : pd.DataFrame
        Must contain columns: timestamp (parseable), lat (float), lon (float).
        Any additional columns are preserved in the output.
    tolerance_m : float
        RDP perpendicular-distance tolerance in metres. Points within this
        distance of the line segment connecting their neighbours are dropped.
    min_stop_seconds : float
        Force-retain any point where the gap since the previous point exceeds
        this duration. Preserves dwell events that RDP would collapse.
    velocity_threshold_ms : float
        Force-retain any point where instantaneous velocity (m/s) falls below
        this value. Preserves stop entries/exits and slow manoeuvres.

    Returns
    -------
    pd.DataFrame
        Chronologically sorted, index-reset subset of the input rows.
    """
    if df.empty or len(df) < 3:
        return df.copy()

    # 1. Chronological sort and timestamp normalisation
    df = df.sort_values("timestamp").reset_index(drop=True)
    df["timestamp"] = pd.to_datetime(df["timestamp"])

    # 2. Project to a local UTM zone so tolerance_m is in true metres.
    #    EPSG:4326 degrees must NEVER be used as a distance unit.
    gdf = gpd.GeoDataFrame(
        df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326"
    )
    lon_center = df["lon"].mean()
    lat_center = df["lat"].mean()
    utm_zone = int((lon_center + 180) / 6) + 1
    epsg_utm = 32600 + utm_zone if lat_center >= 0 else 32700 + utm_zone
    gdf_metric = gdf.to_crs(epsg=epsg_utm)

    # 3. RDP simplification via Shapely on projected coordinates.
    #    preserve_topology=True keeps endpoints and prevents ring collapse.
    metric_coords = np.column_stack(
        (gdf_metric.geometry.x.to_numpy(), gdf_metric.geometry.y.to_numpy())
    )
    line = LineString(metric_coords)
    simplified = line.simplify(tolerance=tolerance_m, preserve_topology=True)

    # 4. Map simplified vertices back to original row indices.
    rdp_coords = np.array(simplified.coords)
    rdp_indices: list[int] = []
    for pt in rdp_coords:
        dists = np.linalg.norm(metric_coords - pt, axis=1)
        rdp_indices.append(int(np.argmin(dists)))
    rdp_indices = sorted(set(rdp_indices))

    # 5. Temporal and kinematic masks (operate on projected geometry).
    time_diff_s = df["timestamp"].diff().dt.total_seconds()
    dist_m = gdf_metric.geometry.distance(gdf_metric.geometry.shift(1))
    velocity_ms = dist_m / time_diff_s.replace(0.0, np.nan)

    temporal_mask = (
        (time_diff_s >= min_stop_seconds)   # long gap → force-retain
        | (velocity_ms <= velocity_threshold_ms)  # near-stationary → force-retain
        | time_diff_s.isna()                # first row → always retain
    )
    temporal_indices = sorted(temporal_mask[temporal_mask].index.tolist())

    # 6. Union index sets, deduplicate, and slice the original DataFrame.
    final_indices = sorted(set(rdp_indices) | set(temporal_indices))
    return df.iloc[final_indices].reset_index(drop=True)

Validation Block

Run these checks immediately after downsampling a representative batch to confirm the pipeline is behaving correctly before you deploy it to a production workload.

  • Point count: assert len(result) < len(original) * 0.9 — if fewer than 10% of points were removed, your tolerance is too tight; if fewer than 5% remain, it is too loose.
  • Distance drift: abs(haversine_total(result) - haversine_total(original)) / haversine_total(original) < 0.02 — cumulative route length should deviate less than 2% for logistics use cases.
  • No timestamp inversions: assert result["timestamp"].is_monotonic_increasing — the output must remain chronologically ordered.
  • Endpoint preservation: The first and last rows of the original track must always appear in the result (Shapely’s preserve_topology=True guarantees this for the spatial pass; verify the temporal pass does not strip them).
  • Stop-event coverage: For each known dwell location in your ground-truth dataset, confirm at least one retained point falls within tolerance_m of the stop centroid.
PYTHON 
def validate_downsampled_track(
    original: pd.DataFrame,
    result: pd.DataFrame,
    tolerance_m: float,
) -> None:
    """Raise AssertionError with a descriptive message on any validation failure."""
    assert not result.empty, "Result is empty"
    assert len(result) < len(original), "No points were removed"
    assert result["timestamp"].is_monotonic_increasing, "Timestamps are not monotonic"
    assert result.iloc[0]["timestamp"] == original.iloc[0]["timestamp"], (
        "First point not preserved"
    )
    assert result.iloc[-1]["timestamp"] == original.iloc[-1]["timestamp"], (
        "Last point not preserved"
    )
    compression = 1.0 - len(result) / len(original)
    if compression < 0.10:
        raise AssertionError(
            f"Compression ratio {compression:.1%} < 10% — tolerance_m={tolerance_m} may be too tight"
        )

Parameter Calibration by Transport Mode

Tolerance selection depends on the geometric complexity of the route network and the kinematic profile of the asset. The table below gives empirically derived starting points; always validate against your specific dataset before locking values into a pipeline config.

Use Case tolerance_m min_stop_seconds velocity_threshold_ms Typical Compression
Urban pedestrian / micro-mobility 3–5 m 15 s 0.3 m/s 40–65%
Urban fleet / last-mile delivery 8–12 m 30 s 0.5 m/s 60–80%
Highway logistics / long-haul 20–35 m 60 s 0.8 m/s 80–92%
Off-road / agricultural machinery 50–100 m 120 s 0.5 m/s 85–95%
Maritime / vessel AIS data 100–250 m 300 s 0.2 m/s 90–97%

When operating near UTM zone boundaries (every 6° of longitude) or at latitudes above 60°N/S, replace the auto-detected UTM EPSG code with a custom equal-area projection centred on your track’s bounding box. Coordinate reference system mapping covers the full projection-selection workflow and pyproj integration patterns for these edge cases.

Common Mistakes and Gotchas

  • Applying the tolerance directly to EPSG:4326 coordinates. Shapely’s simplify() is unit-agnostic; in latitude/longitude degrees a 10-unit tolerance eliminates entire countries. Always project first.
  • Using iterrows() to map simplified vertices back to indices. At 1 Hz over 8-hour shifts, a single vehicle generates ~28,800 rows. iterrows() is 50–200× slower than the vectorised np.linalg.norm approach shown above.
  • Ignoring the first row’s NaN time delta. diff() produces NaN for the first row. If you forget to include time_diff_s.isna() in the temporal mask, the first coordinate is not force-retained and can be dropped by the RDP pass, severing the track from its origin.
  • Setting preserve_topology=False in Shapely’s simplify(). This allows the algorithm to collapse coordinate rings and produce self-intersecting geometries. Always use preserve_topology=True for trajectory data.
  • Treating velocity_threshold_ms as a hard stop detector. This threshold retains points near stops; it does not replace a proper stop-detection algorithm. For authoritative dwell segmentation, feed the downsampled output into a dedicated stop-detection pipeline — see GPS precision and error handling for how noise levels affect velocity estimates near stationary events.
  • Skipping temporal validation when the track has signal gaps. A cellular-based device entering a tunnel generates a multi-minute gap. Without min_stop_seconds, both endpoints of the gap may be dropped by RDP, creating a phantom straight-line segment that never occurred on the road network.

Frequently Asked Questions

Why does RDP collapse dwell events even when I set a tight tolerance?

RDP is purely geometric. A vehicle stopped at a loading dock for 3 minutes generates hundreds of near-identical coordinates; every one of them lies within a few metres of its neighbours. They all satisfy the tolerance criterion and get legally dropped. The temporal mask (min_stop_seconds) exists precisely to override this: it force-retains one representative point from each dwell region regardless of geometric proximity.

Can I use this function with a movingpandas TrajectoryCollection?

Yes. Apply downsample_gps_track() to each trajectory’s underlying DataFrame before constructing Trajectory objects, or use a collection-level apply pattern. Downsampling after construction is possible but requires re-building the spatial index. Pre-process the raw DataFrames first; it is both faster and safer.

How do I validate turn preservation at intersections?

Compute bearing change at each retained point (np.arctan2(dy, dx) in projected coordinates) and compare the distribution of turn angles against the original. A well-calibrated tolerance preserves the shape of the bearing-change histogram. If sharp-turn bins shrink after simplification, lower tolerance_m or add an angular-threshold mask that force-retains points where the bearing change exceeds a configured limit (e.g., 20°).

What compression ratio is realistic for 1 Hz urban fleet data?

With tolerance_m=10, min_stop_seconds=30, velocity_threshold_ms=0.5, expect 65–80% point reduction on typical urban routes. Highway data at the same frequency compresses further (80–90%) because long straight segments are geometrically trivial. Parking-heavy datasets compress less because stop-point retention forces many points to survive the spatial pass.

Back to Sampling Rate Optimization