Attribute Join or Spatial Join? Choosing How to Combine Two Layers #

Problem statement #

You have ward boundaries and a spreadsheet of population figures. You need them together. There are two ways, and the wrong one is slower, less reliable, and occasionally wrong:

# spatial: match by where things are
joined = gpd.sjoin(wards, population_points, predicate="contains")

# attribute: match by a shared code
joined = wards.merge(population_df, on="ward_code")

Eleven seconds versus four milliseconds. And the spatial version can produce a different answer, because a point on a boundary lands in two wards while a code matches exactly once.

The reverse mistake is just as common. Someone has incident records with a ward_name column and joins on it:

joined = wards.merge(incidents, on="ward_name", how="left")
print(joined["incident_id"].isna().sum())      # 4,118 unmatched

Four thousand unmatched, because the incidents say "Ancoats & Beswick" and the boundaries say "Ancoats and Beswick". The geometry was right there and would have matched every one.

Choosing between them is not a style question. It is a question about where the relationship actually lives — in a shared key, or in the geometry.

Quick answer #

# a shared key — always prefer this
wards = wards.merge(population, on="ward_code", how="left", validate="one_to_one")

# no key — the relationship is where things are
incidents = gpd.sjoin(incidents, wards[["ward_code", "geometry"]],
                      how="left", predicate="within")

# derive the key once, reuse it many times
incidents = gpd.sjoin(incidents, wards[["ward_code", "geometry"]],
                      how="left", predicate="within")
incidents.to_parquet("incidents_with_ward.parquet")     # now merges are free

Use validate= on every merge. It turns a silent cardinality surprise into an exception.

Step-by-step solution #

1. Find out whether a usable key exists #

import pandas as pd

def find_shared_keys(left, right, *, sample=1000):
    """Columns present in both, with how well their values overlap."""
    shared = set(left.columns) & set(right.columns) - {"geometry"}
    if not shared:
        print("no shared column names")
        return []

    results = []
    for col in sorted(shared):
        lv = set(left[col].dropna().astype(str).str.strip())
        rv = set(right[col].dropna().astype(str).str.strip())
        if not lv or not rv:
            continue
        overlap = len(lv & rv)
        results.append({
            "column": col,
            "left_distinct": len(lv),
            "right_distinct": len(rv),
            "shared_values": overlap,
            "left_matched_%": round(100 * overlap / len(lv), 1),
            "left_unique": left[col].is_unique,
            "right_unique": right[col].is_unique,
        })
    df = pd.DataFrame(results).sort_values("left_matched_%", ascending=False)
    print(df.to_string(index=False))
    return df

find_shared_keys(wards, population)

    column  left_distinct  right_distinct  shared_values  left_matched_%  left_unique  right_unique
 ward_code            215             215            215           100.0         True         True
 ward_name            215             215            188            87.4         True         True

ward_code matches 100% and is unique on both sides — a perfect key, use it. ward_name matches 87.4%, which is the classic symptom of ampersands, apostrophes and case differences.

Codes beat names, always. A code is issued by an authority and does not change when someone reformats a spreadsheet.

2. Prefer the attribute join when a key exists #

merged = wards.merge(population, on="ward_code", how="left", validate="one_to_one")
print(f"{len(wards):,} → {len(merged):,}, "
      f"{merged['population'].isna().sum():,} unmatched")

215 → 215, 0 unmatched

validate= is the argument that makes a merge safe:

wards.merge(population, on="ward_code", validate="one_to_one")

MergeError: Merge keys are not unique in right dataset; not a one-to-one merge

That error is worth having. Without it, a duplicated key on the right silently multiplies your ward count, and the map still renders.

indicator=True shows exactly what matched:

merged = wards.merge(population, on="ward_code", how="outer", indicator=True)
print(merged["_merge"].value_counts().to_dict())

{'both': 215, 'left_only': 0, 'right_only': 0}

right_only rows are population figures for wards that do not exist in your boundaries — a real finding, usually a boundary revision.

3. Use a spatial join when the relationship is positional #

Some relationships have no key by nature:

# an incident has coordinates, not a ward code
incidents = gpd.sjoin(incidents, wards[["ward_code", "geometry"]],
                      how="left", predicate="within")

# which parcels does this flood zone affect?
affected = gpd.sjoin(parcels, flood_zones, predicate="intersects")

# how far to the nearest school?
properties = gpd.sjoin_nearest(properties, schools, max_distance=2_000,
                               distance_col="dist_m")

