How to Enrich Many Files Against One Reference Layer

Problem statement

Nine hundred parcel files need a ward name attached. The wards live in one small layer.

for path in sorted(Path("data").rglob("*.gpkg")):
    parcels = gpd.read_file(path)
    wards = gpd.read_file("reference/wards.gpkg")      # ← re-read 900 times
    out = gpd.sjoin(parcels, wards, predicate="within")
    out.to_file(dst / path.name, driver="GPKG")

The wards layer is read 900 times, its spatial index is rebuilt 900 times, and in a parallel run each of sixteen workers holds its own copy. The job spends more time loading the reference than joining against it.

Then the results are wrong in three quiet ways: parcels on a ward boundary produce duplicate rows, parcels outside every ward vanish because the join defaults to inner, and any file in a different CRS silently matches nothing.

Enriching many files against one reference is the most common batch join there is, and it has four decisions that decide whether it is fast and correct.

Quick answer

Load the reference once, reuse it, join left, and check the match rate:

import geopandas as gpd
from functools import lru_cache
from pathlib import Path

TARGET_CRS = 27700

@lru_cache(maxsize=1)
def reference(path: str = "reference/wards.gpkg"):
    """Loaded once per process, index built once."""
    ref = gpd.read_file(path).to_crs(TARGET_CRS)
    ref = ref[["ward_name", "geometry"]]        # only what the join needs
    ref.sindex                                   # build the index now, not per call
    return ref

def enrich(path: Path, dst: Path):
    parcels = gpd.read_file(path)
    if parcels.crs is None:
        raise ValueError("no CRS")
    parcels = parcels.to_crs(TARGET_CRS)

    joined = gpd.sjoin(parcels, reference(), predicate="within", how="left")

    # a boundary parcel can match two wards — resolve deterministically
    joined = (joined.sort_values(["index_right"])
                    .loc[lambda d: ~d.index.duplicated(keep="first")]
                    .drop(columns="index_right"))

    matched = joined["ward_name"].notna().sum()
    joined.to_file(dst / path.name, driver="GPKG")
    return {"rows": len(joined), "matched": int(matched),
            "rate": round(matched / max(len(joined), 1), 4)}

Decision	Wrong choice	Right choice
loading	read the reference per file	load once per process, cache it
CRS	assume they match	reproject both to one target, explicitly
join type	`how="inner"` (the default)	`how="left"`, so non-matches stay visible
duplicates	ignore them	resolve deterministically, and count them

before:  900 files, 41 min, reference read 900×
after :  900 files,  6 min, reference read 1× per worker

Where the time actually goes

The join is cheap. Reading and indexing the reference 900 times is not.

Step-by-step solution

The match rate at the end is what turns a silent failure into a visible one.

1. Load the reference once — and once per worker

from functools import lru_cache

@lru_cache(maxsize=1)
def reference(path="reference/wards.gpkg"):
    ref = gpd.read_file(path).to_crs(TARGET_CRS)
    ref.sindex                     # force the index build here
    return ref

lru_cache gives one copy per process. In a ProcessPoolExecutor each worker builds its own — sixteen copies of a small reference layer, which is fine, and vastly better than 900.

Touching .sindex matters. GeoPandas builds the spatial index lazily on first use, so without this the index is built inside the first join and the cost lands unpredictably.

For a large reference in a parallel run, sixteen copies may be too much memory. Two options:

# a) initialise once per worker at pool start
def init_worker(path):
    global REF
    REF = gpd.read_file(path).to_crs(TARGET_CRS)
    REF.sindex

with ProcessPoolExecutor(8, initializer=init_worker, initargs=("reference/wards.gpkg",)) as pool:
    ...

# b) push the join into a database and skip the copies entirely
#    see: spatial SQL queries from Python

2. Trim the reference to what the join needs

ref = gpd.read_file(path)[["ward_name", "ward_code", "geometry"]]

A reference with 60 attribute columns brings all 60 into every output file. That inflates every result, collides with existing column names, and makes the join measurably slower. Select before joining, always.

Simplifying the reference geometry is tempting and usually wrong — it moves boundaries, so features near an edge join to a different ward. Only do it if the tolerance is far below the accuracy you need, and say so in a comment.

