Spatial Join vs Attribute Join in GeoPandas

Choosing between merge and sjoin is a decision about what relates your two tables — a shared key or a shared location — and getting it wrong produces either an empty result or a needlessly expensive one. This guide draws the line precisely and shows which is correct for each case; it is for anyone combining two datasets in GeoPandas. It sits under Spatial Joins & Merging in Geospatial Data Ingestion & Processing Workflows.

Why This Approach / What Goes Wrong

An attribute join and a spatial join answer different questions and run on different machinery. DataFrame.merge is a hash join on a key column: it hashes the join key on both sides and matches equal values in roughly linear time, and it never inspects geometry. gpd.sjoin is a spatial join: it builds an R-tree index over one layer's bounding boxes, queries candidate pairs, then evaluates a geometric predicate (within, intersects, contains) on each candidate. If both datasets already share a reliable key — a parcel_id, a census GEOID, an ISO region code — the attribute merge is exact and cheap, and involving geometry only adds cost and ambiguity.

The mistakes go both ways. Reaching for sjoin when a key exists pays for an index build and a predicate pass you did not need, and it introduces boundary ambiguity: a point sitting exactly on a shared edge can match two adjacent polygons, duplicating the left row. Conversely, trying to merge two spatial layers that have no common key returns nothing, because there is no equal value to hash against. The distinction is orthogonal to join direction — inner versus left is a separate choice covered in Performing Left Joins with GeoPandas sjoin.

There is one failure mode unique to the spatial path: the predicate runs on raw coordinates, so if the two layers disagree on their coordinate reference system every predicate evaluates false and you get an all-NaN result — silently, on older GeoPandas. Aligning projections first, ideally against the definitions in Coordinate Systems with PyProj, is not optional. A merge has no such trap because it never looks at coordinates.

The deeper difference is what each mechanism lets you promise about the output. An attribute join has a cardinality you can state in advance and have the library enforce: if parcel_id is unique in the assessment table, the join is one-to-one and pandas will raise if that turns out to be false. A spatial join has a cardinality that is a property of the geometry itself and is therefore unknowable until the predicate has run. Two zoning polygons that overlap by a hand's width will duplicate every parcel in that strip; a district layer with a hairline gap along a river will silently produce unmatched rows. You cannot assert this away, only measure it afterwards — which is why the verification step below counts rows for the merge and match rates for the sjoin. Treating the two as interchangeable enrichment calls is how a population total ends up 4 % high with no error anywhere in the log.

There is also a repeatability difference that matters for anything audited. A key join is deterministic across machines and library versions: the same two CSV files produce the same matched set forever. A spatial join depends on the GEOS build underneath Shapely, on the floating-point result of the reprojection that preceded it, and on which side of a boundary a coordinate lands after rounding. Those inputs are stable in practice but not guaranteed, so a pipeline that reruns a spatial join every night and compares counts to yesterday will occasionally see a single row flip. Keeping the geometric step as far upstream as possible — resolve membership once, store the resulting key, then join on the key thereafter — removes that entire class of drift.

Prerequisites

• geopandas>=0.14
• pandas>=2.0
• shapely>=2.0

conda install -c conda-forge "geopandas=0.14.*" "pandas=2.0.*" "shapely=2.0.*"

Three version boundaries change the answers on this page. GeoPandas 0.10 renamed the sjoin predicate argument from op to predicate; code written against the old name still runs on 0.10–0.13 with a FutureWarning and stops working after that, so a snippet copied from an old answer will fail in a way that looks like an API removal rather than a rename. GeoPandas 0.14 raises a ValueError when the two operands carry different CRSs, where earlier releases only emitted a UserWarning and went on to produce the all-NaN result described above — if you inherit a pipeline that "used to work" and now raises, the join was almost certainly wrong before, not broken now. On the attribute side, pandas 2.0 turned on copy-on-write semantics for many operations, which changes whether merge output shares memory with its inputs; it does not change which rows match, but it does change whether an in-place edit to the result leaks back into parcels. Print geopandas.__version__ and pandas.__version__ at the top of any pipeline whose join behaviour you intend to reason about.

Step-by-Step Implementation

1. When a shared key exists, use an attribute merge. The geometry column rides along untouched; rows match purely on the key. Call merge on the GeoDataFrame (the left side) so the result stays a GeoDataFrame.

import geopandas as gpd
import pandas as pd

