Spatial Joins & Merging in Python Pipelines

Spatial joins and merging are the stage of Geospatial Data Ingestion & Processing Workflows where separately ingested layers are stitched into a single analysis-ready table by geometric relationship rather than by a shared key. A point layer of sensors gains its administrative district; parcels inherit the flood zone they fall inside; fragmented polygons collapse into one authoritative boundary. This stage runs directly after Coordinate Reference System Transformations — because a join across mismatched projections silently returns nothing — and it leans on the vectorized predicate engine documented under Shapely Geometry Operations. This guide covers CRS-safe joins with GeoPandas, predicate selection, cardinality control, topology-safe geometric merging, and scaling the same operation to hundreds of millions of rows.

Architecture & Data Structures

A spatial join has two inputs — a left and a right GeoDataFrame — and one relationship, the predicate. geopandas.sjoin walks the left geometries, uses the right layer's spatial index to find candidate matches, tests the predicate on each candidate, and emits one output row per matching pair. The result is a GeoDataFrame that keeps the left geometry column and appends the right layer's attribute columns, together with an index_right column recording which right feature matched. Nothing about the coordinates changes; only attributes are transferred.

The engine underneath is an R-tree, exposed as GeoDataFrame.sindex. GeoPandas builds it lazily on first access and caches it, which is why the naive O(n²) pairwise comparison never actually runs — each left geometry is tested only against right features whose bounding boxes overlap. The bounding-box pass is a fast filter; the predicate is the exact refinement.

import geopandas as gpd

# Two layers, already normalized to one projected CRS (see CRS section below)
sensors = gpd.read_parquet("sensors_utm.parquet")      # points, EPSG:25832
districts = gpd.read_parquet("districts_utm.parquet")  # polygons, EPSG:25832

# One output row per (sensor, district) pair where the predicate holds.
# index_right carries the matched district's index; predicate is explicit.
tagged = gpd.sjoin(sensors, districts, how="inner", predicate="within")
print(tagged.columns)   # sensor cols + district cols + 'index_right'

The how parameter decides which unmatched rows survive: inner keeps only matched pairs, left keeps every left feature (unmatched ones get NaN in the right columns), and right keeps every right feature. The distinction between transferring attributes this way and a plain key-based merge is important enough to have its own deep dive in Spatial Join vs Attribute Join in GeoPandas.

Underneath the convenience wrapper is an operation you can call directly, and doing so once is the fastest way to understand what the join actually costs. sindex.query takes an array of geometries and a predicate and returns a two-row integer array: the first row holds positional offsets into the queried layer, the second holds positional offsets into the indexed one. That array is the join, in its rawest form — everything sjoin adds on top is index bookkeeping, column concatenation, and the how logic that reinstates unmatched rows.

import numpy as np
import geopandas as gpd

# The raw pair array: nothing but integer offsets, no attribute copying yet
pairs = districts.sindex.query(sensors.geometry, predicate="contains")
left_pos, right_pos = pairs          # shape (2, n_matches)

print(pairs.shape[1], "matching pairs from",
      len(sensors), "x", len(districts), "candidate combinations")

# Reconstruct the join by position — .iloc, never .loc, because these are offsets
matched = sensors.iloc[left_pos].copy()
matched["district_id"] = districts["district_id"].to_numpy()[right_pos]

Two properties of that array explain most surprises. It is positional, so feeding the values to .loc on a frame with a non-default index quietly returns the wrong rows — a bug that only appears once someone upstream stops calling reset_index. And its length is the true size of the operation: the pair count, not the input row count, determines how much memory the join needs, because each pair becomes an output row carrying a copy of both sides' attributes. A hundred thousand points against a hundred thousand overlapping polygons is a small join; a hundred thousand against a hundred densely overlapping ones can be a much larger one.

The predicate is evaluated inside the index query, not afterwards. GEOS tests bounding boxes to build the candidate set and then applies the real predicate to each candidate before the pair is emitted, which means an expensive predicate on a layer with few box overlaps is cheap, while a cheap predicate on a layer of long, thin, heavily overlapping envelopes — road centrelines, river networks — is not. Geometry shape, more than row count, drives the cost.

Environment Configuration & Dependency Resolution

