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.

import time
import geopandas as gpd
from pathlib import Path

dst = Path("data/naive")
dst.mkdir(parents=True, exist_ok=True)
start = time.perf_counter()
for path in sorted(Path("data/parcels").rglob("*.gpkg")):
    parcels = gpd.read_file(path)
    wards = gpd.read_file("reference/wards.gpkg")      # ← re-read for every file
    out = gpd.sjoin(parcels, wards, predicate="within")
    out.to_file(dst / path.name, driver="GPKG")
print(f"{len(list(dst.glob('*.gpkg')))} files in {time.perf_counter() - start:.1f} s")

start = time.perf_counter()
for _ in range(52):
    gpd.read_file("reference/wards.gpkg").sindex
print(f"reading and indexing the 32 wards 52 times: {time.perf_counter() - start:.2f} s")

52 files in 1.3 s
reading and indexing the 32 wards 52 times: 0.11 s

The wards layer is read once per file — 52 times here — and its spatial index rebuilt each time. With 32 generalised wards that costs almost nothing, as the second timing shows; with a reference of a few hundred thousand polygons, the loading dominates and grows with every file. Measure it before deciding it matters.

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)}

import time

dst = Path("data/enriched")
dst.mkdir(parents=True, exist_ok=True)
start = time.perf_counter()
results = [{"file": path.name, **enrich(path, dst)} for path in sorted(Path("data/parcels").glob("*.gpkg"))]
rows, matched = sum(r["rows"] for r in results), sum(r["matched"] for r in results)
print(f"{len(results)} files in {time.perf_counter() - start:.1f} s: {rows:,} parcels, {matched / rows:.2%} matched")
print("lowest:", [(r["file"], r["rows"], float(r["rate"])) for r in sorted(results, key=lambda r: r["rate"])[:4]])

52 files in 1.1 s: 126,051 parcels, 94.25% matched
lowest: [('sq_378_392.gpkg', 1, 0.0), ('sq_384_386.gpkg', 8, 0.0), ('sq_386_388.gpkg', 1, 0.0), ('zz_relabelled.gpkg', 4987, 0.0)]

Four files matched nothing. Three are squares on the city's edge holding one to eight parcels, whose parcels sit outside the generalised ward boundaries — expected, and harmless. The fourth is 4,987 parcels at 0%: the square saved with a Web Mercator label, reprojected faithfully into the wrong place. Sort by rate alone and it hides among the edge squares; the row count is what gives it away.

Where the time actually goes

Load the reference once, trim it, one CRS, a deliberate predicate

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("reference/wards.gpkg")[["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

parcels = gpd.read_file("data/parcels/sq_384_398.gpkg")   # the city centre
wards = reference()
for predicate in ("within", "covered_by", "covers", "intersects"):
    joined = gpd.sjoin(parcels, wards, predicate=predicate)
    print(f"{predicate:10} {joined.index.nunique():5} of {len(parcels)} parcels matched, "
          f"{joined.index.duplicated().sum()} extra rows from parcels matching two wards")

within      2479 of 2515 parcels matched, 0 extra rows from parcels matching two wards
covered_by  2479 of 2515 parcels matched, 0 extra rows from parcels matching two wards
covers         0 of 2515 parcels matched, 0 extra rows from parcels matching two wards
intersects  2515 of 2515 parcels matched, 40 extra rows from parcels matching two wards

sjoin asks whether each left geometry has the predicate against a right one, so "covers" asks whether a parcel covers a ward — never, here. "Inside or exactly on the edge" is "covered_by". The first version of this block used "covers" for that, and would have tagged nothing.

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"]]
tagged = parcels.join(tags)
print(f"representative point: {tagged['ward_name'].notna().sum()} of {len(parcels)} tagged, "
      f"{tagged.index.duplicated().sum()} duplicates")

representative point: 2515 of 2515 tagged, 0 duplicates

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")
print(len(joined), "rows,", joined["ward_name"].isna().sum(), "with no ward")

2515 rows, 36 with no ward

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:

joined = gpd.sjoin(parcels, wards, predicate="intersects", how="left")
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")])
print(len(joined), "rows after resolving")

  76 rows from 36 multi-matched features
2515 rows after resolving

Or keep both matches, collapsed into one row:

multi = gpd.sjoin(parcels, wards, predicate="intersects", how="left")
tags = (multi.groupby(multi.index)["ward_name"]
             .agg(lambda s: ";".join(sorted(s.dropna().unique()))))
print(tags[tags.str.contains(";")].head(3).to_string())

99     Miles Platting & Newton Heath;Piccadilly
103     Harpurhey;Miles Platting & Newton Heath
171    Miles Platting & Newton Heath;Piccadilly

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)}

