Geometric Intersections & Overlays in Python: Precision Spatial Pipelines

Overlays cut two polygon layers against each other to produce a new one — the set-theoretic core of Spatial Analysis & Advanced Query Techniques. Where a spatial join attaches attributes and leaves geometry intact, an overlay rebuilds geometry at every boundary crossing, so it is the operation you reach for when zoning must be split by flood risk, land use clipped to a watershed, or parcels differenced against a right-of-way. This guide covers the predicate-and-operation model, the topology hygiene that keeps output clean, and the scaling path onto DuckDB Spatial Analytics; the mode-by-mode recipe lives in Computing Overlay Union & Difference with GeoPandas, and distance-filtered variants in Proximity & Buffer Analysis.

The four set-theoretic overlay modes One pair of overlapping polygon layers, A in cyan and B in violet, drawn four times. Each panel shades the region that survives a different overlay mode: intersection keeps only the shared lens, union keeps everything, difference keeps A minus the overlap, and symmetric difference keeps everything except the overlap. One pair of layers, four overlay modes A B Intersection A ∩ B A B Union A ∪ B A B Difference A − B A B Symmetric diff. A △ B
Each panel keeps the shaded pieces of the same A/B pair: intersection, union, difference, and symmetric_difference — with identity (A carrying B's attributes) as the fifth mode. Validate geometry and share one metric CRS before any cut.

Architecture & Data Structures

An overlay is a constructive operation over two GeoPandas DataFrames: geopandas.overlay(left, right, how=...) computes the planar intersection graph of both geometry columns, then reassembles the pieces according to the how mode and carries attributes from whichever inputs contributed to each piece. The five modes — intersection, union, difference, symmetric_difference, and identity — map directly onto set theory, and the engine underneath is GEOS, the same C library that powers Shapely geometry operations. Because GEOS assumes a planar (Cartesian) coordinate space, the correctness of every mode depends on two invariants the API does not check for you: both layers must be valid and both must share a single projected CRS.

import geopandas as gpd

# Two polygon layers: administrative zoning and FEMA flood hazard areas
zones = gpd.read_parquet("zoning_zones.parquet")
floodplains = gpd.read_parquet("fema_floodplains.parquet")

# The overlay API surface — one call, five set-theoretic modes
zoning_at_risk = gpd.overlay(zones, floodplains, how="intersection")
print(zoning_at_risk.columns.tolist())   # attributes from BOTH inputs, suffixed on clash
print(zoning_at_risk.geom_type.unique()) # ['Polygon'] once keep_geom_type is respected

The output is a new GeoDataFrame whose rows are the intersection pieces, not the original features — a single zone that straddles three flood zones becomes three rows. That row multiplication is the mental model to hold onto: overlay is a cut, and everything downstream (area accounting, attribute reconciliation, deduplication) follows from it.

Underneath the single overlay() call are three distinct phases, and knowing which one is misbehaving is most of the debugging. First a bounding-box join over the spatial index produces candidate pairs. Second, GEOS computes an exact constructive intersection for each surviving pair. Third — and only for union, symmetric_difference and identity — GeoPandas computes the leftovers: the parts of each input not consumed by any intersection, obtained by differencing each feature against the union of everything it met. That third phase is why union is routinely five to ten times slower than intersection on the same pair of layers, and why a union can exhaust memory on inputs an intersection handles comfortably. If a mode is slow, check whether you actually need the leftovers.

The API has three neighbouring operations that people reach for interchangeably and should not. gpd.overlay is the table-to-table cut described here. GeoSeries.intersection is element-wise: it aligns two series on the DataFrame index and intersects row i with row i, which is what you want for a pre-paired frame from a spatial join and never what you want for two independent layers — mismatched indexes yield rows of None rather than an error. gpd.clip is a mask operation: it cuts one layer to the outline of another and keeps only the left layer's attributes, so it is the right call for "restrict this to the study area" and the wrong one when you need to know which right-hand feature each piece fell into.

import geopandas as gpd

parcels = gpd.read_parquet("parcels.parquet")
hazard_zones = gpd.read_parquet("fema_floodplains.parquet")
study_area = gpd.read_file("city_boundary.gpkg")

# Three operations, three different results — pick deliberately
cut = gpd.overlay(parcels, hazard_zones, how="intersection")  # every pair, both attribute sets
masked = gpd.clip(parcels, study_area)                        # trimmed to an outline, left attrs only

# Element-wise: only meaningful once the two frames are already paired row-for-row
pairs = gpd.sjoin(parcels, hazard_zones, predicate="intersects")
paired = pairs.geometry.intersection(
    hazard_zones.geometry.loc[pairs["index_right"]].reset_index(drop=True).set_axis(pairs.index)
)

Two API changes in GeoPandas 1.0 affect overlay code written against 0.13 or earlier. GeoSeries.unary_union was renamed union_all() — the old name still works but warns, and the new one accepts a method argument, where method="coverage" is dramatically faster on layers that are already an edge-matched partition of space with no overlaps. And keep_geom_type now defaults to True in overlay, so code that silently relied on point and line fragments coming through will find them gone; that is almost always the correct new behaviour, but it changes row counts on upgrade. Both belong in an upgrade checklist alongside the vectorization changes described in Shapely 1.x vs Shapely 2 Vectorization.

Environment Configuration & Dependency Resolution

Overlay semantics changed materially with Shapely 2.0, which rewired the constructive operations onto vectorized GEOS calls. Pin a version floor of shapely>=2.0 so make_valid and the vectorized overlay engine are available, and keep the GEOS/GDAL/PROJ stack aligned by installing from a single channel — mixing pip and conda GEOS builds silently shifts predicate and overlay results between versions.

# conda-forge keeps GEOS/GDAL/PROJ consistent across the whole stack
conda create -n overlays -c conda-forge \
  "python=3.12" "geopandas>=1.0" "shapely>=2.0" \
  "pyproj>=3.6" "duckdb>=1.1" "pyarrow>=14"
conda activate overlays
# Fail loudly if the vectorized constructive layer is unavailable
import shapely
assert tuple(int(p) for p in shapely.__version__.split(".")[:2]) >= (2, 0), (
    "Shapely 2.0+ is required for vectorized make_valid and overlay."
)

The GEOS version underneath matters as much as the Shapely version, because the overlay algorithm itself was replaced. GEOS 3.9 introduced OverlayNG, a rewrite that made constructive operations robust by snapping intersections onto a computed precision grid instead of failing outright. In practice that turned a large class of hard TopologyException crashes into results — usually the right ones, occasionally results with a vertex nudged by a nanometre. Two consequences follow. First, an overlay that raised on an old stack and now succeeds is not necessarily proof that the input got cleaner; it may just be that OverlayNG papered over the same defect. Second, output is only bit-for-bit reproducible across machines if the GEOS versions match, which is why a pinned environment matters for anything whose numbers get published.

import shapely

print(shapely.geos_version)          # (3, 12, 1)
print(shapely.geos_capi_version_string)
assert shapely.geos_version >= (3, 9, 0), "OverlayNG requires GEOS 3.9+"

If you need determinism stronger than "same version, same answer", impose the precision model yourself rather than letting GEOS choose one. shapely.set_precision snaps coordinates onto an explicit grid before the cut, so every machine sees identical input and the overlay becomes reproducible by construction — the technique is developed further in Snapping and Simplifying Polygons Without Gaps.

For datasets past a comfortable in-memory size, add DuckDB's spatial extension (covered below) rather than a heavier server dependency — it ships as a single wheel with no system libraries to reconcile.

Vectorized Operations & Core Workflow

A production overlay is a fixed sequence: standardize CRS, validate topology, build a spatial index, cut, reconcile attributes, and export. Each stage is vectorized — the C engine processes the whole geometry array at once, so never loop rows when a native method exists. The block below runs end to end on two polygon layers.

The production overlay pipeline Seven left-to-right stages, each labelled with the GeoPandas call that performs it: align CRS with to_crs, validate with make_valid, build a spatial index with sindex, prefilter candidate pairs with sjoin, cut with overlay, reconcile attributes with dissolve, and export with to_parquet. A fixed sequence: two raw layers to one web-ready overlay 1 Align CRS to_crs() 2 Validate make_valid() 3 Index sindex 4 Prefilter sjoin() 5 Cut overlay() 6 Reconcile dissolve() 7 Export to_parquet()
The overlay() cut (stage 5) is the expensive step; stages 1–4 make it correct and cheap, and stages 6–7 reconcile and persist the result.
import geopandas as gpd
import pandas as pd
from shapely import make_valid

# 1. Ingest and align both layers to ONE projected CRS (metres, not degrees)
zones = gpd.read_parquet("zoning_zones.parquet")
floodplains = gpd.read_parquet("fema_floodplains.parquet")

if not zones.crs.equals(floodplains.crs):
    floodplains = floodplains.to_crs(zones.crs)
if zones.crs.is_geographic:
    metric = zones.estimate_utm_crs()          # auto-pick the right UTM zone
    zones = zones.to_crs(metric)
    floodplains = floodplains.to_crs(metric)

# 2. Validate and drop degenerate geometry before the cut
for layer in (zones, floodplains):
    invalid = ~layer.geometry.is_valid
    layer.loc[invalid, "geometry"] = layer.loc[invalid, "geometry"].apply(make_valid)
zones = zones[zones.geometry.is_valid & ~zones.geometry.is_empty]
floodplains = floodplains[floodplains.geometry.is_valid & ~floodplains.geometry.is_empty]

# 3. Cut: intersection keeps only the shared area, attributes from both inputs
overlay = gpd.overlay(zones, floodplains, how="intersection", keep_geom_type=True)

# 4. Reconcile attributes — highest flood risk wins per zone
overlay["risk_level"] = pd.Categorical(
    overlay["risk_level"], categories=["Low", "Medium", "High", "Extreme"], ordered=True
)
reconciled = (
    overlay.dissolve(by="zone_id", aggfunc={"risk_level": "max"})
    .reset_index()
)

# 5. Export to a columnar format that carries the CRS with the file
reconciled.to_parquet("zoning_flood_overlay.parquet", compression="zstd")

Using dissolve rather than a plain groupby keeps the result a GeoDataFrame and merges the per-zone pieces back into one geometry, so the CRS and geometry column survive the aggregation intact.

Geometry / Data Processing Details

Real-world layers arrive with self-intersections, ring-orientation errors, and duplicate vertices that turn a union into a GEOS TopologyException — or, worse, silently corrupt the result before it raises. Validate at ingestion, not after failure. Prefer make_valid() over the legacy buffer(0) trick: it preserves geometry type and is deterministic. The deeper repair strategies live in Topology Validation & Repair and Fixing Self-Intersecting Polygons Programmatically.

Two data-processing concerns dominate once geometry is valid: sliver control and candidate pruning. Slivers are the hair-thin polygons produced where two boundaries almost coincide; drop them by area after the cut. Candidate pruning matters because a naive overlay is O(n×m) — every left feature tested against every right feature. GeoPandas builds an R-tree via sindex (backed by Shapely 2.0's STRtree) to pre-filter with a bounding-box join before the expensive exact intersection, the same indexing discipline behind Nearest Neighbor & KD-Tree Search.

Sliver polygons before and after an area threshold On the left, a zoning block and a floodplain block digitised from different sources have boundaries that almost but not exactly coincide. The raw intersection returns three hair-thin polygons along that seam, each well under a square metre, drawn in red between the two blocks. On the right, the same cut after filtering the result by area: the seam remains but the sliver rows are gone, leaving only the substantive pieces. Where slivers come from, and how they leave Raw cut gpd.overlay(how="intersection") zones floodplains 3 sliver rows, each under 1 m² After the area threshold result[result.area > 1.0] zones floodplains seam kept, sliver rows discarded A bounding-box sjoin prunes the candidate pairs first; an area threshold removes what the cut leaves behind.
Near-coincident boundaries from two different digitising sources survive the cut as hair-thin polygons — filtering the output by area is the deterministic way to drop them.
import geopandas as gpd

# Prune candidate pairs with a bounding-box sjoin before the exact cut.
# A small positive buffer captures near-touching features without generating slivers.
zones_probe = gpd.GeoDataFrame(
    zones[["zone_id"]], geometry=zones.buffer(0.5), crs=zones.crs
)
candidates = gpd.sjoin(zones_probe, floodplains, how="inner", predicate="intersects")

# Restrict the overlay to the touched features only, then drop sliver artefacts
touched = zones[zones["zone_id"].isin(candidates["zone_id"])]
result = gpd.overlay(touched, floodplains, how="intersection", keep_geom_type=True)
result = result[result.geometry.area > 1.0]   # discard sub-1 m² slivers

For layers too large for a comfortable GeoDataFrame, push the exact intersection into DuckDB's spatial extension, which runs ST_Intersects over columnar data in C++ and reads GeoParquet directly. DuckDB treats geometry as WKB/WKT, so hand it well-known text and reconstruct on the way out.

import duckdb

con = duckdb.connect()
con.install_extension("spatial"); con.load_extension("spatial")

# SQL-native intersection over GeoParquet — no server, no full materialization
con.sql("""
    SELECT z.zone_id, f.risk_level, f.elevation
    FROM 'zoning_zones.parquet' z
    JOIN 'fema_floodplains.parquet' f
      ON ST_Intersects(z.geom, f.geom)
""").df()

Diagnosing slivers properly. A raw area threshold is a blunt instrument: it deletes hair-thin seam artefacts and also deletes the genuine 0.8 m² triangle of floodplain that clips a corner of a real parcel. The distinguishing property of a sliver is not that it is small — it is that it is thin, which means its area is tiny relative to its perimeter. The Polsby-Popper compactness ratio, 4π·area / perimeter², is 1 for a circle, around 0.79 for a square, and collapses toward zero as a polygon degenerates into a ribbon. Filtering on compactness and area removes seam artefacts while keeping small-but-chunky pieces:

import numpy as np

result["compactness"] = (
    4 * np.pi * result.geometry.area / result.geometry.length.pow(2)
)
slivers = result[(result["compactness"] < 0.02) & (result.geometry.area < 5.0)]
clean = result.drop(index=slivers.index)

print(f"{len(slivers)} sliver rows removed, {slivers.geometry.area.sum():.2f} m² total")
# 1,284 sliver rows removed, 41.06 m² total

Printing the total discarded area is the part worth keeping. If the slivers you dropped account for 41 m² out of 184 km², you have removed digitising noise; if they account for several hectares, the two layers disagree substantively and the right response is to reconcile the sources, not to delete the evidence.

A deterministic alternative to filtering. Slivers exist because two boundaries that represent the same real-world line were digitised to slightly different coordinates. Snapping both layers onto a common grid coarser than the disagreement makes the boundaries coincide exactly, so the slivers are never generated in the first place. Choose the grid from the data's true positional accuracy — 0.01 m for survey-grade cadastre, 1 m for a hand-digitised hazard layer — and never finer than the least accurate input:

from shapely import set_precision

# Snap both layers to a 0.05 m grid so near-coincident boundaries become identical
zones["geometry"] = set_precision(zones.geometry.values, grid_size=0.05)
floodplains["geometry"] = set_precision(floodplains.geometry.values, grid_size=0.05)

clean_cut = gpd.overlay(zones, floodplains, how="intersection", keep_geom_type=True)

Snapping can collapse a genuinely narrow feature — a 3 cm alley on a 5 cm grid — into an empty geometry, so re-run the validity and emptiness filter afterwards. It is the more honest fix of the two, because it states the positional tolerance explicitly instead of inferring it from the artefacts.

What the cost is actually driven by. Overlay time scales with total vertex count far more than with row count: a hundred thousand rectangles cut faster than five thousand coastline polygons. Measure the input the way the engine sees it before you decide whether an overlay is feasible.

import shapely

for name, layer in (("zones", zones), ("floodplains", floodplains)):
    verts = shapely.get_num_coordinates(layer.geometry.values)
    print(f"{name}: {len(layer):>7,} rows  {verts.sum():>10,} vertices  "
          f"max {verts.max():,} on one feature")
# zones:       4,812 rows    1,204,318 vertices  max 41,209 on one feature
# floodplains:   937 rows   11,882,004 vertices  max 2,904,551 on one feature

A single feature holding three million vertices — a merged coastline or a dissolved national boundary — will dominate the runtime on its own, and simplifying only that feature is usually a bigger win than simplifying the whole layer. Past roughly ten million vertices per side, a single-process overlay stops being interactive: expect tens of minutes and a memory footprint several times the input, because OverlayNG builds an edge graph over both inputs simultaneously. The escape hatches, in increasing order of effort, are simplifying to the tolerance the output actually needs, tiling the cut by bounding box and concatenating the pieces, pushing the exact intersection into DuckDB Spatial Analytics, or partitioning it across workers with Scaling with Dask-GeoPandas.

Tiling deserves a note because it is easy to get subtly wrong. Cutting by tile and concatenating works only if each feature is assigned to exactly one tile; clipping both layers to overlapping tiles double-counts everything in the overlap region and inflates every area you compute afterwards. Clip the left layer to disjoint tiles, take the right layer as everything intersecting that tile, and the pieces reassemble without duplication:

import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import box

minx, miny, maxx, maxy = zones.total_bounds
edges_x = np.linspace(minx, maxx, 5)      # 4 x 4 = 16 disjoint tiles
edges_y = np.linspace(miny, maxy, 5)

pieces = []
for i in range(len(edges_x) - 1):
    for j in range(len(edges_y) - 1):
        tile = box(edges_x[i], edges_y[j], edges_x[i + 1], edges_y[j + 1])
        left = gpd.clip(zones, tile)                       # disjoint: no double counting
        if left.empty:
            continue
        right = floodplains[floodplains.intersects(tile)]   # NOT clipped — full features
        pieces.append(gpd.overlay(left, right, how="intersection", keep_geom_type=True))

tiled = gpd.GeoDataFrame(pd.concat(pieces, ignore_index=True), crs=zones.crs)

The result is geometrically identical to a single-shot overlay, but peak memory is bounded by the largest tile rather than the whole dataset. Features that straddle a tile edge come back as two adjacent rows, so dissolve on the feature id afterwards if downstream code expects one row per input feature.

CRS Alignment & Projection Pipeline

Overlay is a planar operation, so a shared, projected CRS is not optional — it is the difference between a correct cut and silent nonsense. Areas computed from an overlay in unprojected WGS84 (EPSG:4326) are square degrees, and the distortion grows with latitude, so a "shared area" figure at 60° N is meaningless. Align both layers to a metric CRS first; estimate_utm_crs() picks the right UTM zone automatically, and the full mechanics — axis order, always_xy, EPSG versus PROJ strings — live in Coordinate Systems with PyProj.

import geopandas as gpd

def align_for_overlay(left: gpd.GeoDataFrame, right: gpd.GeoDataFrame) -> tuple:
    """Reproject both layers to one metric CRS before any overlay."""
    if not left.crs.equals(right.crs):
        right = right.to_crs(left.crs)
    if left.crs.is_geographic:
        metric = left.estimate_utm_crs()
        left, right = left.to_crs(metric), right.to_crs(metric)
    assert left.crs.axis_info[0].unit_name == "metre", "Overlay CRS must be metric."
    return left, right

estimate_utm_crs() is the right default and the wrong answer for one specific case: extents that span more than a UTM zone or two. UTM is a set of 6°-wide strips, and scale error grows away from each strip's central meridian — tolerable inside one zone, but by the time a national dataset stretches across four zones, forcing it all into one of them inflates areas at the edges by well over a percent. That is invisible on a map and fatal to an area-based statistic. For anything wider than roughly two zones, use an equal-area projection built for the region instead: ETRS89-LAEA (EPSG:3035) for Europe, an Albers Equal Area for the conterminous United States, or a custom Lambert Azimuthal Equal Area centred on the data.

import geopandas as gpd

def metric_crs_for(layer: gpd.GeoDataFrame) -> str:
    """UTM inside one or two zones; an equal-area CRS for anything wider."""
    minx, _, maxx, _ = layer.to_crs("EPSG:4326").total_bounds
    zones_spanned = int(maxx // 6) - int(minx // 6) + 1
    if zones_spanned <= 2:
        return layer.estimate_utm_crs().to_string()
    return "EPSG:3035"           # ETRS89-LAEA — swap for your region's equal-area CRS

print(metric_crs_for(zones))     # EPSG:32632 for a Turin-sized extent

Two further alignment traps show up specifically in overlay work rather than in reprojection generally. The first is the antimeridian: a layer crossing 180° longitude comes back from a naive reprojection with polygons that wrap the entire globe, and any overlay against them produces nonsense over the whole extent. Split those features at the antimeridian before projecting. The second is subtler — two layers that both report EPSG:4326 but were realized against different datums, say WGS84 versus a national ETRS realization. They align to within a metre or two, which is enough to produce a continuous ribbon of slivers along every shared boundary and not enough to look wrong. If a supposedly edge-matched pair produces slivers everywhere, suspect the datum before the digitising, and check the transformation pipeline as described in Coordinate Reference System Transformations.

Production Export & Integration

from shapely import set_precision

# Web-ready export: reproject to WGS84, snap to a 6-decimal grid, write GeoJSON
web_ready = reconciled.to_crs("EPSG:4326").copy()
web_ready["geometry"] = web_ready.geometry.apply(lambda g: set_precision(g, grid_size=1e-6))
web_ready.to_file("overlay_dashboard.geojson", driver="GeoJSON")

# Optional vector tiles (run externally):
# tippecanoe -o overlay.mbtiles -z 14 -Z 0 --drop-densest-as-needed overlay_dashboard.geojson

Area-weighted apportionment: the analysis most overlays exist for. Cutting the layers is rarely the deliverable. The question behind it is usually "how much of this falls inside that" — how much population lives in the flood zone, how much assessed value sits in a proposed right-of-way, how much cropland falls in each watershed. The overlay supplies the geometry; apportionment turns it into a number, on the assumption that the quantity is spread evenly across each source feature. That assumption is the whole method, and it is worth stating out loud: it is defensible for cropland across a field and shaky for population across a census tract that is half industrial park.

import geopandas as gpd

# census_tracts carries a population count; floodplains delimits the hazard
tracts = gpd.read_parquet("census_tracts.parquet")
floodplains = gpd.read_parquet("fema_floodplains.parquet")
tracts, floodplains = align_for_overlay(tracts, floodplains)

tracts["tract_area_m2"] = tracts.geometry.area          # record BEFORE the cut

flooded = gpd.overlay(tracts, floodplains, how="intersection", keep_geom_type=True)
flooded["piece_area_m2"] = flooded.geometry.area
flooded["area_share"] = flooded["piece_area_m2"] / flooded["tract_area_m2"]
flooded["population_at_risk"] = flooded["population"] * flooded["area_share"]

by_risk = flooded.groupby("risk_level", observed=True)["population_at_risk"].sum().round(0)
print(by_risk)
# risk_level
# High       18,442.0
# Medium     51,097.0
# Extreme     3,806.0

Compute the denominator before the overlay, never after — after the cut, tracts.geometry.area no longer refers to a whole tract, and dividing by the piece area gives every piece a share of 1.0. Two sanity checks catch nearly every apportionment bug: the per-feature shares must sum to at most 1.0 (more than that means the right layer self-overlaps and the same ground is being counted twice), and the apportioned total must not exceed the source total. Where the even-spread assumption is too weak to accept, the honest upgrade is to weight by a raster of built-up area or night-time lights instead of by bare area, which is the sampling problem covered in Zonal Statistics & Raster Sampling.

# The two checks that catch double counting and leakage
shares = flooded.groupby("tract_id", observed=True)["area_share"].sum()
assert (shares <= 1.0 + 1e-9).all(), "Right layer self-overlaps — dissolve it first"
assert flooded["population_at_risk"].sum() <= tracts["population"].sum() + 1

If the first assertion fires, the fix is upstream: dissolve the right-hand layer on its classification column so its features no longer overlap each other, then re-cut. Overlapping hazard polygons are extremely common in published hazard data, where a 100-year and a 500-year zone are stored as two overlapping features rather than as nested rings.

Windows / Platform Edge Cases & Debugging

Frequently Asked Questions

When should I use overlay instead of a spatial join? Use a spatial join when the geometry you want to keep already exists and you only need to attach attributes from another layer — "which district is each sensor in" leaves the sensor points untouched. Use overlay when you need geometry that does not exist yet, because the answer is a portion of a feature: the part of a parcel inside the flood zone, the part of a district outside the service area. If your next step is to measure the area or length of the shared part, you need an overlay; the distinction is worked through from the join side in Spatial Join vs Attribute Join in GeoPandas.

Which mode do I want when I only care about one layer's footprint? identity. It keeps every part of the left layer — inside and outside the right one — and tags the overlapping parts with the right layer's attributes, so the total area is unchanged and every original feature is still represented. intersection silently drops features that never met the right layer, which is what makes it dangerous for accounting: the rows that vanish are exactly the ones with zero exposure, and they still belong in the denominator.

Should I clean the data or just repair it during the overlay? Clean upstream, every time. make_valid inside the overlay script fixes the run in front of you and hides a defect that will reappear in every other job reading the same layer, usually differently. Repair once at ingestion, write the repaired layer to GeoParquet, and let every downstream job consume the clean copy — the workflow argued for in Topology Validation & Repair. Keep a validity assertion in the overlay script as a tripwire, not as the repair.

How big is too big for a GeoPandas overlay? Count vertices, not rows. Under about a million vertices per side, an overlay is interactive on a laptop. Around ten million per side you are looking at tens of minutes and a memory footprint several times the input. Past that, tile the cut, simplify to the tolerance the deliverable actually needs, or move the exact intersection into DuckDB or Dask-GeoPandas. A single pathological feature with millions of vertices can push you over that line on its own, so check the per-feature maximum before blaming the row count.

Is a sliver always something to delete? No — and deciding by size alone is how real slivers of geometry get thrown away. A sliver is thin, not merely small, so filter on a compactness ratio together with area, and always print the total area you discarded. Sub-square-metre noise along a seam is safe to drop; several hectares means the two layers genuinely disagree about where a boundary is, and the fix belongs in the source data rather than in a where clause.

Can I overlay a polygon layer against lines or points? gpd.overlay expects polygonal inputs on both sides and will drop or complain about mixed types once keep_geom_type is honoured. For lines cut by polygons — a road network split at district boundaries — use gpd.clip per polygon or GeoSeries.intersection on paired frames, which return line pieces cleanly. For points, no cutting is needed at all: a spatial join with predicate="within" gives you the same answer with far less work.

Related