Spatial joins depend on the Shapely 2.0 vectorized predicate path and on GeoPandas 0.14+, where predicate (not the removed op) is the argument name and sjoin_nearest is available. Version skew between GeoPandas, Shapely, and the GEOS that Shapely binds is the usual cause of "predicate returns different matches on another machine." Pin the whole stack from conda-forge so GEOS, GDAL, and PROJ are built together.

conda install -c conda-forge \
  "geopandas>=0.14" "shapely>=2.0" "pyproj>=3.6" \
  "duckdb>=0.10" "dask-geopandas>=0.3"

Confirm the resolved libraries at startup — a Shapely still linked against GEOS 3.8 will silently fall back to the slow, non-vectorized predicate loop.

import geopandas, shapely
print("geopandas:", geopandas.__version__)
print("shapely:", shapely.__version__)
print("GEOS:", shapely.geos_version)   # (3, 12, x) on a current conda-forge build

Four version boundaries change what this page's code does rather than merely how fast it runs, and all four sit inside the range of GeoPandas releases still in production use. GeoPandas 0.10 introduced the frame method form (left.sjoin(right)) alongside the module function and began deprecating op= in favour of predicate=. GeoPandas 0.13 removed the optional pygeos and rtree index backends entirely and made Shapely 2.0 mandatory, so the spatial index is now always a Shapely STRtree and the old advice about choosing an index engine no longer applies to anything. GeoPandas 0.14 folded sindex.query_bulk into sindex.query, which now accepts either a single geometry or an array — code written against query_bulk still runs but emits a deprecation warning. And the dwithin predicate, which answers "within d units of" without building a buffer, requires GEOS 3.10 or newer underneath Shapely; on an older GEOS it raises rather than degrading, so gate on shapely.geos_version before shipping it.

The practical consequence is that a "works on my machine" spatial join is nearly always an environment report, not a code report. When a join returns different matches on two machines, print the four version numbers above from both before touching the code — a GEOS difference changes boundary-case predicate results at the last bit of floating-point precision, and that is enough to move a point that sits exactly on a shared edge from one polygon to the other.

Vectorized Operations & Core Workflow

The canonical workflow is: load both layers, force them onto one projected CRS, run the join with an explicit predicate, then collapse any one-to-many expansion the join produced. The example below tags each air-quality sensor with the district it sits inside and never assumes the two files arrived in the same projection.

import geopandas as gpd

TARGET_EPSG = 25832   # ETRS89 / UTM 32N — metric, correct for the study area

sensors = gpd.read_file("sensors.geojson")          # likely EPSG:4326
districts = gpd.read_file("districts.gpkg")         # may be a national grid

# 1. Bring BOTH operands onto the identical projected CRS before joining.
sensors = sensors.to_crs(epsg=TARGET_EPSG)
districts = districts.to_crs(epsg=TARGET_EPSG)

# 2. Join: keep every sensor, attach the district it falls within.
tagged = gpd.sjoin(
    sensors,
    districts[["district_id", "district_name", "geometry"]],
    how="left",
    predicate="within",
)

# 3. A sensor on a shared boundary can match two districts (1:N). Keep the
#    first deterministic match so the left cardinality is preserved.
tagged = tagged[~tagged.index.duplicated(keep="first")]

print(f"{tagged['district_id'].isna().sum()} sensors fell outside all districts")

Selecting only the columns you need from the right layer (district_id, district_name, geometry) keeps the output narrow and avoids column-name collisions. When names do clash, control the disambiguating suffixes with lsuffix and rsuffix rather than accepting the defaults.

Keeping the first match is only defensible when the duplicates are an artefact — a point on a border that genuinely belongs to one district. When the multiplicity is real information, throwing it away is the bug. A parcel that overlaps three zoning categories is in three zoning categories, and the honest reduction is an aggregation rather than a choice: collect the matched values into a list, count them, or pivot them into indicator columns. All three keep one row per left feature without pretending the ambiguity does not exist.

import geopandas as gpd

overlapping = gpd.sjoin(parcels, zones[["zone_code", "geometry"]],
                        how="left", predicate="intersects")

# One row per parcel, with the full match set preserved rather than truncated
summary = (
    overlapping.groupby(level=0)
    .agg(zone_codes=("zone_code", lambda s: sorted(s.dropna().unique())),
         zone_count=("zone_code", "nunique"))
)

