Spatial Epi Public Health GIS

Querying GeoParquet With DuckDB Spatial

This guide, part of Cloud-Native Geospatial Storage for Surveillance, covers the day-to-day use of the archive: an epidemiologist asking a question in SQL without copying a partition, extracting a file, or waiting for a cluster.

Problem Context & Constraints

An in-process analytical engine reading Parquet directly gives an analyst something a database usually does not: the ability to query a terabyte archive from a laptop session, over the network, with no server to provision and no copy to govern. The constraint is that this only works when the query is written so the engine can eliminate data before reading it.

Three patterns determine whether that happens.

Partition predicates must be on partition columns, literally. A filter on a derived expression — extracting a month from a timestamp column rather than filtering the month partition key — defeats directory elimination, and the engine opens every partition.

Spatial predicates need a bounding-box path. A spatial function applied to the decoded geometry cannot be evaluated until the geometry is read. A predicate expressed against the covering bounding-box columns, or written so the engine can derive one, eliminates row groups before decoding anything.

Joins against a large reference layer should be the small side. Joining a billion-row case archive to a boundary table works well when the boundary table is broadcast; it works badly when the planner chooses otherwise because statistics are missing.

The fourth constraint is governance rather than performance: a query interface over protected data still needs suppression. An analyst aggregating to blocks and exporting the result has produced a small-cell release, and the query layer is the natural place to prevent that.

Two Queries With the Same Answer and Very Different Costs Two query plans side by side over the same archive. The first filters on the partition key directly and uses a bounding-box predicate, so the engine eliminates all but one directory and then all but a few row groups, decoding geometry for a small fraction of rows. The second wraps the date column in a function and applies a spatial predicate to the decoded geometry, so no directory or row group can be eliminated and every row is decoded before being discarded. Both return the same result. The planner can only skip what the predicate lets it skip Predicate on the partition key 48 directories → 1 opened 420 row groups → 6 read geometry decoded for 0.2% of rows 1.4 seconds Predicate on a derived expression 48 directories → 48 opened 420 row groups → 420 read geometry decoded for every row 6 minutes Identical results; the difference is entirely in how the filter was written

Prerequisites

  • python 3.11
  • duckdb 1.0.0 with the spatial extension, pyarrow 16.1.0, pandas 2.2.2
  • A partitioned archive written as in the sibling guides, with covering bounding-box columns
  • Boundary layers stored per vintage
  • The disclosure policy’s minimum cell size, for the suppression wrapper

Step-by-Step Solution

Open the archive as a view, write predicates that eliminate, join boundaries by vintage, and wrap aggregation in a suppression step so a small cell cannot leave the session.

# Query patterns over a partitioned GeoParquet surveillance archive.
# Pinned: duckdb==1.0.0, pyarrow==16.1.0, pandas==2.2.2
import logging
import os

import duckdb
import pandas as pd

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("storage.query")

MIN_CELL = 11        # from the disclosure policy, not chosen here


def connect(archive_root: str, boundary_root: str) -> duckdb.DuckDBPyConnection:
    """A session with the archive and the boundary vintages registered as views."""
    con = duckdb.connect()
    con.execute("INSTALL spatial; LOAD spatial;")
    con.execute("SET preserve_insertion_order=false;")   # lower memory on big scans
    con.execute(f"""
        CREATE VIEW cases AS
        SELECT * FROM read_parquet('{os.path.join(archive_root, "**", "*.parquet")}',
                                   hive_partitioning=true, union_by_name=true)
    """)
    con.execute(f"""
        CREATE VIEW boundaries AS
        SELECT * FROM read_parquet('{os.path.join(boundary_root, "**", "*.parquet")}',
                                   hive_partitioning=true)
    """)
    log.info("session ready: archive=%s boundaries=%s", archive_root, boundary_root)
    return con


def cases_in_bbox(con, event_month: str, minx: float, miny: float,
                  maxx: float, maxy: float) -> pd.DataFrame:
    """Selective read: partition predicate first, then the bbox columns."""
    sql = """
        SELECT geoid, event_date, case_count
        FROM cases
        WHERE event_month = ?                      -- eliminates directories
          AND bbox.xmax >= ? AND bbox.xmin <= ?    -- eliminates row groups
          AND bbox.ymax >= ? AND bbox.ymin <= ?
    """
    df = con.execute(sql, [event_month, minx, maxx, miny, maxy]).df()
    log.info("bbox query month=%s rows=%d", event_month, len(df))
    return df