3. Force both sides into one CRS, and fail if you cannot

def conform(gdf, crs, *, label=""):
    if gdf.crs is None:
        raise ValueError(f"{label}: no CRS — refusing to guess")
    return gdf if gdf.crs.to_epsg() == crs else gdf.to_crs(crs)

A join between mismatched CRS is the classic silent failure: GeoPandas warns once, the predicate is evaluated on incomparable coordinates, and the result is zero matches. Zero matches looks exactly like "this file covers a region with no wards", so nobody investigates. See spatial join returns empty results.

Reproject the files, not the reference — the reference is loaded once and should stay in the target CRS.

4. Choose the predicate deliberately

gpd.sjoin(parcels, wards, predicate="within")     # strictly inside
gpd.sjoin(parcels, wards, predicate="covers")     # inside or exactly on the edge
gpd.sjoin(parcels, wards, predicate="intersects") # any overlap at all

For polygons-in-polygons, within requires the parcel to be entirely inside one ward — so a parcel straddling a boundary matches nothing at all. That is usually not what people mean.

The robust choice for polygon enrichment is to join on a representative point:

probes = parcels.copy()
probes["geometry"] = parcels.geometry.representative_point()
tags = gpd.sjoin(probes, wards, predicate="within", how="left")[["ward_name"]]
parcels = parcels.join(tags)

One point is inside exactly one ward, so this gives one row per parcel by construction — no duplicates, no unmatched straddlers. It is also faster, because point-in-polygon is cheaper than polygon-in-polygon. See spatial predicates explained.

5. Join left, and resolve duplicates on purpose

joined = gpd.sjoin(parcels, wards, predicate="within", how="left")

how="inner" — the default — drops every parcel that matched nothing. The file gets smaller, the job reports success, and the missing parcels are invisible. how="left" keeps them with a null ward, which is both honest and countable.

Then deal with multi-matches:

dupes = joined.index.duplicated(keep=False)
if dupes.any():
    print(f"  {dupes.sum()} rows from {joined.index[dupes].nunique()} multi-matched features")

joined = (joined.sort_values(["index_right"])          # deterministic, not file order
                .loc[lambda d: ~d.index.duplicated(keep="first")])

Or keep both matches, collapsed into one row:

tags = (joined.groupby(joined.index)["ward_name"]
              .agg(lambda s: ";".join(sorted(s.dropna().unique()))))

6. Check the match rate — this is the real test

def match_report(joined, key="ward_name"):
    total = len(joined)
    matched = int(joined[key].notna().sum())
    return {"rows": total, "matched": matched,
            "rate": round(matched / max(total, 1), 4)}

rate = result["rate"]
if rate < 0.95:
    print(f"  ⚠ {path.name}: only {rate:.1%} matched")

A file at 100% and a file at 0% both look fine in isolation; a fleet of files where one is at 3% is obviously a CRS problem or an out-of-area delivery. Aggregating the rate across the batch is what surfaces it:

import pandas as pd
summary = pd.DataFrame(results)
print(summary["rate"].describe())
suspect = summary[summary["rate"] < 0.9]
print(f"{len(suspect)} files below 90%:")
print(suspect.nsmallest(5, "rate")[["file", "rows", "matched", "rate"]])

Code examples

Example 1: the complete enrichment job

import geopandas as gpd
import pandas as pd
from functools import lru_cache
from pathlib import Path
from concurrent.futures import ProcessPoolExecutor, as_completed

TARGET_CRS = 27700
REF_PATH = "reference/wards.gpkg"
REF_COLS = ["ward_name", "ward_code"]

@lru_cache(maxsize=1)
def reference():
    ref = gpd.read_file(REF_PATH)
    ref = conform(ref, TARGET_CRS, label="reference")[REF_COLS + ["geometry"]]
    ref.sindex
    return ref