parcels_tagged = parcels.join(summary)
parcels_tagged["is_split_zoning"] = parcels_tagged["zone_count"] > 1

print(parcels_tagged["zone_count"].value_counts().sort_index())
# 0      412      parcels in no zone at all
# 1    18693      unambiguous
# 2      1104     genuinely split between two zones

Grouping on level=0 rather than a column is deliberate: the left index survives a spatial join intact, so it is the only key guaranteed to identify a left feature even when the layer has no id column of its own. The is_split_zoning flag is worth carrying downstream — it converts a silent data-quality issue into a queryable column, and it tells a reviewer which rows a later area-weighted decision was applied to.

Geometry & Data Processing Details

Choosing the predicate. The predicate defines the topological test, and picking the wrong one is the most common correctness bug in a join. within and contains are strict interior/boundary relationships (a point exactly on a polygon edge is not within it under some geometries); intersects is the permissive "share any point" test; covers/covered_by behave like contains/within but include boundary touches; crosses and touches are rarer edge relationships. For point-in-polygon tagging, within is usually right but leaks features that land precisely on shared edges — switch to intersects if boundary sensors must be retained (accepting that they then match both neighbours).

# Point-in-polygon: which flood zone does each parcel centroid fall in?
parcels["rep_point"] = parcels.representative_point()   # guaranteed inside
lookup = parcels.set_geometry("rep_point")
zoned = gpd.sjoin(lookup, flood_zones, how="left", predicate="within")

Using representative_point() instead of centroid matters for concave or multipart polygons, whose centroid can fall outside the polygon and silently miss its own zone.

The predicates are not a loose vocabulary; they are named cases of the DE-9IM intersection model, and the distinctions that look pedantic are exactly the ones that decide borderline features. The model describes the relationship between two geometries as a matrix of intersections between their interiors, boundaries and exteriors, and each predicate is a pattern over that matrix. contains requires the contained geometry's interior to meet the container's interior and forbids any part of it from touching the container's exterior — but it also requires at least one interior point in common, which is why a polygon does not contain a line lying entirely along its own boundary. covers drops that last requirement, so it is true whenever nothing sticks out. For point-in-polygon work the difference shows up on exactly one class of feature: a point sitting on a shared edge is covered_by both neighbours, within neither, and intersects both.

from shapely.geometry import Point, Polygon

block = Polygon([(0, 0), (100, 0), (100, 100), (0, 100)])
corner_sensor = Point(100, 50)      # exactly on the eastern boundary

print(block.contains(corner_sensor))    # False — boundary is not the interior
print(block.covers(corner_sensor))      # True  — nothing lies outside the block
print(block.touches(corner_sensor))     # True  — boundaries meet, interiors do not
print(block.intersects(corner_sensor))  # True  — the permissive catch-all

Which one is correct depends on whether the boundary belongs to the feature in your domain, and that is a data-governance question, not a geometry one. Cadastral parcels usually share a boundary that belongs to both neighbours; administrative reporting usually needs each point counted exactly once. When the second rule applies and points genuinely land on edges, within plus a deterministic tie-break beats intersects plus deduplication, because the tie-break is then visible in the code instead of hidden in row ordering.

Cardinality. A spatial join is a relational join: if a parcel overlaps three zoning polygons under intersects, the output contains three rows for that parcel. Decide the collapse rule explicitly — keep the largest-overlap match, aggregate the matched attributes into a list, or deduplicate on the left index as shown above. Ignoring this is how downstream counts silently inflate.

Nearest-neighbour joins. When features do not overlap at all — matching each sensor to the closest road segment — use sjoin_nearest, which pairs every left feature with its nearest right feature and can return the gap in a distance_col. Because the search radius is expressed in CRS units, this operation is only meaningful on a projected (metric) CRS; the deeper trade-off against a raw KD-tree lives in Nearest Neighbor & KD-Tree Search.

# Attach the closest road within 250 m; distance in METRES (projected CRS)
near_road = gpd.sjoin_nearest(
    sensors, roads, how="left", max_distance=250, distance_col="dist_m",
)

