Spatial Storage Formats

A spatial storage format is the on-disk encoding — columns, geometry serialisation, partition layout, and embedded metadata — that determines how fast and how correctly a trajectory dataset can be scanned, filtered, and joined by an analytical engine.

Prerequisites Checklist

Before designing a storage layout for trajectory data, confirm your environment and inputs. This topic builds on the modelling and cleaning stages in Spatiotemporal Data Foundations & Structures, and expects trajectories that have already been cleaned and structured.

• Python ≥ 3.10 with geopandas ≥ 1.0, pyarrow ≥ 15, and shapely ≥ 2.0 for geometry handling.
• duckdb ≥ 0.10 with the spatial extension installable at runtime, for querying without a separate database server.
• h3 ≥ 4.0 (optional) if you partition by hexagonal cell — the same discrete grid used in discrete global grid systems.
• Input schema: each row must have a geometry column (point or linestring), a ts (UTC datetime64) column, a traj_id identifier, and any per-point attributes (speed_ms, _qa_flag). Timestamps must be timezone-aware — timezone-naive columns silently break time-partition pruning.
• A declared CRS. The GeoDataFrame must carry an explicit .crs. Geometry stored in a metric projected CRS (a UTM zone) keeps distance and area predicates correct; store the CRS in metadata so downstream readers reproject deliberately rather than by accident.
• A trajectory object contract. Ideally the frame originates from a Trajectory wrapper defined in trajectory object design patterns so identifiers and column names are consistent.

Storage Problem Taxonomy

Most storage bugs are silent: the write succeeds, the file opens, and the corruption only surfaces as wrong query results or a 50× slowdown weeks later. Classifying the failure mode tells you which part of the layout to fix.

Deterministic Pipeline Overview

The five stages below are the standard order for turning an in-memory trajectory into a queryable analytical dataset. Each stage is deterministic — given the same input frame and parameters, the byte layout and metadata are reproducible, which matters for cache invalidation and data lineage.

1. Build a validated GeoDataFrame — assemble the trajectory with an explicit CRS and UTC-aware timestamps; reject frames with a missing or ambiguous CRS before proceeding.
2. Choose a geometry encoding — pick WKB for interoperability or GeoArrow for analytical scans, and pin the choice so every file in the dataset is homogeneous.
3. Partition by space and time — derive partition columns (an H3 cell and a date) so that a typical query touches a small fraction of the directories.
4. Write partitioned GeoParquet — write with embedded CRS, cleaning version, and source metadata, and a row-group size tuned to the workload; one row group should hold one coherent partition slice.
5. Query with DuckDB — read the dataset through the spatial extension with Hive partitioning, letting predicate pushdown on the partition columns and row-group statistics skip irrelevant data.

Implementation Walkthrough

The functions below implement stages 1–4 (write) and stage 5 (query). Geometry is stored in a metric projected CRS (EPSG:32633 by default) so that any spatial predicate — distance, area, containment — stays metrically correct; raw EPSG:4326 is never used for computation, only recorded when the source is geographic.

from __future__ import annotations

import logging
from pathlib import Path

import geopandas as gpd
import pandas as pd

logger = logging.getLogger(__name__)

# Dataset-level contract: every file in one dataset must share these.
GEOMETRY_ENCODING = "WKB"     # "WKB" (interop) or "geoarrow" (analytical)
DEFAULT_METRIC_EPSG = 32633   # UTM zone 33N — adjust to your data's extent


