Spatial Join Cardinality Explained: One-to-One, One-to-Many and What to Do
On this page
Problem statement
You join Leeds' 293,729 registered land titles to its 33 electoral wards and get more rows back than you put in:
joined = gpd.sjoin(parcels, wards, predicate="intersects")
print(f"{len(parcels):,} → {len(joined):,}")293,729 → 295,766
Two thousand extra rows. Nothing errored. Every row is correct — each of those titles really does intersect two wards. And every total computed from that frame is now wrong, because those titles are counted twice.
Or the opposite:
joined = gpd.sjoin(parcels, wards, how="inner", predicate="within")
print(f"{len(joined):,}")290,842
Nearly three thousand titles vanished. They were not deleted; they simply matched nothing — a title that crosses a ward boundary is within neither ward — 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
sjoin alone guarantees none of them.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()) |
The cheap alternative for a clean coverage: join by representative point instead, which matches at most one polygon.
Why rows multiply, and getting one row per input
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)
toy_wards = gpd.GeoDataFrame({"ward": ["west", "east"]}, geometry=[west, east], crs=27700)
toy_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), # inside east, its left edge on the boundary
box(30, 2, 32, 4)], # outside both
crs=27700)
joined = gpd.sjoin(toy_parcels, toy_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, six rows out. straddling intersects both wards, which is true — it does. touching lies entirely inside east, but its left edge sits on the boundary, and a shared edge is an intersection — so it matches west too, with zero area in common. 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}
report = join_report(parcels, gpd.sjoin(parcels, wards, how="left", predicate="intersects"))join
input rows 293,729
output rows 295,784 (+2,055)
distinct inputs 293,729
duplicated 2,055 ⚠ some inputs matched more than once
unmatched 18
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 title touching several wards is plausible for a long road or a railway, and one touching dozens 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):>9,} rows "
f"{j.index.duplicated().sum():>7,} duplicated "
f"{j['index_right'].isna().sum():>8,} unmatched")intersects 295,784 rows 2,055 duplicated 18 unmatched
within 293,729 rows 0 duplicated 2,887 unmatched
contains 293,729 rows 0 duplicated 293,729 unmatched
overlaps 295,783 rows 2,054 duplicated 290,860 unmatched
touches 293,729 rows 0 duplicated 293,728 unmatched
intersects duplicates the 2,055 titles that touch more than one ward; within never duplicates, because a title cannot lie wholly inside two non-overlapping wards, and pays for it with 2,887 straddling titles unmatched. The other three answer different questions — contains is backwards here, overlaps finds the straddlers and touches the boundary contacts — and none of them assigns. What each predicate asks is spatial predicates explained; for cardinality, the point is that the choice moves rows between "duplicated" and "unmatched", and no predicate makes both zero.
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())
by_point = gpd.sjoin(points[["geometry"]], wards[["ward_name", "geometry"]],
how="left", predicate="within")
print(f"{len(by_point):,} rows, {by_point['ward_name'].isna().sum():,} unmatched")293,729 rows, 402 unmatched
One point-in-polygon test per title, 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 few unmatched are titles whose point falls outside the district boundary altogether — land registered to Leeds titles that crosses into a neighbouring authority.
The caveat: a point landing exactly on a shared boundary is within 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 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[left.index.get_indexer(matched.index)]
right_geoms = right.geometry.values[right.index.get_indexer(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
by_overlap = assign_by_largest_overlap(parcels, wards, "ward_name")
differ = (by_overlap["ward_name"] != by_point["ward_name"].reindex(by_overlap.index)).sum()
print(f"largest overlap and representative point disagree for {differ:,} titles")largest overlap and representative point disagree for 482 titles
(c) Keep every match, in a list — when losing the others is unacceptable:
joined = gpd.sjoin(parcels, wards[["ward_name", "geometry"]], how="left", predicate="intersects")
lists = joined.groupby(level=0)["ward_name"].agg(lambda s: sorted(v for v in s if isinstance(v, str)))
n_wards = lists.str.len()
print(n_wards.value_counts().sort_index().to_dict()){0: 18, 1: 291707, 2: 1957, 3: 43, 4: 4}
(d) Split the geometry — when a title genuinely belongs to both in proportion:
split = gpd.overlay(parcels, wards[["ward_name", "geometry"]], how="intersection")
split["area_m2"] = split.geometry.area
print(f"{len(parcels):,} titles → {len(split):,} title-ward pieces")293,729 titles → 295,765 title-ward pieces
overlay cuts each title at every ward boundary, so a straddling title becomes several 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. Here are the crimes West Yorkshire Police recorded in Leeds in June 2026, counted per ward:
joined = gpd.sjoin(wards, incidents, how="left", predicate="contains")
wrong = joined.groupby(level=0).size() # ← counts rows
right = joined.groupby(level=0)["index_right"].count() # ← counts matches
print(f"size(): {wrong.sum():,} in total, smallest {wrong.min()}")
print(f"count(): {right.sum():,} in total, smallest {right.min()}")
print(f"crimes: {len(incidents):,}")size(): 9,224 in total, smallest 72
count(): 9,224 in total, smallest 72
crimes: 9,224
size() counts rows, so a ward with no incidents would get 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. Every Leeds ward had crimes in June, so the two agree here — which is exactly how the size() version survives testing and then miscounts the day a smaller area or a rarer category has an empty zone.
Reconcile the total against the number of points. If it exceeds it, points are being counted in more than one ward — meaning the wards overlap, or the predicate is intersects and points sit on boundaries. If it falls short, some points lie outside every ward.
Leeds parcels against wards, with cardinality enforced
The data used on this page
parcels is HM Land Registry's INSPIRE Index Polygons for Leeds — the extent of every registered freehold title. wards is the 33 Leeds electoral wards (December 2024, full resolution) from the ONS. incidents is every street-level crime West Yorkshire Police published for June 2026 inside those wards, from data.police.uk, and zones is Leeds' 488 LSOAs with their Census 2021 population, used in Example 3. Land Registry data © Crown copyright and database right 2026, Open Government Licence v3.0.
import time
import geopandas as gpd
import numpy as np
import pandas as pd
import pyogrio
import requests
# HM Land Registry's download sets a cookie and redirects to itself, so use a session
reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
f.write(reply.content)
pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"}) # the file's schema URL is unreachable
parcels = gpd.read_file("zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"])
ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
wards = gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
"?where=LAD24NM%3D%27Leeds%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
wards = (wards.set_crs(27700, allow_override=True).rename(columns={"WD24CD": "ward_code", "WD24NM": "ward_name"})
.sort_values("ward_code").reset_index(drop=True))
# June 2026 street-level crime, fetched in tiles (the API returns at most 10,000 per request)
west, south, east, north = wards.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):
for attempt in range(5): # the API answers bursts with 429
r = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
if r.status_code != 429:
break
time.sleep(2 ** attempt)
r.raise_for_status()
crimes += [{"id": c["id"], "lon": float(c["location"]["longitude"]),
"lat": float(c["location"]["latitude"])} for c in r.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes[["id"]], crs=4326,
geometry=gpd.points_from_xy(crimes["lon"], crimes["lat"])).to_crs(27700)
incidents = incidents[incidents.within(wards.union_all())].reset_index(drop=True)
# Leeds LSOAs with their Census 2021 population, for Example 3
LSOA = f"{ONS}/Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query"
features, offset = [], 0
while True:
page = requests.get(LSOA, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD",
"outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
features += page["features"]
if not page.get("properties", {}).get("exceededTransferLimit"):
break
offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700)
population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
f"&geography={','.join(zones['LSOA21CD'])}&c2021_restype_3=0&measures=20100"
"&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = zones["LSOA21CD"].map(population)
print(f"{len(parcels):,} titles, {len(wards)} wards, {len(incidents):,} crimes, "
f"{len(zones)} LSOAs with {zones['population'].sum():,} residents")293,729 titles, 33 wards, 9,224 crimes, 488 LSOAs with 811,957 residents
Example 1: a join wrapper that enforces cardinality
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()]
li = left.index.get_indexer(matched.index)
ri = right.index.get_indexer(matched["index_right"].astype(int))
if strategy == "largest_overlap" and len(matched):
overlap = shapely.area(shapely.intersection(left.geometry.values[li], right.geometry.values[ri]))
keep = matched.assign(_overlap=overlap).sort_values("_overlap", ascending=False).groupby(level=0).head(1)
elif strategy == "nearest_centroid" and len(matched):
d = shapely.distance(shapely.centroid(left.geometry.values[li]),
shapely.centroid(right.geometry.values[ri]))
keep = matched.assign(_dist=d).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
assigned = sjoin_one_to_one(parcels, wards, ["ward_code", "ward_name"], strategy="largest_overlap") raw join produced 295,784 rows for 293,729 inputs — reduced to one each by largest_overlap
293,729 rows out (== 293,729 in), 18 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
def diagnose_cardinality(left, right, *, predicate="intersects", top=2):
"""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("\n distribution of matches per input:")
print(" ", per_input.value_counts().sort_index().to_dict())
print("\n worst offenders:")
for idx, n in multi.head(top).items():
geom = left.geometry.loc[idx]
print(f" index {idx}: {n} matches, {geom.geom_type}, area {geom.area:,.1f}")
for r in joined.loc[[idx], "index_right"].dropna().astype(int):
shared = geom.intersection(right.geometry.loc[r])
print(f" → right {r}: shared {shared.geom_type}, area {shared.area:,.4f}")
kinds = {}
for idx in multi.index:
for r in joined.loc[[idx], "index_right"].dropna().astype(int):
shared = left.geometry.loc[idx].intersection(right.geometry.loc[r])
kind = "area" if shared.area > 0 else "boundary only"
kinds[kind] = kinds.get(kind, 0) + 1
print(f"\n what the multi-matched inputs share with each ward: {kinds}")
return joined
diagnosed = diagnose_cardinality(parcels.head(20_000), wards)predicate 'intersects': 20,000 → 20,171 rows
matched exactly once 19,831
matched more than one 169
matched nothing 2
distribution of matches per input:
{1: 19831, 2: 167, 3: 2}
worst offenders:
index 11480: 3 matches, Polygon, area 92,286.4
→ right 25: shared Polygon, area 29,686.7800
→ right 8: shared Polygon, area 62,581.7979
→ right 26: shared MultiPolygon, area 17.8056
index 12158: 3 matches, Polygon, area 509,882.0
→ right 25: shared MultiPolygon, area 26,519.7329
→ right 8: shared Polygon, area 483,348.1257
→ right 26: shared Polygon, area 14.1821
what the multi-matched inputs share with each ward: {'area': 339, 'boundary only': 1}
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.
On the Leeds titles the answer is the less comfortable one: almost every extra match shares real area. Titles genuinely cross ward lines — a ward boundary drawn down a road or a stream cuts through the land either side — and many of the crossings are slivers, like the 14 and 18 m² pieces above. Changing the predicate would not remove them; deciding by area would. That is what largest_overlap does, or a threshold on the share of each title that falls in the ward.
Knowing whether the extra matches are real overlaps or zero-area contacts turns a vague "the join gave too many rows" into a specific fix: filter by intersection area or take the largest overlap for real overlaps; use within or a representative point for boundary contact.
Example 3: making the totals reconcile
When counts must add up, overlay is the honest tool. Leeds' LSOAs do not nest exactly inside its wards, so a ward population built from LSOAs has to split some of them:
def apportion(left, right, left_key, right_key, value_col):
"""Split left features at right boundaries and share value_col among the pieces by area."""
pieces = gpd.overlay(left[[left_key, value_col, "geometry"]], right[[right_key, "geometry"]],
how="intersection")
pieces["piece_area"] = pieces.geometry.area
share = pieces["piece_area"] / pieces.groupby(left_key)["piece_area"].transform("sum")
pieces[f"{value_col}_apportioned"] = pieces[value_col] * share
print(f"{len(left):,} features → {len(pieces):,} pieces, "
f"{(pieces.groupby(left_key).size() > 1).sum():,} of them split")
print(f"area in {left.geometry.area.sum() / 1e6:,.3f} km² out {pieces['piece_area'].sum() / 1e6:,.3f} km²")
print(f"{value_col} in {left[value_col].sum():,} out {pieces[f'{value_col}_apportioned'].sum():,.0f}")
return pieces
pieces = apportion(zones, wards, "LSOA21CD", "ward_code", "population")
by_ward = pieces.groupby("ward_code")["population_apportioned"].sum().round()
print(by_ward.describe().round(0).to_string())488 features → 913 pieces, 298 of them split
area in 551.765 km² out 551.326 km²
population in 811,957 out 811,957
count 33.0
mean 24605.0
std 3149.0
min 19105.0
25% 22806.0
50% 23818.0
75% 26450.0
max 31268.0
The population out equals the population in: every LSOA's people are shared among its pieces, and the shares sum to one. That is the property a join cannot give you — a straddling feature becomes pieces whose values sum to the original, so every ward total adds up and the grand total matches the input. The areas need not match: the LSOA boundaries here are the ONS's generalised version and the wards are full resolution, so their edges disagree by metres all along the district boundary, and the slivers that fall outside every ward drop out of the overlay.
The cost is real. Geometry is modified, row count grows, and the pieces are no longer the original features — an LSOA code now appears more than once. 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 an LSOA's population by area is a reasonable first approximation; weighting by address points is better, and splitting its postcode is meaningless.
Cardinality comes from geometry
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 two levers that change the count are independent, and the steps measured both: the predicate moves rows between duplicated and unmatched (step 3), and the join type decides whether the unmatched are kept as rows with nulls or silently dropped. Default to a left join, because unmatched rows are usually a finding — 18 Leeds titles lie entirely outside the city's wards — and an inner join makes that finding invisible.
Nearest ties, overlapping right layers and counting rows
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.
index_rightis only called that when the right frame's index has no name. Give the index a name and GeoPandas uses the name instead, so code that readsjoined["index_right"]breaks.- 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.
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.