Computing Overlay Union & Difference with GeoPandas

geopandas.overlay implements the set-theoretic operations — intersection, union, difference, symmetric difference, and identity — that cut two polygon layers against each other into a new one. This guide runs each mode on realistic zoning-versus-floodplain layers and shows the topology hygiene that keeps the output free of slivers and TopologyExceptions. It is for anyone overlaying zoning against hazard areas, land use against parcels, or administrative boundaries against a watershed. It sits under Geometric Intersections & Overlays in Spatial Analysis & Advanced Query Techniques.

Why This Approach / What Goes Wrong

overlay is not the same as a spatial join. A spatial join attaches attributes wherever geometries relate but leaves the original geometries intact; an overlay rebuilds geometry, cutting features against each other and emitting new polygons at every boundary crossing. A single zone that straddles three flood zones comes back as three rows — that row multiplication is the mental model everything downstream (area accounting, attribute reconciliation, deduplication) follows from.

Two things break overlays, and neither is checked for you by the API. The first is invalid input geometry: self-intersecting rings and bad orientation turn a union into a GEOS TopologyException, or — worse — silently corrupt the result before it raises. The second is an unprojected coordinate system: because the engine underneath is GEOS (the same C library behind Shapely geometry operations), every mode assumes a planar Cartesian space, so an overlay computed in EPSG:4326 produces areas in square degrees that grow more distorted with latitude. Validate first with the tooling from Topology Validation & Repair, work in a projected metric CRS, and overlay becomes the precise instrument it is meant to be.

Prerequisites

• geopandas>=0.14 — the overlay API and GeoParquet I/O
• shapely>=2.0 — the vectorized GEOS overlay engine and make_valid
• pyproj>=3.6 — CRS alignment and estimate_utm_crs()

conda install -c conda-forge "geopandas>=0.14" "shapely>=2.0" "pyproj>=3.6"

Installing the whole GEOS/GDAL/PROJ stack from a single conda-forge channel keeps the C libraries aligned; mixing a pip GEOS build with a conda one can shift overlay results subtly between versions.

Step-by-Step Implementation

1. Load two polygon layers and align both to one projected CRS.

GeoPandas will not reproject inside overlay, so the alignment is your responsibility. estimate_utm_crs() picks the correct UTM zone automatically; never overlay in a geographic CRS, and avoid Web Mercator (EPSG:3857) for the cut because its area distortion away from the equator makes every measured result unreliable. The full axis-order and always_xy mechanics live in Coordinate Systems with PyProj.

import geopandas as gpd

# Administrative zoning and FEMA flood hazard areas as GeoDataFrames
zoning = gpd.read_file("zoning.gpkg")
floodplain = gpd.read_file("floodplain.gpkg")

# Bring both onto ONE metric CRS so areas come out in m², not degrees
if not zoning.crs.equals(floodplain.crs):
    floodplain = floodplain.to_crs(zoning.crs)
if zoning.crs.is_geographic:
    metric = zoning.estimate_utm_crs()      # auto-pick the right UTM zone
    zoning = zoning.to_crs(metric)
    floodplain = floodplain.to_crs(metric)

2. Validate both layers before overlaying.

Repair invalid geometry with make_valid() rather than the legacy buffer(0) trick — it preserves geometry type and is deterministic. Deeper repair strategies are covered in Fixing Self-Intersecting Polygons Programmatically.

from shapely import make_valid

for layer in (zoning, floodplain):
    invalid = ~layer.geometry.is_valid
    if invalid.any():
        layer.loc[invalid, "geometry"] = layer.loc[invalid, "geometry"].apply(make_valid)

# Drop degenerate geometry that would poison the cut
zoning = zoning[zoning.geometry.is_valid & ~zoning.geometry.is_empty]
floodplain = floodplain[floodplain.geometry.is_valid & ~floodplain.geometry.is_empty]

3. Intersection — the area in both layers (zoning and flood risk).

Intersection keeps only the shared area and carries attributes from both inputs. Because it emits new polygons, area is meaningful to compute immediately after.

zoning_at_risk = gpd.overlay(zoning, floodplain, how="intersection", keep_geom_type=True)
zoning_at_risk["risk_area_m2"] = zoning_at_risk.geometry.area

4. Difference — zoning outside the floodplain.

Difference subtracts the second layer from the first, so zoning − floodplain is the buildable land clear of the hazard. Order matters: overlay(a, b, "difference") keeps the part of a not covered by b.

zoning_safe = gpd.overlay(zoning, floodplain, how="difference", keep_geom_type=True)

5. Union — every distinct piece from both, attributes preserved.

Union yields all intersection pieces plus the non-overlapping remainder of each input, with each output row tagged by whichever inputs contributed to it (unmatched attributes come back as NaN).

combined = gpd.overlay(zoning, floodplain, how="union", keep_geom_type=True)

6. Symmetric difference and identity — the two remaining modes.

Symmetric difference keeps area in exactly one layer but not both; identity cuts the first layer by the second while staying bounded to the first layer's extent (an intersection plus the leftover difference of the left input).