def write_trajectory_geoparquet(
    gdf: gpd.GeoDataFrame,
    dataset_root: str | Path,
    *,
    cleaning_version: str,
    source: str,
    metric_epsg: int = DEFAULT_METRIC_EPSG,
    h3_resolution: int = 7,
    row_group_size: int = 512_000,
    geometry_encoding: str = GEOMETRY_ENCODING,
) -> Path:
    """
    Write a trajectory GeoDataFrame as a space/time-partitioned GeoParquet dataset.

    Parameters
    ----------
    gdf : GeoDataFrame with a geometry column, a UTC-aware ``ts`` column,
          a ``traj_id`` column, and an explicit ``.crs``.
    dataset_root : Directory root of the partitioned dataset.
    cleaning_version : Version tag of the cleaning pipeline that produced ``gdf``.
    source : Provenance string (e.g. feed name) embedded in file metadata.
    metric_epsg : Metric projected CRS the geometry is stored in. NEVER 4326 —
                  distances/areas in degrees are not metric.
    h3_resolution : H3 resolution for the spatial partition key.
    row_group_size : Rows per Parquet row group; size to one coherent slice.
    geometry_encoding : "WKB" or "geoarrow"; recorded in dataset metadata.

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

    Raises
    ------
    ValueError
        If required columns are missing, the frame is empty, the CRS is
        undefined, or timestamps are timezone-naive.
    """
    required = {"geometry", "ts", "traj_id"}
    missing = required - set(gdf.columns)
    if missing:
        raise ValueError(f"GeoDataFrame missing columns: {missing}")

    if gdf.empty:
        raise ValueError("Refusing to write an empty GeoDataFrame.")

    if gdf.crs is None:
        raise ValueError("GeoDataFrame has no CRS; set one before writing.")

    ts = pd.to_datetime(gdf["ts"])
    if ts.dt.tz is None:
        raise ValueError("`ts` is timezone-naive; convert to UTC before writing.")

    out = gdf.copy()
    out["ts"] = ts

    # ── Store geometry in a metric projected CRS ──────────────────────────────
    # NEVER compute distances/areas on EPSG:4326. Reproject to the metric CRS
    # so any spatial predicate applied later stays metrically correct.
    if out.crs.to_epsg() != metric_epsg:
        out = out.to_crs(epsg=metric_epsg)

    # ── Derive partition columns (space + time) ───────────────────────────────
    # H3 keys are computed from geographic lon/lat, so index off a WGS84 view;
    # the stored geometry stays in the metric CRS.
    geo = out.geometry.to_crs(epsg=4326).representative_point()
    try:
        import h3

        out["h3_cell"] = [
            h3.latlng_to_cell(pt.y, pt.x, h3_resolution) for pt in geo
        ]
    except ImportError:
        logger.warning("h3 not installed; falling back to a coarse geohash grid.")
        out["h3_cell"] = (
            geo.y.round(1).astype(str) + "_" + geo.x.round(1).astype(str)
        )
    out["date"] = out["ts"].dt.strftime("%Y-%m-%d")

    # ── Write partitioned GeoParquet with embedded metadata ───────────────────
    # GeoDataFrame.to_parquet serialises the CRS into the file-level ``geo``
    # metadata; the extra kwargs travel as key/value metadata for lineage.
    root = Path(dataset_root)
    root.mkdir(parents=True, exist_ok=True)

    out.to_parquet(
        root,
        partition_cols=["h3_cell", "date"],
        geometry_encoding=geometry_encoding,
        write_covering_bbox=True,        # per-row bbox → stronger spatial pushdown
        row_group_size=row_group_size,
        schema_version="1.1.0",
        index=False,
    )

    # Sidecar lineage file — human-auditable, mirrors what is embedded per file.
    (root / "_dataset_meta.json").write_text(
        pd.Series(
            {
                "cleaning_version": cleaning_version,
                "source": source,
                "metric_epsg": metric_epsg,
                "geometry_encoding": geometry_encoding,
                "h3_resolution": h3_resolution,
                "row_group_size": row_group_size,
            }
        ).to_json(indent=2)
    )

    logger.info("Wrote %d rows to %s", len(out), root)
    return root

The read side uses DuckDB so that no database server is required and the partition columns become first-class query columns:

import duckdb


def query_trajectory_dataset(
    dataset_root: str,
    *,
    h3_cells: list[str] | None = None,
    date_from: str | None = None,
    date_to: str | None = None,
) -> "duckdb.DuckDBPyRelation":
    """
    Query a partitioned GeoParquet trajectory dataset with predicate pushdown.

    Parameters
    ----------
    dataset_root : Root directory written by ``write_trajectory_geoparquet``.
    h3_cells : Optional spatial filter on the ``h3_cell`` partition column.
    date_from, date_to : Optional inclusive date bounds (YYYY-MM-DD).

    Returns
    -------
    A lazy DuckDB relation; call ``.df()`` to materialise a DataFrame.
    """
    con = duckdb.connect()
    con.install_extension("spatial")
    con.load_extension("spatial")

    # hive_partitioning exposes h3_cell/date as columns; a WHERE clause on them
    # prunes whole directories before any Parquet file is opened.
    filters: list[str] = []
    params: list[object] = []
    if h3_cells:
        placeholders = ", ".join(["?"] * len(h3_cells))
        filters.append(f"h3_cell IN ({placeholders})")
        params.extend(h3_cells)
    if date_from:
        filters.append("date >= ?")
        params.append(date_from)
    if date_to:
        filters.append("date <= ?")
        params.append(date_to)

    where = f"WHERE {' AND '.join(filters)}" if filters else ""
    sql = f"""
        SELECT traj_id, ts, h3_cell, date, geometry
        FROM read_parquet('{dataset_root}/**/*.parquet', hive_partitioning = true)
        {where}
        ORDER BY traj_id, ts
    """
    return con.execute(sql, params)

Geometric and Mathematical Grounding

Columnar layout. Parquet stores each column contiguously rather than row-by-row. A query that reads three of thirty columns touches only those three column chunks, so I/O scales with columns read, not table width. Within a column, values are dictionary- or run-length-encoded and then compressed, which is highly effective for trajectory attributes such as traj_id and h3_cell that repeat across many consecutive rows.

