Union Overlay of Layers in PyQGIS

Intersection keeps only where two layers overlap; difference keeps only where they do not. Union keeps everything: it splits both layers along each other's boundaries and returns every resulting piece, each carrying the attributes of whichever input features cover it. The result is a single layer that partitions the combined area into zones with a consistent combination of attributes — the basis for cross-tabulations such as land use by zoning designation, or for comparing two classifications of the same area.

This recipe belongs to Vector Data Manipulation. It runs native:union on two layers, classifies the pieces by origin, recomputes areas, builds a cross-tabulation, uses self-union to count overlaps within one layer, and handles the slivers and NULLs that union inevitably produces.

Union keeps every pieceTwo overlapping polygons, A and B. Union returns three kinds of piece: parts only in A, with A's attributes and NULLs for B; parts only in B, with B's attributes and NULLs for A; and the overlap, with attributes from both. Together the pieces cover exactly the combined area of A and B with no overlaps.Everything, split by everythingA onlyB onlybothA onlyA attributes, B fields NULLbothattributes from A and BB onlyB attributes, A fields NULL

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series.
  • Two polygon layers in the same projected CRS, with valid geometries.
  • Distinct field names, or a prefix to tell them apart — union combines both attribute tables.

Run the union

The algorithm takes an input and an overlay layer and returns pieces with fields from both.

import processing
from qgis.core import QgsProject

landuse = QgsProject.instance().mapLayersByName("landuse_2026")[0]
zoning = QgsProject.instance().mapLayersByName("zoning_plan")[0]

pieces = processing.run("native:union", {
    "INPUT": landuse,
    "OVERLAY": zoning,
    "OVERLAY_FIELDS_PREFIX": "zone_",
    "OUTPUT": "memory:landuse_u_zoning",
})["OUTPUT"]

print(pieces.featureCount(), "pieces;", pieces.fields().names())
QgsProject.instance().addMapLayer(pieces)

Breakdown: The prefix renames every overlay field — designation becomes zone_designation — so fields with the same name in both layers do not collide and the origin of every column is obvious. The output covers the combined extent of both layers. Unlike intersection, there is no field selection parameter, so very wide inputs produce very wide outputs; drop columns afterwards or refactor the inputs first with rename and reorder fields.

Classify pieces by origin

Every piece came from A only, B only, or both. Recording which turns the union into something easy to style and filter.

from qgis.core import QgsField, edit
from qgis.PyQt.QtCore import QVariant

def is_null(v):
    return v is None or (hasattr(v, "isNull") and v.isNull())

prov = pieces.dataProvider()
prov.addAttributes([QgsField("origin", QVariant.String, len=6)])
pieces.updateFields()
idx = pieces.fields().indexOf("origin")

changes = {}
for f in pieces.getFeatures():
    in_a = not is_null(f["landuse_id"])
    in_b = not is_null(f["zone_zone_id"])
    changes[f.id()] = {idx: "both" if in_a and in_b else "A" if in_a else "B"}
prov.changeAttributeValues(changes)

Breakdown: A key field from each input — here each layer's own id — is non-NULL exactly when the piece is covered by that layer. Classifying on keys rather than on optional attributes avoids misclassifying pieces where an attribute is legitimately empty. Pieces covered only by the zoning plan show where the land-use survey has gaps; pieces covered only by land use show areas outside any zone. Both are often more interesting than the overlap.

Recompute areas

As with every overlay, area fields inherited from the inputs describe the original features, not the pieces.

Area on pieces, not on originalsAfter union, each piece inherits area fields from the land-use polygon and the zone that cover it, describing those whole polygons. A new field computed from the piece's own geometry is the only area that can be summed. Summing it over all pieces equals the area of the combined footprint.Sum only what you measured on the pieceinheritedlanduse_area_m2zone_area_m2whole originalsrecomputedpiece_m2 = area($geometry)sums correctly

measured = processing.run("native:fieldcalculator", {
    "INPUT": pieces, "FIELD_NAME": "piece_m2", "FIELD_TYPE": 0,
    "FIELD_LENGTH": 14, "FIELD_PRECISION": 1, "FORMULA": "area($geometry)",
    "OUTPUT": "memory:union_measured"})["OUTPUT"]
total = sum(f["piece_m2"] for f in measured.getFeatures())
print(f"combined footprint: {total / 1e6:.2f} km²")

Breakdown: The sum of piece areas is the area of the combined footprint of both layers, with every location counted once — a useful check, since it should equal the area of a dissolve of both layers together. Inherited area fields are harmless as long as nobody sums them; dropping them avoids temptation.

Cross-tabulate the result

The classic use of a union is a matrix: how much of each land-use class lies in each zoning designation, including land in no zone and zones with no surveyed land use.

from collections import defaultdict

matrix = defaultdict(float)
for f in measured.getFeatures():
    lu = f["landuse_class"] if not is_null(f["landuse_class"]) else "(no survey)"
    zn = f["zone_designation"] if not is_null(f["zone_designation"]) else "(no zone)"
    matrix[(lu, zn)] += f["piece_m2"] / 1e4

zones = sorted({z for _, z in matrix})
print("land use".ljust(18) + "".join(z[:12].rjust(13) for z in zones))
for lu in sorted({l for l, _ in matrix}):
    print(lu[:18].ljust(18) + "".join(f"{matrix.get((lu, z), 0):13.1f}" for z in zones))

