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.

Controlling the Variables

A benchmark is a controlled experiment, and most spatial benchmarks are uncontrolled ones with a number at the end. Five variables account for nearly all the variance between published results.

Five variables that decide a spatial benchmarkCache state, index presence, scope, data shape and repetition account for most of the variance between published comparisons.Cache state, stated and equalthe largest single factor, in both directionsIndexes present on both sidesverify with an explain plan, not by assumptionSame scope, end to endincluding read and writeReal data shapevertex distribution and selectivity, not synthetic uniformityRepeated runs, distribution reportedmedian and spread, first run discarded
Control these five and the number means something. Leave any uncontrolled and it is a number with an experiment-shaped hole behind it.

Cache state is the largest. A PostGIS query against a warm buffer cache can be an order of magnitude faster than the same query cold, and a GeoPandas run against a file in the OS page cache skips the read entirely. State both explicitly and measure both conditions, because production sees both.

Index presence comes second, and it is embarrassing how often one side has an index and the other does not. Verify with an explain plan on the database side and by checking that sindex was actually consulted on the Python side.

Scope decides what is being compared. If one timing includes reading and writing and the other does not, the comparison is meaningless. Time the whole operation the pipeline will actually perform.

Data shape matters more than data size. Vertex-count distribution, clustering of the overlay layer, and the selectivity of the join all move the result substantially, and they differ between synthetic test data and real parcel fabric.

Repetition turns a number into a measurement. Run each condition several times, discard the first, and report the median and the spread rather than a single figure.

A crossover curve says more than a single timingRelative wall time of the two engines across dataset sizes, showing where the in-database approach overtakes the in-process one on one particular fabric.10k parcels2.1× in-process time100k parcels1.3× in-process time1M parcels0.7× in-process time5M parcels0.3× in-process timeAbove 1.0 the in-process path is faster; below it the database is. Re-measure after any version change on either side.
The crossover is a property of your data and hardware. Measuring the curve costs a little more than measuring one point and answers far more.

The crossover point between the two engines is the number most people want from a benchmark, and it is genuinely useful — provided it is understood as a property of your data and hardware rather than of the libraries. Measuring it means running the same comparison at several dataset sizes rather than one, which costs a little more time and produces a curve instead of a single verdict.

Reading the Result Without Over-Concluding

The output of a fair benchmark is narrower than it looks, and treating it as a general verdict is where benchmarks do damage.

What a fair benchmark does and does not tell youThe scope of a benchmark result: valid for this data, hardware, predicate and cache condition, and not generalisable beyond them.It tells youIt does not tell youEngine choiceWhich was faster on this workloadWhich is faster in generalHardwareBehaviour on the machine you measuredHow it scales on different storageOperationsThroughputDeployment, concurrency or reproducibility costDurabilityThe state of things todayWhether it still holds after a version bump
The most useful artefact is the re-runnable script, not the number it printed.

What a benchmark tells you is how these two implementations performed on this data, on this hardware, under these cache conditions, for this predicate. That is genuinely useful for the decision in front of you. It does not tell you which engine is faster in general, and it will not transfer to a different county, a different predicate or a different machine without re-measuring.

Two further cautions. A benchmark measures throughput and rarely measures the things that dominate operations — deployment complexity, concurrency behaviour, reproducibility — which are covered in the parent comparison and often outweigh a factor of two in speed. And a benchmark run once becomes stale: library versions change, data grows, and the conclusion drawn from it quietly stops being true. If the decision matters enough to benchmark, it matters enough to re-run when either side changes materially.

The most useful artefact from a benchmarking exercise is usually not the timing at all: it is the script that produced it. A committed, re-runnable comparison lets the next person repeat the measurement in minutes instead of re-arguing the question from memory.

Part of: PostGIS vs GeoPandas for setback batch processing

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.