38 min read 8422 words

Environmental Constraint & Exclusion Screening

Exclusion screening is the stage where a study area stops being land and becomes buildable area, and it is part of the core energy-GIS data and spatial fundamentals pipeline. The failure mode it addresses is quiet and expensive: a screening script loads a dozen constraint layers, subtracts each one in turn from the study polygon, sums what is left, and reports a buildable-area figure that is wrong by more than the margin the project was won on. Nothing raises. The map looks right. The error only surfaces when an environmental consultant redoes the overlay with the same layers and gets a different number.

Three things produce that divergence, and all three are arithmetic rather than ecological. Constraint layers overlap, so subtracting them one at a time removes the shared area more than once. Most national constraint products are rasters on their own grids, so mixing them with vector parcels without an explicit resolution decision quantises the boundary to a cell size nobody chose. And area is not preserved by most projections, so a figure computed in the frame the data happened to arrive in is not the figure a permit reviewer will compute in an equal-area one. This page builds the screening stage that removes all three: layers are harmonised into one metric frame, classified by the legal force they actually carry, unioned once, subtracted once, and reported with a per-layer accounting that reconciles back to the gross area.

Four constraint layers, one union, and the 588 hectares between them On the left, four stacked translucent layers over one study-area outline: a wetlands layer covering 11 percent, a floodplain covering 9 percent, a habitat corridor covering 14 percent and a slope mask covering 21 percent, drawn so their heavy mutual overlap along the drainage is visible. On the right, two bars: the naive sum of the four coverages at 55 percent of the study area, and the union at 41 percent. The difference of 14 percentage points, or 588 hectares, is annotated as the double-counted overlap. 4 200 ha study area · four constraint layers wetlands — 11% floodplain — 9% habitat corridor — 14% slope > 15° — 21% sum of the four coverages 55% — 1 890 ha buildable union of the four layers 41% — 2 478 ha buildable the 14-point gap is overlap counted twice — 588 ha, enough to change the turbine count

Why subtracting constraints one at a time fails

Consider a 4,200-hectare study area with four constraint layers: a National Wetlands Inventory polygon set, a FEMA 100-year floodplain, a state-designated habitat corridor, and a slope mask derived from a 10-metre DEM. Each layer covers a different fraction of the study area — 11, 9, 14 and 21 percent respectively — and a naive sum of those fractions is 55 percent, implying 1,890 hectares of buildable land. The true excluded area is 41 percent, because wetlands sit inside floodplains almost by definition and the habitat corridor follows the same drainage, so the layers overlap heavily. The correct answer is 2,478 hectares, and the naive one under-states buildable land by 588 hectares — enough to change how many turbines fit, or whether the project clears its minimum size.

The inverse error is just as common. A script that subtracts layers sequentially from a shrinking remainder — study.difference(wetlands).difference(floodplain).difference(habitat) — gets the overlap arithmetic right but pays for it in geometry: each difference operation produces a more complex polygon than the last, and by the fourth subtraction the result carries tens of thousands of vertices, slivers along every shared edge, and enough topological noise that the subsequent area calculation is sensitive to the order the layers were applied in. Unioning the constraints first, then subtracting once, produces the same answer with a fraction of the vertices and no order dependence.

The third failure is the unit one. A constraint layer downloaded in EPSG:4326 and subtracted from a parcel layer already reprojected to a UTM zone will silently produce an empty intersection, because the two coordinate ranges do not overlap — the same defect covered in detail under coordinate reference systems for energy projects. When it does not fail outright, it fails softly: the areas come out in square degrees, which are plausible-looking numbers with no physical meaning.

The most consequential modelling decision on this page is not geometric. Constraint layers differ in what they legally do, and collapsing that difference into a single exclusion mask throws away the information a developer most needs. Three classes are worth separating.

Hard exclusions remove land unconditionally: open water, designated wilderness, existing structures, and slopes above the crane specification. No permit makes these buildable, so they belong in the mask that produces the buildable-area headline figure.