The distinction is not about which data you have but about where the relationship is stored. A ward code in an incident record is someone's earlier spatial join, written down. If it is absent, or you do not trust it, the geometry still knows.

4. Measure the cost difference #

import time

def timed(label, fn):
    t0 = time.perf_counter()
    result = fn()
    print(f"{label:<32} {time.perf_counter() - t0:>8.3f} s   {len(result):>10,} rows")
    return result

timed("merge on ward_code",
      lambda: parcels.merge(ward_attrs, on="ward_code", how="left"))
timed("sjoin within",
      lambda: gpd.sjoin(parcels, wards[["ward_code", "geometry"]],
                        how="left", predicate="within"))
timed("sjoin intersects",
      lambda: gpd.sjoin(parcels, wards[["ward_code", "geometry"]],
                        how="left", predicate="intersects"))

merge on ward_code                  0.412 s    4,012,884 rows
sjoin within                       11.204 s    3,998,204 rows
sjoin intersects                   14.882 s    4,318,552 rows

Twenty-seven times slower, and note the row counts differ across all three. The merge is exact and one-to-one; within drops straddling parcels; intersects duplicates them.

The cost gap comes from what each does. A merge builds a hash table on the key and looks each row up — O(n). A spatial join builds an R-tree, queries it per feature, then runs an exact geometry predicate on every candidate. Even with the index, it is fundamentally more work.

5. Use both: derive the key once, then merge #

The best pattern for repeated work is a spatial join once, stored, and attribute joins thereafter:

# expensive, run once when the data arrives
incidents = gpd.sjoin(incidents, wards[["ward_code", "geometry"]],
                      how="left", predicate="within")
incidents = incidents.drop(columns="index_right")
incidents.to_parquet("incidents_with_ward.parquet")

# cheap, run whenever
incidents = gpd.read_parquet("incidents_with_ward.parquet")
by_ward = incidents.groupby("ward_code").size()
wards = wards.merge(by_ward.rename("incidents"), on="ward_code", how="left")

This is exactly what a ward_code column in a supplier's incident extract is — their spatial join, precomputed. Which is also why you should check it rather than trust it:

def verify_key_against_geometry(gdf, areas, key, *, sample=2000, seed=0):
    """Does the recorded key agree with where the feature actually is?"""
    import numpy as np
    rng = np.random.default_rng(seed)
    idx = rng.choice(len(gdf), min(sample, len(gdf)), replace=False)
    probe = gdf.iloc[idx][[key, "geometry"]].copy()

    checked = gpd.sjoin(probe, areas[[key, "geometry"]],
                        how="left", predicate="within", rsuffix="_geom")
    checked = checked[~checked.index.duplicated(keep="first")]

    recorded = checked[key].astype(str)
    actual = checked[f"{key}_geom"].astype(str)
    disagree = (recorded != actual) & actual.notna() & (actual != "nan")

    print(f"{len(checked):,} sampled")
    print(f"  agree      {(~disagree).sum():,}")
    print(f"  disagree   {disagree.sum():,} ({100 * disagree.mean():.2f}%)")
    print(f"  no match   {(actual == 'nan').sum():,}")
    if disagree.any():
        print(checked.loc[disagree, [key, f"{key}_geom"]].head(5).to_string())
    return disagree.sum()

verify_key_against_geometry(incidents, wards, "ward_code")

2,000 sampled
  agree      1,982
  disagree      12 (0.60%)
  no match       6

Twelve disagreements out of 2,000 means the supplier's codes were computed against a different boundary vintage, or with a different rule. Six with no match are outside every ward. Both are findings; neither is visible without checking.

6. Handle the case where the key is nearly right #

Names almost match, and fixing them is often better than falling back to geometry:

def normalise_name(s):
    return (s.astype(str).str.lower()
             .str.replace(r"\s*&\s*", " and ", regex=True)
             .str.replace(r"[’']", "", regex=True)
             .str.replace(r"[^a-z0-9 ]", " ", regex=True)
             .str.replace(r"\s+", " ", regex=True)
             .str.strip())

wards["_key"] = normalise_name(wards["ward_name"])
incidents["_key"] = normalise_name(incidents["ward_name"])

before = wards["ward_name"].isin(incidents["ward_name"]).sum()
after = wards["_key"].isin(incidents["_key"]).sum()
print(f"exact match: {before}/{len(wards)} → normalised: {after}/{len(wards)}")

exact match: 188/215 → normalised: 213/215

Two still fail, which is where a spatial join earns its place — or fuzzy matching, if there is no geometry to fall back on.

Normalising is cheap and repeatable, and it is worth doing before concluding that a key is unusable.

