Find Gaps and Overlaps Between Polygons in PyQGIS

Many polygon layers are meant to be coverages: parcels, land-use zones, school catchments, census areas. In a coverage every point of the area belongs to exactly one polygon — no two polygons overlap and there are no holes between them. Real data breaks that rule constantly. Two neighbouring parcels digitized separately overlap by a few centimetres along their shared edge; a zoning update leaves a thin sliver of nothing between two zones; an imported layer double-counts a building plot. Every one of those defects distorts area totals, spatial joins and overlays.

This recipe belongs to Data Quality & Topology Validation. It finds overlaps between pairs of polygons, finds gaps inside the covered area, filters out numerical noise with sensible tolerances, and writes both to layers you can review.

Overlap, gap and sliverFour neighbouring polygons. Two of them overlap along their shared edge, shown as a highlighted lens shape counted twice in any area total. Between two others is a thin gap, an unassigned strip inside the covered area. A tiny sliver of overlap only millimetres wide is numerical noise from separate digitizing and is filtered by an area tolerance.What a broken coverage looks likeoverlapcounted twice in area totalsgapbelongs to no polygonsliver < tolerancenoise: filter by area

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series.
  • A polygon layer in a projected CRS so that areas and tolerances are in metres. If the layer is geographic, reproject first — see choosing a projected CRS for analysis.
  • Valid geometries. Run a validity report first; overlay operations on invalid polygons give unreliable results.

Find overlaps with a spatial index

An overlap is the intersection of two polygons that has area. Testing every pair would be quadratic; a spatial index reduces the candidates to neighbours whose bounding boxes touch.

from qgis.core import (QgsProject, QgsSpatialIndex, QgsFeatureRequest,
                       QgsVectorLayer, QgsFeature, QgsWkbTypes)

parcels = QgsProject.instance().mapLayersByName("parcels")[0]
MIN_AREA = 0.5     # m² — smaller overlaps are digitizing noise

geoms = {f.id(): f.geometry() for f in parcels.getFeatures()}
keys = {f.id(): f["parcel_id"] for f in parcels.getFeatures(
    QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry))}
index = QgsSpatialIndex(parcels.getFeatures(),
                        flags=QgsSpatialIndex.FlagStoreFeatureGeometries)

overlaps = QgsVectorLayer(
    f"MultiPolygon?crs={parcels.crs().authid()}"
    "&field=a:string&field=b:string&field=area_m2:double", "overlaps", "memory")
out = []
for fid, geom in geoms.items():
    engine = QgsGeometry.createGeometryEngine(geom.constGet())
    engine.prepareGeometry()
    for other in index.intersects(geom.boundingBox()):
        if other <= fid:
            continue                       # each pair once
        g2 = geoms[other]
        if not engine.intersects(g2.constGet()):
            continue
        inter = geom.intersection(g2)
        if inter.type() != QgsWkbTypes.PolygonGeometry or inter.area() < MIN_AREA:
            continue
        f = QgsFeature(overlaps.fields())
        f.setAttributes([keys[fid], keys[other], round(inter.area(), 2)])
        inter.convertToMultiType()
        f.setGeometry(inter)
        out.append(f)
overlaps.dataProvider().addFeatures(out)
QgsProject.instance().addMapLayer(overlaps)
print(len(out), "overlaps above", MIN_AREA, "m²")

Breakdown: other <= fid ensures each pair is tested once, halving the work and avoiding duplicate output. A prepared geometry engine makes repeated intersects tests against the same polygon much faster, which matters because each polygon is tested against all its neighbours. Neighbours that only share an edge intersect, but the intersection is a line, not an area; checking the result's type filters those out. The area threshold separates real overlaps from floating-point noise along shared edges. Add from qgis.core import QgsGeometry if running the snippet on its own.

Find gaps by dissolving and differencing

A gap is space inside the covered area that no polygon claims. The trick is to define "inside": dissolve everything into one shape, then holes in the dissolved shape are gaps — but so is the space between islands of coverage, and the outer boundary is not.