Permittable constraints remove land unless a permit is obtained: wetlands under a Section 404 individual permit, floodplain development under a local ordinance, some habitat overlaps under an incidental-take permit. These belong in a second mask, and the honest output is two numbers — buildable without permits, and buildable with a defined permitting path — rather than one.

Advisory layers carry no legal force at all but predict opposition or cost: viewshed sensitivity, prime farmland classifications in states without a farmland statute, informal habitat mapping. These belong in the scoring stage described under building a site suitability scoring pipeline, not in the exclusion mask, because a weighted score can express “expensive” while a mask can only express “impossible”.

Hard, permittable, advisory — three classes, three destinations Three columns. The hard column lists open water, designated wilderness, existing structures and slopes above the crane specification, and routes to the buildable-area headline figure. The permittable column lists Section 404 wetlands, floodplain development and incidental take habitat overlaps, and routes to a second buildable figure conditional on a permitting path. The advisory column lists viewshed sensitivity, prime farmland outside farmland-statute states and informal habitat mapping, and routes to the weighted suitability score rather than to any mask. Classify before you union — the class decides the destination hard exclusion open water designated wilderness existing structures slope above crane spec buildable-area headline permittable Section 404 wetlands floodplain development incidental-take habitat second figure, with a permitting path advisory viewshed sensitivity prime farmland (no statute) informal habitat mapping weighted score, never a mask A layer with no class is not a constraint — it is a shapefile.

Prerequisites and data requirements

The workflow assumes Python 3.11+ with geopandas>=0.14, shapely>=2.0, rasterio>=1.3 and pyproj>=3.6. Inputs are a study-area polygon, a set of vector constraint layers in any GDAL-readable format, and optionally raster constraint masks such as a slope threshold or a land cover classification.

Three input requirements are non-negotiable. Every layer must arrive with a declared CRS or be explicitly tagged at read time; the working frame must be projected and metric, with an equal-area frame such as EPSG:5070 used for any reported hectare figure; and every constraint layer must carry a constraint_class attribute drawn from the three classes above, assigned at ingestion rather than inferred later. A layer without a class is not a constraint, it is a shapefile.

Geometry validity matters more here than in most stages, because union_all on invalid input either raises or — worse — returns a valid-looking result whose area is wrong. Repair before you union, using the predicates described in validating geometry topology with Shapely 2.

Core implementation: one union, one subtraction, one accounting

The function below harmonises the layers, unions each class separately, subtracts once, and returns both the buildable geometry and a per-layer accounting that reconciles to the gross area. The per-layer figures are deliberately computed against the study area rather than against each other: they overlap, they sum to more than the total exclusion, and that is the honest way to report them.

python
from dataclasses import dataclass, field

import geopandas as gpd
from shapely.ops import unary_union

EQUAL_AREA_EPSG = 5070  # NAD83 / CONUS Albers — hectares are only defensible here


@dataclass
class ExclusionResult:
    buildable: gpd.GeoDataFrame
    gross_ha: float
    hard_ha: float
    permittable_ha: float
    per_layer_ha: dict = field(default_factory=dict)


