Spatiotemporal Query Indexing

Spatiotemporal query indexing is the set of physical-layout decisions that let a question of the form “every fix inside this polygon between these two timestamps” touch a small fraction of a trajectory archive instead of all of it.

It is a separate discipline from spatial indexing because trajectory archives have an unusual shape: they are shallow in space and extremely deep in time. A city fleet covers the same few hundred square kilometres every single day, so after three years the archive contains a thousand near-identical spatial layers stacked on top of each other. A spatial index navigates that stack perfectly and still returns all thousand layers, because nothing in it knows about time.

Prerequisites

• A settled row shape. Whether the archive stores one row per fix or one row per trip changes every decision below; see spatial storage formats.
• A metric CRS, consistently applied. Index keys computed in one CRS and queried in another silently return nothing; see coordinate reference system mapping.
• A monotone UTC timestamp column. Time partitioning against a local-time column produces partitions that overlap twice a year.
• Python stack. duckdb >= 0.10, pyarrow >= 14, geopandas >= 0.14, shapely >= 2.0, and psycopg/PostGIS 3.4 if you use the database path.

The Shape of the Problem

Mechanism Taxonomy

The last row is under-used and often the highest-value addition to a trajectory archive. A summary table with one row per trip — entity, start and end time, bounding box, fix count — is small enough to keep entirely in memory and answers a large class of queries (“which trips could possibly intersect this polygon in this window?”) without touching the fix-level data at all. The fix-level scan then runs against a candidate list of hundreds rather than a partition of millions.

Deterministic Pipeline Overview

Implementation Walkthrough

The example below builds a partitioned, curve-clustered archive and queries it with DuckDB. The write path is where all the performance is decided; the read path is then ordinary SQL.

import duckdb
import numpy as np
import pandas as pd


def hilbert_key(x: np.ndarray, y: np.ndarray, order: int = 16) -> np.ndarray:
    """
    Map projected metric coordinates to a 1-D Hilbert index.

    Nearby points in the plane receive nearby keys, so sorting on the result
    lays spatially adjacent rows out contiguously on disk.

    Parameters
    ----------
    x, y : np.ndarray
        PROJECTED metric coordinates. Passing degrees produces a valid-looking
        key with the wrong locality.
    order : int
        Bits per axis; 16 gives a 65 536 x 65 536 grid over the extent.

    Raises
    ------
    ValueError
        If the arrays differ in length or are empty.
    """
    if x.shape != y.shape:
        raise ValueError("x and y must have the same shape.")
    if x.size == 0:
        raise ValueError("Empty coordinate arrays.")

    side = 1 << order
    # Normalise to the integer grid over the data extent.
    xs = np.clip(((x - x.min()) / max(np.ptp(x), 1e-9) * (side - 1)).astype(np.int64), 0, side - 1)
    ys = np.clip(((y - y.min()) / max(np.ptp(y), 1e-9) * (side - 1)).astype(np.int64), 0, side - 1)

    rx = np.zeros_like(xs)
    ry = np.zeros_like(ys)
    d = np.zeros_like(xs)
    s = side // 2
    while s > 0:
        rx = ((xs & s) > 0).astype(np.int64)
        ry = ((ys & s) > 0).astype(np.int64)
        d += s * s * ((3 * rx) ^ ry)
        # Rotate the quadrant so the curve stays continuous.
        swap = ry == 0
        flip = swap & (rx == 1)
        xs_f, ys_f = xs.copy(), ys.copy()
        xs = np.where(flip, s - 1 - xs_f, xs_f)
        ys = np.where(flip, s - 1 - ys_f, ys_f)
        xs_s, ys_s = xs.copy(), ys.copy()
        xs = np.where(swap, ys_s, xs_s)
        ys = np.where(swap, xs_s, ys_s)
        s //= 2
    return d


def write_partitioned(df: pd.DataFrame, out_dir: str) -> None:
    """
    Write a fix-level frame as day-partitioned, Hilbert-clustered Parquet.

    The two decisions that matter both happen here: the partition column,
    and the sort order inside each partition.
    """
    required = {"t", "x", "y", "entity_id"}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if df.empty:
        raise ValueError("Input DataFrame is empty.")

    df = df.copy()
    df["day"] = df["t"].dt.strftime("%Y-%m-%d")
    df["h_key"] = hilbert_key(df["x"].to_numpy(), df["y"].to_numpy())
    # Sorting inside the partition is what makes row-group statistics useful.
    df = df.sort_values(["day", "h_key", "t"])

    con = duckdb.connect()
    con.register("fixes", df)
    con.execute(
        f"""
        COPY (SELECT * EXCLUDE (h_key) FROM fixes)
        TO '{out_dir}'
        (FORMAT PARQUET, PARTITION_BY (day),
         ROW_GROUP_SIZE 400000, OVERWRITE_OR_IGNORE 1)
        """
    )


def query_window(out_dir: str, bbox: tuple, t0: str, t1: str) -> pd.DataFrame:
    """Space-and-time range query; both predicates prune, in that order."""
    minx, miny, maxx, maxy = bbox
    con = duckdb.connect()
    return con.execute(
        f"""
        SELECT entity_id, t, x, y
        FROM read_parquet('{out_dir}/**/*.parquet', hive_partitioning = 1)
        WHERE day BETWEEN strftime(TIMESTAMP '{t0}', '%Y-%m-%d')
                      AND strftime(TIMESTAMP '{t1}', '%Y-%m-%d')
          AND t BETWEEN TIMESTAMP '{t0}' AND TIMESTAMP '{t1}'
          AND x BETWEEN {minx} AND {maxx}
          AND y BETWEEN {miny} AND {maxy}
        ORDER BY entity_id, t
        """
    ).df()

