Resampling high-frequency telemetry natively in Polars

Downsampling raw GPS telemetry in Polars comes down to one primitive: group_by_dynamic, which bins a sorted timestamp axis into fixed-width windows and aggregates each window in a single parallel pass. Raw fleet and wearable devices log at 5–50 Hz, which is both noisier and far larger than any downstream analysis needs, so the first stage of most mobility pipelines collapses that stream to a stable 1 s or 5 s cadence. Polars does this faster than pandas, spills to disk through lazy frames when the data exceeds memory, and keeps every device separated through a single group_by argument — but only if you normalize the time axis first and defer all speed and distance work until after a projection to a metric coordinate system.

Why resampling comes first

High-frequency telemetry is rarely uniform. Devices drop samples during power-saving states, buffer bursts over flaky cellular links, and emit duplicate epochs when a firmware clock resets. Feeding that irregular stream directly into speed derivation or stay-point detection produces phantom accelerations and unstable window statistics, because every derivative divides by an inconsistent dt. Resampling to a regular grid is the deterministic fix, and it is the entry point to the wider sampling rate optimization discipline that trades resolution against storage and compute.

Polars is the right engine for this at scale. Its Rust query engine is multi-threaded and columnar, and its lazy API reorders and prunes work before any bytes are read — the same reasons it outperforms pandas on the large-file preprocessing that precedes trajectory modeling in Spatiotemporal Data Foundations & Structures. Where pandas materializes the full frame in RAM and resamples on one core, Polars streams the file in chunks and parallelizes the windowed aggregation across every core.

Core resampling pipeline

The pipeline has four stages. The order matters: you cannot bin an unsorted axis, and you must not derive motion before projecting.

Polars Resampling Pipeline Four sequential stages: Normalize Time Axis, group_by_dynamic Binning, Project to Metric CRS, and Lazy Streaming Scale-out, connected by arrows. 1. Normalize Time Axis 2. Bin with group_by_dynamic 3. Project to Metric CRS 4. Lazy Scale-out
  1. Normalize the time axis. Cast the timestamp to a UTC Datetime, sort by device_id then time, and resolve duplicate epochs. group_by_dynamic assumes a monotonic index within each group.
  2. Bin with group_by_dynamic. Choose a window (every="1s") and aggregate each bin — a mean or median position, a point count, and any sensor channels you carry forward.
  3. Project to a metric CRS. Only after binning do you transform coordinates to a UTM zone and derive distance and speed, because kinematics in degree-space are not metric.
  4. Scale with lazy frames. Swap read_parquet for scan_parquet and collect with the streaming engine so the same code runs out-of-core on files larger than memory.

Production-ready Polars implementation

The function below downsamples multi-device telemetry to a fixed cadence. It handles the empty frame, unsorted and duplicate timestamps, and the single-device versus multi-device cases through one group_by argument. Distance and speed are computed only after projecting each binned point to a metric UTM CRS derived from the data centroid — never from raw longitude and latitude. The projection uses pyproj with always_xy=True to guarantee lon-before-lat ordering; the same UTM zone selection logic is covered in CRS transformation best practices.

PYTHON 
from __future__ import annotations

import polars as pl
from pyproj import Transformer