def screen_exclusions(
    study_area: gpd.GeoDataFrame,
    constraints: dict[str, gpd.GeoDataFrame],
    classes: dict[str, str],
    *,
    working_epsg: int,
) -> ExclusionResult:
    """Subtract classified constraint layers from a study area, once.

    `constraints` maps a layer name to its GeoDataFrame; `classes` maps the same
    names to 'hard', 'permittable' or 'advisory'. Advisory layers are measured and
    reported but never subtracted.
    """
    study = study_area.to_crs(working_epsg)
    study_geom = unary_union(study.geometry.values)

    per_layer, hard_parts, permittable_parts = {}, [], []
    for name, gdf in constraints.items():
        cls = classes[name]
        layer = gdf.to_crs(working_epsg)
        # Repair before union: an invalid ring makes the union area meaningless.
        layer["geometry"] = layer.geometry.make_valid()
        clipped = layer.geometry.intersection(study_geom)
        clipped = clipped[~clipped.is_empty]
        if clipped.empty:
            per_layer[name] = 0.0
            continue
        merged = unary_union(clipped.values)
        per_layer[name] = _hectares(merged, working_epsg)
        if cls == "hard":
            hard_parts.append(merged)
        elif cls == "permittable":
            permittable_parts.append(merged)

    hard = unary_union(hard_parts) if hard_parts else None
    permittable = unary_union(permittable_parts) if permittable_parts else None

    buildable_geom = study_geom.difference(hard) if hard is not None else study_geom
    strict_geom = (
        buildable_geom.difference(permittable) if permittable is not None else buildable_geom
    )

    buildable = gpd.GeoDataFrame(
        {"scenario": ["hard_only", "hard_and_permittable"]},
        geometry=[buildable_geom, strict_geom],
        crs=study.crs,
    )
    return ExclusionResult(
        buildable=buildable,
        gross_ha=_hectares(study_geom, working_epsg),
        hard_ha=_hectares(hard, working_epsg) if hard is not None else 0.0,
        permittable_ha=_hectares(permittable, working_epsg) if permittable is not None else 0.0,
        per_layer_ha=per_layer,
    )


def _hectares(geom, working_epsg: int) -> float:
    """Area in hectares, measured in an equal-area frame regardless of the working one."""
    s = gpd.GeoSeries([geom], crs=working_epsg).to_crs(EQUAL_AREA_EPSG)
    return float(s.area.iloc[0]) / 10_000.0

Two details in that function carry most of its value. The intersection with the study area before unioning keeps the constraint geometries bounded — a national wetlands layer unioned in full is an expensive way to compute nothing — and the separate _hectares helper reprojects to an equal-area frame for every measurement, so the working frame can be chosen for distance operations without corrupting the area figures.

Raster constraints: the resolution decision nobody makes explicitly

Slope masks, land cover and many habitat products are rasters, and combining them with vector parcels forces a choice: vectorise the raster, or rasterise the vectors. The choice is usually made by accident, and it changes the answer.

Vectorising a 30-metre slope mask produces polygon boundaries that step in 30-metre increments, so a parcel edge that runs diagonally across the grid gains a staircase whose area error is proportional to the cell size and to the length of the boundary. On a 4,200-hectare study area with 40 kilometres of constraint boundary, a 30-metre quantisation moves roughly 60 hectares — 1.4 percent — and the direction of the error depends on whether cells are included when their centre or any part is covered.

Rasterising the vectors instead makes the whole computation cell-based, which is fast and internally consistent, at the cost of quantising the parcel boundary too. That is usually the right trade for a screening pass over thousands of parcels and the wrong one for the final buildable-area figure on a shortlisted site, where the vector boundary is the legal object and the raster is the approximation.

Cell size, boundary staircase and the hectares it moves On the left, a diagonal constraint boundary drawn over a coarse grid with the cells that would be marked excluded shaded, producing a visible staircase against the true line. In the middle, the same boundary over a grid three times finer, with a much closer staircase. On the right, a table of the area moved on 40 kilometres of constraint boundary: about 60 hectares at 30 metre cells, 20 hectares at 10 metre cells and 2 hectares at 1 metre cells, with a note that centre-based and any-part-covered rules push the error in opposite directions. The same boundary at two cell sizes 30 m cells 10 m cells red line: the true constraint boundary · shaded: cells marked excluded area moved on 40 km of boundary 30 m cells → ≈ 60 ha 10 m cells → ≈ 20 ha 1 m cells → ≈ 2 ha centre-covered under-excludes; any-part-covered over-excludes — state the rule with the figure Rasterise for screening, vectorise for the shortlisted site, and record which produced the number.

The practical rule: rasterise for screening, vectorise for the final site, and record which was used alongside the figure. A buildable-area number without its resolution provenance cannot be reconciled against anyone else’s.

