Test Whether Two Geometries Intersect in PyQGIS

intersects() returns a boolean and looks like the simplest method in the geometry API. It is also the one most often used where a different predicate was meant, because its definition is broader than the everyday word suggests: two shapes that merely share a boundary point intersect, even though their interiors never meet. Spatial joins that double-count along shared edges almost always trace back to that single fact.

This page is a focused recipe within Geometry Operations and Spatial Predicates in PyQGIS. It covers the basic test, the boundary cases that make it surprising, how to choose between intersects, contains, within and overlaps, and how to make the test fast over a whole layer.

Prerequisites

• QGIS 3.34 LTR (bundled Python 3.12) or newer.
• Two layers, or two geometries, in the same CRS. Predicates compare raw coordinates and do not reproject.
• Valid geometry — GEOS predicates can raise on self-intersecting input.

The basic test

Both geometries must be QgsGeometry objects, and the method is symmetric: a.intersects(b) and b.intersects(a) always agree.

from qgis.core import QgsProject

parcels = QgsProject.instance().mapLayersByName("parcels")[0]
zones = QgsProject.instance().mapLayersByName("flood_zones")[0]

zone_geometry = next(zones.getFeatures()).geometry()

hits = [
    feature["parcel_id"]
    for feature in parcels.getFeatures()
    if feature.geometry().intersects(zone_geometry)
]
print(len(hits), "parcel(s) touch the flood zone")

Breakdown: feature.geometry() returns the geometry handle; calling intersects() on it does the full GEOS test. A list comprehension is fine here because the predicate is cheap, but note that every feature is read from the provider — the prefiltering section below is what makes this scale.

If the two layers are in different CRSs the answer will be silently wrong rather than an error, because the coordinate ranges simply never overlap. Guard for it:

if parcels.crs() != zones.crs():
    raise ValueError(
        f"CRS mismatch: {parcels.crs().authid()} vs {zones.crs().authid()}"
    )

Breakdown: Comparing QgsCoordinateReferenceSystem objects with != is reliable. Raising rather than reprojecting silently is the right default in a script someone else will run — a mismatch usually means the wrong layer was picked, not that a transform is needed.

Choose the right predicate

intersects() answers "do these share any point at all". Most real questions are narrower, and picking the specific predicate makes the intent readable and the result correct.

for feature in parcels.getFeatures():
    geometry = feature.geometry()
    if geometry.within(zone_geometry):
        category = "wholly inside"
    elif geometry.intersects(zone_geometry):
        category = "partly inside"
    else:
        category = "outside"
    print(feature["parcel_id"], category)

Breakdown: Testing within() first and intersects() second encodes the specificity — every geometry that is within another also intersects it, so the reverse order would classify everything as "partly inside". This three-way split is the shape most real classification code takes.

Prefilter with bounding boxes

An exact predicate is a real geometric computation. A bounding-box comparison is four number comparisons. Running the cheap test first discards the overwhelming majority of non-matches for almost nothing.

from qgis.core import QgsFeatureRequest

zone_box = zone_geometry.boundingBox()
request = QgsFeatureRequest().setFilterRect(zone_box)

hits = [
    f["parcel_id"]
    for f in parcels.getFeatures(request)
    if f.geometry().intersects(zone_geometry)
]

Breakdown: setFilterRect() pushes the box test down to the provider, so features outside it are never read at all — a much bigger saving than filtering in Python. The exact intersects() still runs on what survives, because a bounding box is a coarse approximation: a long diagonal line has a large box that overlaps shapes the line itself misses entirely. For repeated queries against many geometries, promote this to a QgsSpatialIndex — see Build and Use a Spatial Index in PyQGIS.

Test many pairs efficiently

When both sides are layers, prepare the query geometry once. GEOS can build an internal index for a geometry that will be tested repeatedly, and QgsGeometry exposes it through a prepared form used automatically when you call the predicate in a tight loop against the same object.

from qgis.core import QgsFeatureRequest, QgsGeometry

zone_geometries = [f.geometry() for f in zones.getFeatures()]
merged = QgsGeometry.unaryUnion(zone_geometries)
box = merged.boundingBox()

