Intersect Two Vector Layers in PyQGIS

Intersection is the overlay that answers "how much of A lies in B": how many hectares of each soil type fall inside each field, which length of each road runs through a flood zone, which part of each parcel is covered by a protection designation. The operation cuts features of one layer by the features of another and keeps only the overlapping pieces, each carrying attributes from both sides. It is one of the most used geoprocessing steps, and one whose output is easy to misread when areas are not recomputed or slivers are not filtered.

This recipe belongs to Vector Data Manipulation. It runs native:intersection, controls which attributes are carried across, recomputes measurements on the pieces, filters slivers, deals with invalid input, and summarises the result into the table the question asked for.

What an intersection returnsTwo input layers, fields and soil zones, overlap. The intersection returns only the overlapping pieces: each piece is a part of one field that lies in one soil zone, carrying the field's attributes and the soil zone's attributes. Parts of fields outside every soil zone are dropped, as are parts of soil zones outside every field.Only the overlap, with attributes from bothfieldssoil zones→pieceattributesfield_idcropsoil_type

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series.
  • Two layers in the same projected CRS. Processing reprojects the overlay on the fly when CRSs differ, but measuring areas in a geographic CRS gives meaningless numbers, so reproject first if needed.
  • Valid geometries in both layers; a validity report takes seconds.

Run the intersection

The algorithm takes an input layer, an overlay layer and optional lists of fields to keep from each.

import processing
from qgis.core import QgsProject

fields_layer = QgsProject.instance().mapLayersByName("fields")[0]
soils = QgsProject.instance().mapLayersByName("soil_zones")[0]

pieces = processing.run("native:intersection", {
    "INPUT": fields_layer,
    "OVERLAY": soils,
    "INPUT_FIELDS": ["field_id", "crop"],
    "OVERLAY_FIELDS": ["soil_type", "drainage"],
    "OVERLAY_FIELDS_PREFIX": "",
    "OUTPUT": "memory:fields_x_soils",
})["OUTPUT"]

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

Breakdown: Each output feature is the geometric intersection of one input feature with one overlay feature, so a field lying across three soil zones becomes three pieces. Listing fields explicitly keeps the output narrow and readable; empty lists keep all fields from that side. OVERLAY_FIELDS_PREFIX adds a prefix such as soil_ to overlay field names, which avoids clashes when both layers have a field called name or id. The output geometry type follows the input: intersecting polygons with polygons gives polygons, lines with polygons gives lines.

Recompute measurements on the pieces

Attributes copied from the inputs describe the whole original features. An area_ha field from the fields layer still holds the whole field's area on every piece — the most common mistake in overlay analysis.

Copied areas are wrong after cuttingA field of 12 hectares is cut into three pieces of 5, 4 and 3 hectares. If the input's area field is copied, every piece claims 12 hectares and a sum gives 36. Recomputing area on each piece gives 5, 4 and 3, which sum correctly to 12.Recompute after every overlaycopied fieldpiece 1: 12 hapiece 2: 12 hapiece 3: 12 hasum: 36 ha ✗recomputed $areapiece 1: 5 hapiece 2: 4 hapiece 3: 3 hasum: 12 ha ✓

measured = processing.run("native:fieldcalculator", {
    "INPUT": pieces, "FIELD_NAME": "piece_ha", "FIELD_TYPE": 0,
    "FIELD_LENGTH": 12, "FIELD_PRECISION": 4,
    "FORMULA": "$area / 10000",
    "OUTPUT": "memory:fields_x_soils_measured"})["OUTPUT"]

total_pieces = sum(f["piece_ha"] for f in measured.getFeatures())
total_fields = sum(f.geometry().area() for f in fields_layer.getFeatures()) / 1e4
print(f"pieces {total_pieces:,.1f} ha of {total_fields:,.1f} ha of fields")

Breakdown: $area measures each piece using the project's ellipsoid settings; area($geometry) gives planar area in layer units, which is what you want when the layer is in an equal-area or local projected CRS. Comparing the total of the pieces with the total of the inputs is a quick sanity check: the pieces can only add up to less than or equal to the fields' area, the difference being the parts outside every soil zone. For line layers, recompute $length the same way.

Filter slivers

Where two layers were digitized independently, their boundaries rarely coincide exactly. Intersection then produces slivers — thin pieces along shared boundaries that are artefacts, not real overlaps.

import math

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

slivers = [f.id() for f in measured.getFeatures()
           if f["piece_ha"] < 0.01 or thinness(f.geometry()) < 0.05]
print(len(slivers), "probable slivers out of", measured.featureCount())

from qgis.core import edit
with edit(measured):
    measured.deleteFeatures(slivers)

Breakdown: Very small pieces (here under 100 m²) and very thin ones (low isoperimetric quotient) are almost always boundary mismatch. Choose thresholds from the data's accuracy rather than by eye; reporting how many pieces and how much area were removed keeps the cleanup honest. The better long-term fix is to make the layers agree first by snapping one to the other, after which slivers mostly disappear.

Handle invalid input

Overlays are where invalid geometries fail loudly — or worse, quietly lose features. Processing's invalid-feature setting decides what happens; setting it explicitly in a script makes the behaviour predictable.

from qgis.core import QgsProcessingContext, QgsFeatureRequest, QgsProcessingFeedback