python
import numpy as np
import rasterio
from rasterio.features import geometry_mask


def slope_exclusion_ha(
    dem_path: str,
    study_geom,
    *,
    max_slope_deg: float = 15.0,
) -> tuple[float, float]:
    """Excluded hectares from a slope threshold, plus the cell-size uncertainty band."""
    with rasterio.open(dem_path) as src:
        window = src.window(*study_geom.bounds)
        dem = src.read(1, window=window, masked=True).astype("float32")
        transform = src.window_transform(window)
        cell_m = abs(transform.a)

        dzdy, dzdx = np.gradient(dem, cell_m)
        slope_deg = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))

        inside = ~geometry_mask(
            [study_geom], out_shape=dem.shape, transform=transform, invert=False
        )
        steep = (slope_deg > max_slope_deg) & inside & ~dem.mask

        cells = int(steep.sum())
        ha = cells * (cell_m ** 2) / 10_000.0
        # Boundary cells are the uncertainty: half a cell along the excluded perimeter.
        perimeter_cells = int(
            (steep ^ np.roll(steep, 1, axis=0)).sum() + (steep ^ np.roll(steep, 1, axis=1)).sum()
        )
        band_ha = perimeter_cells * 0.5 * (cell_m ** 2) / 10_000.0
    return ha, band_ha

Returning the uncertainty band alongside the figure is what makes a raster-derived exclusion defensible. A slope exclusion of “412 hectares” invites a challenge; “412 hectares ± 18 from a 10-metre grid” answers it in advance.

Error handling and edge cases

A constraint layer that does not intersect the study area. This is normal — a national layer clipped to a county usually contributes nothing — and must not be an error. What must be an error is a layer that intersects nothing across every study area in a batch, which almost always means a CRS or extent problem rather than an absence of constraints.

A study area entirely excluded. Report it as zero buildable hectares with the binding layer named, never as a dropped row. A parcel that vanishes from the output is indistinguishable from one that was never submitted, and that distinction is exactly what a landowner will ask about.

Multipart and nested geometry. Wetland layers routinely contain polygons with holes, and a hole inside a wetland is not buildable land unless the hole is genuinely upland. Preserve the rings as authored, repair rather than simplify, and let difference handle the topology — manual ring manipulation is where most of the subtle area errors in this stage come from.

Layers with different vintages. A 2019 wetlands delineation and a 2024 floodplain revision describe different moments. Record the vintage per layer in the accounting output; a reviewer comparing against current data needs to know which layer is stale, and the answer is rarely the one they assume.

Performance and scalability

The expensive operation is the union, and its cost is driven by vertex count rather than by feature count. Three levers matter, in order. Clip each constraint layer to the study area before unioning — this is the single largest win, and it is why the implementation above does it first. Simplify only the layers whose boundary precision does not carry legal weight, and never the ones that do; a simplify(tolerance=1.0) on an advisory viewshed layer is free, and the same call on a wetland delineation is a compliance problem. Finally, use the spatial index to skip constraint features whose bounding box misses the study area entirely, which on a national layer removes the large majority before any geometry is touched.

For portfolio runs, the shape of the parallelism is per study area rather than per layer: each site is independent, the constraint layers are read-only, and a worker pool over sites scales linearly until the object store rather than the CPU becomes the limit. Reading the constraint layers once into each worker and reusing them across that worker’s sites avoids re-reading a national layer for every parcel.

Validation and audit trail

The accounting output is the deliverable, not a by-product. For every screening run, record the gross study area, the per-layer excluded area, the unioned hard and permittable areas, the resulting buildable figures for both scenarios, and the sum of the per-layer figures alongside the union — because the gap between those two numbers is the overlap, and a reviewer will want it.