Code examples #

Example 1: a joiner that picks the method and reports #

import time
import geopandas as gpd
import pandas as pd

def join_layers(left, right, right_cols, *, key=None, predicate="within",
                how="left", validate="many_to_one", verbose=True):
    """Attribute-join on `key` if it is usable, otherwise spatial-join."""
    t0 = time.perf_counter()

    usable_key = None
    if key and key in left.columns and key in right.columns:
        lv = set(left[key].dropna().astype(str).str.strip())
        rv = set(right[key].dropna().astype(str).str.strip())
        coverage = len(lv & rv) / max(len(lv), 1)
        if verbose:
            print(f"key '{key}': {100 * coverage:.1f}% of left values found on the right"
                  f"{', unique on the right' if right[key].is_unique else ''}")
        if coverage >= 0.98 and right[key].is_unique:
            usable_key = key
        elif verbose:
            print(f"  → not usable as a key; falling back to a spatial join")

    if usable_key:
        out = left.merge(right[[usable_key, *right_cols]], on=usable_key,
                         how=how, validate=validate)
        method = f"attribute join on '{usable_key}'"
    else:
        joined = gpd.sjoin(left, right[[*right_cols, right.geometry.name]],
                           how=how, predicate=predicate)
        dup = int(joined.index.duplicated().sum())
        joined = joined[~joined.index.duplicated(keep="first")]
        out = joined.drop(columns=[c for c in ("index_right",) if c in joined])
        method = f"spatial join, predicate='{predicate}'"
        if verbose and dup:
            print(f"  {dup:,} features matched more than one — kept the first")

    if verbose:
        unmatched = int(out[right_cols[0]].isna().sum())
        print(f"{method}")
        print(f"  {len(left):,} in → {len(out):,} out "
              f"({len(out) - len(left):+,})")
        print(f"  {unmatched:,} unmatched ({100 * unmatched / max(len(out), 1):.2f}%)")
        print(f"  {time.perf_counter() - t0:.3f} s")
    return out

wards = join_layers(wards, population, ["population", "households"],
                    key="ward_code")
incidents = join_layers(incidents, wards, ["ward_code", "ward_name"],
                        key="ward_code", predicate="within")

key 'ward_code': 100.0% of left values found on the right, unique on the right
attribute join on 'ward_code'
  215 in → 215 out (+0)
  0 unmatched (0.00%)
  0.004 s

key 'ward_code': 0.0% of left values found on the right
  → not usable as a key; falling back to a spatial join
spatial join, predicate='within'
  8,436 in → 8,436 out (+0)
  112 unmatched (1.33%)
  11.204 s

The 98% coverage threshold is a judgement, and it is deliberately high: a key matching 90% of rows is not a key, it is a partially-correct lookup, and falling back to geometry gives a complete answer instead of a mostly-complete one.

Requiring uniqueness on the right is the other half. A non-unique right key would fan the merge out, which is exactly the surprise validate= exists to prevent — so the function refuses that path rather than raising later.

Printing the timing makes the cost visible: 0.004 s against 11.2 s is the argument for deriving the key once and storing it.

Example 2: deriving and caching a spatial key #

from pathlib import Path
import hashlib
import geopandas as gpd

def spatial_key(points, areas, key_col, *, cache_dir="cache",
                predicate="within", force=False):
    """Assign an area key to each point, caching the expensive join."""
    cache_dir = Path(cache_dir)
    cache_dir.mkdir(exist_ok=True)

    fingerprint = hashlib.sha256(
        f"{len(points)}|{tuple(round(v, 3) for v in points.total_bounds)}|"
        f"{len(areas)}|{tuple(round(v, 3) for v in areas.total_bounds)}|"
        f"{key_col}|{predicate}".encode()).hexdigest()[:16]
    cached = cache_dir / f"spatial_key_{fingerprint}.parquet"

    if cached.exists() and not force:
        keys = gpd.pd.read_parquet(cached)
        print(f"  cache hit: {cached.name}")
        out = points.copy()
        out[key_col] = keys[key_col].reindex(out.index).values
        return out

    areas_p = areas.to_crs(points.crs) if areas.crs != points.crs else areas
    joined = gpd.sjoin(points[[points.geometry.name]],
                       areas_p[[key_col, areas_p.geometry.name]],
                       how="left", predicate=predicate)
    dup = int(joined.index.duplicated().sum())
    joined = joined[~joined.index.duplicated(keep="first")].reindex(points.index)

    out = points.copy()
    out[key_col] = joined[key_col].values
    gpd.pd.DataFrame({key_col: out[key_col]}).to_parquet(cached)

    print(f"  computed and cached: {cached.name}")
    print(f"  {dup:,} boundary duplicates resolved, "
          f"{int(out[key_col].isna().sum()):,} unmatched")
    return out