# parcels carry geometry + parcel_id; a non-spatial assessment table shares parcel_id
parcels = gpd.read_file("parcels.gpkg")
assessments = pd.read_csv("assessments.csv")    # columns: parcel_id, assessed_value

parcels_valued = parcels.merge(assessments, on="parcel_id", how="left")
# Geometry untouched; rows matched purely on the key — no CRS, no predicate

2. Normalise the key dtype before merging. The single most common empty-merge cause is a key that is an integer on one side and a zero-padded string on the other (5 vs "05"). Coerce both explicitly.

parcels["parcel_id"] = parcels["parcel_id"].astype(str).str.strip()
assessments["parcel_id"] = assessments["parcel_id"].astype(str).str.strip().str.zfill(6)

3. When the only relationship is location, use sjoin. Here air-quality sensors (points) gain the district polygon each falls within. Both layers must share a CRS, and for anything metric — nearest joins, distance filters — that CRS should be a local projected system, never geographic EPSG:4326 or Web Mercator EPSG:3857, whose degrees and distorted metres corrupt distances.

import geopandas as gpd

sensors = gpd.read_file("sensors.gpkg")
districts = gpd.read_file("districts.gpkg")

# Align the CRS FIRST — sjoin evaluates its predicate on raw coordinates
sensors = sensors.to_crs(districts.crs)

sensors_in_district = gpd.sjoin(
    sensors,
    districts[["district_name", "geometry"]],
    how="left",
    predicate="within",   # strict containment; avoids boundary double-matches
)

Name the predicate explicitly. The default is intersects, which counts a point on a shared boundary as matching both neighbouring polygons; within (point strictly inside) removes that ambiguity for point-in-polygon enrichment.

4. Combine both when the enrichment is two-stage — spatially assign a region, then merge regional attributes by the resulting key. This is the common pattern where geometry establishes membership and a plain key carries the payload.

import pandas as pd

region_stats = pd.read_csv("region_stats.csv")   # district_name, avg_income
enriched = sensors_in_district.merge(region_stats, on="district_name", how="left")
# sjoin produced district_name from geometry; merge attaches its attributes by key

5. Declare the cardinality you expect on the merge. validate= turns an assumption into an enforced contract, and indicator= labels every output row with which side it came from. Together they convert the two most expensive silent failures — a duplicated right key that inflates counts, and a systematically unmatched key that quietly nulls a column — into an exception and a printed tally.

parcels_valued = parcels.merge(
    assessments,
    on="parcel_id",
    how="left",
    validate="one_to_one",     # raises MergeError if either side has duplicate keys
    indicator="key_source",    # 'both' or 'left_only' per row
)

print(parcels_valued["key_source"].value_counts())
# both         12841
# left_only      117

validate accepts "one_to_one", "one_to_many", "many_to_one" and "many_to_many". The one worth reaching for by default is "many_to_one": it allows the left layer to repeat a key (several parcels in the same tax district) while forbidding duplicates on the right, which is the shape of almost every real enrichment. There is no equivalent for sjoin — the closest you get is comparing len(result) to len(left) after the fact.

6. Use a composite key when no single column identifies a row. United States census geography is the standard example: tract numbers repeat between counties and county codes repeat between states, so only the tuple is unique. Build the composite explicitly rather than trusting a partial key, and use left_on/right_on when the two sides spell the columns differently.

key_cols = ["state_fips", "county_fips", "tract_code"]
for frame in (census_tracts, acs_estimates):
    for col in key_cols:
        frame[col] = frame[col].astype(str).str.strip()
census_tracts["geoid"] = census_tracts[key_cols].agg("".join, axis=1)
acs_estimates["geoid"] = acs_estimates[key_cols].agg("".join, axis=1)

tracts_enriched = census_tracts.merge(
    acs_estimates.drop(columns=key_cols),
    on="geoid",
    how="left",
    validate="one_to_one",
)

Collapsing the tuple into one string column before the join is worth the extra line: a multi-column merge compares each column independently, so a single mis-typed component fails the whole match with no indication of which one was wrong, whereas a concatenated geoid can be inspected, length-checked (str.len().eq(11).all()), and diffed against the other side with set arithmetic.