context = QgsProcessingContext()
context.setInvalidGeometryCheck(QgsFeatureRequest.GeometryAbortOnInvalid)
feedback = QgsProcessingFeedback()
try:
    processing.run("native:intersection", {
        "INPUT": fields_layer, "OVERLAY": soils, "OUTPUT": "memory:"},
        context=context, feedback=feedback)
except Exception as err:
    print("invalid geometry stopped the run:", err)

Breakdown: Aborting on invalid geometry is the safe setting for production scripts: a run either succeeds on clean data or stops with a clear message. The alternatives — skipping invalid features or ignoring the check — produce a result that silently lacks some features or contains wrong pieces. Fix the input with native:fixgeometries and run again, as described in fixing invalid geometries.

Intersect with several overlays at once

Questions often involve more than two layers: the parts of each field that are on clay soil, inside a water-protection zone and below 200 m elevation. Chaining pairwise intersections works, but native:multiintersection overlays several layers in one call and keeps attributes from all of them.

protection = QgsProject.instance().mapLayersByName("water_protection")[0]
lowland = QgsProject.instance().mapLayersByName("below_200m")[0]

combined = processing.run("native:multiintersection", {
    "INPUT": fields_layer,
    "OVERLAYS": [soils, protection, lowland],
    "OVERLAY_FIELDS_PREFIX": "",
    "OUTPUT": "memory:fields_all_overlays"})["OUTPUT"]
print(combined.featureCount(), "pieces inside all three overlays")

Breakdown: Each output piece lies inside the input feature and at least one feature of every overlay, so only the area satisfying all conditions remains — an "and" across layers. Attributes from every overlay are appended, which is why a prefix or a careful field selection beforehand matters when overlays share field names. One call is easier to read than three chained ones and avoids writing intermediate layers, though the measurement and sliver steps above still apply to the result. For "or" logic — inside any of several zones — dissolve the overlays into one layer first and intersect with that.

Summarise the result

The intersection is rarely the final answer. The question — hectares of each soil type per crop, or per field — is a group-by over the pieces.

From pieces to a summary tableThe pieces layer has one row per field and soil combination with recomputed area. Grouping by crop and soil type and summing area gives a compact table of hectares, for example wheat on clay 120 hectares, wheat on loam 340 hectares. The table answers the question; the pieces layer remains for mapping.Group the pieces to answer the questionpiecesfield × soilpiece_hagroup bycrop, soil_typesum piece_hasummarywheat · clay 120 hawheat · loam 340 ha

from collections import defaultdict

table = defaultdict(float)
for f in measured.getFeatures():
    table[(f["crop"], f["soil_type"])] += f["piece_ha"]

for (crop, soil), ha in sorted(table.items()):
    print(f"{crop:<12} {soil:<10} {ha:10.1f} ha")

Breakdown: Summing recomputed piece areas by the attribute combination of interest produces the table directly. For larger tables or further analysis, the pieces can go to pandas as in analysing an attribute table with pandas. Keep the pieces layer: it lets readers check any number in the table by looking at the map.

Intersect a layer with itself

Overlaps within one layer — conflicting designations, double-digitized parcels — are found by intersecting a layer with itself and discarding each feature's intersection with itself.

self_overlaps = processing.run("native:intersection", {
    "INPUT": soils, "OVERLAY": soils,
    "INPUT_FIELDS": ["zone_id"], "OVERLAY_FIELDS": ["zone_id"],
    "OVERLAY_FIELDS_PREFIX": "other_", "OUTPUT": "memory:"})["OUTPUT"]
real = [f for f in self_overlaps.getFeatures() if f["zone_id"] != f["other_zone_id"]]
print(len(real) // 2, "overlapping pairs of soil zones")

Breakdown: Every feature intersects itself completely, so those pieces are dropped by comparing ids; each genuine overlap appears twice, once from each side, hence the halving. This is a quick check rather than a full topology audit — finding gaps and overlaps handles tolerances and reporting properly.

QGIS version compatibility

native:intersection with field selection and prefix works on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. Since 3.20 it can use multiple overlay layers through native:multiintersection. On QGIS 4, QgsFeatureRequest.GeometryAbortOnInvalid is Qgis.InvalidGeometryCheck.AbortOnInvalid.

Troubleshooting

  • Totals are far too large. Area fields were copied from the inputs; recompute on the pieces.
  • The output is empty. The layers do not overlap, often because one is in a different CRS with a wrong definition.
  • Thousands of tiny pieces. Boundaries do not coincide; filter slivers or snap first.
  • The run stops with a geometry error. Fix invalid geometries before overlaying.

Conclusion

Intersect with explicit field lists and a prefix to avoid name clashes, recompute area or length on every piece, filter slivers with area and thinness thresholds tied to accuracy, fail fast on invalid geometry, and group the pieces into the summary the question asks for.

Frequently Asked Questions

What is the difference between intersection and clip? Clip keeps only the input's attributes and cuts it by the overlay's outline; intersection keeps attributes from both and splits by each overlay feature.

Can I intersect points with polygons? Yes — the output is the points inside polygons, with polygon attributes. A spatial join is usually simpler.

Is the result different if I swap input and overlay? Geometry is the same; field order and the output geometry type follow the input.

Why does the output have more features than either input? Each input feature is split once per overlay feature it overlaps, so a field crossing four soil zones becomes four pieces. Feature counts grow with how finely the two layers cut each other.

How do I keep the parts outside the overlay too? Use union, covered in union overlay of layers.