Merging geometries. Joining transfers attributes; merging fuses geometry. dissolve() groups by an attribute and unions each group's geometries into one — the standard way to roll parcels up to a neighbourhood. For lower-level topological consolidation, shapely.ops.unary_union merges a raw geometry collection. Both are sensitive to slivers where boundaries almost, but not exactly, coincide; snapping to a small tolerance before the union removes the hairline gaps. Overlay-style set operations (union, difference, intersection between whole layers) are a separate concern covered in Geometric Intersections & Overlays.

import geopandas as gpd
from shapely import set_precision
from shapely.ops import unary_union

def dissolve_safe(gdf: gpd.GeoDataFrame, group_col: str, grid: float = 0.001):
    """Dissolve by attribute, snapping to a precision grid to kill slivers."""
    rows = []
    for key, group in gdf.groupby(group_col):
        # set_precision snaps coordinates onto a grid (metres) before union,
        # collapsing near-coincident boundaries that would leave slivers.
        snapped = [set_precision(g, grid) for g in group.geometry]
        rows.append({group_col: key, "geometry": unary_union(snapped)})
    out = gpd.GeoDataFrame(rows, crs=gdf.crs)
    return out[out.geometry.is_valid].reset_index(drop=True)

Any input geometry that is invalid before the union (self-intersecting rings, unclosed shells) will corrupt the result, so run the merge only on layers that have passed Topology Validation & Repair.

CRS Alignment & Projection Pipeline

A spatial join makes no sense across two coordinate systems, and GeoPandas enforces this: joining layers with different CRSs raises a ValueError, and joining a layer whose CRS is None proceeds on raw numbers that almost never align. The rule is absolute — normalize both operands to one CRS before the join, and make that CRS projected whenever distance or area enters the logic.

from pyproj import CRS

TARGET = CRS.from_epsg(25832)   # metric, project-wide target

# Fail loudly instead of joining mislabelled or mismatched coordinates.
for name, layer in {"sensors": sensors, "districts": districts}.items():
    if layer.crs is None:
        raise ValueError(f"{name} has no CRS — backfill it before joining")

sensors = sensors.to_crs(TARGET)
districts = districts.to_crs(TARGET)
assert sensors.crs.equals(districts.crs), "operands must share one CRS"

Reprojection is not a neutral operation for a join, and this is the part that surprises people who have already got the CRS right. Coordinates move, but the edges between them are still interpreted as straight lines in whatever plane they now live in. A polygon boundary that ran straight in a national grid becomes a slightly curved true path after reprojection, yet is still stored and tested as a straight segment between the transformed endpoints — so a point that sat a few centimetres inside the original polygon can fall outside the reprojected one. On administrative boundaries with vertices kilometres apart, the discrepancy is metres. Densify long edges before reprojecting when the join has to be defensible at the boundary, and reproject both layers from their native definitions rather than chaining transformations through an intermediate CRS, since each hop repeats the error.

import geopandas as gpd

# Add vertices every 100 m so reprojection follows the true boundary path
districts = districts.copy()
districts["geometry"] = districts.segmentize(100)   # units of the SOURCE CRS
districts = districts.to_crs(TARGET)

The antimeridian is the other geometry-level failure, and it does not announce itself. A layer spanning ±180° in a geographic CRS holds features whose coordinate lists jump from 179.9 to −179.9, which every planar predicate reads as a shape wrapping the entire globe. Every join against such a feature matches almost everything. If your study area touches the seam — the Pacific, eastern Russia, Fiji, New Zealand's outlying islands — either split the geometries at ±180° before joining or work entirely in a projected CRS whose central meridian sits near the data, which moves the discontinuity somewhere harmless.

The subtle trap is sjoin_nearest's max_distance and any buffer built before an intersects join: both are in CRS units. Set max_distance=250 on an EPSG:4326 layer and you have asked for 250 degrees — effectively unbounded. Reproject to a local UTM zone or national grid first, and never run metric joins in EPSG:4326 or EPSG:3857, whose scale distorts with latitude. The full projection mechanics, including always_xy and datum shifts, are in Coordinate Reference System Transformations and Coordinate Systems with PyProj.

Production Export & Integration

