Bbox Spatial Filtering of GeoParquet Inputs

Attach a bounding-box predicate to the read and the Parquet reader will skip every row group whose spatial extent does not intersect your region — reading only the features inside the area you are rebuilding, instead of scanning the whole continental file.

When to use this

Bbox filtering is the mechanism behind regional and incremental rebuilds. You maintain one large national or continental GeoParquet file, but a given tile job only needs Bavaria, or the Bay Area, or the single metro whose data just changed. Rather than convert the entire dataset every time, you read the window that matters and tile only that.

Use it when:

• The source is a large single file (or partitioned dataset) covering far more area than the current build.
• You rebuild sub-regions on their own cadence — a nightly refresh of one province against a whole-country archive.
• You are stitching regional .pmtiles outputs that will later be merged, and each shard reads its own bbox window.

The prerequisite for it to actually save I/O rather than just filter after reading is spatial locality on disk: the row groups must be spatially clustered, so that any one region touches only a few of them. Files written with the H3/S2 spatial sort from the ingestion pipeline satisfy this; a randomly ordered file does not.

Specification detail

GeoParquet 1.1 defines a covering — a companion column holding a per-row bounding box as a struct {xmin, ymin, xmax, ymax}, all float/double. The geo metadata records where it lives:

{
  "primary_column": "geometry",
  "columns": {
    "geometry": {
      "encoding": "WKB",
      "covering": {
        "bbox": {
          "xmin": ["bbox", "xmin"], "ymin": ["bbox", "ymin"],
          "xmax": ["bbox", "xmax"], "ymax": ["bbox", "ymax"]
        }
      }
    }
  }
}

Because the covering is an ordinary numeric column, Parquet writes min/max statistics per row group for each of its four members. A reader evaluating bbox.xmax >= region_xmin AND bbox.xmin <= region_xmax AND ... consults those statistics in the footer and skips any row group whose range cannot contain a match — no geometry is decoded for skipped groups. This is predicate pushdown, and it is what makes the covering column worth its modest storage cost.

Two overlap tests matter:

The efficient pattern is bbox-overlap first (skips row groups), then optionally ST_Intersects as an exact refinement on the survivors. The four inequalities that express rectangle overlap are: xmax >= region.xmin, xmin <= region.xmax, ymax >= region.ymin, ymin <= region.ymax — all four must hold.

Two subtleties are worth internalising. The row-group predicate is conservative: it returns a superset, because a row group whose bbox statistics overlap the region may still contain individual features that fall entirely outside it. That is fine for tiling — Tippecanoe clips to tile boundaries anyway — but if you need an exact clip, add ST_Intersects as a second-stage refinement. And for a partitioned GeoParquet dataset laid out with Hive-style directories (for example .../region=us-west/part-0.parquet), the partition key prunes whole files before any row group is even opened, so a directory layout aligned to your rebuild regions is the coarsest and cheapest form of spatial filtering.

Production command

The DuckDB query below rebuilds one region from a continental file. The WHERE clause on the covering columns is what the reader pushes down to row-group statistics; the result streams straight to NDJSON on stdout, ready to pipe into Tippecanoe:

-- region_bay_area.sql  — window: lon [-122.6,-121.7], lat [37.2,38.1]
INSTALL spatial; LOAD spatial;
COPY (
  SELECT id, class, name, geometry
  FROM read_parquet('data/buildings_north_america.geoparquet')
  WHERE bbox.xmax >= -122.60 AND bbox.xmin <= -121.70
    AND bbox.ymax >=   37.20 AND bbox.ymin <=   38.10
) TO '/dev/stdout' (FORMAT JSON, ARRAY false);

Wired into a shell pipeline, the region never touches disk as GeoJSON:

duckdb < region_bay_area.sql \
  | tippecanoe \
      --output=dist/bay_area.pmtiles \
      --layer=buildings \
      --minimum-zoom=8 \
      --maximum-zoom=15 \
      --drop-densest-as-needed \
      --read-parallel \
      --force \
      -

If you need an exact regional clip rather than the conservative superset, keep the cheap bbox predicate to prune row groups and add ST_Intersects as a second stage that runs only on the survivors:

WHERE bbox.xmax >= -122.60 AND bbox.xmin <= -121.70
  AND bbox.ymax >=   37.20 AND bbox.ymin <=   38.10
  AND ST_Intersects(geometry,
        ST_MakeEnvelope(-122.60, 37.20, -121.70, 38.10))

