PostGIS vs GeoPandas for Validation at Scale

You have a validation workload that runs fine on a sample but you need to know how it behaves when the parcel layer grows from a hundred thousand features to ten million — and whether the GeoPandas job that finishes in seconds today will run out of memory next quarter. This page gives concrete throughput and memory figures at three scales, explains the Generalized Search Tree (GiST) versus Sort-Tile-Recursive tree (STRtree) index behaviour that drives them, and provides a benchmark harness you can run on your own data. It refines the engine-selection method in Comparing Spatial Validation Engines with the specific numbers that decide the PostGIS-versus-GeoPandas call.

Prerequisites

• PostGIS 3.4+ on PostgreSQL 15+, with shared_buffers and work_mem sized for your host. An under-provisioned work_mem forces spatial joins to spill to disk and distorts every PostGIS timing downward.
• GeoPandas 0.14+ with Shapely 2.0+ and pyarrow 14+. Shapely 2.0 vectorised functions and the STRtree are prerequisites; without them the Python side is single-geometry-at-a-time and the comparison is unfair.
• A single projected coordinate reference system (CRS). Fix on one projected CRS (for example EPSG:3857) so area and distance predicates are valid and identical in both engines. Never benchmark one engine in geographic degrees and the other in metres.
• Enough disk for 10M-feature fixtures. A ten-million-polygon GeoPackage plus its PostGIS copy and GiST index can exceed tens of gigabytes. Provision before you generate.
• Gotcha: psql \timing measures round-trip time including network latency, while EXPLAIN ANALYZE measures server-side execution. Use the latter for engine comparison so client latency does not contaminate the number.

Step-by-Step Procedure

Step 1 — Generate reproducible scaled fixtures

import numpy as np
import geopandas as gpd
from shapely.geometry import Polygon

def make_fixture(n: int, seed: int = 42, epsg: int = 3857) -> gpd.GeoDataFrame:
    """Deterministic grid of small square polygons with a zoning attribute."""
    rng = np.random.default_rng(seed)
    side = int(np.ceil(np.sqrt(n)))
    xs, ys = np.meshgrid(np.arange(side), np.arange(side))
    xs, ys = xs.ravel()[:n] * 10.0, ys.ravel()[:n] * 10.0
    geoms = [Polygon([(x, y), (x + 8, y), (x + 8, y + 8), (x, y + 8)]) for x, y in zip(xs, ys)]
    return gpd.GeoDataFrame(
        {"id": np.arange(n), "zoning_code": rng.choice(["R1", "C2", None], size=n)},
        geometry=geoms, crs=f"EPSG:{epsg}",
    )

for n in (100_000, 1_000_000):
    make_fixture(n).to_parquet(f"fixture_{n}.parquet")

Verification: reload each Parquet file and assert len(gdf) == n and gdf.crs.to_epsg() == 3857. A silent CRS drop here would invalidate every downstream distance and area comparison. Generate the 10M fixture on a host with adequate disk; the same function scales but the write is large.

Step 2 — Measure the GeoPandas validation pass

import time, tracemalloc, shapely, pandas as pd

def geopandas_pass(gdf: gpd.GeoDataFrame) -> dict:
    tracemalloc.start()
    t0 = time.perf_counter()

    valid = shapely.is_valid(gdf.geometry.values)          # per-feature
    tree = shapely.STRtree(gdf.geometry.values)            # index BUILD (timed)
    t_build = time.perf_counter()
    pairs = tree.query(gdf.geometry.values, predicate="intersects")  # set-level
    t_query = time.perf_counter()

    _, peak = tracemalloc.get_traced_memory()
    tracemalloc.stop()
    return {
        "valid_pass_s": round(t_build - t0, 3),
        "strtree_build_s": round(t_query - t_build, 3),
        "peak_mb": round(peak / 1e6, 1),
        "invalid": int((~valid).sum()),
        "intersect_pairs": int(pairs.shape[1]),
    }

print(geopandas_pass(gpd.read_parquet("fixture_100000.parquet")))

Verification: run twice and confirm peak_mb is stable while strtree_build_s may fall on the second run as allocator caches warm. Peak memory should scale roughly linearly with feature count — a super-linear jump signals coordinate-array copies you should eliminate.

Step 3 — Measure the equivalent PostGIS pass

-- Ensure the index exists and stats are fresh before timing
CREATE INDEX IF NOT EXISTS fx_gix ON fixture USING GIST (geometry);
ANALYZE fixture;

-- Per-feature validity + set-level self-intersection, server-side timing
EXPLAIN (ANALYZE, BUFFERS)
SELECT a.id
FROM fixture a
JOIN fixture b
  ON a.id < b.id
 AND ST_Intersects(a.geometry, b.geometry);

Read Execution Time from the plan, and confirm the join uses an Index Scan on fx_gix rather than a Seq Scan. Compare shared hit (warm buffer cache) against shared read (cold disk reads) to label the run.

Verification: the first run after ANALYZE is typically cold with high shared read; a second immediate run should show mostly shared hit and a lower execution time. Record both — the warm number is your steady-state throughput, the cold number your worst case.

Step 4 — Locate the crossover

import pandas as pd

# Populate with your measured numbers per scale
rows = [
    {"features": 100_000, "gp_peak_mb": 190,  "gp_time_s": 0.4,  "pg_time_s": 0.6},
    {"features": 1_000_000, "gp_peak_mb": 2100, "gp_time_s": 5.1, "pg_time_s": 3.9},
    {"features": 10_000_000, "gp_peak_mb": None, "gp_time_s": None, "pg_time_s": 41.0},
]
df = pd.DataFrame(rows)
df["gp_viable"] = df["gp_peak_mb"].notna()
print(df)
# Crossover = smallest feature count where PostGIS time < GeoPandas time
# or GeoPandas exhausts memory (gp_viable == False)