The reconciliation every exclusion report needs to survive review A table of the screening accounting. Gross study area: 4,200 hectares. Per-layer excluded area: wetlands 462, floodplain 378, habitat corridor 588 and slope 882 hectares, summing to 2,310. Union of all constraints: 1,722 hectares, with the 588 hectare difference labelled as overlap. Two buildable results follow: 2,478 hectares with hard exclusions removed, and 2,142 hectares with permittable constraints removed as well. Three CI assertions are listed beneath: buildable never exceeds gross, the union never exceeds gross, and the per-layer sum is never less than the union. The numbers a reviewer will ask for, published together gross study area 4 200 ha − wetlands 462 ha − floodplain 378 ha − habitat corridor 588 ha − slope > 15° 882 ha per-layer sum 2 310 ha union of constraints 1 722 ha overlap (sum − union) 588 ha buildable · hard only 2 478 ha buildable · hard + permittable 2 142 ha CI assertions buildable ≤ gross union ≤ gross per-layer sum ≥ union Each assertion catches one bug: a subtraction in the wrong frame, an unclipped national layer, a layer measured in another CRS Publish the overlap explicitly — it is the first question asked

Three assertions belong in CI. Buildable area must never exceed gross area, which catches a subtraction performed in the wrong frame. The union of the constraint layers must never exceed the gross area either, which catches an unclipped layer. And the sum of the per-layer areas must be greater than or equal to the union area, which is a tautology when the arithmetic is right and fails loudly when a layer was measured in one frame and unioned in another.

python
def assert_exclusion_integrity(result: ExclusionResult) -> None:
    """CI gate: the accounting has to reconcile before the figure is published."""
    hard_only = result.buildable.loc[
        result.buildable["scenario"] == "hard_only", "geometry"
    ].iloc[0]
    strict = result.buildable.loc[
        result.buildable["scenario"] == "hard_and_permittable", "geometry"
    ].iloc[0]

    assert result.gross_ha > 0, "empty study area"
    assert _hectares(hard_only, EQUAL_AREA_EPSG) <= result.gross_ha * 1.0001, (
        "buildable exceeds gross — the subtraction ran in a non-equal-area frame"
    )
    assert strict.area <= hard_only.area * 1.0001, (
        "the permittable scenario is larger than the hard-only one — masks were swapped"
    )
    assert sum(result.per_layer_ha.values()) >= result.hard_ha * 0.9999, (
        "per-layer areas sum to less than their union — a layer was measured in another CRS"
    )

Frequently asked questions

Should the exclusion mask be built once per region or once per project?

Once per project, from region-wide source layers. The geometry is shared, but the classification is not: whether a wetland is a hard exclusion or a permittable one depends on the project’s permitting strategy, and whether a slope is excluded depends on the crane specification. Sharing the resolved mask between projects silently imposes one project’s assumptions on another.

How should overlapping constraint layers be reported?

Both ways, with the overlap named. Report the per-layer area — which sums to more than the total — and the unioned area, and give the difference a label. Reviewers ask the overlap question every time, and a report that pre-empts it avoids a rebuild.

Is a 30-metre national land-cover product good enough for exclusion screening?

For screening, yes; for a final buildable-area figure, no. The quantisation error scales with the length of the constraint boundary, so it is small on compact parcels and material on long, irregular ones. Use the national product to rank sites, and a site-specific delineation for the shortlist.

What about constraints that are not spatial at all?

Encode them as attributes on the study area rather than forcing them into geometry. A parcel whose title carries a conservation easement is fully excluded without any polygon being involved, and representing that as a geometric mask makes it invisible in a spatial audit. The exclusion pipeline should accept both geometric and attribute-based exclusions and report them in the same accounting.

How do I handle a constraint layer that is itself derived from another?

Record the derivation and do not double-count. A hydric-soils layer and a wetlands delineation overlap by construction because one informs the other, so treating them as independent constraints inflates the apparent constraint coverage without changing the union. The accounting handles this correctly as long as both are unioned rather than summed — which is another reason the union figure, not the sum, is the headline.

Explore this section

Back to Core Energy-GIS Data & Spatial Fundamentals