Breakdown: Replacing NULLs with explicit labels such as "(no zone)" keeps uncovered areas in the table instead of silently dropping them. The printed matrix — hectares of each land-use class per zoning designation — shows, for example, how much residential land use falls in zones designated for industry. For a spreadsheet, build the same table as a pandas pivot as in analysing an attribute table with pandas.

Self-union: count overlaps within one layer

Running union on a single layer, with no overlay, splits it wherever its own features overlap. Each resulting piece appears once per feature covering it, which makes counting overlap depth straightforward.

Overlap depth from a self-unionThree overlapping buffer polygons in one layer. A self-union splits them into pieces at every boundary. Pieces with identical geometry are grouped; the size of each group is the number of features covering that piece: one where a single buffer lies, two where two overlap, three in the centre where all three overlap.How many features cover each place?3112self-uniongroup identical piecescount = overlap depth

self_u = processing.run("native:union", {
    "INPUT": QgsProject.instance().mapLayersByName("school_catchments")[0],
    "OVERLAY": None, "OUTPUT": "memory:"})["OUTPUT"]

groups = defaultdict(int)
for f in self_u.getFeatures():
    groups[f.geometry().asWkb().toHex().data()] += 1
print(max(groups.values()), "catchments overlap at the most crowded place")

Breakdown: With no overlay, union resolves the layer against itself, producing duplicate pieces where features overlap — one per covering feature. Grouping pieces by identical geometry (here by their WKB) and counting gives the overlap depth for each piece: how many school catchments cover each address area, how many licence areas claim each square metre. Rounding coordinates before hashing guards against floating-point differences.

Compare two versions of a classification

Union is also the tool for comparing two classifications of the same area — this year's land-use survey against last year's, a contractor's habitat map against your own. The union's pieces carry both classes, so agreement is simply the share of area where they match.

old = QgsProject.instance().mapLayersByName("landuse_2025")[0]
cmp = processing.run("native:union", {
    "INPUT": landuse, "OVERLAY": old, "OVERLAY_FIELDS_PREFIX": "old_",
    "OUTPUT": "memory:"})["OUTPUT"]

same = changed = 0.0
transitions = defaultdict(float)
for f in cmp.getFeatures():
    a = f.geometry().area()
    new_cls, old_cls = f["landuse_class"], f["old_landuse_class"]
    if is_null(new_cls) or is_null(old_cls):
        continue
    if new_cls == old_cls:
        same += a
    else:
        changed += a
        transitions[(old_cls, new_cls)] += a / 1e4
print(f"agreement {same / (same + changed):.1%}")
for (o, n), ha in sorted(transitions.items(), key=lambda kv: -kv[1])[:5]:
    print(f"{o} → {n}: {ha:.1f} ha")

Breakdown: Pieces covered by both versions are compared class by class; pieces covered by only one are skipped, since they are coverage changes rather than reclassifications. The agreement percentage summarises how stable the classification is, and the largest transitions — farmland to residential, forest to clear-cut — are usually the real story. Boundary noise between the two versions inflates the changed area with slivers, so filter slivers before computing agreement, as described below.

Deal with slivers and NULL attributes

Union generates the most slivers of any overlay, because every boundary of both layers cuts every piece. Thin pieces along near-coincident boundaries are noise and inflate the "A only" and "B only" categories.

import math

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

noise = [f.id() for f in measured.getFeatures()
         if f["piece_m2"] < 50 or thinness(f.geometry()) < 0.05]
with edit(measured):
    measured.deleteFeatures(noise)
print(len(noise), "sliver pieces removed")

Breakdown: Dropping slivers removes small areas from the totals; for most cross-tabulations that is the right trade, but report how much area was removed. An alternative that loses nothing is to merge each sliver into its largest neighbour with native:eliminateselectedpolygons, which keeps the total area intact while removing the noise from the table. Snapping the layers to each other before the union is the cleanest fix, as covered in snapping geometries.

QGIS version compatibility

native:union with an optional overlay and field prefix works on QGIS 3.34 LTR, 3.40 LTR and QGIS 4; the prefix parameter has existed since 3.8. native:eliminateselectedpolygons is available on all current releases. On QGIS 4, QgsField takes QMetaType.Type.QString, and NULL attributes arrive as None.

Troubleshooting

  • Field names collide or are renamed oddly. Set OVERLAY_FIELDS_PREFIX.
  • The matrix is missing uncovered areas. NULL attributes were skipped; label them explicitly.
  • Piece counts are huge. Many slivers along near-coincident boundaries; filter, eliminate or snap.
  • Union fails with a topology error. Fix invalid geometries in both inputs first.

Conclusion

Run native:union with an overlay prefix, classify pieces as A only, B only or both from each layer's key, recompute piece areas, cross-tabulate with explicit labels for uncovered areas, use self-union to count overlap depth, and remove or eliminate slivers before reporting.

Frequently Asked Questions

How is union different from merge? Merge stacks features without splitting them, so overlaps remain; union splits at every boundary so no two pieces overlap.

Can I union more than two layers? Chain unions, or merge all layers and self-union the result.

Does union work with lines? It is designed for polygons; for lines, use split with lines or node the network instead.

Why are there pieces with all attributes NULL? Only in self-union with no overlay are all pieces from the one layer; in a two-layer union every piece belongs to at least one input. All-NULL pieces usually come from invalid input geometry.