Verification: the gp_viable == False row marks the hard ceiling; the first row where pg_time_s < gp_time_s marks the soft crossover. If both land at the same scale, your decision is clean; if they diverge, the band between them is where a hybrid split is worth considering.

Interpreting Results

The numbers in Step 4 are illustrative, but the shape is consistent across real workloads. At 100k features GeoPandas usually wins on wall-clock time because the whole dataset fits in cache and the STRtree build is cheap, while PostGIS pays a small round-trip and planning overhead. At 1M features the two converge: GeoPandas memory climbs into the gigabytes and the STRtree build becomes a visible fraction of total time, while PostGIS holds steady because it never materialises the table. At 10M features GeoPandas frequently cannot complete a naive in-memory self-join at all, whereas PostGIS finishes because it streams pages through the GiST index.

The decisive variable is memory, not raw speed. GeoPandas peak memory grows linearly with both feature count and vertices per geometry; PostGIS memory is bounded by work_mem and shared_buffers regardless of table size. That is why the crossover is better described as a ceiling than as a speed line — GeoPandas does not merely get slower, it stops fitting. This is the same memory-versus-SQL trade-off developed conceptually in Comparing Spatial Validation Engines, made concrete here with figures.

Index behaviour explains the rest. The GiST index is built once and persists on disk, so a repeated set-level check pays the build cost zero additional times; the STRtree is rebuilt on every process start, so each GeoPandas run re-pays it. For one-shot validation the STRtree cost is part of the total and often acceptable; for a scheduled job that runs hourly, re-building the tree every run is pure waste that the database avoids.

There is also a serialisation cost that the raw predicate timings hide. When PostGIS runs a validation query it can return a compact result — a count, or a list of failing feature identifiers — without ever shipping geometries to the client. When GeoPandas validates, the geometries are already in the Python process, so there is no wire cost, but getting them there in the first place required reading and parsing the whole file. The honest comparison therefore depends on where the data starts: if it already lives in the database, the PostGIS pass avoids a multi-gigabyte read that the GeoPandas pass must pay; if it arrives as files that never touch a database, the GeoPandas pass avoids a bulk load and index build that the PostGIS pass must pay. At 10M features these ingestion costs frequently dwarf the predicate evaluation itself, which is why a benchmark that times only the check and ignores the round-trip can point you at the wrong engine.

Tuning shifts the crossover but does not remove it. Raising work_mem lets PostGIS keep larger hash and sort structures in memory and can turn a disk-spilling join into an in-memory one, sharpening its advantage at the 1M and 10M scales. On the Python side, switching from row-wise object access to Shapely 2.0 array functions and reading via pyarrow-backed Parquet reduces both time and peak memory, pushing the GeoPandas ceiling higher. Neither change alters the fundamental shape: linear memory growth for GeoPandas against bounded memory for PostGIS. Treat the reported crossover as specific to a host and a configuration, and re-measure whenever either changes materially.

Gotchas & Edge Cases

Vertex count dominates memory more than feature count. A million simple squares is light; a million multipolygon coastlines with thousands of vertices each can exhaust RAM at a tenth of that count. Always report average vertices per geometry alongside feature count, or your crossover figure will not transfer to a different dataset.

A missing ANALYZE makes PostGIS look terrible. Immediately after a bulk load the planner has no statistics and may choose a sequential scan over the GiST index, inflating times by an order of magnitude. Run ANALYZE before every PostGIS benchmark and confirm the index scan in EXPLAIN output.

GeoPandas .sindex is not the same as building an STRtree yourself. The GeoDataFrame.sindex accessor lazily builds and caches an STRtree on first access. If your benchmark touches .sindex before the timed block, you have already paid the build cost and will under-report it. Decide deliberately whether index build is inside or outside your measured window and be consistent across engines.

Self-joins explode without the a.id < b.id guard. Omitting the ordering predicate doubles every intersection pair and includes each feature intersecting itself, inflating both result counts and time. This is a frequent source of “PostGIS is mysteriously slow” reports that are actually a broken query.

The 10M GeoPandas run may swap rather than error. On a host with large swap, GeoPandas can thrash for many minutes instead of raising MemoryError, producing a misleading “it completed” result at absurd wall-clock cost. Cap the run with a timeout and treat swapping as a failure, not a slow success.

When to Escalate

If your workload lives past the GeoPandas ceiling but you want to keep the Python API and rule structure, partition the job with Scaling GeoPandas Validation with Dask; Dask-GeoPandas runs per-feature checks across partitions that each fit in memory. If the deciding factor is set-level predicate throughput rather than Python-native logic, the micro-level evidence for choosing the database is in Benchmarking Shapely vs PostGIS Predicate Checks, which isolates individual predicates from the full-pipeline noise measured here. When neither engine wins cleanly across your scale range, return to the decision matrix in Comparing Spatial Validation Engines and design a hybrid split at the crossover point.

Frequently Asked Questions

Related

• Comparing Spatial Validation Engines — the decision matrix and hybrid-split method this page supplies numbers for
• Benchmarking Shapely vs PostGIS Predicate Checks — predicate-level micro-benchmarks behind these full-pipeline figures
• Scaling GeoPandas Validation with Dask — partitioned Python execution past the memory ceiling

Back to Comparing Spatial Validation Engines