incidents = spatial_key(incidents, wards, "ward_code")
incidents = spatial_key(incidents, wards, "ward_code")     # second call is instant

  computed and cached: spatial_key_8f2c41a9b0e14d7c.parquet
  187 boundary duplicates resolved, 112 unmatched
  cache hit: spatial_key_8f2c41a9b0e14d7c.parquet

The fingerprint covers row counts, bounds, the key column and the predicate, so a change to either layer invalidates the cache. It is not a content hash — two different layers with identical counts and bounds would collide — but it catches the realistic case of data being updated, and it costs microseconds rather than hashing gigabytes.

Caching only the key column keeps the file tiny: 8,436 strings instead of 8,436 geometries. The join it replaces took eleven seconds.

Example 3: reconciling two layers that should match and do not #

import pandas as pd
import geopandas as gpd

def reconcile(left, right, key, *, areas=None, sample=500):
    """Why do these two layers not join cleanly?"""
    lv = left[key].dropna().astype(str).str.strip()
    rv = right[key].dropna().astype(str).str.strip()
    ls, rs = set(lv), set(rv)

    print(f"key '{key}'")
    print(f"  left  {len(lv):,} rows, {len(ls):,} distinct, "
          f"{'unique' if left[key].is_unique else 'NOT unique'}")
    print(f"  right {len(rv):,} rows, {len(rs):,} distinct, "
          f"{'unique' if right[key].is_unique else 'NOT unique'}")
    print(f"  in both      {len(ls & rs):,}")
    print(f"  left only    {len(ls - rs):,}  {sorted(ls - rs)[:4]}")
    print(f"  right only   {len(rs - ls):,}  {sorted(rs - ls)[:4]}")

    # near-misses that normalisation would fix
    def norm(s):
        return (s.lower().replace("&", "and").replace("'", "")
                 .replace("-", " ").replace("  ", " ").strip())
    ln = {norm(v): v for v in ls - rs}
    rn = {norm(v): v for v in rs - ls}
    fixable = set(ln) & set(rn)
    if fixable:
        print(f"\n  {len(fixable):,} would match after normalising:")
        for k in sorted(fixable)[:5]:
            print(f"    {ln[k]!r} ↔ {rn[k]!r}")

    # would geometry do better?
    if areas is not None and isinstance(left, gpd.GeoDataFrame):
        probe = left.head(sample)
        by_geom = gpd.sjoin(probe[[probe.geometry.name]],
                            areas[[key, areas.geometry.name]],
                            how="left", predicate="within")
        by_geom = by_geom[~by_geom.index.duplicated(keep="first")]
        geom_rate = 100 * by_geom[key].notna().mean()
        key_rate = 100 * probe[key].astype(str).isin(rs).mean()
        print(f"\n  on {len(probe):,} sampled rows:")
        print(f"    key match      {key_rate:.1f}%")
        print(f"    geometry match {geom_rate:.1f}%")
        if geom_rate > key_rate + 5:
            print(f"    → the geometry is more reliable; use a spatial join")
    return ls, rs

reconcile(incidents, wards, "ward_name", areas=wards)

key 'ward_name'
  left  8,436 rows, 214 distinct, NOT unique
  right 215 rows, 215 distinct, unique
  in both      188
  left only    26  ["Ancoats & Beswick", "Cheetham's", ...]
  right only   27  ["Ancoats and Beswick", "Cheethams", ...]

  26 would match after normalising:
    "Ancoats & Beswick" ↔ "Ancoats and Beswick"
    "Cheetham's" ↔ "Cheethams"

  on 500 sampled rows:
    key match      87.2%
    geometry match 98.8%
    → the geometry is more reliable; use a spatial join

Two conclusions from one function. Twenty-six of the twenty-seven mismatches are ampersands and apostrophes, so normalising recovers nearly all of them. And the geometry matches 98.8% against the key's 87.2%, which settles the method question directly.

The NOT unique on the left is expected — many incidents per ward — and it means validate="many_to_one" is the right assertion, not "one_to_one".

Explanation #

The choice between an attribute join and a spatial join is really a question about where the relationship is recorded.

An attribute join reads a relationship someone has already written down. A ward_code on an incident record means somebody decided that incident belongs to that ward — a decision made at some point, by some process, using some boundary vintage. The join looks it up, exactly and cheaply.

