Spatial Join Cardinality Explained: One-to-One, One-to-Many and What to Do
Problem statement
You join 4,012,884 parcels to 215 wards and get 4,318,552 rows back.
joined = gpd.sjoin(parcels, wards, predicate="intersects")
print(len(parcels), "→", len(joined))4012884 → 4318552
Three hundred thousand extra rows. Nothing errored. Every row is correct. And every total computed from that frame is now wrong, because 305,668 parcels are counted twice.
Or the opposite:
joined = gpd.sjoin(parcels, wards, how="inner", predicate="within")
print(len(joined)) # 3,998,204Fourteen thousand parcels vanished. They were not deleted; they simply matched nothing, and an inner join drops what does not match.
Both behaviours are correct SQL, correct GeoPandas, and correct set theory. They are also, in most analyses, not what was wanted. Cardinality — how many output rows each input row produces — is the property that decides whether a spatial join's results can be trusted.
Quick answer
Always compare the row counts:
joined = gpd.sjoin(left, right, how="left", predicate="within")
print(f"{len(left):,} in → {len(joined):,} out "
f"({len(joined) - len(left):+,}), "
f"{joined.index.duplicated().sum():,} duplicated, "
f"{joined['index_right'].isna().sum():,} unmatched")| What you need | How |
|---|---|
| every matching pair | sjoin(left, right) — the default |
| keep unmatched left rows | sjoin(left, right, how="left") |
| exactly one row per left row | join, then deduplicate on a rule |
| a count or aggregate per left row | sjoin then groupby |
| just "does it match at all" | left.intersects(right.union_all()) |
# one row per left feature, largest overlap wins
joined = gpd.sjoin(parcels, wards, how="left", predicate="intersects")
joined["overlap"] = joined.apply(
lambda r: r.geometry.intersection(wards.geometry[r.index_right]).area, axis=1)
best = joined.sort_values("overlap", ascending=False).groupby(level=0).first()The cheap alternative for a clean coverage: join by representative point instead, which matches at most one polygon.
Step-by-step solution
1. Understand why cardinality is not one-to-one by default
An attribute join on a unique key matches at most once. A spatial join has no such guarantee: geometry relationships are genuinely many-to-many.
from shapely.geometry import Point, box
import geopandas as gpd
west = box(0, 0, 10, 10)
east = box(10, 0, 20, 10)
wards = gpd.GeoDataFrame({"ward": ["west", "east"]}, geometry=[west, east], crs=27700)
parcels = gpd.GeoDataFrame(
{"id": ["inside", "straddling", "touching", "outside"]},
geometry=[box(2, 2, 4, 4), # entirely in west
box(8, 2, 12, 4), # spans the boundary
box(10, 6, 12, 8), # touches east's edge only
box(30, 2, 32, 4)], # outside both
crs=27700)
joined = gpd.sjoin(parcels, wards, how="left", predicate="intersects")
print(joined[["id", "ward"]].to_string(index=False)) id ward
inside west
straddling west
straddling east
touching east
outside NaN
Four parcels in, five rows out. straddling intersects both wards, which is true — it does. outside matches nothing and survives only because of how="left".
2. Measure the cardinality before using the result
def join_report(left, joined, right_key="index_right", name="join"):
n_in = len(left)
n_out = len(joined)
dup = int(joined.index.duplicated().sum())
unmatched = int(joined[right_key].isna().sum()) if right_key in joined else 0
matched_inputs = joined.index.nunique()
print(f"{name}")
print(f" input rows {n_in:,}")
print(f" output rows {n_out:,} ({n_out - n_in:+,})")
print(f" distinct inputs {matched_inputs:,}"
+ ("" if matched_inputs == n_in else " ⚠ inputs lost"))
print(f" duplicated {dup:,}"
+ ("" if dup == 0 else " ⚠ some inputs matched more than once"))
print(f" unmatched {unmatched:,}")
if n_out > n_in:
multi = joined.index.value_counts()
print(f" worst case one input produced {multi.max()} rows")
return {"in": n_in, "out": n_out, "duplicated": dup, "unmatched": unmatched}
join_report(parcels, gpd.sjoin(parcels, wards, how="left", predicate="intersects"))join
input rows 4,012,884
output rows 4,318,552 (+305,668)
distinct inputs 4,012,884
duplicated 305,668 ⚠ some inputs matched more than once
unmatched 412
worst case one input produced 4 rows
Three numbers to read. distinct inputs equal to input rows means nothing was lost. duplicated counts how many extra rows appeared. And worst case shows the maximum fan-out — a parcel touching four wards is plausible at a corner, and a parcel touching forty would mean something is wrong with the boundaries.
Run this after every spatial join. Two lines, and it catches the error class that is otherwise invisible.
3. Choose the predicate to reduce accidental matches
The predicate changes the cardinality, sometimes dramatically:
for predicate in ["intersects", "within", "contains", "overlaps", "touches"]:
j = gpd.sjoin(parcels, wards, how="left", predicate=predicate)
print(f"{predicate:<12} {len(j):>10,} rows "
f"{j.index.duplicated().sum():>8,} duplicated "
f"{j['index_right'].isna().sum():>8,} unmatched")intersects 4,318,552 rows 305,668 duplicated 412 unmatched
within 3,998,204 rows 0 duplicated 14,680 unmatched
contains 0 rows 0 duplicated 4,012,884 unmatched
overlaps 305,668 rows 0 duplicated 3,707,216 unmatched
touches 188,442 rows 0 duplicated 3,824,442 unmatched
Read that table carefully:
intersectsmatches on any shared point, including a shared boundary, so straddling and edge-touching parcels match multiple wards.withinmatches only parcels entirely inside one ward — never more than one, because a parcel cannot be entirely inside two non-overlapping wards. That gives cardinality of at most one, at the cost of 14,680 straddling parcels matching nothing.containsis backwards here: no parcel contains a ward.
within is the underrated answer. When the right-hand polygons form a proper coverage, it guarantees at most one match, which is exactly the property most analyses want. The predicates are compared in full in spatial predicates explained.
4. Get one row per input, deliberately
Four strategies, each answering a different question:
(a) By representative point — cheapest, and exact for a clean coverage:
points = parcels.set_geometry(parcels.geometry.representative_point())
joined = gpd.sjoin(points[["geometry"]], wards[["ward", "geometry"]],
how="left", predicate="within")
parcels["ward"] = joined["ward"]One point-in-polygon test per parcel, and a point is inside at most one polygon of a non-overlapping set. representative_point() rather than centroid because a centroid can fall outside a concave shape.
The caveat: a point landing exactly on a shared boundary is within neither (strictly inside neither) or intersects both, depending on the predicate. Rare with real coordinates, routine with snapped data.
(b) By largest overlap — the usual intent for polygons:
import numpy as np
import shapely
def assign_by_largest_overlap(left, right, key):
joined = gpd.sjoin(left, right[[key, "geometry"]], how="left",
predicate="intersects")
matched = joined[joined["index_right"].notna()].copy()
# vectorised intersection area — no apply(axis=1)
left_geoms = left.geometry.values[matched.index.map(left.index.get_loc)]
right_geoms = right.geometry.values[matched["index_right"].astype(int)]
matched["overlap"] = shapely.area(shapely.intersection(left_geoms, right_geoms))
best = (matched.sort_values("overlap", ascending=False)
.groupby(level=0)[key].first())
out = left.copy()
out[key] = best.reindex(out.index)
return out
parcels = assign_by_largest_overlap(parcels, wards, "ward")(c) Keep every match, in a list — when losing the others is unacceptable:
joined = gpd.sjoin(parcels, wards[["ward", "geometry"]], how="left",
predicate="intersects")
lists = joined.groupby(level=0)["ward"].agg(
lambda s: sorted(v for v in s if isinstance(v, str)))
parcels["wards"] = lists
parcels["n_wards"] = lists.str.len()(d) Split the geometry — when a parcel genuinely belongs to both in proportion:
split = gpd.overlay(parcels, wards[["ward", "geometry"]], how="intersection")
split["area_m2"] = split.geometry.area
print(f"{len(parcels):,} parcels → {len(split):,} parcel-ward pieces")overlay cuts each parcel at every ward boundary, so a straddling parcel becomes two pieces with correct areas. Totals then reconcile exactly, at the cost of more rows and modified geometry.
5. Get one row per right-hand feature — the aggregate case
Counting points in polygons is the mirror image, and the errors are different:
joined = gpd.sjoin(wards, incidents, how="left", predicate="contains")
counts = joined.groupby(level=0).size() # ← wrong
counts = joined.groupby(level=0)["index_right"].count() # ← rightsize() counts rows, so a ward with no incidents gets 1 — the row where index_right is NaN. count() on the right-hand column counts non-null values, so an empty ward gets 0. This is the same distinction as SQL's COUNT(*) versus COUNT(column), and it silently inflates every zero in the result.
wards["incidents"] = joined.groupby(level=0)["index_right"].count()
print(f"{(wards['incidents'] == 0).sum()} wards with no incidents")
print(f"total counted: {wards['incidents'].sum():,} of {len(incidents):,}")Reconcile the total. If it exceeds the number of incidents, points are being counted in more than one ward — meaning the wards overlap, or the predicate is intersects and points sit on boundaries.
Code examples
Example 1: a join wrapper that enforces cardinality
import numpy as np
import shapely
import geopandas as gpd
def sjoin_one_to_one(left, right, right_cols, *, predicate="intersects",
strategy="largest_overlap", verbose=True):
"""Spatial join guaranteeing exactly one output row per left row."""
if strategy == "representative_point":
probe = left.set_geometry(left.geometry.representative_point())
joined = gpd.sjoin(probe[[probe.geometry.name]],
right[[*right_cols, right.geometry.name]],
how="left", predicate="within")
joined = joined[~joined.index.duplicated(keep="first")]
else:
joined = gpd.sjoin(left, right[[*right_cols, right.geometry.name]],
how="left", predicate=predicate)
n_raw = len(joined)
matched = joined[joined["index_right"].notna()]
if strategy == "largest_overlap" and len(matched):
li = np.array([left.index.get_loc(i) for i in matched.index])
ri = matched["index_right"].to_numpy().astype(int)
overlap = shapely.area(shapely.intersection(
left.geometry.values[li], right.geometry.values[ri]))
matched = matched.assign(_overlap=overlap)
keep = matched.sort_values("_overlap", ascending=False) \
.groupby(level=0).head(1)
elif strategy == "nearest_centroid" and len(matched):
li = np.array([left.index.get_loc(i) for i in matched.index])
ri = matched["index_right"].to_numpy().astype(int)
d = shapely.distance(shapely.centroid(left.geometry.values[li]),
shapely.centroid(right.geometry.values[ri]))
matched = matched.assign(_dist=d)
keep = matched.sort_values("_dist").groupby(level=0).head(1)
else:
keep = matched.groupby(level=0).head(1)
joined = keep.reindex(left.index)
if verbose:
print(f" raw join produced {n_raw:,} rows for {len(left):,} inputs "
f"— reduced to one each by {strategy}")
out = left.copy()
for col in right_cols:
out[col] = joined[col].reindex(out.index)
if verbose:
unmatched = int(out[right_cols[0]].isna().sum())
print(f" {len(out):,} rows out (== {len(left):,} in), "
f"{unmatched:,} unmatched ({100 * unmatched / len(out):.2f}%)")
assert len(out) == len(left), "cardinality guarantee violated"
return out
parcels = sjoin_one_to_one(parcels, wards, ["ward_code", "ward_name"],
strategy="largest_overlap") raw join produced 4,318,552 rows for 4,012,884 inputs — reduced to one each by largest_overlap
4,012,884 rows out (== 4,012,884 in), 412 unmatched (0.01%)
The assert is the point. A function promising one-to-one should fail loudly if it does not deliver, rather than returning a frame that quietly breaks a total three steps later.
The three strategies answer different questions. largest_overlap is the right default for polygons. representative_point is far cheaper and correct for a clean coverage. nearest_centroid suits cases where overlap area is not meaningful — a point layer against zones, say.
reindex(left.index) restores rows that matched nothing, keeping the guarantee even for unmatched inputs.
Example 2: diagnosing an unexpected row count
import geopandas as gpd
import pandas as pd
def diagnose_cardinality(left, right, *, predicate="intersects", top=5):
"""Explain why a join produced the row count it did."""
joined = gpd.sjoin(left, right, how="left", predicate=predicate)
per_input = joined.index.value_counts()
print(f"predicate '{predicate}': {len(left):,} → {len(joined):,} rows")
print(f" matched exactly once {(per_input == 1).sum():,}")
print(f" matched more than one {(per_input > 1).sum():,}")
print(f" matched nothing {int(joined['index_right'].isna().sum()):,}")
multi = per_input[per_input > 1]
if len(multi):
print(f"\n distribution of matches per input:")
print(" ", per_input.value_counts().sort_index().to_dict())
print(f"\n worst offenders:")
for idx, n in multi.head(top).items():
geom = left.geometry.loc[idx]
matches = joined.loc[[idx]]
print(f" index {idx}: {n} matches, "
f"{geom.geom_type}, area {geom.area:,.1f}")
for _, row in matches.iterrows():
if pd.notna(row["index_right"]):
other = right.geometry.loc[row["index_right"]]
shared = geom.intersection(other)
print(f" → right {row['index_right']}: "
f"shared {shared.geom_type}, area {shared.area:,.4f}")
zero_area = 0
if len(multi):
for idx in multi.head(50).index:
for _, row in joined.loc[[idx]].iterrows():
if pd.notna(row["index_right"]):
shared = left.geometry.loc[idx].intersection(
right.geometry.loc[row["index_right"]])
if shared.area == 0:
zero_area += 1
if zero_area:
print(f"\n → {zero_area} of the sampled extra matches share zero area "
f"(boundary contact only). Try predicate='within', or "
f"predicate='overlaps' to exclude them.")
return joined
diagnose_cardinality(parcels.head(50_000), wards)predicate 'intersects': 50,000 → 53,882 rows
matched exactly once 46,204
matched more than one 3,791
matched nothing 5
distribution of matches per input:
{1: 46204, 2: 3702, 3: 84, 4: 5}
worst offenders:
index 1882: 4 matches, Polygon, area 4,118.2
→ right 41: shared Polygon, area 1,204.8801
→ right 42: shared LineString, area 0.0000
→ right 88: shared Point, area 0.0000
→ right 89: shared MultiPoint, area 0.0000
→ 41 of the sampled extra matches share zero area (boundary contact only). Try
predicate='within', or predicate='overlaps' to exclude them.
The geometry type of the shared region is the diagnosis. A Polygon intersection is a real overlap; a LineString means the two features share only an edge; a Point means they meet at a corner. Under intersects all three count as a match, and the last two are almost never what an analysis wants.
Knowing that the extra matches are zero-area contacts turns a vague "the join gave too many rows" into a specific fix: use overlaps, or filter by intersection area.
Example 3: making the totals reconcile
When counts must add up, overlay is the honest tool:
import geopandas as gpd
def apportion(left, right, right_key, value_col=None):
"""Split left features at right boundaries so totals reconcile exactly."""
pieces = gpd.overlay(left, right[[right_key, "geometry"]], how="intersection")
pieces["piece_area"] = pieces.geometry.area
original = left.geometry.area
total = pieces.groupby(pieces.index.name or "index")["piece_area"].sum() \
if pieces.index.name else None
if value_col:
# apportion an attribute by area share
share = pieces["piece_area"] / pieces.groupby(
pieces[left.index.name or "id"])["piece_area"].transform("sum")
pieces[f"{value_col}_apportioned"] = pieces[value_col] * share
print(f"{len(left):,} features → {len(pieces):,} pieces")
print(f"area in {original.sum() / 1e6:,.3f} km²")
print(f"area out {pieces['piece_area'].sum() / 1e6:,.3f} km²")
print(f"difference {abs(original.sum() - pieces['piece_area'].sum()):,.1f} m²")
return pieces
pieces = apportion(parcels, wards, "ward_code", value_col="population")
by_ward = pieces.groupby("ward_code").agg(
area_km2=("piece_area", lambda s: s.sum() / 1e6),
population=("population_apportioned", "sum"))
print(by_ward.head())4,012,884 features → 4,318,552 pieces
area in 8,412.993 km²
area out 8,412.993 km²
difference 0.4 m²
Area in equals area out to within floating-point noise. That is the property a join cannot give you: a straddling parcel becomes two pieces whose areas sum to the original, so every ward total is exactly right and the grand total matches the input.
The cost is real. Geometry is modified, row count grows, and the pieces are no longer the original features — a parcel id now appears twice. Use it when totals must reconcile; use a one-to-one join when the features must stay whole.
Apportioning an attribute by area share assumes the attribute is uniformly distributed within the feature, which is a modelling assumption worth stating. Splitting a parcel's population by area is reasonable; splitting its address is not.
Explanation
Cardinality is the property that distinguishes a spatial join from an attribute join, and it is a consequence of geometry rather than of any library's design.
An attribute join on a unique key is a function: each left row maps to at most one right row, because the key is unique. That guarantee comes from the data model, and it is why merge rarely surprises anyone.
A spatial join has no key and no uniqueness. ST_Intersects and .intersects() are relations, and a relation between two sets of geometry is many-to-many in general. A parcel can touch four wards; a road can cross twenty; a buffered point can contain hundreds of features. The join reports every true pair, which is the only defensible thing for it to do — any deduplication would require knowing which match you consider authoritative, and only your question determines that.
So the fan-out is not an error to be prevented but a fact to be resolved. The resolution is a decision, and the four strategies encode four different intentions: largest overlap says "the ward containing most of it", representative point says "the ward containing its centre", the list form says "all of them, recorded", and overlay says "split it, because it genuinely belongs to both".
The predicate is the first lever, and it is underused. within cannot match more than once against a non-overlapping coverage, because a shape entirely inside one polygon cannot also be entirely inside a disjoint one. That single fact makes within the cheapest way to get a cardinality guarantee — at the cost of dropping features that straddle boundaries, which then need handling separately. intersects, by contrast, matches on any shared point including a corner, which is why so many "why did my join duplicate rows" cases turn out to be zero-area boundary contacts.
The join type controls loss rather than duplication, and the two are independent. how="inner" drops unmatched left rows; how="left" keeps them with nulls. Neither affects fan-out. The reason to default to left is that unmatched rows are usually a finding — features outside the study area, gaps in the boundaries, bad coordinates — and an inner join makes that finding invisible by construction.
And the aggregate case has its own trap. After a left join, groupby().size() counts rows, so an unmatched group gets 1 rather than 0. groupby()["index_right"].count() counts non-null values and gets it right. This is exactly SQL's COUNT(*) versus COUNT(column) distinction, and it silently turns every zero into a one — which on a choropleth of counts is invisible and on a total is a systematic overstatement.
The practical discipline is one line: compare the row count in to the row count out, every time. More means duplication; fewer means loss; equal means neither, or both cancelling. It is the cheapest check available and it catches the error class that no map, no summary statistic and no visual inspection will reveal.
Edge cases or notes
intersectsincludes boundary contact. A shared edge or a single corner counts as a match.withingives at most one match against a non-overlapping coverage — the cheapest cardinality guarantee.sjoin_nearestcan return several rows on exact distance ties. Deduplicate on the index.groupby().size()counts rows;groupby()[col].count()counts non-nulls. After a left join, use the second.- Overlapping right-hand polygons duplicate regardless of predicate. Check with a self-join first.
overlay(how="intersection")makes totals reconcile at the cost of splitting geometry and growing the row count.representative_point()beatscentroid— a centroid can fall outside a concave polygon.- A point exactly on a shared boundary matches both neighbours under
intersectsand neither underwithin. index_rightis dropped by some GeoPandas versions after a join. Capture what you need before it goes.- A self-join needs
left.index != right.index, or every feature matches itself.
Internal links
- How to perform a spatial join in Python (GeoPandas) — the operation itself
- GeoPandas spatial join returns duplicate rows — the symptom, and its fixes
- Spatial predicates explained — choosing the relationship
- How to do a spatial join in PostGIS with SQL — the same problem in SQL
- Attribute join or spatial join? — when a key would do
- Nearest-neighbour joins explained — cardinality for distance joins
- Spatial join returns empty results in GeoPandas — the opposite failure
- How to count points in polygons with GeoPandas — the aggregate case
FAQ
Why does my spatial join return more rows than I put in?
Because features match more than one feature on the other side — a parcel straddling two wards intersects both. That is correct; the join reports every true pair.
How do I guarantee one row per input?
Join, then deduplicate on an explicit rule: largest overlap, representative point, or first match. Or use predicate="within", which cannot match twice against a non-overlapping coverage.
Why did rows disappear from my join?
how="inner" drops unmatched left rows. Use how="left" so features that match nothing survive with nulls — they are usually a finding worth seeing.
Why does a ward with no incidents show a count of 1?
groupby().size() counts rows, including the row where the right-hand side is null. Use groupby()["index_right"].count(), which counts non-null values.
Which predicate should I use?
within when the right-hand polygons form a coverage and you want at most one match. intersects when any shared point counts, accepting the fan-out. overlaps to exclude zero-area boundary contact.
How do I make my totals reconcile exactly?
gpd.overlay(left, right, how="intersection"), which splits features at boundaries so areas sum correctly. The cost is modified geometry and more rows.
What is the single most useful check?
Compare the input row count with the output row count after every join. More means duplication, fewer means loss, and both are invisible on a map.