13 min read 2906 words

Building Multi-Layer Exclusion Masks with GeoPandas Overlay

The scenario is a screening script that runs, produces a buildable-area figure, and disagrees with the environmental consultant’s figure by 14 percent. No exception was raised, both analysts used the same layers, and neither can immediately say why. This page fixes that class of disagreement, and it is the mechanical half of environmental constraint and exclusion screening.

Root-cause analysis

Three mechanisms produce a divergent buildable-area figure, and a real disagreement usually involves more than one.

  1. Sequential subtraction with a sum. Computing each layer’s area, summing those areas, and subtracting the sum from the gross area double-counts every overlap. Wetlands sit inside floodplains, habitat corridors follow drainage, and the overlap on a typical study area is a quarter to a third of the total constrained area.
  2. Chained differences. study.difference(a).difference(b).difference(c) gets the arithmetic right and the geometry wrong: each operation adds vertices and slivers along shared edges, and by the fourth layer the result carries enough topological noise that its area depends on the order the layers were applied in.
  3. A frame that does not preserve area. A hectare figure computed in a conformal frame is latitude-dependent, so two analysts working in different UTM zones or in Web Mercator will legitimately disagree — the same distortion covered under projection and CRS quick reference.
Sum, chain, or union — three combination strategies compared Three rows over the same four constraint layers. The first, sum then subtract, produces 1,890 buildable hectares and is marked as double-counting every overlap. The second, chained difference, produces about 2,478 hectares but with 41,000 output vertices and an answer that changes slightly with layer order. The third, union then difference, produces 2,478 hectares with 4,100 vertices and no order dependence. Each row carries its vertex count and whether its answer depends on layer ordering. Same four layers, three combination strategies sum of layer areas, then subtract double-counts every overlap 1 890 ha — chained .difference() per layer answer drifts with layer order 2 478 ha 41 000 vertices union_all(), then one difference order-independent 2 478 ha 4 100 vertices The middle row is the dangerous one: it is arithmetically right, so nobody looks at it again — and its geometry degrades until the area depends on the order the shapefiles happened to be listed in.

Pre-flight validation

Before any overlay runs, three properties have to hold: every layer carries a declared CRS, every geometry is valid, and every layer has been clipped to the study area. The third is not just a performance measure — an unclipped national layer unioned in full can produce a geometry whose area exceeds the study area, which then makes the accounting assertions meaningless.

python
import geopandas as gpd


def preflight_constraint_layers(
    study: gpd.GeoDataFrame,
    layers: dict[str, gpd.GeoDataFrame],
    *,
    working_epsg: int,
) -> dict[str, dict]:
    """Surface the three faults that make an overlay disagree, before it runs."""
    report: dict[str, dict] = {}
    study_area = study.to_crs(working_epsg).union_all()

    for name, gdf in layers.items():
        if gdf.crs is None:
            raise ValueError(f"{name}: undeclared CRS — set_crs() the true source frame first")
        projected = gdf.to_crs(working_epsg)
        invalid = int((~projected.is_valid).sum())
        empty = int(projected.geometry.is_empty.sum())
        intersects = int(projected.geometry.intersects(study_area).sum())
        report[name] = {
            "features": len(projected),
            "invalid": invalid,
            "empty": empty,
            "intersecting_study_area": intersects,
            "source_epsg": gdf.crs.to_epsg(),
        }
        if intersects == 0:
            # Not fatal on its own — a national layer may genuinely miss this county.
            report[name]["warning"] = "no features intersect the study area"
    return report

A layer reporting zero intersecting features is normal in isolation and alarming in aggregate: when every layer reports zero, the study area is almost certainly in a different frame from the constraints, not in an unconstrained paradise.

Fix implementation

The correct shape is: clip, repair, union per class, subtract once. GeoPandas overlay with how="difference" does the subtraction, and union_all() does the union — the important part is that the union happens before the difference, exactly once.

