Clipping Solar Parcels to County Setback Boundaries in GeoPandas #

You run gpd.clip on a parcel layer, subtract the county setback, and the reported buildable area comes back as 0.0007 — or the whole overlay dies with a GEOSException: TopologyException, or every parcel silently returns an empty envelope. This is the failure this page exists to eliminate: turning a parcel and a county setback ordinance into a defensible buildable_area_ha figure that a solar developer can size a project against. It is a downstream step of Regulatory Boundary Mapping — the compliance mask tells you which land is excluded, and this operation measures how much of a given parcel survives once the setback rule is applied. Get the projection or the geometry handling wrong and the buildable acreage feeding the entire financial model is wrong, usually without an exception to warn you.

The operation itself is small. The buildable geometry is the parcel eroded inward by the setback distance $s$, with hard exclusions $E$ (roads, wetlands, occupied structures) removed:

where $P\ominus s$ is the morphological erosion — a negative buffer — and division by $10^4$ converts square metres to hectares. Every production failure lives in the data $P$ and $E$ consume, and in the CRS the area is measured in, not in that formula.

Root-cause analysis #

Four compounding causes account for nearly every broken setback calculation, and each maps to a distinct fix below.

1. clip / overlay run on a geographic CRS. If parcels_gdf.crs is EPSG:4326, .area and .buffer(-30) operate in degrees, not metres. A 30 m setback becomes a 30-degree buffer that erases the parcel to empty, and areas print as tiny decimal degree² values that are physically meaningless. This is the single most common cause and it never raises — it just returns nonsense.
2. Invalid parcel geometry aborts the overlay. County parcel fabrics are full of self-intersecting rings, duplicate vertices, and bowtie polygons. gpd.overlay and gpd.clip push these straight into the GEOS engine, which throws a TopologyException mid-run and kills the batch — or, worse, returns a partially corrupted result.
3. Negative-buffer collapse. geometry.buffer(-s) on a parcel narrower than $2s$ returns an empty polygon; on an L-shaped or pinched parcel it returns a multipart polygon of disconnected slivers. Feed either downstream and you get empty overlays, exploded row counts, or buildable fragments too small to site a single table on.
4. One global setback for every county. Setback ordinances vary by jurisdiction — property-line, occupied-dwelling, and road setbacks differ from county to county. Hard-coding a single setback_m applies the wrong rule to most parcels in a multi-county portfolio, and the error is invisible in the output geometry.

Pre-flight validation #

Surface the root cause before the buffer and overlay run. Refusing a geographic CRS and a missing per-county setback column up front converts three of the four silent failures into a fast, precise error — which is exactly what you want in a CI/CD screening job that processes thousands of parcels unattended.

import geopandas as gpd


def preflight_setback_inputs(
    parcels_gdf: gpd.GeoDataFrame,
    exclusions_gdf: gpd.GeoDataFrame,
) -> None:
    """Raise on the exact root cause before any buffer or overlay is evaluated."""
    # Cause 1: geographic CRS makes buffer() and area() return degrees, not metres
    for name, gdf in (("parcels", parcels_gdf), ("exclusions", exclusions_gdf)):
        if gdf.crs is None:
            raise ValueError(f"{name} has no CRS; refusing to buffer or measure blindly.")
        if gdf.crs.is_geographic:
            raise ValueError(
                f"{name} is geographic ({gdf.crs.to_epsg()}); reproject to a metric CRS "
                "(EPSG:5070 or a local UTM zone) before applying a metre setback."
            )
    # Cause 4: a per-county rule must be attached per parcel, not hard-coded
    if "setback_m" not in parcels_gdf.columns:
        raise ValueError(
            "parcels missing per-county 'setback_m'; join the county rule table first."
        )
    if parcels_gdf["setback_m"].isna().any():
        raise ValueError("null setback_m present; every parcel needs a county setback rule.")
    # Cause 2: report invalid rings so the repair step is not a surprise
    invalid = int((~parcels_gdf.geometry.is_valid).sum())
    if invalid:
        print(f"[preflight] {invalid} parcels have invalid geometry; make_valid will repair them.")

Fix implementation #

The corrected function reprojects to a metric CRS, repairs geometry, applies the per-parcel setback as a negative buffer, erases hard exclusions with overlay(how="difference"), then cleans up the collapse artefacts before measuring. Parameter choices are justified for solar siting: EPSG:5070 (CONUS Albers equal-area) keeps hectare totals comparable across a multi-county portfolio; a single project can instead use its local UTM zone. The min_fragment_m2 floor discards slivers too small to hold a module table, and keep_geom_type=True stops the difference from leaking stray lines or points into the result.

import geopandas as gpd
import pandas as pd