def enrich_one(path: Path, src_root: Path, dst_root: Path) -> dict:
    gdf = gpd.read_file(path)
    if gdf.empty:
        return {"file": str(path), "status": "empty", "rows": 0}
    gdf = conform(gdf, TARGET_CRS, label=path.name)

    probes = gdf.copy()
    probes["geometry"] = gdf.geometry.representative_point()
    tags = gpd.sjoin(probes, reference(), predicate="within", how="left")

    # a probe on a shared edge can still match twice — collapse deterministically
    tags = tags.sort_values("index_right").loc[lambda d: ~d.index.duplicated(keep="first")]
    out = gdf.join(tags[REF_COLS])

    dst = dst_root / path.relative_to(src_root)
    dst.parent.mkdir(parents=True, exist_ok=True)
    tmp = dst.with_suffix(".tmp.gpkg")
    out.to_file(tmp, driver="GPKG")
    tmp.replace(dst)

    matched = int(out[REF_COLS[0]].notna().sum())
    return {"file": str(path), "status": "ok", "rows": len(out),
            "matched": matched, "rate": round(matched / len(out), 4)}

def run(src: Path, dst: Path, workers=8):
    files = sorted(src.rglob("*.gpkg"))
    results = []
    with ProcessPoolExecutor(workers) as pool:
        futures = {pool.submit(enrich_one, p, src, dst): p for p in files}
        for fut in as_completed(futures):
            path = futures[fut]
            try:
                results.append(fut.result())
            except Exception as exc:
                results.append({"file": str(path), "status": "failed",
                                "error": f"{type(exc).__name__}: {exc}"})

    df = pd.DataFrame(results)
    ok = df[df["status"] == "ok"]
    print(f"{len(ok)}/{len(files)} ok · {ok['rows'].sum():,} rows · "
          f"overall match {ok['matched'].sum()/max(ok['rows'].sum(),1):.1%}")
    low = ok[ok["rate"] < 0.9]
    if len(low):
        print(f"⚠ {len(low)} files below 90% match:")
        print(low.nsmallest(5, "rate")[["file", "rows", "matched", "rate"]].to_string(index=False))
    return df

Example 2: enriching from an attribute join instead of a spatial one

Not every reference join is spatial. When the files carry a code, a plain merge is far cheaper — and has its own duplicate trap.

def enrich_by_code(gdf, lookup: pd.DataFrame, on="ward_code"):
    assert lookup[on].is_unique, (
        f"lookup has {lookup[on].duplicated().sum()} duplicate {on} values — "
        f"the merge would multiply rows"
    )
    before = len(gdf)
    out = gdf.merge(lookup, on=on, how="left", validate="many_to_one")
    assert len(out) == before, f"merge changed row count: {before} → {len(out)}"
    return out

validate="many_to_one" makes pandas raise if the lookup is not unique, which is the single most valuable argument in merge and almost nobody passes it. Without it a duplicated lookup key silently multiplies every matching row. See my cleaned layer has more rows.

Example 3: when the reference is too big to copy per worker

# push the join into PostGIS: the reference stays put, the index is persistent
from sqlalchemy import create_engine

def enrich_via_postgis(path, engine):
    gdf = conform(gpd.read_file(path), TARGET_CRS, label=path.name)
    gdf["_row"] = range(len(gdf))
    gdf[["_row", "geometry"]].to_postgis("tmp_probe", engine, if_exists="replace",
                                         index=False)
    tags = pd.read_sql("""
        SELECT p._row, w.ward_name, w.ward_code
        FROM tmp_probe p
        LEFT JOIN wards w ON ST_Within(ST_PointOnSurface(p.geom), w.geom)
    """, engine)
    return gdf.merge(tags, on="_row", how="left").drop(columns="_row")

The reference is indexed once, in the database, and no worker holds a copy. Worth it when the reference is hundreds of megabytes or when several jobs share it. See how to run spatial SQL queries from Python.

Explanation

Triage rows pairing each silent enrichment failure with its cause and detection. — All three produce a file that opens fine and is wrong.

The performance story is simple and worth stating plainly: the join is not the expensive part. Reading a reference layer and building its R-tree costs a fixed amount, and doing it once per file multiplies that fixed cost by the number of files. With 900 files and a two-second load, that is thirty minutes of pure repetition before any joining happens.

lru_cache at module level solves it because processes are the caching unit. Each worker pays the cost once and reuses the loaded, indexed frame for every file it handles. The same reasoning applies to database connections, model objects and compiled regexes — anything expensive to build and cheap to reuse.

