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.

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.

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.

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 Reducing Memory Bottlenecks in Geospatial Workflows