for result in results:
    rate = result["rate"]
    if rate < 0.95:
        print(f"  ⚠ {result['file']}: only {rate:.1%} matched")

  ⚠ sq_378_382.gpkg: only 33.3% matched
  ⚠ sq_378_384.gpkg: only 53.1% matched
  ⚠ sq_378_392.gpkg: only 0.0% matched
  ⚠ sq_380_382.gpkg: only 72.0% matched
  ⚠ sq_380_384.gpkg: only 85.2% matched
  ⚠ sq_382_396.gpkg: only 92.0% matched
  ⚠ sq_382_400.gpkg: only 92.1% matched
  ⚠ sq_382_402.gpkg: only 91.6% matched
  ⚠ sq_382_404.gpkg: only 14.3% matched
  ⚠ sq_384_384.gpkg: only 92.1% matched
  ⚠ sq_384_386.gpkg: only 0.0% matched
  ⚠ sq_386_388.gpkg: only 0.0% matched
  ⚠ sq_386_390.gpkg: only 79.0% matched
  ⚠ sq_386_404.gpkg: only 86.4% matched
  ⚠ sq_388_392.gpkg: only 84.0% matched
  ⚠ sq_388_400.gpkg: only 94.9% matched
  ⚠ sq_390_394.gpkg: only 79.7% matched
  ⚠ sq_390_396.gpkg: only 68.4% matched
  ⚠ zz_relabelled.gpkg: only 0.0% 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"]])

count    52.000000
mean      0.839192
std       0.294548
min       0.000000
25%       0.860825
50%       0.976700
75%       0.989075
max       0.998400
Name: rate, dtype: float64
14 files below 90%:
                  file  rows  matched    rate
4      sq_378_392.gpkg     1        0  0.0000
24     sq_384_386.gpkg     8        0  0.0000
34     sq_386_388.gpkg     1        0  0.0000
51  zz_relabelled.gpkg  4987        0  0.0000
22     sq_382_404.gpkg    21        3  0.1429

Enriching Manchester parcel tiles, spatially and by key

The data used on this page

Real files and a real reference layer for Manchester. The files are HM Land Registry's INSPIRE Index Polygons for Manchester City Council (Crown copyright, OGL), split into the 2 km grid squares they would arrive in from a tiling job — 51 GeoPackages in data/parcels/, 121,064 freehold parcels in all. The reference is the city's 32 wards from the ONS Wards (December 2024) generalised boundaries, as reference/wards.gpkg with ward_name and ward_code. One more file is damaged on purpose in a way that happens: data/parcels/zz_relabelled.gpkg is one of the squares whose British National Grid coordinates were saved with a Web Mercator label.

import math
import os
from pathlib import Path

import geopandas as gpd
import requests
from shapely.geometry import box

WARDS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
         "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query"
         "?where=LAD24NM%3D%27Manchester%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
Path("reference").mkdir(exist_ok=True)
if not Path("reference/wards.gpkg").exists():
    wards = gpd.read_file(WARDS).rename(columns={"WD24NM": "ward_name", "WD24CD": "ward_code"})
    wards[["ward_name", "ward_code", "geometry"]].to_file("reference/wards.gpkg")

