Erase Features with Difference in PyQGIS

Many analyses start by taking things away. Developable land is parcels minus protected areas, minus flood zones, minus existing buildings. Usable habitat is forest minus a road buffer. A service area is the 10-minute catchment minus the lake in the middle. The difference operation — "erase" in some other software — cuts the parts of one layer that are covered by another and keeps the rest, with the input's attributes intact.

This recipe belongs to Vector Data Manipulation. It runs native:difference, chains several exclusions, keeps measurements honest, handles multipart results and slivers, uses symmetric difference to find change between two versions of a layer, and keeps large overlays fast.

Difference keeps what is not coveredA parcel polygon overlapped by a protected area. Difference removes the overlapping part and returns the remainder of the parcel with the parcel's attributes; the protected area's attributes are not added. If the protected area splits the parcel, the remainder becomes a multipart polygon. Features entirely covered disappear from the output.Input minus overlay, input's attributes onlyparcelprotected→remainder: multipart, parcel attributes

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series.
  • An input layer and one or more overlay layers in the same projected CRS.
  • Valid geometries; difference is an overlay and fails or misbehaves on invalid input.

Erase one layer with another

The algorithm takes the layer to cut and the layer to cut away.

import processing
from qgis.core import QgsProject

parcels = QgsProject.instance().mapLayersByName("parcels")[0]
protected = QgsProject.instance().mapLayersByName("protected_areas")[0]

remaining = processing.run("native:difference", {
    "INPUT": parcels,
    "OVERLAY": protected,
    "OUTPUT": "memory:parcels_unprotected",
})["OUTPUT"]

print(parcels.featureCount(), "parcels →", remaining.featureCount(), "with unprotected land")
QgsProject.instance().addMapLayer(remaining)

Breakdown: Each input feature has every overlapping overlay feature subtracted from it. The output keeps the input's attributes only — nothing from the overlay is added, which is the main practical difference from intersection. Parcels wholly inside a protected area vanish from the output, so the feature count can drop; parcels split by a protected corridor become multipart features. All overlay features are subtracted regardless of their attributes; filter the overlay first if only some designations should count.

Chain several exclusions

Typical suitability analyses subtract several constraint layers in turn. Running difference repeatedly works, but dissolving the constraints into one overlay first is simpler and often much faster.

Chained or merged constraintsChained: parcels minus protected areas, then minus flood zones, then minus buildings, three overlay passes with intermediate layers. Merged: protected areas, flood zones and buffered buildings merged and dissolved into one constraint layer, then a single difference. The result is the same; the merged route is easier to inspect and usually faster.Three passes, or onechained− protected− flood zones− buildingsthree overlaysmergedmerge + dissolveconstraints → one layer− onceeasier to inspect

constraints = [protected,
               QgsProject.instance().mapLayersByName("flood_zones")[0],
               processing.run("native:buffer", {
                   "INPUT": QgsProject.instance().mapLayersByName("buildings")[0],
                   "DISTANCE": 10, "SEGMENTS": 8, "DISSOLVE": False,
                   "OUTPUT": "memory:"})["OUTPUT"]]

merged = processing.run("native:mergevectorlayers", {
    "LAYERS": constraints, "CRS": parcels.crs(), "OUTPUT": "memory:"})["OUTPUT"]
mask = processing.run("native:dissolve", {
    "INPUT": merged, "FIELD": [], "OUTPUT": "memory:constraints"})["OUTPUT"]

developable = processing.run("native:difference", {
    "INPUT": parcels, "OVERLAY": mask, "OUTPUT": "memory:developable"})["OUTPUT"]

Breakdown: Merging brings the constraint layers into one, reprojecting to the parcels' CRS; dissolving fuses overlapping constraints so the difference subtracts each place once. The buildings are buffered first to exclude a margin around them. A single constraint layer is also a useful product in itself: styled on the map it shows exactly what was excluded and why, which is the first thing anyone reviewing a suitability result asks. Dissolving a very large constraint set can be slow; for national data, dissolve per region.

Keep measurements honest

Like every overlay, difference leaves stale measurement fields behind. A parcel's area_m2 still holds the full parcel area on the remainder.

measured = processing.run("native:fieldcalculator", {
    "INPUT": developable, "FIELD_NAME": "free_m2", "FIELD_TYPE": 0,
    "FIELD_LENGTH": 14, "FIELD_PRECISION": 1,
    "FORMULA": "area($geometry)", "OUTPUT": "memory:developable_measured"})["OUTPUT"]

share = processing.run("native:fieldcalculator", {
    "INPUT": measured, "FIELD_NAME": "free_share", "FIELD_TYPE": 0,
    "FIELD_LENGTH": 6, "FIELD_PRECISION": 3,
    "FORMULA": '"free_m2" / "area_m2"', "OUTPUT": "memory:developable_final"})["OUTPUT"]

Breakdown: Recomputing area on the remainder gives the free area per parcel; dividing by the original area — which the stale field still conveniently holds — gives the share of each parcel left after constraints. That share is often the most useful single number in the result, because it ranks parcels by how constrained they are. area($geometry) is planar in layer units; use $area for ellipsoidal measurement.

Clean up slivers and fragments

Where constraint boundaries almost coincide with parcel boundaries, difference leaves thin slivers and tiny fragments that are not usable land. Removing them makes counts and maps meaningful.

from qgis.core import edit, QgsGeometry

MIN_PART_M2 = 200
with edit(share):
    for f in share.getFeatures():
        g = f.geometry()
        parts = [QgsGeometry(p.clone()) for p in g.constParts()]
        keep = [p for p in parts if p.area() >= MIN_PART_M2]
        if not keep:
            share.deleteFeature(f.id())
        elif len(keep) < len(parts):
            new = QgsGeometry.collectGeometry(keep)
            share.changeGeometry(f.id(), new)

