Benchmarking Spatial Joins: PostGIS vs GeoPandas

Opinions about which engine runs spatial joins faster multiply endlessly, but a ten-minute micro-benchmark on your own data settles the argument for good. This guide builds a controlled experiment that runs the same proximity and intersection joins over identical parcels and streets, once through PostGIS with ST_DWithin and ST_Intersects, once through GeoPandas with gpd.sjoin, then measures wall-clock time and peak memory so you can interpret the numbers honestly. It expands the decision framework in the parent PostGIS vs GeoPandas comparison with reproducible evidence.

Prerequisites

Step-by-step

Step 1: Standardize the inputs in a metric CRS

A fair benchmark demands that both engines see byte-identical geometry in the same projection. Reproject once in GeoPandas, write GeoParquet for the Python side, and load the same frame into PostGIS so neither engine gets a head start.

import geopandas as gpd

CRS = "EPSG:32617"  # UTM 17N, metric, so ST_DWithin meters are real meters
parcels = gpd.read_file("parcels.gpkg").to_crs(CRS)
streets = gpd.read_file("streets.gpkg").to_crs(CRS)

parcels.to_parquet("bench_parcels.parquet")   # GeoPandas side reads this
streets.to_parquet("bench_streets.parquet")
print(len(parcels), "parcels;", len(streets), "streets")

Step 2: Load the same features into PostGIS and index them

Push the identical frames into PostGIS via GeoPandas’ to_postgis, then build GiST indexes. Without the indexes, ST_DWithin and ST_Intersects fall back to sequential scans and the comparison is meaningless.

from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://user:pass@localhost:5432/bench")
parcels.to_postgis("parcels", engine, if_exists="replace", index=False)
streets.to_postgis("streets", engine, if_exists="replace", index=False)

with engine.begin() as con:
    con.execute(text("CREATE INDEX ON parcels USING GIST (geometry)"))
    con.execute(text("CREATE INDEX ON streets USING GIST (geometry)"))
    con.execute(text("ANALYZE parcels"))   # let the planner see the index
    con.execute(text("ANALYZE streets"))

Step 3: Time the PostGIS proximity and intersection joins

Run a proximity join (ST_DWithin, 10 m) and a strict ST_Intersects join server-side. Wrap each in a timer and read EXPLAIN ANALYZE once to confirm the planner used the index rather than a sequential scan.

import time

DWITHIN = text("""
    SELECT p.ctid FROM parcels p JOIN streets s
    ON ST_DWithin(p.geometry, s.geometry, 10.0)
""")
INTERSECTS = text("""
    SELECT p.ctid FROM parcels p JOIN streets s
    ON ST_Intersects(p.geometry, s.geometry)
""")

def time_sql(stmt):
    with engine.connect() as con:
        t0 = time.perf_counter()
        n = con.execute(stmt).rowcount
        return time.perf_counter() - t0, n

pg_dwithin_s, pg_dwithin_n = time_sql(DWITHIN)
pg_intersects_s, pg_intersects_n = time_sql(INTERSECTS)

Step 4: Time the GeoPandas joins and sample memory

Mirror both joins in GeoPandas. Proximity has no direct sjoin predicate, so buffer the parcels by 10 m in the metric CRS first, then intersect. Track peak memory with tracemalloc so the comparison covers cost, not just speed.

import time, tracemalloc

def time_py(fn):
    tracemalloc.start()
    t0 = time.perf_counter()
    result = fn()
    elapsed = time.perf_counter() - t0
    peak_mb = tracemalloc.get_traced_memory()[1] / 1e6
    tracemalloc.stop()
    return elapsed, len(result), peak_mb

def gpd_dwithin():
    probe = parcels.copy()
    probe["geometry"] = probe.buffer(10.0)          # 10 m proximity zone
    return gpd.sjoin(probe, streets, predicate="intersects")

def gpd_intersects():
    return gpd.sjoin(parcels, streets, predicate="intersects")

gp_dwithin_s, gp_dwithin_n, gp_dwithin_mb = time_py(gpd_dwithin)
gp_intersects_s, gp_intersects_n, gp_intersects_mb = time_py(gpd_intersects)

Verification

The benchmark is only trustworthy if both engines return the same matches. Compare row counts per join; small differences usually mean a boundary-touch tie or a missing index, not random noise.

print(f"DWithin   PostGIS={pg_dwithin_n:>7}  GeoPandas={gp_dwithin_n:>7}")
print(f"Intersect PostGIS={pg_intersects_n:>7}  GeoPandas={gp_intersects_n:>7}")
assert abs(pg_intersects_n - gp_intersects_n) <= 0.01 * pg_intersects_n
print(f"Times (s): pg_dw={pg_dwithin_s:.3f} gp_dw={gp_dwithin_s:.3f} "
      f"gp_peak={gp_dwithin_mb:.0f}MB")

Interpret the results with the workload in mind. PostGIS typically wins as parcel counts climb, because the persistent GiST index amortizes across runs and execution stays server-side. GeoPandas often wins at small-to-medium scale where its in-memory R-tree and vectorized Shapely 2.0 operations avoid client-server round trips. Whichever leads on your data, feed that conclusion back into the parent PostGIS vs GeoPandas comparison, and if the winning path is GeoPandas at scale, shard it with parallelizing parcel validation with Dask-GeoPandas.

Common Pitfalls

  • Benchmarking a cold cache. The first query pays for disk reads and index warm-up. Run each join two or three times and report the median, or you will credit PostGIS with a penalty that vanishes on the second call.
  • Forgetting the index, then blaming the engine. A missing GiST index turns ST_Intersects into a full cross-product scan. Confirm the plan with EXPLAIN ANALYZE before trusting any timing.
  • Comparing a geographic CRS against a projected one. If one engine runs in EPSG:4326 and the other in a metric CRS, the joins are not equivalent and the proximity distances differ. Pin both to the same projected CRS.

Frequently Asked Questions

Why does the proximity join need a buffer in GeoPandas but not in PostGIS?

PostGIS exposes a dedicated distance predicate, so ST_DWithin tests proximity directly against the index. GeoPandas’ sjoin only supports topological predicates like intersects and within, so you emulate a distance test by buffering one layer by the threshold and intersecting. Buffer in a metric CRS so the radius is expressed in real meters.

How many features do I need before PostGIS clearly wins?

There is no universal crossover; it depends on geometry complexity, hardware, and index freshness. On typical parcel and street data the balance often tips toward PostGIS somewhere in the low millions of features, but the point of this micro-benchmark is to find your own crossover rather than trust a rule of thumb.

Is tracemalloc an accurate way to compare memory?

It captures Python-level allocations well, which is what dominates the GeoPandas path, but it does not see memory held inside the PostGIS server process. Treat it as a fair gauge of client-side footprint for the Python engine and use database server metrics separately if you need the PostGIS side’s memory.

Does row-count parity guarantee the joins are equivalent?

It is a strong signal but not a proof, since two different match sets could share a count by coincidence. For high-stakes benchmarks, compare the actual matched ID pairs, not just totals, and investigate any boundary-touch cases where the two engines legitimately disagree.

Conclusion

A reproducible micro-benchmark replaces engine folklore with numbers from your own parcels. By standardizing the CRS, indexing both stores, timing identical proximity and intersection joins, and checking match parity, you learn exactly where the PostGIS-versus-GeoPandas crossover falls for your workload and can size the rest of the setback batch pipeline with confidence.