def counts_by_area(con, event_month: str, vintage: str,
                   min_cell: int = MIN_CELL) -> pd.DataFrame:
    """Aggregate to areas of a given boundary vintage, suppressing small cells.

    Suppression happens INSIDE the query, so a small cell never materialises in
    the analyst's session and cannot be exported by accident.
    """
    sql = """
        WITH b AS (
            SELECT geoid, geom FROM boundaries WHERE vintage = ?
        ),
        joined AS (
            SELECT b.geoid, c.case_count
            FROM cases c
            JOIN b ON ST_Within(c.geom, b.geom)
            WHERE c.event_month = ?
        ),
        agg AS (
            SELECT geoid, SUM(case_count) AS cases, COUNT(*) AS records
            FROM joined GROUP BY geoid
        )
        SELECT geoid,
               CASE WHEN cases < ? THEN NULL ELSE cases END AS cases,
               cases < ? AS suppressed
        FROM agg
        ORDER BY geoid
    """
    df = con.execute(sql, [vintage, event_month, min_cell, min_cell]).df()
    n_sup = int(df["suppressed"].sum())
    log.info("counts_by_area month=%s vintage=%s areas=%d suppressed=%d (min_cell=%d)",
             event_month, vintage, len(df), n_sup, min_cell)
    if n_sup and df["suppressed"].mean() > 0.5:
        log.warning("more than half of areas suppressed — the geography is too fine "
                    "for this month's counts")
    return df


def explain(con, sql: str, params=()) -> str:
    """Read the plan before running an expensive query."""
    plan = con.execute("EXPLAIN " + sql, list(params)).fetchall()
    text = "\n".join(r[1] for r in plan)
    log.info("plan:\n%s", text[:800])
    return text

Putting suppression inside the query rather than after it is the design choice that matters most. An analyst who receives a result set with small cells already present has them in memory, in a notebook, and in whatever they export next; a result set where the suppression happened in the engine never had them. It costs a CASE expression and removes an entire class of accidental disclosure.

The bounding-box predicate is worth writing explicitly even when the engine could infer one, because inference depends on version and on how the geometry predicate is phrased. Writing the box filter alongside the exact spatial predicate gives the planner an eliminable condition and keeps the precise one for correctness.

Where Suppression Happens Decides What Exists Two paths for the same aggregation. In the first, the engine applies the small-cell rule during aggregation, so the analyst's session receives a result in which suppressed areas are null and no small count ever existed outside the protected store. In the second, the raw aggregate including small counts is returned to the session and suppressed afterwards in the analyst's own code, so the small counts existed in memory, in the notebook output and in anything exported before the suppression line ran. A small cell that never leaves the engine cannot be exported by accident Suppress inside the query aggregate and apply the rule session receives nulls nothing small ever existed here Suppress after the result aggregate and return everything small counts in memory and output one forgotten line exports them The cost of the safe path is a CASE expression and a helper function so analysts do not write the aggregation by hand

Validation & Edge Cases

Predicates that defeat elimination. Wrapping a partition column in a cast or a function, comparing it to a non-literal, or filtering on a column derived in a subquery all prevent directory pruning. Read the plan rather than guessing; the row-group count in the plan tells you immediately.

Geometry stored in a different coordinate system than the query box. A bounding box in geographic coordinates against a projected archive silently matches nothing or everything. Assert the system on the view and transform the box, not the archive.

Union by name hiding a schema error. Reading partitions with differing schemas by name is convenient and will silently null a renamed column. Compare the union schema against the expected one and fail on a surprise.

Aggregations that suppress most areas. If half the areas are suppressed, the geography is too fine for the counts, and the analyst should aggregate to a coarser level rather than publishing a mostly empty map. The warning in the code exists to prompt that.

Complementary suppression left undone. Nulling one small cell in a row whose total is published lets the reader recover it by subtraction. The query layer can suppress the primary cells; it cannot decide which additional cells must go without knowing what else is being released. Where a result will be published rather than used internally, route it through the disclosure review process rather than treating the query’s suppression as sufficient.

Results cached by a notebook. An analytical notebook keeps outputs in its saved file, so a result set that was correctly suppressed still ends up stored wherever the notebook is. That is usually inside the protected environment and occasionally is not. Treat notebook files as data artifacts subject to the same location rules as any export.

Memory on a wide scan. A query that selects every column including geometry over many partitions will exhaust memory. Project only the columns needed, which in a columnar format is free.

INFO storage.query bbox query month=2026-08 rows=1204
INFO storage.query counts_by_area month=2026-08 vintage=2020 areas=1842 suppressed=619 (min_cell=11)
WARNING storage.query more than half of areas suppressed — the geography is too fine for this month's counts

The coordinate reference system mismatch in that list produces a result that is empty or total rather than wrong, which makes it easy to recognise once you know the pattern:

A Degree Box Against a Metre Archive An archive stored in a projected coordinate system whose coordinates run into the millions of metres. A bounding-box query written in longitude and latitude uses values between minus one hundred twenty and fifty, a tiny box near the projected origin where no data lies, so the query returns zero rows without error. Transforming the box into the archive's system before querying produces a box covering the intended county and returns the expected rows. Zero rows with no error usually means the box is in the wrong units projected archive, coordinates in metres degree box near origin: 0 rows transformed box 1 204 rows transform the query box into the archive's system, never the archive into the box's

Compliance Notes

Record the minimum cell size the query layer enforces and its policy source, the helper functions analysts are expected to use rather than writing raw aggregations, and the boundary vintage each query joined against.

A query interface over protected data is an access path and belongs in the system’s access control model. Two practices keep it defensible. Log queries with the run or session identity, the partitions touched and the row counts returned — not the results — so unusual access patterns are visible. And provide the suppression-aware helpers prominently enough that the safe path is also the easy one, because an analyst who has to write the aggregation by hand will eventually write it without the rule.