Once a layer exceeds a few million rows, the in-memory GeoPandas join becomes the bottleneck, and the operation is better pushed down into an engine built for it. DuckDB spatial evaluates ST_Intersects in C++ and reads GeoParquet directly, so the join runs where the data lives. Both operands must already be in the same CRS — DuckDB's spatial functions are CRS-agnostic and will happily join misaligned coordinates.

import duckdb

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;")

# Read GeoParquet natively — no WKT round-trip needed on recent DuckDB.
result = con.execute("""
    SELECT s.sensor_id, d.district_name
    FROM read_parquet('sensors_utm.parquet') s
    JOIN read_parquet('districts_utm.parquet') d
      ON ST_Intersects(s.geometry, d.geometry)
""").df()

For datasets that outgrow one node, Scaling with Dask-GeoPandas partitions both layers and runs the join per partition; spatially shuffling first with a Hilbert curve keeps matching features co-located so each partition compares only nearby data — the detailed recipe is Parallel Spatial Joins with Dask-GeoPandas.

import dask_geopandas as dgpd

parcels = dgpd.read_parquet("s3://bucket/parcels_utm/")
zones = dgpd.read_parquet("s3://bucket/zones_utm/")

parcels = parcels.spatial_shuffle(how="hilbert")   # co-locate nearby features
joined = dgpd.sjoin(parcels, zones, how="left", predicate="intersects")
joined.to_parquet("s3://bucket/parcels_zoned/")     # partitioned GeoParquet out

When the join result is destined for a database, load it via PostGIS integration and let the server run the join behind a GiST index; the SRID stored on each geometry column must match the CRS you reprojected to. Pushing the predicate into SQL is worth it precisely when the result is small relative to the inputs — a count or an aggregate per polygon — because then nothing crosses the wire but the answer.

-- Runs against the GiST index; only the aggregate returns to Python
SELECT d.district_id,
       count(*)            AS sensor_count,
       avg(s.reading)      AS mean_reading
FROM   sensors  s
JOIN   districts d
  ON   ST_Intersects(d.geom, s.geom)
GROUP  BY d.district_id;

Knowing when to make that move matters more than knowing how. The in-memory join has three distinct scaling limits and they bite in a predictable order. Index construction is the cheapest and most forgiving: an STRtree over a few million bounding boxes builds in seconds and costs tens of megabytes, so the tree is almost never the problem. Predicate evaluation is next, and it scales with total vertex count rather than row count — ten thousand coastline polygons with a thousand vertices each cost far more than a million rectangles. The wall is the output: because every matching pair materialises a full row, a join whose average multiplicity is three turns a five-million-row left layer into fifteen million rows, each carrying both layers' attribute columns, and the concatenation step needs that entire result in memory at once alongside the inputs.

That ordering suggests the cheap interventions before the architectural one. Trim the right layer's columns to the two or three you actually need, so each output row is small. Simplify heavy geometries when the analysis tolerance allows it. Chunk the left layer and write each chunk's result straight to partitioned GeoParquet rather than concatenating in memory, which caps peak usage at one chunk regardless of total size. Only when those run out — typically past ten million left rows or a genuinely global extent — is a distributed or database-backed join the right answer rather than an expensive detour.

Production checklist

1. Reproject both operands to one projected CRS, then assert they are equal.
2. Pass predicate explicitly — never rely on a default.
3. Decide and apply a cardinality-collapse rule immediately after every join.
4. Express max_distance / buffers in metres on a metric CRS, never degrees.
5. Validate geometry before any dissolve/unary_union; snap to a precision grid to prevent slivers.
6. Export intermediates as GeoParquet (CRS preserved, predicate pushdown) rather than Shapefile.

Windows & Platform Edge Cases & Debugging