python
import geopandas as gpd

EQUAL_AREA_EPSG = 5070


def build_exclusion_mask(
    study: gpd.GeoDataFrame,
    layers: dict[str, gpd.GeoDataFrame],
    *,
    working_epsg: int,
    buffer_m: dict[str, float] | None = None,
) -> tuple[gpd.GeoDataFrame, dict[str, float]]:
    """Return the buildable remainder plus a per-layer hectare accounting."""
    buffer_m = buffer_m or {}
    study_p = study.to_crs(working_epsg)
    study_geom = study_p.union_all()

    per_layer: dict[str, float] = {}
    pieces = []
    for name, gdf in layers.items():
        layer = gdf.to_crs(working_epsg)
        layer["geometry"] = layer.geometry.make_valid()
        # Clip first: this is the single largest cost saving in the whole function.
        clipped = layer.clip(study_geom)
        if clipped.empty:
            per_layer[name] = 0.0
            continue
        geom = clipped.union_all()
        if buffer_m.get(name):
            geom = geom.buffer(buffer_m[name])   # working room, in metres
        per_layer[name] = _ha(geom, working_epsg)
        pieces.append(geom)

    if not pieces:
        return study_p.assign(exclusion_ha=0.0), per_layer

    excluded = gpd.GeoSeries(pieces, crs=study_p.crs).union_all()
    buildable = gpd.GeoDataFrame(
        {"excluded_ha": [_ha(excluded, working_epsg)]},
        geometry=[study_geom.difference(excluded)],
        crs=study_p.crs,
    )
    return buildable, per_layer


def _ha(geom, working_epsg: int) -> float:
    return float(
        gpd.GeoSeries([geom], crs=working_epsg).to_crs(EQUAL_AREA_EPSG).area.iloc[0]
    ) / 10_000.0

Two details carry the correctness. make_valid() runs before any union, because unioning an invalid ring produces an area that is wrong without raising. And _ha reprojects to an equal-area frame for every measurement, so the working frame stays free for distance operations without contaminating the hectare figures.

Vertex growth under chained differences versus one union A chart of output vertex count after each of four constraint layers is applied. The chained difference series rises steeply: 1,100 vertices after the first layer, 6,800 after the second, 19,400 after the third and 41,000 after the fourth. A single flat marker shows the union-then-difference result at 4,100 vertices regardless of layer count. An inset shows the sliver polygons that appear along nearly-coincident boundaries in the chained result and are absent from the unioned one. Vertices in the buildable polygon after each layer 0k 10k 20k 30k 40k layer 1 layer 2 layer 3 layer 4 union then one difference — 4 100 chained — 41 000 the extra vertices are slivers along nearly coincident boundaries they cost area, memory and every downstream predicate that touches them Slivers are not cosmetic: each one is a tiny polygon that a later intersection has to evaluate, and a handful of them is what turns a two-second overlay into a two-minute one.

Fallback routing and performance tuning

  • Clip before you union. A national wetlands layer has millions of vertices; the part inside one county has thousands. Clipping first turns the union from the dominant cost into a rounding error.
  • Simplify only what carries no legal weight. A simplify(tolerance=1.0) on an advisory viewshed layer is free; the same call on a wetland delineation changes a regulated boundary.
  • Use the spatial index implicitly. GeoDataFrame.clip already queries the index, so an explicit sindex.query before it buys nothing — the common mistake is doing neither and calling intersection on the full layer.
  • Union in one call, not in a loop. union_all() on a list is substantially faster than repeated pairwise unions, because it can use a cascaded strategy rather than rebuilding the accumulated geometry each time.
  • Keep the mask, not the difference, when reusing. For a portfolio, computing the exclusion union once and differencing it against each parcel is far cheaper than rebuilding the union per parcel.

Downstream validation

The accounting is what makes two analysts agree. Publish the gross area, each layer’s clipped area, the union area, and the difference between the per-layer sum and the union — that difference is the overlap, and it is the number the disagreement was always about.