request = QgsFeatureRequest().setFilterRect(box)
affected = [
    f["parcel_id"]
    for f in parcels.getFeatures(request)
    if f.geometry().intersects(merged)
]
print(len(affected), "parcel(s) affected by any zone")

Breakdown: Unioning the zones once turns n × m comparisons into n × 1. That is almost always the right move when the question is "does this touch any of them" rather than "which one does it touch". unaryUnion() merges the whole collection in a single optimised pass rather than pairwise. If you do need to know which zone, keep them separate and use a spatial index instead.

Point in polygon: the special case

Point-in-polygon is the most common predicate question of all, and it has one genuine ambiguity: a point lying exactly on the boundary. within() excludes it, intersects() includes it, and which you want depends on whether the polygons tile a surface.

from qgis.core import QgsFeatureRequest, QgsProject, QgsSpatialIndex

points = QgsProject.instance().mapLayersByName("incidents")[0]
districts = QgsProject.instance().mapLayersByName("districts")[0]

index = QgsSpatialIndex(
    districts.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries
)
names = {f.id(): f["district_name"] for f in districts.getFeatures()}

for point in points.getFeatures():
    geometry = point.geometry()
    assigned = None
    for did in index.intersects(geometry.boundingBox()):
        if geometry.within(index.geometry(did)):
            assigned = names[did]
            break
    print(point["incident_id"], assigned or "unassigned")

Breakdown: Using within() rather than intersects() means a point sitting exactly on a shared district boundary is assigned to neither district instead of both — which is what you want when the districts tile a surface and each incident must be counted once. The break stops at the first containing polygon, which is safe precisely because within() cannot match two neighbours. If unassigned boundary points are unacceptable, fall back to intersects() for those and pick deterministically, for example the district with the lowest feature ID, so repeated runs agree.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

The predicate methods on QgsGeometry have been stable throughout 3.x, and their semantics follow the OGC simple-features definitions, so they match PostGIS behaviour as well.

Troubleshooting

• Everything returns false. The two layers are in different CRSs, so their coordinate ranges do not overlap. Compare authid() on both.
• Adjacent parcels are counted twice in a join. They share a boundary, which intersects() treats as a hit. Use overlaps() if a shared edge should not count.
• A geometry appears to intersect itself. within() and intersects() are both true for identical geometries. Exclude self-matches by comparing feature IDs in a self-join.
• A GEOS exception is raised. One of the inputs is invalid. Repair it with makeValid() first.
• The test is slow on large layers. No bounding-box prefilter. Add QgsFeatureRequest().setFilterRect() or build a spatial index.
• Results differ from the GUI's spatial join. The GUI defaults to a specific predicate set; check which predicates were ticked rather than assuming it used intersects.

Conclusion

intersects() is the broadest of the spatial predicates: it is true for a shared corner, a shared edge, a partial overlap and full containment alike. Choosing the narrower predicate that matches your actual question makes the code self-documenting and avoids the double-counting that shared boundaries otherwise cause. For performance, prefilter with a bounding box and union the query side when you only need "any".

Frequently Asked Questions

Do two polygons that share only a border intersect? Yes. intersects() is true whenever the geometries share any point, including a single boundary point. If a shared edge should not count, use overlaps(), which requires the interiors to share area.

Does intersects() reproject the geometries? No. It compares raw coordinates and assumes both are in the same CRS. A mismatch produces False for every pair with no warning, so check crs().authid() on both layers before testing.

Is a bounding-box test enough on its own? No. It is an approximation that produces false positives — a diagonal line's box covers a large area the line never touches. Use it to shorten the candidate list, then run the exact predicate.

What is the difference between contains and within? They are inverses. a.contains(b) and b.within(a) mean the same thing. Both are true when the geometries are identical, because a shape contains itself.

Is disjoint just the opposite of intersects? Yes, exactly — a.disjoint(b) is always the negation of a.intersects(b). Use whichever reads more naturally in the surrounding code.

Related

• Geometry Operations and Spatial Predicates in PyQGIS — the guide this recipe belongs to
• Build and Use a Spatial Index in PyQGIS
• Calculate the Distance Between Features in PyQGIS
• Clip a Vector Layer in PyQGIS