Breakdown: Working part by part removes small fragments from multipart remainders while keeping the substantial parts. Features with nothing above the threshold are deleted outright. The threshold should reflect the purpose — the minimum plot size for development, the minimum patch size for a habitat. Recompute areas afterwards. A negative-then-positive buffer (buffer(-1) then buffer(1)) is an alternative that removes slivers narrower than two metres regardless of area.

Report what was excluded, and why

A suitability result invites the question "why is this parcel not on the list?" Answering it needs the exclusion broken down by constraint. Intersecting the parcels with each constraint separately, before the merged difference, gives that breakdown.

from collections import defaultdict

names = ["protected areas", "flood zones", "building margin"]
excluded = defaultdict(dict)
for name, layer in zip(names, constraints):
    hit = processing.run("native:intersection", {
        "INPUT": parcels, "OVERLAY": layer, "INPUT_FIELDS": ["parcel_id"],
        "OVERLAY_FIELDS": [], "OUTPUT": "memory:"})["OUTPUT"]
    for f in hit.getFeatures():
        excluded[f["parcel_id"]][name] = excluded[f["parcel_id"]].get(name, 0) + f.geometry().area()

for pid, parts in list(excluded.items())[:5]:
    reasons = ", ".join(f"{k} {v:,.0f} m²" for k, v in sorted(parts.items(), key=lambda kv: -kv[1]))
    print(pid, "→", reasons)

Breakdown: One intersection per constraint records how much of each parcel each constraint covers. Because constraints overlap — a flood zone inside a protected area — the per-constraint areas can sum to more than the excluded total; they explain the exclusion rather than partition it. Writing the reasons into a text field on the parcels, or into a separate table joined by parcel_id, lets a reviewer click a parcel and see why it was ruled out. That transparency is what turns a suitability map from a black box into something people trust.

Find change with symmetric difference

Symmetric difference returns the parts of either layer not covered by the other — exactly what changed between two versions of a coverage, such as last year's and this year's forest.

Symmetric difference as changeTwo versions of a forest layer overlap mostly. Symmetric difference returns the areas only in the old version, which were lost, and the areas only in the new version, which were gained, each with attributes from the layer it came from. The unchanged core is dropped.Only in old, or only in newonly in 2025forest lostin bothunchangeddroppedonly in 2026forest gained

old = QgsProject.instance().mapLayersByName("forest_2025")[0]
new = QgsProject.instance().mapLayersByName("forest_2026")[0]

change = processing.run("native:symmetricaldifference", {
    "INPUT": old, "OVERLAY": new, "OVERLAY_FIELDS_PREFIX": "new_",
    "OUTPUT": "memory:forest_change"})["OUTPUT"]

lost = gained = 0.0
for f in change.getFeatures():
    a = f.geometry().area() / 1e4
    if f["new_forest_id"] is None or (hasattr(f["new_forest_id"], "isNull") and f["new_forest_id"].isNull()):
        lost += a
    else:
        gained += a
print(f"lost {lost:.1f} ha, gained {gained:.1f} ha")

Breakdown: Pieces from the old layer carry the old attributes with NULL in the prefixed new fields, and vice versa, which is how the code tells losses from gains. Simple differences in each direction — old minus new, new minus old — give the same result as two layers if that is easier to use. Boundary noise between versions shows up as slivers of fake change, so apply the same sliver filter before reporting change areas.

Performance on large overlays

Difference against a large, detailed overlay can be slow because each input feature is tested against many complex overlay polygons.

from qgis.core import QgsFeatureRequest

simple_mask = processing.run("native:simplifygeometries", {
    "INPUT": mask, "METHOD": 0, "TOLERANCE": 1.0, "OUTPUT": "memory:"})["OUTPUT"]
subdivided = processing.run("native:subdivide", {
    "INPUT": simple_mask, "MAX_NODES": 256, "OUTPUT": "memory:"})["OUTPUT"]
fast = processing.run("native:difference", {
    "INPUT": parcels, "OVERLAY": subdivided, "OUTPUT": "memory:"})["OUTPUT"]

Breakdown: Subdividing a huge dissolved polygon into pieces of at most 256 vertices lets the spatial index find only the nearby pieces for each parcel, and each geometry operation works on small shapes — often an order of magnitude faster on national datasets. Light simplification at a tolerance below the data's accuracy reduces vertex counts further without changing results meaningfully. Both steps are worth trying whenever an overlay takes minutes.

QGIS version compatibility

native:difference, native:symmetricaldifference, native:mergevectorlayers and native:subdivide work on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. Multiple overlay layers in one difference call are supported through native:multidifference from QGIS 3.20. On QGIS 4, NULL attributes arrive as None.

Troubleshooting

  • Features vanished. They were entirely covered by the overlay — expected; count them if it matters.
  • The overlay's attributes are missing. Difference keeps only input attributes; use intersection or union to combine them.
  • Thin slivers everywhere. Boundaries nearly coincide; filter by part area or snap the layers first.
  • The run takes very long. Dissolve, simplify slightly and subdivide the overlay.

Conclusion

Use native:difference to remove covered areas, merge and dissolve several constraints into one overlay before subtracting, recompute area and the remaining share, drop slivers and tiny fragments by part, use symmetric difference to measure change between versions, and subdivide large overlays for speed.

Frequently Asked Questions

Is "erase" the same as difference? Yes. ArcGIS calls it Erase; QGIS calls it Difference.

Can difference be applied to lines? Yes — lines minus polygons gives the line parts outside the polygons, useful for road length outside protected areas.

Does the order of input and overlay matter? Very much: A minus B is not B minus A. Symmetric difference is the order-independent version.

How do I keep only the parts inside instead? Use intersection or clip.