7. Know where each mechanism stops scaling. The two have different cost curves, and the crossover is not where intuition puts it. A hash join is roughly linear in the total number of rows and is essentially free up to tens of millions: joining a 5-million-row parcel table to a 200-row district lookup on district_id takes well under a second and a few hundred megabytes. A spatial join pays three separate costs — building the R-tree over the right layer, querying it once per left geometry, and evaluating the predicate on every candidate pair — and the third term is the one that bites, because it scales with the vertex count of the candidate polygons, not their row count. Joining a million points against 200 simple districts is fast; joining the same million points against 200 coastline-detailed administrative polygons with 50,000 vertices each can take minutes.

import time
import geopandas as gpd
from shapely import get_num_coordinates, simplify

# Vertex count, not row count, predicts spatial-join cost
print(int(get_num_coordinates(districts.geometry.values).sum()))   # 1284402

# Pre-filter with the index once, and thin geometry that is more detailed
# than the question needs (10 m tolerance in a metric CRS)
districts_lite = districts.copy()
districts_lite["geometry"] = simplify(districts.geometry.values, tolerance=10)

t0 = time.perf_counter()
tagged = gpd.sjoin(sensors, districts_lite, how="left", predicate="within")
print(f"sjoin: {time.perf_counter() - t0:.2f}s over {len(sensors):,} points")
# sjoin: 3.41s over 1,000,000 points

Simplifying the join geometry is safe only when the tolerance is far below the precision of the question being asked, and it must be done on a copy — the thinned polygons are a join accelerator, not the layer you export. When the left layer itself runs to hundreds of millions of rows the in-memory path is finished regardless of vertex count; push the predicate into an engine that streams, either DuckDB spatial over GeoParquet or Parallel Spatial Joins with Dask-GeoPandas. A key join at that size needs no special handling at all, which is itself a reason to prefer it when both are available.

Verification

Confirm the attribute join preserved row count and the spatial join matched plausibly. A left merge on a unique key must not change the number of rows; a wildly low spatial-match rate almost always means the CRS was misaligned.

# Attribute merge: a left join on a unique key must not change the parcel count
assert len(parcels_valued) == len(parcels), "Key merge duplicated rows — non-unique key on the right"
print("Unmatched assessments:", parcels_valued["assessed_value"].isna().sum())
# Unmatched assessments: 0

# Spatial join: most sensors should land in some district
matched = sensors_in_district["district_name"].notna().mean()
print(f"Sensors matched to a district: {matched:.1%}")
# Sensors matched to a district: 98.7%
assert matched > 0.5, "Few matches — almost certainly a CRS mismatch between inputs"

The spatial side needs one more check that the merge does not: whether the join expanded the left layer. Because sjoin emits one row per matching pair, a duplicated left index is the signature of overlapping right geometries, and it is the difference between "98.7 % matched" being reassuring and being meaningless.

dupes = sensors_in_district.index.duplicated().sum()
print(f"Left rows duplicated by multiple matches: {dupes}")
# Left rows duplicated by multiple matches: 0

# Which districts overlap, if the count is non-zero
if dupes:
    offenders = (sensors_in_district.index.value_counts()
                 .loc[lambda s: s > 1].index)
    print(sensors_in_district.loc[offenders, "district_name"].unique())

A useful cross-check when both routes are available is to run them against each other. If the district layer also carries the district_id that the sensors table already stores, join both ways and compare: the spatial result should agree with the key result on every row, and the rows where it disagrees are exactly the records whose stored key is stale or whose coordinates are wrong. That disagreement set is far more valuable than either join on its own, and it is the cheapest data-quality report you will ever write.

# A plain lookup table: the district attributes without their geometry
districts_lookup = districts.drop(columns="geometry")

by_key = sensors.merge(districts_lookup, on="district_id", how="left")
disagree = (by_key["district_name"].values
            != sensors_in_district["district_name"].reindex(by_key.index).values)
print(f"Rows where the stored key and the geometry disagree: {disagree.sum()}")
# Rows where the stored key and the geometry disagree: 23

Edge Cases & Debugging

