Geostatistical Modeling Spatial Statistics with Python

Using GeoParquet for Large Spatial Datasets

TL;DR: Write with gdf.to_parquet("data.parquet", geometry_encoding="WKB", write_covering_bbox=True) and read a subset with gpd.read_parquet("data.parquet", bbox=(minx, miny, maxx, maxy)). For attribute pruning add filters=[("landuse", "==", "forest")]. For datasets that do not fit in memory, write a partitioned PyArrow dataset with partition_cols and let bounding-box plus predicate pushdown skip whole files and row groups before anything is decoded.

Why This Matters

Shapefile and GeoJSON force a full parse of every feature before you can touch a single record, they store geometry row by row, and Shapefile additionally caps files at 2 GB and field names at ten characters. For the multi-million-feature layers common in geostatistical work — parcel databases, sensor archives, gridded prediction points — that model collapses. GeoParquet stores geometry and attributes columnarly, compressed, with per-row-group statistics, so a reader can prune by spatial extent and by attribute value before decoding. This page applies the lazy, read-only-what-you-need philosophy of the Memory-Efficient Processing guide to on-disk storage, and it underpins every ingestion step in the wider Python Workflows for Spatial Modeling & Regression stack.

The payoff is concrete: a spatial query over a national point layer that took minutes to scan as GeoJSON becomes a sub-second read that touches only the row groups intersecting your window. Combined with partitioning, it lets you keep terabyte-scale archives on disk or object storage and pull working slices into RAM on demand.

Environment

bash
pip install \
  geopandas==0.14.4 \
  pyarrow==15.0.2 \
  shapely==2.0.4 \
  pyproj==3.6.1 \
  numpy==1.26.4
python
import numpy as np
import geopandas as gpd
import pyarrow as pa
import pyarrow.parquet as pq
import pyarrow.dataset as ds
from shapely.geometry import box

GeoParquet 1.1 bounding-box support requires geopandas>=0.14 and pyarrow>=13. Confirm versions before relying on bbox= pushdown.

Step-by-Step Implementation

Step 1 — Write a Single GeoParquet File

GeoDataFrame.to_parquet serialises geometry as WKB and records CRS and schema in the file metadata. Sort features spatially first so consecutive rows are geographically close, which makes row-group bounding boxes tight and prunable.

python
gdf = gpd.read_file("sensors.gpkg")                 # any source
gdf = gdf.to_crs(3857)                              # a projected CRS for tidy bboxes

# Spatial sort by Hilbert distance so row groups are compact blocks
gdf = gdf.sort_values(by=gdf.geometry.hilbert_distance()).reset_index(drop=True)

gdf.to_parquet(
    "sensors.parquet",
    geometry_encoding="WKB",
    write_covering_bbox=True,     # per-row-group bbox statistics for pushdown
    row_group_size=100_000,       # features per independently prunable block
    compression="zstd",
)

write_covering_bbox=True adds the covering-bbox column that makes spatial pushdown possible. Without it, read_parquet(bbox=...) must fall back to reading everything.

Why row-group pruning needs a spatial sort Each panel plots the along-x extent of six row groups as horizontal bars against a shared axis, with the query bounding box drawn as a vertical dashed band. On the left the file was written in arbitrary row order, so every bar spans almost the whole extent and crosses the band, and all six groups are decoded. On the right the file was sorted by Hilbert distance, so the six bars are short and non-overlapping and only the third crosses the band, leaving one group decoded and five skipped. Why row-group pruning needs a spatial sort Arbitrary row order every row group's bbox spans the whole extent RG 0 RG 1 RG 2 RG 3 RG 4 RG 5 decode decode decode decode decode decode x min query bbox x max 6 of 6 row groups decoded 600,000 features touched — a full scan Sorted by hilbert_distance() each row group covers one compact patch RG 0 RG 1 RG 2 RG 3 RG 4 RG 5 skip skip decode skip skip skip x min query bbox x max 1 of 6 row groups decoded 100,000 features touched — one row_group_size block

Step 2 — Partition a Larger-than-Memory Dataset

When the source itself does not fit in RAM, write a partitioned dataset: many files under a directory, split by a categorical or spatial-tile column. Readers skip whole partitions whose value cannot match the query.