def resample_telemetry(
    df: pl.DataFrame,
    every: str = "1s",
    time_col: str = "timestamp",
    lon_col: str = "lon",
    lat_col: str = "lat",
    device_col: str | None = "device_id",
) -> pl.DataFrame:
    """
    Downsample high-frequency GPS telemetry to a fixed cadence in Polars.

    Bins each device's track into fixed-width windows with group_by_dynamic,
    then derives distance and speed in a metric projected CRS (UTM zone
    chosen from the data centroid) — NOT from raw lon/lat degrees.

    Parameters
    ----------
    df : pl.DataFrame
        Must contain time_col, lon_col, lat_col (and device_col if given).
    every : str
        Target window width as a Polars duration string ("1s", "5s", "1m").
    device_col : str | None
        Column separating independent tracks. None for a single device.

    Returns
    -------
    pl.DataFrame
        One row per (device, window): binned lon/lat, point count,
        x_m/y_m metric coordinates, step distance, and speed_ms.

    Raises
    ------
    ValueError
        If required columns are missing.
    """
    required = {time_col, lon_col, lat_col}
    if device_col is not None:
        required.add(device_col)
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")

    # Edge case: empty input — return an empty frame with the output schema.
    if df.height == 0:
        return df.clear()

    # ── Stage 1: normalize the time axis ──────────────────────────────
    # Cast to UTC Datetime; sort within device; drop duplicate epochs.
    # group_by_dynamic REQUIRES a monotonic index per group.
    sort_keys = [device_col, time_col] if device_col else [time_col]
    df = (
        df.with_columns(pl.col(time_col).cast(pl.Datetime("us", "UTC")))
        .sort(sort_keys)
        .unique(subset=sort_keys, keep="first", maintain_order=True)
    )

    # ── Stage 2: bin with group_by_dynamic ────────────────────────────
    # group_by isolates each device so no window spans two vehicles.
    agg = (
        df.group_by_dynamic(
            index_column=time_col,
            every=every,
            group_by=device_col if device_col else None,
        )
        .agg(
            pl.col(lon_col).mean().alias(lon_col),
            pl.col(lat_col).mean().alias(lat_col),
            pl.len().alias("n_points"),
        )
        .sort(sort_keys)
    )

    # Edge case: a single point produces one bin — no motion to derive.
    if agg.height < 2:
        return agg.with_columns(
            x_m=pl.lit(None, dtype=pl.Float64),
            y_m=pl.lit(None, dtype=pl.Float64),
            step_m=pl.lit(None, dtype=pl.Float64),
            speed_ms=pl.lit(None, dtype=pl.Float64),
        )

    # ── Stage 3: project to a metric CRS, THEN derive kinematics ──────
    # NEVER compute distance/speed on EPSG:4326 degrees. Derive the UTM
    # zone from the centroid so points sit in the low-distortion corridor.
    c_lon = agg[lon_col].mean()
    c_lat = agg[lat_col].mean()
    utm_zone = int((c_lon + 180) // 6) + 1
    utm_epsg = (32600 if c_lat >= 0 else 32700) + utm_zone
    to_utm = Transformer.from_crs("EPSG:4326", f"EPSG:{utm_epsg}", always_xy=True)
    x, y = to_utm.transform(agg[lon_col].to_numpy(), agg[lat_col].to_numpy())

    over = device_col if device_col else pl.lit(1)
    agg = agg.with_columns(
        x_m=pl.Series("x_m", x),
        y_m=pl.Series("y_m", y),
    ).with_columns(
        # Metre deltas and elapsed seconds within each device track.
        dx=(pl.col("x_m") - pl.col("x_m").shift(1)).over(over),
        dy=(pl.col("y_m") - pl.col("y_m").shift(1)).over(over),
        dt_s=(pl.col(time_col) - pl.col(time_col).shift(1))
        .over(over)
        .dt.total_nanoseconds()
        / 1e9,
    ).with_columns(
        step_m=(pl.col("dx") ** 2 + pl.col("dy") ** 2).sqrt(),
    ).with_columns(
        # Guard the first row of each track (dt is null) against div-by-zero.
        speed_ms=pl.when(pl.col("dt_s") > 0)
        .then(pl.col("step_m") / pl.col("dt_s"))
        .otherwise(None),
    )

    return agg.drop(["dx", "dy", "dt_s"])

For files larger than memory, keep the identical logic but enter through a LazyFrame. scan_parquet never reads the whole file, and collect(streaming=True) executes the query in chunks:

PYTHON 
import polars as pl

# scan_parquet builds a lazy plan; nothing is read until collect().
lazy = (
    pl.scan_parquet("s3://fleet/telemetry/*.parquet")
    .with_columns(pl.col("timestamp").cast(pl.Datetime("us", "UTC")))
    .sort(["device_id", "timestamp"])
    .group_by_dynamic(
        index_column="timestamp",
        every="5s",
        group_by="device_id",
    )
    .agg(
        pl.col("lon").mean(),
        pl.col("lat").mean(),
        pl.len().alias("n_points"),
    )
)

# The streaming engine processes the scan in chunks — out-of-core safe.
downsampled = lazy.collect(streaming=True)

To go the other direction — increasing cadence rather than reducing it — use upsample, which inserts empty rows on a regular grid that you then fill. This is the mechanism behind downsampling high-frequency GPS tracks without losing path integrity, where a resample-then-refill step re-establishes a uniform axis:

PYTHON 
import polars as pl

# upsample fills GAPS on a fixed grid; forward-fill then carries the
# last known position. Interpolation of coordinates must still happen
# in a metric CRS if you derive distance from the result.
regular = (
    df.sort("timestamp")
    .upsample(time_column="timestamp", every="1s", group_by="device_id")
    .with_columns(pl.col(["lon", "lat"]).forward_fill().over("device_id"))
)

Validation block

Run these checks after resampling. Each targets a distinct failure: a broken time grid, a silently dropped device, or an implausible derived speed that signals a degree-space calculation slipped through.

PYTHON 
import polars as pl


def validate_resampled(
    raw: pl.DataFrame,
    out: pl.DataFrame,
    every_seconds: int = 1,
    device_col: str = "device_id",
    time_col: str = "timestamp",
    max_speed_ms: float = 70.0,
) -> None:
    """
    Sanity-check a resampled telemetry frame.
    Raises AssertionError on failure; prints a summary on success.
    """
    assert out.height > 0, "Resampled frame is empty"

    # 1. No device is silently dropped by the windowing.
    raw_devices = set(raw[device_col].unique().to_list())
    out_devices = set(out[device_col].unique().to_list())
    missing = raw_devices - out_devices
    assert not missing, f"Devices lost during resampling: {missing}"

    # 2. Timestamps are aligned to the target grid within each device.
    grid_ns = every_seconds * 1_000_000_000
    remainder = (
        out.select(
            (pl.col(time_col).dt.epoch("ns") % grid_ns).alias("r")
        )["r"]
    )
    assert (remainder == 0).all(), "Timestamps not aligned to the resample grid"

    # 3. Derived speed is physically plausible — a huge value usually means
    #    kinematics were computed on lon/lat degrees, not a metric CRS.
    if "speed_ms" in out.columns:
        p999 = out["speed_ms"].drop_nulls().quantile(0.999) or 0.0
        assert p999 < max_speed_ms, (
            f"p99.9 speed {p999:.1f} m/s exceeds {max_speed_ms} m/s — "
            "check that distance was computed in a projected CRS"
        )

    reduction = 100 * (1 - out.height / raw.height)
    print(
        f"Validation passed. Rows: {raw.height} -> {out.height} "
        f"({reduction:.1f}% reduction) across {len(out_devices)} devices."
    )

Common mistakes and gotchas

  • Binning an unsorted axis. group_by_dynamic trusts that the index is ascending within each group. If you skip the sort, or sort only globally without partitioning by device_id, windows silently span the wrong rows and per-device tracks bleed into one another. Always sort(["device_id", "timestamp"]) immediately before binning.

  • Leaving duplicate epochs in place. Two rows sharing one timestamp break the monotonic assumption and can produce a dt_s of zero, which turns into an infinite speed. Resolve collisions with unique(subset=[...]) or aggregate them before the resample rather than after.

  • Computing speed on longitude and latitude. A step of 0.001° means something very different in metres at the equator than at 60° latitude. Deriving speed = degrees / seconds yields a number that is neither metric nor comparable across regions. Project each binned point to a metric CRS first — the implementation above derives the UTM zone from the centroid, matching the discipline in Spatiotemporal Data Foundations & Structures.

  • Calling collect() without streaming on huge files. A plain collect() materializes the entire result in memory and defeats the purpose of scan_parquet. For multi-gigabyte logs use collect(streaming=True) so the engine processes the scan in bounded-memory chunks.

  • Mixing shift across device boundaries. A bare pl.col("x_m").shift(1) pulls the last row of the previous device into the first row of the next, injecting a phantom giant step. Every windowed derivative must be wrapped in .over("device_id") so shifts reset at each track boundary.

  • Reaching for pandas out of habit. For a few thousand rows either library is fine, but on the multi-million-row logs typical of fleet telemetry, pandas resampling runs single-threaded and must hold the whole frame in RAM. Polars parallelizes the same aggregation and streams from disk, which is why it is the default preprocessing engine for large sampling rate optimization jobs.

Back to Sampling Rate Optimization

Related in this section