Computing Boolean Operations on CAD Footprints with Shapely

To union, difference or intersect CAD-derived footprints reliably, repair every ring before the operation, shift the coordinates to a local origin, use a single collection-wide call rather than a pairwise fold, and handle the multi-part results a CAD source will produce. Boolean algorithms assume valid input and are not obliged to detect that they did not get it. This page is part of Choosing a Geometry Engine for Python Pipelines.

Why CAD Input Breaks Boolean Operations

Geometry that came from a drawing arrives with three properties the algorithms dislike.

Rings are frequently not simple. A closed polyline drafted by a person can cross itself, double back, or repeat a vertex, and none of that is visible at drawing scale. A self-intersecting ring has no well-defined interior, so an operation on it has no well-defined answer.

Coordinates are large. Full projected easting and northing values put the arithmetic at seven significant figures before the decimal point, where the numerical margins that make the predicates robust are proportionally much smaller. Shifting to a local origin costs one subtraction and moves the whole computation into a range where it behaves.

Vertices are nearly coincident rather than coincident. Two footprints that share an edge in the drawing usually share it to within a micron rather than exactly, which produces slivers — zero-width polygons a few nanometres across — in the output. They are valid geometry and they are noise.

Production-Ready Script

# shapely>=2.0, numpy>=1.24, Python 3.9+
from __future__ import annotations

from dataclasses import dataclass
import numpy as np
from shapely.geometry import Polygon, MultiPolygon
from shapely.geometry.base import BaseGeometry
from shapely.ops import unary_union
from shapely.validation import make_valid, explain_validity
from shapely.affinity import translate


@dataclass(frozen=True)
class MergeReport:
    inputs: int
    repaired: int
    dropped: int
    parts_out: int


def _clean(poly: Polygon) -> BaseGeometry | None:
    """Repair a ring, or report it unusable. Never returns invalid geometry."""
    if poly.is_empty or poly.area <= 0:
        return None
    if poly.is_valid:
        return poly
    fixed = make_valid(poly)
    if fixed.is_empty:
        return None
    return fixed


def merge_footprints(
    polygons: list[Polygon], *, sliver_area: float = 1e-6
) -> tuple[BaseGeometry, MergeReport]:
    """Union a collection of CAD footprints, robustly."""
    if not polygons:
        raise ValueError("nothing to merge")

    # A common local origin keeps the arithmetic at small magnitudes.
    all_coords = np.vstack([np.array(p.exterior.coords) for p in polygons if not p.is_empty])
    ox, oy = all_coords.mean(axis=0)

    cleaned, repaired, dropped = [], 0, 0
    for p in polygons:
        shifted = translate(p, xoff=-ox, yoff=-oy)
        was_valid = shifted.is_valid
        fixed = _clean(shifted)
        if fixed is None:
            dropped += 1
            continue
        if not was_valid:
            repaired += 1
        cleaned.append(fixed)

    if not cleaned:
        raise ValueError("every input polygon was empty or unrepairable")

    merged = unary_union(cleaned)          # one pass, not a pairwise fold

    # Drop slivers produced by near-coincident edges.
    if isinstance(merged, MultiPolygon):
        keep = [g for g in merged.geoms if g.area > sliver_area]
        merged = MultiPolygon(keep) if len(keep) > 1 else (keep[0] if keep else merged)

    parts = len(merged.geoms) if hasattr(merged, "geoms") else 1
    return translate(merged, xoff=ox, yoff=oy), MergeReport(
        inputs=len(polygons), repaired=repaired, dropped=dropped, parts_out=parts)


def difference_with_report(a: Polygon, b: Polygon) -> BaseGeometry:
    """Difference with the post-condition a boolean will not check for you."""
    for name, g in (("a", a), ("b", b)):
        if not g.is_valid:
            raise ValueError(f"{name} is invalid: {explain_validity(g)}")
    result = a.difference(b)
    if result.area > a.area + 1e-9:
        raise ValueError("difference increased the area — check ring orientation")
    return result


if __name__ == "__main__":
    merged, report = merge_footprints([...])
    print(report, merged.geom_type)

Key implementation notes:

• Validation happens before every operation, and repair is counted. A run that repaired 300 of 400 inputs is telling you something about the source drawing that a silent repair would hide.
• The origin shift is applied to the inputs and reversed on the output, so the caller sees no difference except robustness.
• unary_union merges the whole collection at once. The pairwise fold is the common first implementation and is markedly slower and less robust at scale.
• Sliver filtering uses an area threshold rather than a buffer trick, because a zero-buffer round trip is itself a source of new invalidity.
• The area post-condition on difference is three lines and catches a class of orientation bug the operation will not report.

Compatibility Matrix

Fallback Strategies

1. TopologyException from GEOS. Almost always invalid input that slipped past validation, or coordinates large enough to lose precision. Confirm the origin shift is applied and that every input passed _clean.

2. The union result is a GeometryCollection. Inputs touch along edges rather than overlapping, producing lines alongside polygons. Filter to polygons, and treat the collection as a signal that the drawing’s edges are coincident-ish rather than shared.

3. Slivers survive the filter. The threshold is too small for the coordinate scale. Set it relative to the smallest meaningful feature area rather than to an absolute constant.

4. Everything merges into one shape. Footprints at different elevations, unioned in plan. This is Shapely working correctly on the wrong question — see the parent page for when a planar engine is the wrong choice.

5. The result has fewer parts than expected. Near-coincident edges bridged features that are separate in reality. Snap the inputs to a deliberate tolerance grid before merging, so the bridging is a decision rather than an accident of precision.

FAQ

Related Pages

• Choosing a Geometry Engine for Python Pipelines — parent reference comparing the planar, mesh and solid engines
• Converting CAD Polylines to GeoJSON with Python — where the footprints these operations consume come from
• Repairing Non-Manifold Meshes with trimesh — the three-dimensional counterpart to ring repair