A spatial join computes the relationship from geometry. Nobody wrote it down; the join derives it from where the features are. That is more expensive and it is authoritative, because geometry does not have typos, does not drift between versions, and does not depend on someone's spreadsheet formatting.

So the trade is cost against authority, and it usually resolves in favour of the key when a good key exists. A hash-table lookup is O(n); a spatial join builds an index, queries it per feature and runs exact geometry predicates on the candidates. Twenty-seven times slower is typical, and the gap widens with size.

But the key is only worth using if it is a good key, and the test is not whether a column exists but whether it matches. A code issued by a statistical agency and present on both sides at 100% coverage is a good key. A name column matching 87% is not a key at all — it is a lookup that will silently drop an eighth of your data, and the eighth it drops will not be random. It will be the wards with ampersands, which correlate with nothing except punctuation, which means the loss looks like noise.

Names are unreliable in specific, predictable ways. Ampersands versus "and", apostrophes present or absent, hyphens versus spaces, capitalisation, trailing whitespace, and genuinely different naming conventions between organisations. Normalising fixes most of it and is worth trying before abandoning the key. What it cannot fix is a genuine disagreement about names, which is where geometry wins outright.

The hybrid pattern is usually best, and it is worth naming: do the spatial join once, store the derived key, and merge thereafter. That is exactly what a supplier's ward_code column is — their spatial join, cached. Recognising that reframes the question: you are not choosing between geometry and a key, you are choosing between your spatial join and someone else's, and the reason to prefer your own is that you know which boundaries and which rule it used.

Which is also why a supplied key deserves verification rather than trust. Checking a sample against the geometry takes seconds and answers whether the two agree. When they do, use the key and enjoy the speed. When they disagree by half a percent, you have learned that the supplier used a different boundary vintage — and that is a finding worth having before it shows up as an unexplained discrepancy in a report.

Edge cases or notes #

• validate= on every merge. "one_to_one", "many_to_one" — it converts a silent fan-out into an exception.
• indicator=True adds a _merge column showing left_only, right_only or both.
• Codes beat names. A code is issued and stable; a name is formatted by whoever last touched the file.
• Normalise before giving up on a name key — ampersands, apostrophes, case and whitespace explain most mismatches.
• A merge on a non-unique right key multiplies rows. That is what validate catches.
• Whitespace is invisible and fatal. .str.strip() on both sides before comparing.
• Type mismatches fail silently: "01" does not match 1. Cast both sides.
• A spatial join can duplicate rows; a merge on a unique key cannot. See cardinality.
• Cache the derived key, not the whole joined frame — it is a fraction of the size.
• A supplied key is someone else's spatial join. Verify it against a sample of your geometry.

Internal links #

• How to join attribute data to a GeoDataFrame in Python — the attribute join in practice
• How to perform a spatial join in Python (GeoPandas) — the spatial join in practice
• Spatial join cardinality explained — why a spatial join changes the row count
• GeoPandas merge returns NaN or no matches — diagnosing a failed attribute join
• How to fuzzy-match place names when joining spatial data — when normalising is not enough
• Nearest-neighbour joins explained — the third kind of join
• How to standardise attribute values against a controlled vocabulary — making a name key usable
• Spatial predicates explained — choosing the spatial relationship

FAQ #

Which join should I use? #

An attribute join whenever both layers share a reliable code — it is exact, one-to-one, and around twenty-seven times faster. A spatial join when the relationship exists only in the geometry.

How do I know if a key is good enough? #

Check the overlap of distinct values and whether the key is unique on the right. Above about 98% coverage with uniqueness, use it. Below that it is a partial lookup, not a key.

Why does my merge on names lose rows? #

Ampersands, apostrophes, hyphens, case and whitespace. Normalise both sides before joining — it typically recovers almost all of the mismatches.

What does validate= do? #

It asserts the cardinality you expect and raises if the data disagrees — "one_to_one", "many_to_one" and so on. Without it, a duplicated key silently multiplies rows.

Should I trust a supplied ward code? #

Verify it. That column is someone else's spatial join, computed against a boundary vintage you cannot see. Checking a sample against your geometry takes seconds.

Can I use both methods? #

That is usually the best pattern: spatial-join once to derive the key, store it, and merge thereafter. You pay the expensive join once instead of on every run.

Why is a spatial join so much slower? #

A merge is a hash lookup per row. A spatial join builds an R-tree, queries it per feature and runs an exact geometry predicate on every candidate — fundamentally more work, even with the index.