The bbox inequalities do the row-group skipping; ST_Intersects decodes geometry only for the rows those groups contain, so the exact test never touches the pruned data. The pyarrow dataset API expresses the coarse predicate declaratively. A pyarrow.dataset filter over the covering fields is evaluated against row-group statistics before any batch is materialised:

import pyarrow.dataset as ds
import pyarrow.compute as pc

dataset = ds.dataset("data/buildings_north_america.geoparquet", format="parquet")
xmin, ymin, xmax, ymax = -122.60, 37.20, -121.70, 38.10

bbox_overlap = (
    (pc.field("bbox", "xmax") >= xmin) & (pc.field("bbox", "xmin") <= xmax) &
    (pc.field("bbox", "ymax") >= ymin) & (pc.field("bbox", "ymin") <= ymax)
)

for batch in dataset.to_batches(
    filter=bbox_overlap,
    columns=["id", "class", "name", "geometry"],   # projection stacks here
    batch_size=200_000,
):
    emit_features(batch)   # decode WKB, write NDJSON to Tippecanoe stdin

Interaction effects

With column projection. The columns=[...] argument in the pyarrow example and the narrow SELECT list in the DuckDB query are column projection pushdown, applied in the same call. Row filtering and column projection compose cleanly: the reader first decides which row groups survive the bbox predicate, then reads only the projected columns for those groups. Together they can turn a whole-continent read into a few megabytes.

With NDJSON streaming. The TO '/dev/stdout' export is the front half of streaming GeoParquet to Tippecanoe’s stdin as NDJSON. Filtering upstream means fewer features cross the pipe, so the NDJSON serializer and the encoder both do proportionally less work — the region is the only thing that ever gets encoded.

With the spatial sort and slippy-map grid. Bbox pushdown only skips row groups when they are spatially clustered, which is exactly what the H3/S2 sort in the ingestion pipeline produces. When you choose the region rectangle, remember it is expressed in the file’s CRS and is conceptually the same lon/lat window that maps to a block of slippy-map tile coordinates; sizing the bbox to your target tile pyramid keeps the read and the tile output aligned.

Performance impact

I/O saved is roughly proportional to the area ratio between your region and the file’s total extent — but only if row groups are spatially sorted. Pulling one metro (~0.5% of a continental extent) from a spatially sorted file commonly reads 95%+ fewer bytes, because only the handful of row groups covering that metro are touched. The same query against an unsorted file where features are shuffled reads nearly everything, because almost every row group’s bbox statistics span the whole continent and none can be skipped.

Wall-clock for a metro rebuild off a sorted continental file typically drops from minutes to a few seconds of read, with the remaining time dominated by Tippecanoe encoding the (now small) feature set.

Common mistakes

Unsorted row groups, so every group intersects the region. Symptom: the query returns the right features but reads nearly the whole file — no speedup. Diagnose by inspecting per-row-group bbox ranges; if every group’s xmin/xmax spans most of the continent, the data is not spatially clustered. Fix: rewrite the file sorted by an H3 or S2 cell as in the ingestion pipeline, so each row group is geographically compact and its statistics become selective.

A lon/lat box against projected data. If the file is stored in EPSG:3857 (Web Mercator metres) but you filter with xmin=-122.6, the predicate matches almost nothing — those degree values are a sub-metre sliver near the origin in metre units. Symptom: an empty or near-empty result for a region you know has data. Fix: express the rectangle in the file’s actual CRS (check crs in the geo metadata) — reproject your degree bounds to 3857, or filter the file back in EPSG:4326 first.

No bbox covering column, so DuckDB does a full scan. Older GeoParquet (pre-1.1) files have no covering column; WHERE bbox.xmin ... then raises Referenced column "bbox" not found, and falling back to ST_Intersects(geometry, ...) still decodes every geometry because there are no row-group statistics to skip on. Fix: add a covering column when rewriting the file, or precompute per-row bbox columns with ST_XMin(geometry), ST_YMin(geometry), etc., and sort by them so future reads can prune.

Parent: GeoParquet Input Processing

Related

• Column Projection Pushdown for GeoParquet — the column-wise counterpart that stacks with bbox filtering in the same read call.
• Streaming GeoParquet to Tippecanoe stdin as NDJSON — pipe the filtered region straight into the encoder with no intermediate file.
• Tile Coordinate Systems and the Slippy Map Grid — how a lon/lat window maps to the tile pyramid you are rebuilding.