Choosing Between In-Database and In-Process Spatial Joins
The benchmark answers which engine is faster; it does not answer which one to use. A compliance pipeline has requirements a throughput number cannot express — the result has to be reproducible from a recorded state, the logic has to be testable in isolation, and the run has to be defensible a year later. This guide turns the choice into a decision you can make for a specific workload by working through the four properties that actually decide it, and shows the hybrid arrangement that most mature pipelines settle on. It is the decision behind PostGIS vs GeoPandas for setback batch processing.
Prerequisites
Step-by-step
Step 1: Characterise the join before comparing engines
Most of the decision is determined by the shape of the join rather than by the engine. Three numbers describe it: the sizes of both inputs and the selectivity — how many rows survive.
def join_profile(left_n, right_n, expected_out, geometry_complexity="simple"):
"""The three numbers that decide more than the engine choice does."""
return {
"left_rows": left_n,
"right_rows": right_n,
"selectivity": expected_out / max(left_n, 1),
"complexity": geometry_complexity, # vertex counts drive the exact-test cost
"shape": ("small-small" if max(left_n, right_n) < 50_000 else
"large-small" if min(left_n, right_n) < 50_000 else "large-large"),
}A large-small join — a county’s parcels against a few dozen districts — is fast everywhere, and the choice should be made on other grounds entirely. A large-large join with low selectivity is where the engines genuinely diverge, and it is also where the index quality matters more than the engine. Profiling first stops the discussion being about the wrong join.
Step 2: Apply the four properties that actually decide it
Speed is one of four, and it is rarely the binding one.
The data volume decides whether in-process is possible at all. GeoPandas holds the working set in memory, so a join whose inputs exceed available RAM either needs tiling or needs a database. Below a few million parcels on a normal machine this constraint simply does not bind.
Reproducibility favours in-process work against immutable snapshots. A run that queries a live database queries whatever the database contained at that moment, and reconstructing that state months later requires either point-in-time recovery or a discipline of versioned tables that most deployments do not have. A run that reads a hashed GeoParquet file is reproducible by construction — which is why the ingest in ingesting ArcGIS feature services writes a snapshot rather than evaluating live.
Testability favours in-process decisively. Rule logic expressed in Python is unit-testable against fixtures with no service running; the same logic expressed in SQL needs a database to test at all, and the fixtures become schema migrations. For a pipeline whose correctness has to be demonstrated, that difference compounds.
Operational fit is the one that most often decides in practice and is least often discussed. A team that already runs PostGIS, whose data lives there, and whose analysts query it directly, will maintain a PostGIS pipeline better than a Python one — and a maintained pipeline beats a theoretically superior abandoned one.
Step 3: Push filtering down, keep judgement up
The arrangement that resolves most of the tension is not to choose but to split along a specific line: the database does selection and the process does adjudication.
import geopandas as gpd
from sqlalchemy import create_engine
def candidates_for_tile(engine, tile_wkt, max_reach_ft):
"""The database finds what is near; it does not decide what is compliant."""
sql = """
SELECT p.parcel_id, p.zoning_district, p.geom
FROM parcels p
WHERE ST_Intersects(p.geom, ST_Buffer(ST_GeomFromText(:wkt, 2263), :reach))
"""
return gpd.read_postgis(sql, engine, geom_col="geom",
params={"wkt": tile_wkt, "reach": max_reach_ft})This gets the property that matters from each side. The spatial index does the work it is best at — reducing a county to a tile — while the rules stay in testable Python, versioned with the rule pack rather than with a schema. It also composes directly with the tiling in chunking county-scale runs by spatial tile.
Step 4: Where the join stays in the database, make it reproducible
If the join genuinely belongs in SQL — very large inputs, or a deployment where that is the only maintained path — the reproducibility gap has to be closed deliberately.
-- Evaluate against a stated snapshot, not against "now".
CREATE TABLE parcels_snapshot_2026_08_08 AS
SELECT * FROM parcels;
CREATE INDEX ON parcels_snapshot_2026_08_08 USING GIST (geom);
-- Record what the run actually read.
INSERT INTO run_manifest (run_id, source_table, row_count, content_hash)
SELECT :run_id, 'parcels_snapshot_2026_08_08', count(*),
md5(string_agg(md5(parcel_id || ST_AsBinary(geom)), '' ORDER BY parcel_id))
FROM parcels_snapshot_2026_08_08;Hashing the snapshot is what makes the SQL path defensible. Without it, “the database” is the answer to what a run read, and that is not an answer — the version and the vintage are what a reviewer needs.
Step 5: Decide per join, not per pipeline
The choice is not global, and treating it as global is what produces both the all-SQL pipeline nobody can test and the all-Python one that cannot handle its largest join.
def route_join(profile, has_maintained_postgis: bool, needs_replay: bool):
if profile["shape"] == "large-large" and profile["selectivity"] < 0.01:
return "database (filter) + process (adjudicate)"
if needs_replay and not has_maintained_postgis:
return "in-process against a snapshot"
if profile["left_rows"] > 5_000_000:
return "database, with a hashed snapshot table"
return "in-process — the simplest thing that is fast enough"The default in the last line is deliberate. In-process against a snapshot is the arrangement with the fewest moving parts and the best testability, and the burden of proof belongs on moving away from it.
Verification
Whichever path is chosen, the two must agree — and that agreement is worth asserting rather than assuming, because the engines differ in edge-case behaviour.
sql_result = gpd.read_postgis(JOIN_SQL, engine, geom_col="geom").sort_values("parcel_id")
py_result = parcels.sjoin(districts, predicate="intersects").sort_values("parcel_id")
assert set(sql_result["parcel_id"]) == set(py_result["parcel_id"]), \
"engines disagree on which parcels match"
area_delta = (sql_result.set_index("parcel_id").area
- py_result.set_index("parcel_id").area).abs().max()
assert area_delta < 1e-6, f"geometry differs between engines by {area_delta}"Two differences are worth knowing about in advance. ST_Intersects and Shapely’s intersects agree on valid geometry and can disagree on invalid geometry, since PostGIS and GEOS versions may differ — which is a reason to repair before joining rather than to prefer one engine. And PostGIS applies its index-backed bounding-box filter before the exact predicate, so an invalid geometry can be excluded by the index in one engine and reach the exact test in the other.
Common Pitfalls
- Choosing on a synthetic benchmark. The shape of your join and the complexity of your geometry dominate; a generic benchmark measures neither.
- Putting rule logic in SQL because the join is there. Selection and adjudication are separable, and the second is where testability is worth the most.
- Querying live tables from a compliance run. The result is unreproducible, and nobody discovers this until a determination is challenged.
- Assuming the index is being used. A join against a geometry column with no GIST index, or with a stale one, degrades to a sequential scan and looks like an engine problem.
- Mixing CRS across the boundary. A geometry sent to PostGIS in the wrong SRID is compared against the wrong coordinates, and PostGIS will do it without complaint if the SRIDs are both declared.
- Ignoring transfer cost. A “fast” database join that returns a million geometries over the wire can be slower end-to-end than an in-process join that never serialises anything.
Frequently Asked Questions
Is PostGIS always faster on large joins?
On large-large joins with a well-maintained GIST index, usually — it has a mature planner and does not materialise intermediate results. The advantage narrows considerably once GeoPandas is reading spatially-sorted GeoParquet with bounding-box pushdown, because much of the database’s edge comes from indexed selection rather than from the join itself, and that part is now available in-process.
What about DuckDB with its spatial extension?
It occupies an interesting middle position: columnar and fast like a database, embedded and snapshot-friendly like in-process work, and it reads GeoParquet directly. For read-only analytical joins over files it is often the best of both, and it does not carry the reproducibility problem that comes with a mutable server. The trade-off is a less mature spatial planner than PostGIS on the hardest joins.
Does the answer change for the nearest-neighbour case?
Yes, and it is the strongest case for the database. PostGIS’s index-backed <-> operator with a LIMIT is genuinely hard to match in-process, because it walks the index in distance order rather than searching a radius and expanding. Measuring distance to the nearest protected feature covers both implementations.
How does this interact with parallelism?
In-process parallelism is straightforward and shares nothing — tiles run in separate workers against separate file reads. Database parallelism concentrates load on one server, which becomes the bottleneck at the point the pipeline most wants to scale. That asymmetry favours in-process for county-scale batch work and matters little for interactive queries.
Should the two paths be maintained in parallel?
Only where the verification above runs as a test. Two implementations that are not continuously compared will diverge, and the divergence is discovered when someone notices two different answers to the same question. If only one can be maintained, keep the in-process one and use the database for selection.
What decides it when everything else is equal?
Who will maintain it. A pipeline is read and modified far more often than it is written, and the engine the team already understands wins on every dimension that shows up after the first month.
Related
Part of: PostGIS vs GeoPandas for setback batch processing
- Benchmarking spatial joins PostGIS vs GeoPandas — the measurement this decision uses.
- Optimizing spatial joins for 100k parcel datasets — making the in-process path fast.
- Chunking county-scale runs by spatial tile — the partitioning the hybrid arrangement uses.
- Converting parcel shapefiles to GeoParquet for batch runs — the snapshot that makes in-process reproducible.
- Tracking data lineage across geospatial ETL steps — recording which path a run took.