• sjoin returns zero matches with no error. The operands are in different CRSs (or one is None), so bounding boxes never overlap. Reproject both to the identical CRS and re-run.
• ValueError: 'left_df' and 'right_df' should have the same crs. GeoPandas is protecting you — call .to_crs() on one operand instead of suppressing the check.
• Output row count is larger than the left layer. A one-to-many match under intersects; deduplicate on the left index or aggregate the matches, as in the core workflow.
• AttributeError: 'GeoDataFrame' object has no attribute 'sindex' / slow joins on Windows. GEOS is not linked or is an old version, usually from mixing pip and conda wheels; reinstall the whole stack from a single conda-forge environment so Shapely 2.0's vectorized predicates are available.
• sjoin_nearest distances look absurd (hundreds of thousands). The layer is still geographic (degrees); reproject to a UTM zone before measuring.
• TypeError: sjoin() got an unexpected keyword argument 'op'. op was removed; the argument is predicate in GeoPandas 0.14+.
• The same join returns a different match count on two machines. GEOS versions differ in last-bit floating-point behaviour on boundary cases. Print shapely.geos_version on both, pin the environment, and stop relying on intersects for features that sit exactly on shared edges.
• predicate="dwithin" raises instead of running. The linked GEOS predates 3.10. Either upgrade the stack or fall back to sjoin_nearest with max_distance, which answers a related question using a different code path.
• Memory spikes at the end of a join that ran fine for minutes. The index query finished and the result concatenation began. The output is larger than either input; chunk the left layer and stream each chunk's result to disk instead of building one frame.
• .loc on the joined frame returns several rows. A left feature matched more than once, so the index is no longer unique. Collapse or aggregate before any label-based lookup — the join is not the bug, the assumption of uniqueness is.
• A join against a layer read from a Shapefile matches nothing on Windows only. The .prj was missing in that copy, so the layer arrived as crs=None and never reprojected. Assert crs is not None at ingestion, as described in Shapefile & GeoJSON Parsing.

Frequently Asked Questions

What is the difference between a spatial join and merge? merge matches rows on equal key values; a spatial join matches on a geometric predicate such as within or intersects. Use merge when both tables share an id column and a spatial join when the only relationship is location. The two are compared in detail in Spatial Join vs Attribute Join in GeoPandas.

Which predicate should I use for point-in-polygon tagging? within for the strict "point inside polygon" case, using representative_point() so a concave polygon's own points are never missed. Switch to intersects only if points that land exactly on shared boundaries must be retained — accepting that those points then match both adjacent polygons.

Why does my join produce more rows than I started with? A left feature matched several right features (one-to-many). Collapse the expansion right after the join by deduplicating on the left index or aggregating the matched attributes; the join itself is behaving correctly. See Performing Left Joins with GeoPandas sjoin.

Do both layers really need the same CRS? Yes. GeoPandas raises a ValueError on mismatched CRSs and gives wrong answers when one is None. Reproject both to one projected CRS before joining, and make it metric whenever max_distance, buffers, or areas are involved.

When should I move the join out of GeoPandas? Around the low millions of rows the in-memory join becomes the bottleneck. Push it into DuckDB spatial for single-node scale, or partition it with Dask-GeoPandas for multi-node workloads.

Should I build the spatial index myself before joining? No — sjoin builds and caches it on first use, and touching sindex beforehand only moves the same work earlier. The one case where explicit index use pays is when you need the raw pair array rather than a joined frame: counting matches, testing a hypothesis about multiplicity, or joining with a custom rule that how= cannot express. Then sindex.query gives you the integer offsets directly, with no attribute copying at all.

Does a spatial join work correctly in a geographic CRS? For containment questions, yes: within and intersects on longitude-latitude coordinates give the right answer everywhere except across the antimeridian and at the poles, because the topological relationship does not depend on the units. For anything involving distance or area it does not — max_distance, dwithin, buffers and largest-overlap tie-breaks are all meaningless in degrees, since a degree of longitude shrinks from about 111 km at the equator to nothing at the poles. Project first whenever a number, rather than a yes or no, comes out of the operation.

How do I join two point layers? Not with sjoin and an overlap predicate — two independently measured points essentially never share coordinates exactly, so intersects returns nothing. Use sjoin_nearest with a max_distance ceiling on a projected CRS, which pairs each left point with its closest right point and lets you reject matches that are too far apart to be credible. For very large point-to-point workloads a KD-tree beats the R-tree path; the measured comparison is in sjoin_nearest vs cKDTree Performance.

Which side of the join should the smaller layer go on? Semantics decide first: how="left" protects the left layer, so the layer that must keep all its rows goes there. Within that constraint the tree is built over the right layer and queried once per left geometry, so a large left layer against a small right one — many sensors, few districts — is both the natural phrasing and the cheaper arrangement. Reversing the operands to chase performance changes which geometry column the result carries, which is rarely a trade worth making.