either_not_both = gpd.overlay(zoning, floodplain, how="symmetric_difference", keep_geom_type=True)
zoning_tagged = gpd.overlay(zoning, floodplain, how="identity", keep_geom_type=True)

7. Reconcile attributes and export.

Overlay multiplies rows, so collapse the pieces back with dissolve (which, unlike a plain groupby, keeps the result a GeoDataFrame and preserves the CRS) and persist to GeoParquet, which carries the CRS in the file metadata.

import pandas as pd

zoning_at_risk["risk_level"] = pd.Categorical(
    zoning_at_risk["risk_level"], categories=["Low", "Medium", "High", "Extreme"], ordered=True
)
reconciled = zoning_at_risk.dissolve(by="zone_id", aggfunc={"risk_level": "max"}).reset_index()
reconciled.to_parquet("zoning_flood_overlay.parquet", compression="zstd")

8. Run all five modes side by side before committing to one.

The fastest way to choose a mode is to look at what each one returns on your actual layers, because the row counts and area totals make the semantics concrete in a way the definitions do not. On a modest pair of layers this takes seconds and settles arguments permanently.

import geopandas as gpd

zoning_area = zoning.geometry.area.sum()
for mode in ("intersection", "union", "difference", "symmetric_difference", "identity"):
    out = gpd.overlay(zoning, floodplain, how=mode, keep_geom_type=True)
    print(f"{mode:<22} rows={len(out):>6,}  "
          f"area={out.geometry.area.sum():>15,.0f} m²  "
          f"cols={len(out.columns)}")

# intersection             rows= 1,842  area=     41,776,940 m²  cols=9
# union                    rows= 4,905  area=    271,118,203 m²  cols=9
# difference               rows= 2,463  area=    142,525,611 m²  cols=5
# symmetric_difference     rows= 3,063  area=    229,341,263 m²  cols=9
# identity                 rows= 4,305  area=    184,302,551 m²  cols=9

Three facts fall straight out of that table. identity returns exactly the left layer's area — it re-cuts zoning without adding or removing ground, which is why it is the right mode for anything that has to reconcile against the original totals. difference returns fewer columns because nothing from the right layer survives; it is the only mode that does not merge the two attribute sets. And union is the only mode whose area exceeds the left layer's, because it also carries the parts of the floodplain that lie outside any zone.

Note the runtime as well as the numbers. union and symmetric_difference are consistently the slowest, because on top of the intersection work they must also compute the leftovers of both inputs. If union is taking minutes and you only ever read the intersected rows, you asked for a more expensive mode than the question needed.

9. Resolve attribute collisions and multipart geometry before writing.

Two mechanical problems show up in the output frame. When both inputs carry a column of the same name, overlay suffixes them _1 and _2 — informative in a notebook, useless in a delivered file — so rename before the cut rather than after, while you still know which layer each column came from. And modes that keep leftovers routinely emit MultiPolygon rows where a single input feature was cut into disconnected pieces; if the consumer expects one geometry per row, explode them and keep the original index so the parts stay traceable.

import geopandas as gpd

# Rename before the cut so no column names collide in the first place
zoning_r = zoning.rename(columns={"name": "zone_name", "updated": "zone_updated"})
flood_r = floodplain.rename(columns={"name": "flood_name", "updated": "flood_updated"})

combined = gpd.overlay(zoning_r, flood_r, how="union", keep_geom_type=True)

# Split multipart rows into single polygons, keeping a trace of the source row
parts = combined.explode(index_parts=True).reset_index(names=["source_row", "part"])
print(parts["part"].max())        # 6  -> one feature was cut into seven pieces

# NaN marks "this layer did not reach here" — make that explicit for consumers
parts["in_floodplain"] = parts["risk_level"].notna()
parts["risk_level"] = parts["risk_level"].fillna("None")

The NaN handling matters more than it looks. In union and symmetric_difference output, a missing right-hand attribute is not missing data — it is a positive statement that the piece lies outside the right layer. Leaving it as NaN means every downstream groupby silently drops those rows, so a summary of flood risk by zone quietly omits all the land that is not at risk. Fill it with an explicit category, and keep the boolean flag alongside it.

Verification

Check that areas reconcile: intersection plus difference should equal the original zoning area to floating-point tolerance. A gap means invalid input or a CRS mismatch distorted the cut.

total_zoning = zoning.geometry.area.sum()
recombined = zoning_at_risk.geometry.area.sum() + zoning_safe.geometry.area.sum()

print(f"Original zoning area:   {total_zoning:,.0f} m²")     # 184,302,551 m²
print(f"Intersection + diff:    {recombined:,.0f} m²")        # 184,302,551 m²

assert abs(total_zoning - recombined) / total_zoning < 1e-6, "Area leak — check validity"
assert zoning_at_risk.geometry.is_valid.all()
assert set(zoning_at_risk.geom_type.unique()) <= {"Polygon", "MultiPolygon"}

Two further identities hold across the modes and catch different bugs. identity must return the left layer's area exactly, since it only re-cuts what was already there — a shortfall means features were dropped by keep_geom_type or by a validity filter. And symmetric_difference plus twice the intersection must equal the union, because the union counts the shared lens once while the other two account for it zero and twice respectively:

identity_out = gpd.overlay(zoning, floodplain, how="identity", keep_geom_type=True)
union_out = gpd.overlay(zoning, floodplain, how="union", keep_geom_type=True)
symdiff_out = gpd.overlay(zoning, floodplain, how="symmetric_difference", keep_geom_type=True)

a_identity = identity_out.geometry.area.sum()
a_union = union_out.geometry.area.sum()
a_symdiff = symdiff_out.geometry.area.sum()
a_inter = zoning_at_risk.geometry.area.sum()

assert abs(a_identity - total_zoning) / total_zoning < 1e-9, "identity lost area"
assert abs((a_symdiff + a_inter) - a_union) / a_union < 1e-9, "union does not reconcile"
print(f"identity {a_identity:,.0f} m² == zoning {total_zoning:,.0f} m²")
# identity 184,302,551 m² == zoning 184,302,551 m²

Run these on a small subset first — a couple of hundred features clipped to one district — so a failure points at a definite pair of layers rather than at forty minutes of wasted compute. If the intersection-plus-difference identity holds on the subset but fails on the full run, the culprit is almost always a handful of pathological features, and zoning[~zoning.geometry.is_valid] names them.

Edge Cases & Debugging

• TopologyException: Input geom 0 is invalid during overlay. Self-intersecting rings; run make_valid on both layers before the cut (step 2).
• Sliver polygons in the output. Hair-thin artifacts where two boundaries almost coincide; drop them by area after the cut with result = result[result.geometry.area > 1.0].
• Areas look astronomical or near zero. You overlaid in EPSG:4326; reproject both layers to a metric CRS first, as in step 1.
• keep_geom_type warning, then missing rows. The cut produced points or lines at tangent boundaries; pass keep_geom_type=True to retain only the polygonal pieces.
• Attribute columns come back as col_1 / col_2. Overlapping column names between inputs get suffixed; rename before overlaying for clarity.
• Memory blows up on a large union. Overlay cost scales with vertex count; simplify() with a domain-appropriate tolerance first, tile the operation by bounding box, or push the exact intersection into an engine built for scale — see Querying GeoParquet with DuckDB Spatial or Scaling with Dask-GeoPandas.
• difference returned the whole left layer unchanged. The arguments are the wrong way round, or the two layers do not actually overlap — check zoning.intersects(floodplain.union_all()).any() before blaming the mode.
• The result is empty and no error was raised. The layers are in different CRSs, so their coordinates never coincide; overlay warns rather than raising, and a warning is easy to miss in a notebook. Assert zoning.crs.equals(floodplain.crs) in the script.
• union output has NaN in the risk column and a groupby lost rows. Those NaNs mean "outside the right layer", not "unknown"; fill them with an explicit category before aggregating, as in step 9.
• Rows multiplied far beyond expectation. The right layer self-overlaps — published hazard data often stores 100-year and 500-year zones as two overlapping polygons — so dissolve it on its classification column first.
• MultiPolygon rows where single polygons were expected. A feature was cut into disconnected pieces; explode(index_parts=True) splits them while keeping the source row traceable.
• unary_union warns on GeoPandas 1.0. It was renamed union_all(); on layers that already tile space without overlaps, union_all(method="coverage") is much faster.
• Different area totals on two machines. Mismatched GEOS builds; pin the stack from one channel and snap both inputs to an explicit precision grid with shapely.set_precision so the cut stops depending on library internals.

Frequently Asked Questions

Which mode should I use to answer "how much of X is inside Y"? intersection gives you the overlapping pieces to measure, but use identity if the answer has to reconcile against X's total. intersection drops every feature of X that never met Y, so the features with zero exposure disappear from the result entirely — and they belong in the denominator. Run identity, measure the tagged pieces, and the totals add up without a second pass.

Why is union so much slower than intersection on the same layers? Because it does strictly more work. Both modes compute the exact intersections, but union then has to derive the leftovers: for each feature, the part not consumed by any intersection, obtained by differencing it against everything it met. That second phase runs over both layers and dominates the runtime. If you never read the rows where one side is NaN, you paid for leftovers you discarded.

Do I need keep_geom_type=True if I am already passing polygons? Yes, because the output can contain non-polygons even when both inputs are clean. Two polygons that touch along an edge or at a corner intersect in a line or a point, and those degenerate pieces come through as geometry-collection members. keep_geom_type=True is the default in GeoPandas 1.0 and keeps only the polygonal parts, which is what any area calculation downstream assumes.

Should I use overlay or clip to restrict a layer to a study area? clip — it is the operation designed for a mask. It cuts the left layer to the outline of the right one, keeps only the left layer's attributes, and does not multiply rows against every mask feature. Use overlay(how="intersection") only when you need to know which right-hand feature each piece fell into and want that feature's attributes carried across.

Can I overlay more than two layers at once? Not in one call — chain the cuts, and expect the row count to multiply at each step. Three layers with a hundred features each can produce thousands of pieces, most of them tiny, so filter and dissolve between steps rather than at the end. Ordering matters for cost too: cut against the most selective layer first so later stages run on fewer features, and put the highest-vertex layer last.