Gaps are holes in the dissolved coverageStep one: dissolve all polygons into one multipolygon. Step two: take its interior rings, the holes, each of which is an area surrounded by coverage but claimed by no polygon. Step three: drop holes smaller than the tolerance and keep the rest as gap polygons. The outside of the coverage is never a gap.Dissolve, take holes, filter1. dissolvenative:dissolveall polygons → one2. interior ringsholes of the resultsurrounded, unclaimed3. filterarea ≥ tolerancekeep as gap polygonsthe outside of the coverage is never a gap

import processing
from qgis.core import QgsGeometry, QgsPolygon, QgsLineString

dissolved = processing.run("native:dissolve", {
    "INPUT": parcels, "FIELD": [], "OUTPUT": "memory:"})["OUTPUT"]
cover = next(dissolved.getFeatures()).geometry()

gaps = QgsVectorLayer(
    f"Polygon?crs={parcels.crs().authid()}&field=area_m2:double", "gaps", "memory")
found = []
for part in cover.constParts():
    for i in range(part.numInteriorRings()):
        ring = part.interiorRing(i).clone()
        hole = QgsGeometry(QgsPolygon(ring))
        if hole.area() >= MIN_AREA:
            f = QgsFeature(gaps.fields())
            f.setAttributes([round(hole.area(), 2)])
            f.setGeometry(hole)
            found.append(f)
gaps.dataProvider().addFeatures(found)
QgsProject.instance().addMapLayer(gaps)
print(len(found), "gaps above", MIN_AREA, "m²")

Breakdown: Dissolving without a field merges every polygon into one multipolygon. Each part's interior rings are holes the coverage surrounds but no polygon fills — exactly the definition of a gap. Converting a ring into a polygon makes its area measurable and gives you something to draw. Holes that are intentional — a lake inside a municipal boundary layer, a courtyard inside a building footprint — appear here too; if the dataset has legitimate holes, compare against a layer of known holes or check the attributes of surrounding polygons before treating every result as an error.

Gaps between separate islands

Interior rings only catch gaps that are fully enclosed. A strip of missing coverage that reaches the edge of the area — between two parishes along a county boundary, say — is not a hole in the dissolve. For those, compare against the boundary the coverage should fill.

county = QgsProject.instance().mapLayersByName("county_boundary")[0]
expected = next(county.getFeatures()).geometry()

edge_gaps = expected.difference(cover)
pieces = [QgsGeometry(p.clone()) for p in edge_gaps.constParts()]
real = [p for p in pieces if p.area() >= MIN_AREA]
print(len(real), "uncovered areas inside the county boundary")

Breakdown: Subtracting the dissolved coverage from the expected extent leaves everything the coverage fails to fill, enclosed or not. This needs an authoritative outline — the county boundary, the project area, the municipal extent — but when one exists it is the more complete test. Both layers must be in the same CRS; if they are not, transform one first. Thin pieces along the boundary itself usually mean the coverage and the outline were digitized at different scales, which is a separate conversation from internal gaps.

Tell slivers from real problems

Area alone is a blunt filter: a long, thin overlap of 0.4 m² along a 200 m shared edge is noise, while a compact 0.4 m² overlap may be a duplicated corner point that matters. A thinness measure separates the two.

Thinness separates slivers from blobsTwo overlaps with the same area. A long thin sliver along a shared edge has a large perimeter for its area and a thinness ratio near zero; it is digitizing noise. A compact blob has a ratio near one and is likely a real error such as a duplicated or misplaced vertex. The ratio is four pi times area divided by perimeter squared.Same area, different meaningsliver: ratio ≈ 0.03 — noiseblob: ratio ≈ 0.9 — checkthinness = 4π · area / perimeter²

import math

def thinness(geom):
    p = geom.length()
    return 4 * math.pi * geom.area() / (p * p) if p else 0

for f in overlaps.getFeatures():
    t = thinness(f.geometry())
    kind = "sliver" if t < 0.1 else "compact"
    print(f["a"], f["b"], f["area_m2"], round(t, 3), kind)