Predicate pushdown and statistics. Each row group carries per-column min/max statistics in the file footer. When a query filters date >= '2026-07-01', the reader compares the predicate against each row group’s date statistics and skips any group whose range cannot match — no bytes of that group are read. Pushdown is only as good as the statistics are tight, which is why data should be sorted by the partition and filter columns before writing: a row group that spans the whole date range prunes nothing. The per-row bounding box written with write_covering_bbox=True extends the same idea to geometry, letting a spatial-range query skip row groups whose bbox does not intersect the query window.

Row-group sizing. A row group is the unit of both parallelism and pruning. If N is the total row count and g the row-group size, the file has ⌈N / g⌉ groups. Small g inflates footer metadata and makes each statistics comparison cover too little data to be worth it; large g means a matching predicate still drags in g rows even when only a handful qualify. The practical sweet spot keeps each row group at 128–512 MB uncompressed so that footer overhead stays negligible while pruning granularity remains useful.

Partition cardinality. Partitioning by h3_cell and date creates one directory per (cell, day) pair. If a dataset covers C cells over D days, the directory count approaches C × D. Choosing too fine an H3 resolution explodes C and produces the small-files problem, where per-file open latency dominates. The discrete global grid systems topic covers how resolution maps to cell edge length so you can match partition granularity to query selectivity.

Calibration & Parameter Tuning

Partition scheme and row-group size are the two knobs that most affect scan cost. Tune them to the query pattern, not to the data volume alone.

Order the partition columns most-selective-first for your dominant query. A dashboard that always filters a small map viewport should lead with h3_cell; a reporting job that scans whole days across a region should lead with date. When queries mix both, DuckDB still prunes on either column, but the leading column governs directory-tree depth and therefore the number of files opened.

Match the H3 partition resolution to the same grid used downstream in aggregating GPS pings into H3 hexagon grids so that a stored partition maps one-to-one onto an aggregation cell — no re-indexing at query time.

Integration & Compatibility

This storage layer sits between cleaning and analytics, so it touches several other stages:

• Trajectory Object Design Patterns: the Trajectory object is the natural input; serialise its underlying GeoDataFrame with the writer above and rehydrate it after a DuckDB read so column contracts stay stable.
• Discrete Global Grid Systems: H3 cell IDs used as partition keys come from the same grid as spatial aggregation, so a partition prune and an aggregation bucket share one key space.
• Sampling Rate Optimization: write downsampled tracks rather than raw high-frequency logs — fewer rows per partition means tighter row groups and cheaper scans, without losing path integrity.
• GPS Precision & Error Handling: carry the _qa_flag and cleaning-version columns into storage so a query can filter to only records that passed QA.
• Exporting Smoothed Trajectories to GeoParquet: the step-by-step deep-dive that takes a single drift-corrected track through this writer and verifies the round trip.

FAQ

Should I encode geometry as WKB or GeoArrow in Parquet?

Use WKB when interoperability across GDAL, QGIS, PostGIS, and older readers matters — it is the GeoParquet default and universally understood. Use GeoArrow native encoding for analytical workloads that read columns into Arrow or DuckDB, because GeoArrow stores coordinates in contiguous typed arrays and avoids the per-row cost of decoding WKB blobs. Record which encoding you wrote in the file’s geo metadata so consumers do not have to guess.

How large should a Parquet row group be for trajectory data?

Aim for 128–512 MB of uncompressed data per row group, roughly 0.5–2 million typical trajectory rows. Smaller groups increase metadata overhead and weaken statistics; larger groups coarsen predicate pushdown. Size each row group to hold one coherent slice of a partition — for example one device-day — so column statistics stay tight enough to prune.

Why did my CRS disappear after writing GeoParquet?

CRS lives in the file-level geo metadata, not in the column values. A plain pyarrow.parquet.write_table on a raw WKB column emits no geo metadata and drops the CRS. Use GeoDataFrame.to_parquet, which serialises the CRS, and always round-trip read the file to assert the recovered CRS matches what you intended.

Do I need to partition small trajectory datasets?

No. Partitioning helps when most queries touch a small slice of a large dataset. Under a few hundred megabytes, a single well-sorted file with good row-group statistics already lets DuckDB skip irrelevant groups. Over-partitioning produces thousands of tiny files, and per-file open cost then dominates.

Can DuckDB read a partitioned GeoParquet dataset directly?

Yes. Install and load the spatial extension, then point read_parquet at a glob such as dataset/**/*.parquet with hive_partitioning enabled. DuckDB exposes the partition columns as real columns, and a WHERE clause on them prunes whole directories before any file opens. The extension also decodes WKB and GeoArrow so ST_ functions work on the geometry column in the same query.

Related

• Trajectory Object Design Patterns
• Exporting Smoothed Trajectories to GeoParquet
• Discrete Global Grid Systems
• Sampling Rate Optimization
• GPS Precision & Error Handling

Back to Spatiotemporal Data Foundations & Structures