The subtlety in query_window is the redundant-looking day predicate. Without it the engine has to open every partition to evaluate t; with it, partition pruning happens on the directory names before a single byte of Parquet is read. Expressing the same constraint twice — once against the partition column and once against the real timestamp — is the standard idiom, and omitting the first half is the most common reason a partitioned archive performs like an unpartitioned one.

Calibration and Parameter Tuning

Partition granularity is the one worth revisiting as the archive grows. Daily partitions are right at city scale; at national scale they produce partitions too large to prune usefully and want a secondary partition on a coarse spatial cell. The signal to watch is the selectivity ratio, not the runtime — hardware improves runtime and hides a layout that is quietly getting worse.

Integration and Compatibility

The index layout is downstream of the write path and upstream of everything else, which makes it awkward to change later. Two decisions in particular are effectively permanent once an archive has years in it: the partition key and the row shape. Both can be changed only by rewriting the whole archive, so it is worth spending an afternoon on the query shapes before the first byte lands.

Two integrations matter in practice. Discrete global grid systems supply an alternative clustering key: an H3 or S2 cell id is a perfectly good sort column, and an S2 id is itself a Hilbert index, which is why range scans on it are spatial queries. And origin-destination flow matrices and kernel density surfaces are both aggregate consumers whose runtime is dominated by how much of the archive they have to read — improving the layout speeds them up more than optimising their own code ever will.

Validation and Testing Patterns

An index is easy to test badly. The usual mistake is to benchmark a warm query on a laptop, see a good number, and ship a layout that degrades quietly for two years. Three checks are worth automating instead.

Selectivity regression. Run a fixed set of representative queries — a district-afternoon, a corridor-week, a single vehicle-month — and record rows scanned and rows returned for each. The ratio is hardware-independent and monotonically reveals a layout drifting away from the query pattern. Fail the build if it doubles.

Cold-cache timing. Measure with the page cache dropped, because a warm benchmark measures memory rather than layout. The gap between warm and cold timing is itself informative: a large gap means the query is I/O-bound and the layout is what matters, a small gap means it is CPU-bound and the predicate is.

Partition-count sanity. Assert that a bounded-time query opens the expected number of partitions. This single check catches the most common silent failure — a predicate written only against the timestamp column, which the planner cannot use for directory pruning — and it catches it at the moment the query is written rather than a year later when the archive is large enough for anyone to notice.

Keep a fixed test archive of a few million rows in the repository. A layout change that improves one query and destroys another is common, and only a fixed corpus makes that visible before deployment.

In This Section

• Indexing trajectories in PostGIS for space-time range queries — GiST, BRIN and partitioning in a database that is also taking writes.
• Using STR-packed R-trees for fast trajectory lookups — bulk-loading an in-process index over a static candidate set.

FAQ

Why is a spatial index alone slow on trajectory data?

Because trajectory archives are deep in time and shallow in space. A city’s fleet revisits the same few hundred square kilometres every day for years, so a purely spatial index returns every year of history for the queried area and filters by time afterwards. The index did its job and still scanned a thousand times more rows than the query needed.

Should time be a partition, an index column, or part of a combined key?

Partition first. Physical partitioning eliminates whole files before any index is consulted, which is a far larger saving than index selectivity. Use a combined key — a curve cell plus a time bucket — for clustering inside the partition, and an index only for the residual. Time as an ordinary indexed column with no partitioning is the worst-performing configuration at scale.

R-tree or space-filling curve?

They solve different halves. An R-tree is a secondary structure that answers arbitrary-shaped range and nearest-neighbour queries but does not change disk layout. A curve turns two or three dimensions into one sortable key, so spatially adjacent rows land adjacently on disk and a range query becomes a few sequential reads. Large systems use both.

Do I need PostGIS, or is DuckDB over Parquet enough?

It depends on writes and concurrency, not on query power. DuckDB over partitioned GeoParquet is excellent for analytical reads over an append-only archive and needs no server, but it is not built for many concurrent writers or row-level updates. PostGIS earns its place when the data is being updated, when many clients query concurrently, or when transactions matter.

What selectivity should I expect from a well-indexed archive?

For one district and one afternoon, under 20 rows scanned per row returned. Ratios in the hundreds mean the clustering does not match the query shape; ratios in the thousands mean time is not pruning at all. Track the ratio rather than the runtime — hardware improves runtime and hides a degrading layout.

Related

• Discrete Global Grid Systems — cell ids as clustering keys, and why an S2 id is already a curve index.
• Spatial Storage Formats — the row-shape and partitioning decisions this page builds on.
• Optimizing Spatial Joins for Trajectory-to-Zone Matching — the in-memory counterpart of the pruning funnel here.
• Origin-Destination Flow Matrices — an aggregate consumer whose cost is dominated by archive layout.
• Exporting Smoothed Trajectories to GeoParquet — the write path that decides whether any of this works.

Back to Spatial Indexing & Density Aggregation