Breakdown: The ratio — the isoperimetric quotient — is 1 for a circle and approaches 0 for long, thin shapes. Slivers along shared edges sit well below 0.1. Keeping both the area and the thinness on the output layer lets you style the results by kind and spend review time on the compact ones. A stricter workflow uses a negative buffer: shapes that vanish when buffered inward by a few centimetres are slivers by definition.

Summarise coverage health in one line

Once the overlap and gap layers exist, a short summary makes the result comparable between runs and between datasets. The useful numbers are relative: what share of the covered area is double-counted, and what share is missing.

cover_area = cover.area()
overlap_area = sum(f["area_m2"] for f in overlaps.getFeatures())
gap_area = sum(f["area_m2"] for f in gaps.getFeatures())

print(f"coverage {cover_area / 1e6:.2f} km² | "
      f"overlaps {len(out)} ({overlap_area / cover_area:.4%}) | "
      f"gaps {len(found)} ({gap_area / cover_area:.4%})")

Breakdown: Expressing overlap and gap areas as fractions of the coverage turns raw counts into a quality measure that scales: three hundred overlaps sound alarming until you see they total 0.002 % of the area, while twelve gaps covering 0.4 % of a cadastre is a serious finding. Recording these figures with a date each time the check runs — in a log table, or in the data quality report — shows whether editing practice is improving. If the coverage should equal a known total, such as the official area of a municipality, compare cover_area with it as well; a large difference points to missing polygons rather than small gaps.

Fix what you found

The right fix depends on which polygon is authoritative. For slivers, snapping both layers to a common tolerance often removes them; for real overlaps, one polygon must yield to the other.

# remove overlaps by giving priority to the polygon with the earlier survey date
with edit(parcels):
    for f in overlaps.getFeatures():
        a = next(parcels.getFeatures(f"\"parcel_id\" = '{f['a']}'"))
        b = next(parcels.getFeatures(f"\"parcel_id\" = '{f['b']}'"))
        loser = b if a["surveyed"] <= b["surveyed"] else a
        winner = a if loser is b else b
        parcels.changeGeometry(loser.id(), loser.geometry().difference(winner.geometry()))

Breakdown: Choosing the winner by a rule — earlier survey, higher accuracy class, official over provisional — keeps the fix explainable. Subtracting the winner from the loser removes the overlap without creating a gap, because the overlapping area stays with the winner. Run the check again afterwards: differencing can leave tiny slivers of its own, which snapping geometries to a layer cleans up. Add from qgis.core import edit when running this on its own.

QGIS version compatibility

The spatial index flag FlagStoreFeatureGeometries and createGeometryEngine exist from QGIS 3.0. native:dissolve is stable across 3.x and 4. On QGIS 4, QgsWkbTypes.PolygonGeometry is Qgis.GeometryType.Polygon, and QgsSpatialIndex.FlagStoreFeatureGeometries is reached through QgsSpatialIndex.Flag.FlagStoreFeatureGeometries.

Troubleshooting

  • Thousands of tiny overlaps. The layer was digitized without snapping; raise MIN_AREA or snap first.
  • No gaps found, but the map shows white strips. The strips reach the outer edge; use the expected-boundary difference.
  • Dissolve is slow or fails. Invalid geometries; fix them first, or dissolve in tiles for very large coverages.
  • Areas look wrong. The layer is in a geographic CRS; reproject to a projected system.

Conclusion

Find overlaps pairwise with a spatial index and a prepared engine, find enclosed gaps as interior rings of the dissolved coverage and edge gaps by differencing against the expected boundary, separate slivers from real defects with area and thinness, and fix by a priority rule before re-checking.

Frequently Asked Questions

Does QGIS have a built-in topology checker? Yes, the Topology Checker core plugin checks rules interactively. It is not exposed as a stable Python API, so scripts use the approach shown here.

What tolerance should I use? Tie it to the data's accuracy: if parcels are accurate to 10 cm, overlaps narrower than that along an edge are noise.

Can I check overlaps between two different layers? Yes. Build the index on one layer and iterate the other; drop the other <= fid test.

Is this faster in PostGIS? For millions of polygons, yes — ST_Intersection with a GiST index runs in the database.