python
gdf["tile"] = (gdf.geometry.x // 100_000).astype("int32").astype(str) + "_" \
            + (gdf.geometry.y // 100_000).astype("int32").astype(str)

table = pa.Table.from_pandas(gdf.to_wkb())          # geometry as WKB columns
pq.write_to_dataset(
    table,
    root_path="sensors_ds",
    partition_cols=["tile"],
    row_group_size=100_000,
)

This lays out sensors_ds/tile=0_0/…, sensors_ds/tile=1_0/… and so on. A query restricted to one tile opens only that subdirectory.

Step 3 — Read with Bounding-Box Pushdown

Pass a bbox tuple to read_parquet. The reader compares your window against each row group’s covering bbox and decodes only intersecting groups.

python
query_window = (400_000, 5_600_000, 450_000, 5_650_000)   # minx, miny, maxx, maxy
subset = gpd.read_parquet("sensors.parquet", bbox=query_window)
print(f"Read {len(subset):,} of {len(gdf):,} features")

Only the row groups overlapping query_window are touched. On a spatially sorted file this typically reads a few percent of the data for a local query.

Step 4 — Add Predicate Pushdown on Attributes

Attribute filters are evaluated against per-row-group min/max statistics, so groups that cannot satisfy the predicate are skipped without decoding. Combine spatial and attribute pushdown for the tightest read.

python
subset = gpd.read_parquet(
    "sensors.parquet",
    bbox=query_window,
    filters=[("status", "==", "active"), ("pm25", ">", 15.0)],
    columns=["geometry", "station_id", "pm25"],   # column projection
)

columns= is the third lever: reading only the columns you model, alongside bbox and predicate pushdown, minimises bytes decoded on all three axes — extent, rows, and columns.

Step 5 — Query a Partitioned Dataset

For the partitioned layout, open it as a PyArrow dataset and let both partition pruning and row-group statistics apply. Convert the result back to a GeoDataFrame.

python
dataset = ds.dataset("sensors_ds", format="parquet", partitioning="hive")
filtered = dataset.to_table(
    filter=(ds.field("tile") == "4_56") & (ds.field("pm25") > 15.0),
    columns=["geometry", "station_id", "pm25"],
)
subset = gpd.GeoDataFrame.from_arrow(filtered).set_crs(3857)
print(subset.shape)

The tile == "4_56" predicate prunes at the directory level; pm25 > 15.0 prunes at the row-group level within the surviving files.

Interpreting the Output

The headline signal is the ratio of features read to features stored: len(subset) / len(gdf). For a local spatial query on a spatially sorted file it should be a small fraction — if it approaches 1.0, pushdown is not engaging. The two usual causes are a missing covering-bbox column (rewrite with write_covering_bbox=True) or an unsorted file where every row group spans the whole extent, making every group intersect every query. Inspect grouping with pq.ParquetFile("sensors.parquet").metadata, which reports num_row_groups and per-group row counts; a single giant row group cannot be pruned at all.

How a bbox query is tested against row-group statistics On the left, sensors.parquet is drawn as three stacked row groups of 100,000 rows each, showing the x and y ranges of the covering bounding box, with the footer metadata block beneath them. On the right, an arrow from each row group carries the reader's verdict for the query window 400,000 to 450,000 east and 5,600,000 to 5,650,000 north: row group 0 is skipped because its maximum x of 395,000 is below the query minimum, row group 1 is decoded because both axes overlap, and row group 2 is skipped because its minimum y of 5,680,000 lies above the query maximum. Which row groups a bbox query actually decodes sensors.parquet · 3 row groups · 300,000 features query bbox = (400,000, 5,600,000, 450,000, 5,650,000) Row group 0 · 100,000 rows bbox x 320,000 – 395,000 y 5,520,000 – 5,610,000 skip — the x ranges are disjoint group max x = 395,000 falls below the query min x = 400,000 no feature here can lie in the window, so nothing is decoded Row group 1 · 100,000 rows bbox x 395,000 – 448,000 y 5,600,000 – 5,660,000 decode — both axes overlap x and y ranges each straddle the window, so the group may contribute 100,000 rows are decoded, then filtered exactly against the bbox Row group 2 · 100,000 rows bbox x 405,000 – 440,000 y 5,680,000 – 5,740,000 skip — the y ranges are disjoint x overlaps, but group min y = 5,680,000 sits above query max y = 5,650,000 both axes must intersect — one disjoint axis is enough to skip Footer metadata — read first num_row_groups = 3 · rows per group column min/max · covering–bbox statistics 1 of 3 groups decoded — 100,000 of 300,000 features read this ratio is the diagnostic: near 1.0 means pushdown never fired a file written as one giant row group reports num_row_groups = 1

For partitioned datasets, list the touched files during a read (PyArrow exposes fragment metadata) to confirm partition pruning skipped the directories you expect. A query that opens every partition means the partition key is not aligned to your query pattern.

Critical Best Practices

Always Write a Covering Bbox and Sort Spatially

Bbox pushdown depends on two things: the covering-bbox column existing, and row groups being geographically compact. Writing without write_covering_bbox=True, or writing in arbitrary row order, leaves you with pruning that never fires. Sort by hilbert_distance() (or GeoHash, or tile id) before writing so each row group covers a small, non-overlapping patch.

Use a Projected CRS for Storage

Bounding boxes in degrees are awkward near the antimeridian and poles, and mixing geographic coordinates with metric query windows silently returns nothing. Store large layers in a projected CRS and document it; downstream distance and buffer work needs metric units anyway, as covered in Reprojecting CRS for Accurate Distance Calculations.

Right-Size Row Groups and Partitions

Row groups near 50,000–200,000 features balance pruning granularity against metadata overhead. Partitions should map to your query pattern: if you always query one region at a time, partition by spatial tile; if you always query one date, partition by date. Over-partitioning creates thousands of tiny files that hurt object-storage listing performance.

Prefer GeoParquet over Shapefile at Every Boundary

Every place a pipeline reads or writes a Shapefile is a place it truncates field names, loses type fidelity, and forces a full parse. Convert legacy Shapefiles to GeoParquet once at ingestion. This dovetails with the profiling and lazy-read patterns in the parent Memory-Efficient Processing guide.

Combine Column, Row, and Extent Projection

The cheapest read narrows on all three axes at once: columns= for the fields you model, filters= for attribute predicates, and bbox= for extent. Reaching for only one leaves bytes on the table. Profile the decoded size, not the file size, to see the effect.

A GeoParquet read is a rectangle of row groups by columns The file is drawn three times as a grid whose rows are five row groups and whose columns are the six stored fields geom, id, pm25, temp, status and date, so each of the thirty cells is one column chunk. A plain read shades all thirty cells; adding bbox=query_window leaves only the two row groups that pass the extent test, twelve chunks; adding columns=[geometry, station_id, pm25] leaves six. A note beneath explains that filters= cuts the remaining rows inside those surviving groups. A read is a rectangle: row groups × columns read_parquet("sensors.parquet") geom id pm25 temp status date 30 of 30 column chunks decoded every field of every row group is read + bbox=query_window geom id pm25 temp status date 12 of 30 column chunks decoded 3 row groups fail the bbox test, never opened + columns=[geometry, station_id, pm25] geom id pm25 temp status date 6 of 30 column chunks decoded 3 of the 6 stored fields are ever touched RG 0 RG 1 RG 2 RG 3 RG 4 filters=[("status", "==", "active"), ("pm25", ">", 15.0)] cuts the third axis Within the two surviving groups, any whose min/max statistics cannot satisfy the predicate is dropped before decoding.

Troubleshooting

Symptom Likely cause Fix
bbox= read returns the whole file No covering-bbox column written Rewrite with write_covering_bbox=True on geopandas>=0.14
Pushdown reads almost everything File not spatially sorted; row groups overlap globally Sort by hilbert_distance() before to_parquet
Empty result for a valid area Query window CRS differs from stored CRS Reproject the bbox to the file’s CRS before querying
filters= has no effect Predicate column not in row-group statistics (e.g. after casting) Ensure the column is a native Arrow type; avoid object dtype
Slow reads from object storage Thousands of tiny partition files Coarsen partition_cols; increase row_group_size
CRS missing after from_arrow Arrow table dropped GeoParquet metadata Re-attach with .set_crs(...) matching the write CRS

Next Steps

GeoParquet is the durable sink for the tiled outputs produced upstream; see Chunked Raster Processing with Dask-GeoPandas for generating partitioned zonal results, and return to the parent Memory-Efficient Processing guide for the downcasting and sparse-weight techniques that pair with columnar storage.

Related:

← Back to Memory-Efficient Processing

Related

Chunked Raster Processing with Dask-GeoPandas Parallel Spatial Workflows with Dask and Joblib Reducing Memory Bottlenecks in Geospatial Workflows