python
def assert_mask_reconciles(gross_ha: float, per_layer: dict[str, float], excluded_ha: float,
                           buildable_ha: float) -> None:
    """CI gate: the three relations that must hold if the arithmetic is right."""
    assert buildable_ha <= gross_ha * 1.0001, "buildable exceeds gross — wrong frame for area"
    assert excluded_ha <= gross_ha * 1.0001, "exclusion exceeds gross — a layer was not clipped"
    assert sum(per_layer.values()) >= excluded_ha * 0.9999, (
        "per-layer sum below the union — a layer was measured in a different CRS"
    )
    assert abs((gross_ha - excluded_ha) - buildable_ha) < max(0.01, gross_ha * 1e-6), (
        "gross − excluded ≠ buildable — the difference and the union disagree"
    )
Four assertions and the specific disagreement each one prevents A four-row table pairing an assertion with the bug it catches. Buildable area less than or equal to gross area catches area measured in a conformal rather than equal-area frame. Excluded area less than or equal to gross catches a constraint layer that was never clipped to the study area. The sum of per-layer areas being at least the union area catches a layer measured in a different CRS from the union. And gross minus excluded equalling buildable catches a difference computed against a stale union. Four assertions, four disagreements they prevent buildable_ha <= gross_ha area computed in a conformal frame excluded_ha <= gross_ha a constraint layer that was never clipped sum(per_layer) >= excluded_ha a layer measured in a different CRS gross − excluded == buildable a difference against a stale union Publish the per-layer sum and the union together — the gap between them is the overlap, and it is the first question asked.

Frequently asked questions

Why does overlay(how="difference") return more rows than the input?

Because the difference can split one polygon into several disjoint parts, and GeoPandas returns them as separate rows unless the geometry is explicitly recombined. That is usually what you want for a buildable-area map and never what you want for a per-parcel accounting — dissolve back to the input key before aggregating, or the same parcel appears several times in the totals.

Should exclusions be buffered before or after the union?

Before, per layer, because different layers need different working-room offsets: a wetland delineation needs the regulatory buffer, a road needs the construction offset, and a viewshed needs none at all. Buffering after the union applies one offset to everything and quietly grows the exclusion by more than any single rule requires.

What is the fastest way to test one parcel against a prepared mask?

Prepare the mask geometry once and test parcels against the prepared version. The mask is reused thousands of times in a portfolio run, which is exactly the case prepared geometry exists for — the preparation cost is amortised after roughly twenty tests.

How should the mask handle holes?

Leave them as authored. A hole inside a wetland polygon is upland the delineation deliberately excluded, and it is buildable unless another layer says otherwise. Filling holes “to clean up the geometry” is a common and expensive tidying instinct: it removes real buildable land and is invisible in the final map.

Does the mask need to be a single polygon?

No, and forcing it to be one is usually a mistake. A study area with several disjoint buildable pockets is genuinely a MultiPolygon, and dissolving it into one geometry with a convex hull or a generous buffer merges pockets that are separated by real constraints. Keep the parts, and report their count and their individual areas — a single 400-hectare pocket and eight 50-hectare pockets are very different projects.

How should the mask be stored between runs?

As a versioned artefact keyed on the layer versions that produced it, in the working CRS, with the per-layer accounting alongside. Recomputing a national mask for every parcel is the most common performance mistake in this stage, and caching the geometry without caching what produced it is the most common correctness one — a mask whose provenance is unknown cannot be reused with confidence.

What tolerance should be used when comparing two analysts’ figures?

Small: a fraction of a percent. Once both figures are computed in an equal-area frame from the same layer versions with the same buffers, they should agree to within floating-point noise. A gap larger than that is not tolerance — it is a difference in inputs or in method, and chasing it down is how the discrepancy in the opening paragraph gets resolved.

Related articles

Back to Environmental Constraint & Exclusion Screening