Exporting smoothed trajectories to GeoParquet

Exporting a smoothed trajectory correctly means four things must survive the write: the coordinates, the coordinate reference system, the UTC time axis, and the provenance of how the data was cleaned. The reliable pattern is to wrap the smoothed points in a GeoDataFrame with an explicit metric CRS, write partitioned GeoParquet with GeoDataFrame.to_parquet so the CRS lands in file metadata, and immediately round-trip read the dataset to prove the CRS and row count came back intact. Skipping that verification is how teams discover, weeks later, that every stored geometry silently defaulted to EPSG:4326.

Why this happens

A smoothed trajectory is the output of a cleaning stage — typically the drift correction described in handling GPS drift in raw trajectory logs — and by that point the coordinates already sit in a metric projected CRS. The failure mode is that Parquet has no native concept of a coordinate system: the CRS is carried only in an optional file-level metadata key. Any writer that treats geometry as an ordinary binary column drops that key, and the file opens later as if it were plain WGS84.

The full treatment of encodings, partition schemes, and row-group sizing lives in the parent spatial storage formats guide, and the broader modelling context in spatiotemporal data foundations and structures. This page narrows to one job: getting a single cleaned track onto disk as GeoParquet with its CRS and lineage guaranteed, ready to be consumed as a trajectory object.

Core export pipeline

The four steps below run in order — you cannot embed a CRS you have not declared, and you cannot verify a round trip you have not written.

Smoothed Trajectory Export Pipeline Four sequential stages: Build GeoDataFrame, Attach Metadata, Write GeoParquet, and Round-trip Verify, connected by arrows. 1. Build GeoDataFrame 2. Attach Metadata 3. Write GeoParquet 4. Round-trip Verify
  1. Build a validated GeoDataFrame. Turn the smoothed x/y (or corrected lon/lat) columns into point geometry, declare the CRS explicitly, and confirm timestamps are UTC-aware.
  2. Attach lineage metadata. Add cleaning_version and source columns so provenance is queryable per row, and derive a date partition column from the UTC timestamp.
  3. Write partitioned GeoParquet. Write with GeoDataFrame.to_parquet so the CRS is embedded, partitioning by date so time-range scans prune directories.
  4. Round-trip read to verify. Read the dataset back and assert row count, CRS equality, and geometry validity before the data is published.

Production-ready Python implementation

The function below takes a smoothed frame — the corrected output of the drift-correction stage — and exports it. It stores geometry in a metric projected CRS (a UTM zone derived from the data, never raw EPSG:4326) so that any later distance or area predicate is correct, and it handles the edge cases that silently corrupt stored datasets: a missing CRS, an empty frame, and mixed geometry types.

PYTHON 
from __future__ import annotations

from pathlib import Path

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point


def export_smoothed_trajectory(
    df: pd.DataFrame,
    dataset_root: str | Path,
    *,
    cleaning_version: str,
    source: str,
    x_col: str = "x_m",
    y_col: str = "y_m",
    ts_col: str = "ts",
    metric_epsg: int = 32633,   # UTM zone 33N — set to the zone the smoother used
) -> Path:
    """
    Export a smoothed/drift-corrected trajectory to partitioned GeoParquet.

    Parameters
    ----------
    df : Smoothed trajectory with projected coordinate columns (``x_col``,
         ``y_col``, in metres), a UTC-aware ``ts_col``, and a ``traj_id`` column.
    dataset_root : Directory root of the partitioned dataset.
    cleaning_version : Tag of the cleaning/smoothing pipeline that produced ``df``.
    source : Provenance string embedded as a column and in file metadata.
    x_col, y_col : Projected coordinate columns in the metric CRS (metres).
    ts_col : UTC-aware timestamp column.
    metric_epsg : EPSG of the metric projected CRS the coordinates are in.
                  NEVER 4326 — x/y here are metres, not degrees.

    Returns
    -------
    Path to the dataset root that was written.

    Raises
    ------
    ValueError
        On missing columns, an empty frame, timezone-naive timestamps, or a
        non-metric target CRS.
    """
    required = {x_col, y_col, ts_col, "traj_id"}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")

    # Edge case: empty frame. A zero-row write yields no partition dirs and an
    # incomplete schema, so fail loudly instead of writing a silent empty set.
    if df.empty:
        raise ValueError("Refusing to export an empty trajectory.")

    if metric_epsg == 4326:
        raise ValueError("metric_epsg must be a projected CRS, not 4326.")

    work = df.copy()

    # Edge case: timezone-naive timestamps break date partitioning.
    work[ts_col] = pd.to_datetime(work[ts_col], utc=True)

    # ── Step 1: build point geometry in the metric CRS ────────────────────────
    # Coordinates are already metres in the projected CRS; build Points directly
    # and declare the CRS so the write can embed it. NEVER label metre values as
    # EPSG:4326 — that mislabels a projected frame as geographic.
    geometry = [Point(xy) for xy in zip(work[x_col], work[y_col])]
    gdf = gpd.GeoDataFrame(work, geometry=geometry, crs=f"EPSG:{metric_epsg}")

    # Edge case: guard against mixed geometry types sneaking in via ``work``.
    geom_types = set(gdf.geometry.geom_type.unique())
    if geom_types != {"Point"}:
        raise ValueError(f"Expected only Point geometry, found: {geom_types}")

    if gdf.crs is None:
        raise ValueError("GeoDataFrame CRS is undefined after construction.")

    # ── Step 2: attach lineage + partition columns ────────────────────────────
    gdf["cleaning_version"] = cleaning_version
    gdf["source"] = source
    gdf["date"] = gdf[ts_col].dt.strftime("%Y-%m-%d")

    # ── Step 3: write partitioned GeoParquet with embedded CRS ────────────────
    # GeoDataFrame.to_parquet serialises the CRS into the file-level ``geo``
    # metadata; partition_cols=["date"] lets time-range scans prune directories.
    root = Path(dataset_root)
    root.mkdir(parents=True, exist_ok=True)
    gdf.to_parquet(
        root,
        partition_cols=["date"],
        write_covering_bbox=True,
        index=False,
    )
    return root