def compute_buildable_area(
    parcels_gdf: gpd.GeoDataFrame,
    exclusions_gdf: gpd.GeoDataFrame,
    target_epsg: int = 5070,          # CONUS Albers equal-area; use a UTM zone for one project
    min_fragment_m2: float = 100.0,   # drop buildable slivers smaller than one table footprint
) -> gpd.GeoDataFrame:
    """Erode each parcel by its county setback, erase hard exclusions, and
    return per-parcel buildable_area_m2 / buildable_area_ha in a metric CRS."""
    preflight_setback_inputs(parcels_gdf, exclusions_gdf)

    # 1. Reproject once so buffer() and area() are in metres (Cause 1)
    parcels = parcels_gdf.to_crs(epsg=target_epsg).copy()
    exclusions = exclusions_gdf.to_crs(epsg=target_epsg).copy()

    # 2. Repair topology before any overlay (Cause 2) — shapely 2.0 / geopandas >= 0.14
    parcels["geometry"] = parcels.geometry.make_valid()
    exclusions["geometry"] = exclusions.geometry.make_valid()

    # 3. Per-parcel negative buffer = county setback from the parcel boundary (Cause 4)
    #    Passing an array-like distance applies each row's own setback_m.
    parcels["geometry"] = parcels.geometry.buffer(-parcels["setback_m"].to_numpy())

    # 4. Negative-buffer collapse cleanup (Cause 3): drop parcels eroded to nothing
    parcels = parcels[~parcels.geometry.is_empty & parcels.geometry.notna()]

    # 5. Erase hard exclusions — difference keeps parcel land NOT covered by exclusions
    buildable = gpd.overlay(parcels, exclusions, how="difference", keep_geom_type=True)

    # 6. Explode multipart remnants and filter slivers, then re-dissolve per parcel
    buildable = buildable.explode(index_parts=False)
    buildable = buildable[buildable.geometry.area >= min_fragment_m2]
    per_parcel = buildable.dissolve(by="parcel_id", as_index=True, aggfunc="first")

    # 7. Report area in both units for siting and audit
    per_parcel["buildable_area_m2"] = per_parcel.geometry.area
    per_parcel["buildable_area_ha"] = per_parcel["buildable_area_m2"] / 1e4
    return per_parcel

Applying the negative buffer as a NumPy array of per-row distances is what keeps the county rule attached to the right parcel — buffer(-30) would silently overwrite every jurisdiction’s ordinance with a single value. Using overlay(how="difference") rather than gpd.clip is deliberate: clip keeps the intersection with a mask, whereas here you want the parcel land that the exclusions do not cover, which is the set difference.

Once you have hectares, the developable capacity follows from a packing density $\rho$ (typically ~0.2–0.4 MW/ha for fixed-tilt utility PV):

so a defensible buildable_area_ha is the load-bearing number the rest of the model rests on.

Fallback routing & performance tuning #

For statewide parcel fabrics or unattended CI/CD screening runs, layer these strategies on top of the core function:

• Reproject once, at the boundary. Transform parcels and exclusions a single time up front, never inside a per-county loop — repeated to_crs is the most common avoidable hotspot in setback batches.
• Pre-buffer exclusions before the difference. Roads and property lines usually arrive as lines or points; buffer them to their own statutory setback first, the same line-to-polygon step covered in Grid Capacity Buffer Analysis, so the difference operates polygon-on-polygon.
• Build the spatial index once. Call exclusions.sindex before the overlay so GEOS prunes non-intersecting pairs, dropping the cost from O(N×M) toward near-linear on large county layers.
• Tune the sliver floor per technology. A single-axis tracker needs a larger contiguous min_fragment_m2 than a rooftop array; expose it as a parameter rather than hard-coding 100.0.
• Chunk by county with dask-geopandas. Partition parcels by county_fips and map the function across workers; each partition already shares one setback rule, so the join and buffer stay cache-friendly.

Downstream validation #

Before a buildable figure feeds a capacity model, gate it with assertions suitable for a CI/CD pipeline. The two invariants that catch the failures above are that no parcel reports a negative buildable area, and that buildable area never exceeds the source parcel area — a setback can only ever remove land, so a result larger than its parent means a CRS or geometry regression slipped through.

def assert_buildable_integrity(
    result_gdf: gpd.GeoDataFrame,
    parcels_gdf: gpd.GeoDataFrame,
    target_epsg: int = 5070,
    area_tol_m2: float = 1.0,
) -> None:
    """CI/CD gate: fail the build if buildable areas are not physically defensible."""
    assert result_gdf.crs.to_epsg() == target_epsg, "result drifted off the metric CRS"
    assert result_gdf.geometry.is_valid.all(), "invalid buildable geometry after overlay"
    assert (result_gdf["buildable_area_m2"] >= 0).all(), "negative buildable area present"

    # buildable must never exceed the original parcel footprint
    src_area = parcels_gdf.to_crs(epsg=target_epsg).set_index("parcel_id").geometry.area
    check = result_gdf["buildable_area_m2"].to_frame().join(src_area.rename("parcel_area_m2"))
    breach = check["buildable_area_m2"] > check["parcel_area_m2"] + area_tol_m2
    assert not breach.any(), f"buildable exceeds parcel area for {list(check.index[breach])}"

Logging the buildable-to-parcel ratio per county as part of the run is what keeps the screening auditable: a reviewer preparing a permitting package can see at a glance that a 640-acre parcel yielded, say, 71% buildable land after setbacks rather than an implausible 103%. The same geometry-integrity discipline is enforced in depth by Spatial Data Quality Validation; pin geopandas, shapely, and pyproj versions in pyproject.toml so a make_valid or overlay behaviour change cannot silently shift the buildable area between runs.

Related #

• Regulatory Boundary Mapping — the parent workflow that produces the exclusion mask this operation clips against.
• Automating US County Boundary Extraction with OSMnx — pull the county polygons that carry each jurisdiction’s setback rule.
• Spatial Data Quality Validation — geometry repair and validity checks the overlay depends on.
• Grid Capacity Buffer Analysis — buffer line and point infrastructure into the polygon exclusions erased here.