• sjoin returns all nulls. The inputs are in different CRSs; reproject one with .to_crs() before joining. A merge never shows this symptom, which is a quick way to confirm the join type was correct.
• merge returns no matches. The key dtypes or formats differ (int vs zero-padded str); normalise both sides before merging.
• Row count explodes after merge. The key is not unique on the right side; drop_duplicates(subset=key) or aggregate the right table before merging.
• Points on borders match two polygons in sjoin. Switch from the intersects default to predicate="within", or de-duplicate on the left index after the join.
• You used sjoin where a key existed. It is slower and adds boundary ambiguity for no benefit — prefer merge whenever a reliable shared key is present.
• Lost geometry after merge. You called merge on the plain DataFrame instead of the GeoDataFrame; put the geometry-bearing frame on the left so the geometry column survives.
• sjoin raises ValueError about the CRS on a pipeline that used to run. GeoPandas 0.14 promoted the old mismatch warning to an error. The join was returning nulls before; add the .to_crs() call rather than pinning back to the older release.
• Both sides have a district_name column and the output has neither. merge appended _x and _y suffixes. Set suffixes=("", "_ref") so the left name survives unchanged, or drop the colliding column from the right frame before the join.
• merge silently turns integer columns into floats. A left join that leaves some rows unmatched has to represent the missing values, and NaN forces float64. Cast back with astype("Int64") — the nullable integer type — once you have decided what an unmatched row means.
• Two sjoin calls in a row fail on index_right. The first join leaves an index_right column behind, and the second refuses to create a second one. Drop it (.drop(columns="index_right")) between joins, or rename it to something meaningful such as district_idx as soon as the first join returns.

Frequently Asked Questions

My two layers share a key, but I do not trust it. Should I use the key or the geometry? Use both, and treat the disagreement as the deliverable. Run the key merge and the spatial join separately, then compare the assigned region column row by row as shown in the verification section. Where they agree, you have corroboration from two independent sources; where they disagree, one of the two is wrong and you now know exactly which records to inspect. Picking one blindly means shipping whichever error happened to be in the source you chose. Once you have reconciled the records and written the corrected key back to storage, switch to the key join permanently — it is cheaper, deterministic, and no longer dependent on the geometry being current.

Is sjoin ever faster than merge for the same enrichment? No, not for an enrichment that both can express. The hash join reads each key once and matches by equality; the spatial join has to build an index, run a bounding-box query per left feature, and then run an exact predicate on every surviving candidate, where the predicate cost grows with the vertex count of the polygons involved. Even on a tiny right layer the spatial route is typically one to two orders of magnitude slower. The only situation where the spatial call wins is the one where the key does not exist at all, in which case the comparison is not a performance question — the merge simply returns nothing.

Which predicate should I use for point-in-polygon work, and does the choice affect the join type decision? For strict point-in-polygon enrichment use within, which requires the point to be in the polygon's interior and therefore cannot match two neighbours across a shared edge. intersects, the default, treats boundary contact as a match and will duplicate any point that lands exactly on a border — rare with real-world floating-point coordinates but common with data that has been snapped to a grid, because snapping is precisely what puts coordinates on the border. The predicate choice is independent of the merge-versus-sjoin decision: it only matters once you have already established there is no shared key. Direction (how="inner" versus how="left") is a third, separate axis, covered in Performing Left Joins with GeoPandas sjoin.

How do I join polygons to polygons when they only partially overlap? A predicate-based join answers a yes-or-no question, so intersects will attach every zoning polygon that touches a parcel at all, including one clipping a corner by a square metre. If you need the dominant overlap rather than all of them, join with intersects and then rank the matches by intersection area, keeping the largest per left row. If you need the geometry of the overlapping pieces themselves — one output feature per intersecting pair, clipped to the shared area — that is not a join at all but an overlay, described in Computing Overlay Union & Difference with GeoPandas.

Do I need to reproject before an attribute merge? No. merge never reads a coordinate, so the CRS of the left frame is irrelevant to whether rows match, and the geometry column passes through untouched with its CRS metadata intact. This is a real practical advantage: an attribute enrichment can run on a layer straight out of storage without a reprojection pass. Reproject when the next operation is metric — a buffer, an area, a distance — and pick a local projected system for it, ideally one derived from the data as in Choosing a UTM Zone Automatically in Python, never Web Mercator.

Should I dissolve before or after the join? After, in almost every case. Dissolving first destroys the per-feature attributes that make the join meaningful and usually produces geometry with many more vertices, which makes any subsequent spatial join slower. Join at the finest granularity you have, then aggregate the enriched result — the grouping and aggregation mechanics are in Dissolving and Aggregating Features by Attribute. The exception is when the right layer is a fragmented version of the regions you actually want, in which case dissolving that side first both fixes the boundaries and reduces the candidate count.