Validation block

Run these checks immediately after the export. Each targets a distinct corruption mode — a dropped CRS, a truncated write, or invalid geometry — so a clean pass gives real confidence before the dataset is handed downstream.

PYTHON 
def verify_geoparquet_export(
    dataset_root: str,
    original: gpd.GeoDataFrame,
    expected_epsg: int = 32633,
) -> gpd.GeoDataFrame:
    """
    Round-trip read a GeoParquet dataset and assert it matches expectations.

    Raises AssertionError on any mismatch; returns the reloaded frame on success.
    """
    reloaded = gpd.read_parquet(dataset_root)

    # 1. CRS survived the write and equals the metric CRS we intended.
    assert reloaded.crs is not None, "CRS was lost on write — not read back."
    assert reloaded.crs.to_epsg() == expected_epsg, (
        f"CRS mismatch: wrote EPSG:{expected_epsg}, read {reloaded.crs.to_epsg()}"
    )

    # 2. No rows were dropped by partitioning or the write.
    assert len(reloaded) == len(original), (
        f"Row count changed: {len(original)} → {len(reloaded)}"
    )

    # 3. Geometry is present, valid, and homogeneous points.
    assert reloaded.geometry.notna().all(), "Null geometry after round trip."
    assert reloaded.geometry.is_valid.all(), "Invalid geometry after round trip."
    assert set(reloaded.geometry.geom_type.unique()) == {"Point"}, (
        "Geometry type changed on round trip."
    )

    # 4. Lineage metadata survived as queryable columns.
    for col in ("cleaning_version", "source", "date"):
        assert col in reloaded.columns, f"Lineage column missing: {col}"

    print(
        f"Verified {len(reloaded)} rows, CRS EPSG:{reloaded.crs.to_epsg()}, "
        f"{reloaded['date'].nunique()} date partition(s)."
    )
    return reloaded

The checklist the code enforces:

  • CRS equality: the recovered EPSG must equal the one written — the single most common silent failure, since a dropped CRS reads back as None or defaults elsewhere in the stack.
  • Row-count contract: a partitioned write that drops rows usually means a null in a partition column; the count check catches it before analytics do.
  • Geometry validity and type: invalid or unexpectedly mixed geometry indicates an axis swap or a stray non-point row leaking in from upstream.

Common mistakes and gotchas

  • Writing with pyarrow.parquet.write_table. It has no geospatial awareness and never emits the geo metadata key, so the CRS is dropped and geometry becomes an opaque binary column. Always write geometry through GeoDataFrame.to_parquet.
  • Mislabelling projected metres as EPSG:4326. After smoothing, x_m/y_m are metres in a UTM zone. Constructing the GeoDataFrame with crs="EPSG:4326" tags a projected frame as geographic, and every later distance query is wrong by orders of magnitude. Declare the actual metric EPSG.
  • Leaving timestamps timezone-naive. A naive ts produces date partitions shifted by the local UTC offset, so a query for one calendar day silently reads the wrong files. Convert to UTC before deriving the partition column.
  • Over-partitioning a single track. Partitioning one short trajectory by date and hour can create dozens of directories each holding a handful of rows — the small-files problem. For a single track, a date-only partition (or none) keeps files large enough to be worth opening. The parent spatial storage formats guide covers when finer partitioning actually pays off.
  • Skipping the round-trip read. A write that returns without error is not proof the CRS or all rows survived. The verification step is cheap relative to re-running an analysis on a silently corrupted dataset — never skip it in a production export.

Back to Spatial Storage Formats