src = Path("data/parcels")
if not src.exists():
    if not os.path.exists("Manchester_City_Council.zip"):
        with requests.Session() as session:                # the service sets a cookie, then redirects
            open("Manchester_City_Council.zip", "wb").write(session.get(
                "https://use-land-property-data.service.gov.uk/datasets/inspire/download/"
                "Manchester_City_Council.zip", timeout=600).content)
    parcels = gpd.read_file("/vsizip/Manchester_City_Council.zip/Land_Registry_Cadastral_Parcels.gml")
    parcels = parcels[["INSPIREID", "geometry"]].rename(columns={"INSPIREID": "inspire_id"})
    parcels = parcels.set_crs(27700, allow_override=True)
    src.mkdir(parents=True)
    x0, y0, x1, y1 = parcels.total_bounds
    centres = parcels.representative_point()
    for i in range(math.floor(x0 / 2000), math.ceil(x1 / 2000)):
        for j in range(math.floor(y0 / 2000), math.ceil(y1 / 2000)):
            inside = centres.within(box(i * 2000, j * 2000, (i + 1) * 2000, (j + 1) * 2000))
            if inside.any():
                parcels[inside].to_file(src / f"sq_{i * 2}_{j * 2}.gpkg")
    # one square delivered with the wrong CRS label: BNG metres, tagged as Web Mercator
    gpd.read_file(sorted(src.glob("sq_*.gpkg"))[40]).set_crs(3857, allow_override=True) \
        .to_file(src / "zz_relabelled.gpkg")

files = sorted(src.glob("*.gpkg"))
print(len(files), "files,", sum(len(gpd.read_file(f)) for f in files), "parcels;",
      len(gpd.read_file("reference/wards.gpkg")), "wards")

52 files, 126051 parcels; 32 wards

Example 1: the complete enrichment job

%%writefile enrich_job.py
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"]

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)

@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

The workers import enrich_one from this module — under Python 3.14's default forkserver start method they cannot use a function defined in a notebook.

from enrich_job import run

df = run(Path("data/parcels"), Path("data/enriched_all"), workers=8)

52/52 ok · 126,051 rows · overall match 95.5%
⚠ 10 files below 90% match:
                           file  rows  matched   rate
   data/parcels/sq_378_392.gpkg     1        0 0.0000
   data/parcels/sq_386_388.gpkg     1        0 0.0000
data/parcels/zz_relabelled.gpkg  4987        0 0.0000
   data/parcels/sq_378_382.gpkg     6        2 0.3333
   data/parcels/sq_382_404.gpkg    21        7 0.3333

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

enriched = gpd.read_file("data/enriched_all/sq_384_398.gpkg")
lookup = gpd.read_file("reference/wards.gpkg")[["ward_code", "ward_name"]].rename(columns={"ward_name": "name_check"})
print(len(enrich_by_code(enriched, lookup)), "rows, as before:", len(enriched))
try:
    enrich_by_code(enriched, pd.concat([lookup, lookup.head(1)]))
except AssertionError as exc:
    print("refused:", exc)

2515 rows, as before: 2515
refused: lookup has 1 duplicate ward_code values — the merge would multiply rows

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"]].rename_geometry("geom").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")

import os, time
if os.system("docker inspect guide-postgis >/dev/null 2>&1") != 0:
    os.system("docker run -d --name guide-postgis -e POSTGRES_USER=user -e POSTGRES_PASSWORD=guide "
              "-e POSTGRES_DB=gis -p 127.0.0.1:55433:5432 postgis/postgis:16-3.4 >/dev/null")
    time.sleep(10)
engine = create_engine("postgresql+psycopg://user:guide@localhost:55433/gis")
reference().rename_geometry("geom").to_postgis("wards", engine, if_exists="replace", index=False)
via_db = enrich_via_postgis(Path("data/parcels/sq_384_398.gpkg"), engine)
print(len(via_db), "parcels,", via_db["ward_name"].notna().sum(), "tagged in PostGIS")

2515 parcels, 2515 tagged in PostGIS

The first version of this function wrote the probes with GeoPandas' default geometry column name, geometry, and then queried p.geom; PostGIS answered column p.geom does not exist. Name the column when you write it.

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.

The join is not the expensive part

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.

Representative points, index_right and suffixes

• 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.