The correctness story has three parts, all silent.

CRS mismatch produces zero matches, which is indistinguishable from "no wards here". The only defence is to conform explicitly and refuse to proceed on a missing CRS, because a guess that happens to be wrong produces plausible output.

Inner join drops unmatched features. The output file is smaller, everything reports success, and the parcels that fell outside every ward are simply gone. how="left" converts that deletion into a null, which a match-rate check can see.

Multi-matching duplicates rows. A parcel exactly on a ward boundary is within both under some predicates, so the output has more rows than the input and every subsequent total is inflated. Joining on a representative point avoids the situation entirely rather than cleaning up after it — one point is in one polygon, so the result is one row per feature by construction.

The match rate is what ties this together operationally. Any single file can legitimately be at 100% or 0%. Across nine hundred files, the distribution is informative: a tight cluster near 100% with one file at 3% is a CRS problem or a delivery from the wrong region, and it is visible in one line of describe(). That is the difference between a batch job that produces output and one that can be trusted.

Edge cases or notes

representative_point() is guaranteed inside the polygon; centroid is not. Use the former for tagging.
sjoin adds an index_right column which must be dropped before writing, or it appears in the output file.
Column-name collisions get _left/_right suffixes silently. Trim the reference before joining.
validate="many_to_one" on merge is the cheapest protection against a duplicated lookup key.
Simplifying the reference moves boundaries and changes which features match near an edge.
A reference in a geographic CRS makes within unreliable near the antimeridian — conform to a projected CRS.
lru_cache keys on the arguments, so calling reference() with different paths caches both; maxsize=1 evicts, which defeats the point if you alternate.
Workers inherit module state on fork but not on spawn (Windows, macOS default). Use initializer= if the load must happen once per worker regardless of platform.

Internal links

How to perform a spatial join in Python (GeoPandas) — the join itself
Spatial predicates explained — why within drops straddling features
Spatial join returns empty results in GeoPandas — the CRS failure this page guards against
GeoPandas spatial join returns duplicate rows — the multi-match case
The anatomy of a batch job: discover, apply, report — where the match-rate report belongs
How to speed up batch GIS jobs with parallel processing — one reference per worker
How to run spatial SQL queries from Python with PostGIS — pushing the join into a database
My cleaned layer has more rows than the original — what an unvalidated merge does

FAQ

Why is the job so slow when the join itself is fast?

Because the reference is being read and indexed once per file. Load it once per process with lru_cache and touch .sindex to build the index up front.

Should I use `within` or a point-in-polygon join?

For tagging polygons with a containing region, join on representative_point(). It gives exactly one row per feature and is faster than polygon-in-polygon.

Why does one file match nothing?

Almost always a CRS mismatch, or a delivery covering a region the reference does not include. Conform both sides explicitly and check the match rate per file.

Inner or left join?

Left. Inner silently deletes features that matched nothing, which is invisible in the output; left keeps them with a null you can count.

How do I stop the reference's columns polluting my output?

Select the columns you need before joining. A reference with 60 attributes brings all 60 into every output file and may collide with existing names.

What if the reference is too large to copy per worker?

Use a pool initializer to load it once per worker, or push the join into PostGIS so the reference and its index stay in one place.

What is the single most useful check?

The match rate per file, aggregated across the batch. One file at 3% among 899 at 100% is obvious in a summary and invisible file by file.

How to Enrich Many Files Against One Reference Layer #

Problem statement #

Quick answer #

Where the time actually goes #

Step-by-step solution #

1. Load the reference once — and once per worker #

2. Trim the reference to what the join needs #

3. Force both sides into one CRS, and fail if you cannot #

4. Choose the predicate deliberately #

5. Join left, and resolve duplicates on purpose #

6. Check the match rate — this is the real test #

Code examples #

Example 1: the complete enrichment job #

Example 2: enriching from an attribute join instead of a spatial one #

Example 3: when the reference is too big to copy per worker #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

Why is the job so slow when the join itself is fast? #

Should I use within or a point-in-polygon join? #

Why does one file match nothing? #

Inner or left join? #

How do I stop the reference's columns polluting my output? #

What if the reference is too large to copy per worker? #

What is the single most useful check? #