Geometry Operations and Spatial Predicates in PyQGIS

Every Processing algorithm that clips, buffers, dissolves, or joins is ultimately a loop over QgsGeometry calls. Working one level down — directly on the geometry objects — is what you do when the algorithm you need does not exist, when you want to avoid writing intermediate layers to disk, or when a per-feature decision has to be made in Python rather than expressed as a parameter dictionary.

This guide sits inside Spatial Data Processing & Automation and complements the algorithm-driven pages around it. Where Vector Data Manipulation shows you how to call native:clip on a whole layer, this one explains what the clip is doing to each geometry, why it sometimes fails, and how to perform the same operation yourself when the layer-level tool is the wrong granularity. It covers the geometry object model, the constructive operations, the spatial predicates and their notorious edge cases, correct distance and area measurement, spatial indexing, and geometry validity.

What a QgsGeometry containsA QgsFeature holds one QgsGeometry. The QgsGeometry is a lightweight handle wrapping a QgsAbstractGeometry, which is concretely a QgsPoint, QgsLineString, QgsPolygon or one of their multi-part counterparts. The handle carries the WKB type and validity state while the wrapped geometry carries the coordinates.A handle, a wrapped geometry, and the concrete types beneath itQgsFeatureattributes + one geometryQgsGeometryvalue-type handlewkbType() · isValid()QgsAbstractGeometryholds the coordinatesconstGet()PointQgsPoint · QgsMultiPointLineQgsLineString · QgsMultiLineStringPolygonQgsPolygon · QgsMultiPolygon

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer. QgsGeometry has been API-stable across 3.x; the few method renames are noted inline.
  • A vector layer loaded in the project, or a path you can open with QgsVectorLayer.
  • A projected CRS for anything involving distance or area. If your data is in EPSG:4326, read Coordinate Reference Systems in PyQGIS first — measuring in degrees is the most common source of wrong answers on this page.
  • Comfort with the Python Console; see QGIS Python Console Basics.

What a QgsGeometry actually holds

QgsGeometry is a thin, copy-on-write handle. The coordinates live in a QgsAbstractGeometry subclass underneath, and the handle carries the WKB type plus a cached validity flag. That indirection matters in practice for two reasons.

First, QgsGeometry behaves like a value, not a reference. Assigning one to another copies it lazily, so mutating a geometry you pulled off a feature does not silently edit the layer. Second, the concrete type is reachable but usually unnecessary: constGet() returns the wrapped abstract geometry when you need vertex-level access, but the convenience methods on the handle cover most work.

from qgis.core import QgsProject, QgsWkbTypes

layer = QgsProject.instance().mapLayersByName("parcels")[0]
feature = next(layer.getFeatures())
geometry = feature.geometry()

print(QgsWkbTypes.displayString(geometry.wkbType()))  # e.g. "Polygon"
print(geometry.isMultipart())                          # False for a single polygon
print(geometry.constGet().vertexCount())               # coordinate count

Breakdown: wkbType() returns a numeric enum; QgsWkbTypes.displayString() turns it into readable text, which is far more useful in a log line. isMultipart() is the check that prevents the single most common geometry bug — a script that assumes one part and silently processes only the first. constGet() gives read-only access to the underlying geometry when you need the vertex count or an individual coordinate.

The single-versus-multi distinction deserves respect. A shapefile of "buildings" may contain both Polygon and MultiPolygon features in the same layer, because the format allows it. Rather than branching everywhere, normalise early:

parts = list(geometry.parts()) if geometry.isMultipart() else [geometry.constGet()]
for part in parts:
    print(part.area())

Breakdown: parts() yields each component of a multi-part geometry. Wrapping the single-part case in a one-element list means the loop body is written once. This is the pattern behind nearly every robust geometry utility in the QGIS codebase.

Constructive operations

Constructive operations take one or two geometries and return a new one. They are the building blocks that layer-level algorithms compose. All of them return a fresh QgsGeometry and leave the input untouched.

Constructive operations on one polygonThe same L-shaped source polygon shown four times in outline. Buffer adds a rounded expansion around it, centroid marks a single interior point, convex hull draws the smallest convex shape enclosing it, and simplify draws a version with fewer vertices.One source polygon, four results — the input is never modifiedbuffer(12)grows by a distancecentroid()may fall outside the shapeconvexHull()smallest convex wrappersimplify(8)drops vertices, keeps shape

from qgis.core import QgsProject

layer = QgsProject.instance().mapLayersByName("parcels")[0]
geometry = next(layer.getFeatures()).geometry()

buffered = geometry.buffer(25.0, 8)      # 25 map units, 8 segments per quarter circle
centroid = geometry.centroid()
hull = geometry.convexHull()
simplified = geometry.simplify(5.0)      # Douglas-Peucker tolerance in map units

print(buffered.area(), geometry.area())

Breakdown: buffer() takes a distance in map units — metres in a projected CRS, degrees in a geographic one, which is why the CRS warning above matters. The second argument controls how many straight segments approximate each quarter of a rounded corner; 8 is a reasonable default, and lowering it produces visibly faceted output. simplify() uses the Douglas-Peucker algorithm and its tolerance is also in map units.

A subtlety worth flagging: centroid() can land outside a concave polygon, which makes it a poor choice for label placement. pointOnSurface() guarantees a point inside the geometry and is what labelling engines actually use. The focused recipes Buffer a Geometry in PyQGIS and Simplify Geometry in PyQGIS go into the parameter trade-offs.

Spatial predicates

Predicates answer yes/no questions about how two geometries relate. PyQGIS exposes them as methods that return a plain bool: intersects(), contains(), within(), touches(), crosses(), overlaps(), disjoint(), and equals(). They are cheap, they are the basis of every spatial join, and their edge-case behaviour is where the bugs live.

Which predicate is true in each configurationFive configurations of a blue shape and a teal shape are shown across the top: disjoint, touching at an edge, partially overlapping, one fully contained in the other, and identical. Beneath each, a column of ticks and crosses records whether intersects, touches, within and overlaps return true.Touching counts as intersecting — but not as overlappingdisjointtouchesoverlapswithinequalsintersects()touches()within()overlaps()✗✗✗✗✓✓✗✗✓✗✗✓✓✗✓✗✓✗✓✗within() is true for identical geometries — a shape contains itself

Three results in that matrix routinely surprise people. Touching counts as intersecting: two parcels sharing a boundary satisfy intersects() even though their interiors never meet, which is why a spatial join on intersects double-counts features along shared edges. Overlapping is stricter than intersecting: overlaps() requires the interiors to share area and neither shape to contain the other. And within() is reflexive: a geometry is within itself, so a self-join on within matches every feature to itself unless you filter on feature ID.

from qgis.core import QgsProject

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

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

for parcel in parcels.getFeatures():
    geometry = parcel.geometry()
    if geometry.within(zone_geometry):
        status = "entirely inside"
    elif geometry.intersects(zone_geometry):
        status = "partially affected"
    else:
        continue
    print(parcel["parcel_id"], status)

Breakdown: Testing within() before intersects() matters — the order encodes the specificity, because every geometry that is within another also intersects it. Skipping the else branch entirely avoids materialising results for the majority of features that are irrelevant. Test Whether Two Geometries Intersect in PyQGIS unpacks the boundary cases with runnable examples.

Overlay operations between two geometries

Predicates answer whether two geometries relate; overlay operations produce the geometry that expresses how. The three that matter are intersection(), combine() (union), and difference(), and at geometry level they are the exact operations the layer-level algorithms in Vector Data Manipulation apply feature by feature.

Intersection, union and difference on one pair of shapesThe same two overlapping circles, labelled A on the left and B on the right, are shown three times. Intersection fills only the lens-shaped region where they overlap. Union fills the whole combined outline. Difference fills the part of A that lies outside B.Same inputs, three different results — order matters only for differencea.intersection(b)ABthe shared region onlya.combine(b)one merged geometrya.difference(b)A minus B — not commutative

from qgis.core import QgsGeometry, QgsProject

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

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

for parcel in parcels.getFeatures():
    geometry = parcel.geometry()
    if not geometry.intersects(flood):
        continue
    inside = geometry.intersection(flood)
    outside = geometry.difference(flood)
    share = inside.area() / geometry.area() if geometry.area() else 0
    print(f"{parcel['parcel_id']}: {share:.1%} in the flood zone, "
          f"{outside.area():.0f} m² outside")

Breakdown: Guarding with intersects() before calling intersection() is not just tidiness — for disjoint inputs the overlay still runs the full GEOS computation before returning an empty geometry, so the cheap predicate is a real saving on large layers. intersection() and difference() both return new geometries; neither modifies the inputs. The if geometry.area() else 0 guard prevents a division error on a degenerate zero-area polygon, which turns up more often than you would expect in digitised data.

combine() is the union, and its most common use is dissolving many geometries into one. Doing that in a loop is quadratic — each step rebuilds a progressively larger polygon — so for more than a handful of parts, collect them and use the static form:

geometries = [f.geometry() for f in parcels.getFeatures()]
dissolved = QgsGeometry.unaryUnion(geometries)
print(dissolved.area(), "m² total, ", len(list(dissolved.parts())), "part(s)")

Breakdown: unaryUnion() unions the whole collection in one optimised GEOS pass instead of n pairwise combines. The result is a multi-polygon whenever the inputs are not all contiguous, which is why counting parts() afterwards is a useful sanity check — a "dissolve" that returns forty parts usually means the inputs had gaps you did not know about.

Working at the vertex level

Sometimes the operation you want has no name in GEOS: snapping a stray node, dropping every third vertex, or rewriting Z values. QgsGeometry exposes vertex-level access for those cases, with an index scheme that runs continuously across all parts and rings.

from qgis.core import QgsGeometry, QgsPointXY, QgsProject

layer = QgsProject.instance().mapLayersByName("survey_lines")[0]
feature = next(layer.getFeatures())
geometry = QgsGeometry(feature.geometry())   # explicit copy — do not edit in place

for index, vertex in enumerate(geometry.vertices()):
    print(index, round(vertex.x(), 2), round(vertex.y(), 2))

geometry.moveVertex(QgsPointXY(0, 0).x(), QgsPointXY(0, 0).y(), 0)
geometry.deleteVertex(1)
print(geometry.asWkt(precision=2))

Breakdown: QgsGeometry(other) makes the copy explicit; without it you are mutating an object that may be shared. vertices() yields QgsPoint objects in traversal order, and the index it produces is the same one moveVertex() and deleteVertex() expect. asWkt(precision=2) is the fastest way to eyeball a geometry in the console — full-precision WKT is unreadable. Bear in mind that deleteVertex() can invalidate a geometry (a triangle with a vertex removed is a line), so re-check isGeosValid() before writing the result back.

Measuring distance and area correctly

geometry.area() and geometry.distance(other) are planar measurements in the layer's own map units. In a projected CRS that is what you want. In EPSG:4326 they return square degrees and degrees, which are not units of area or length and vary with latitude.

For measurements that account for the curvature of the Earth, use QgsDistanceArea, which performs ellipsoidal calculations:

from qgis.core import QgsDistanceArea, QgsProject, QgsUnitTypes

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

calculator = QgsDistanceArea()
calculator.setSourceCrs(layer.crs(), QgsProject.instance().transformContext())
calculator.setEllipsoid(QgsProject.instance().ellipsoid())

for feature in layer.getFeatures():
    square_metres = calculator.measureArea(feature.geometry())
    hectares = calculator.convertAreaMeasurement(
        square_metres, QgsUnitTypes.AreaUnit.AreaHectares
    )
    print(feature["parcel_id"], round(hectares, 3))

Breakdown: setSourceCrs() tells the calculator what the incoming coordinates mean, and setEllipsoid() switches it from planar to ellipsoidal mode — without that second call it silently falls back to planar arithmetic. measureArea() returns square metres regardless of the source CRS, and convertAreaMeasurement() handles the unit conversion so you never hard-code 10000. This is exactly the machinery behind the $area expression variable discussed in Working with QGIS Expressions, and it is what a script to calculate polygon areas uses under the hood.

Making predicate queries fast

The naive spatial join — for every feature in A, test every feature in B — is quadratic. On two layers of 10,000 features that is 100 million predicate evaluations, each involving real geometric work. A QgsSpatialIndex reduces it to a bounding-box lookup followed by a handful of exact tests.

Brute force versus a spatial indexOn the left, a nested loop draws lines from one query shape to all sixteen candidate features, labelled as ten thousand by ten thousand exact tests. On the right, a spatial index bounding box selects three nearby candidates, and only those three receive an exact predicate test.The index does not answer the question — it shortens the candidate listnested loop — every pair10 000 × 10 000 exact testsspatial index — candidates first3 exact tests after a box lookup

from qgis.core import QgsProject, QgsSpatialIndex

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

index = QgsSpatialIndex(zones.getFeatures())
zone_geometries = {f.id(): f.geometry() for f in zones.getFeatures()}

for parcel in parcels.getFeatures():
    geometry = parcel.geometry()
    for zone_id in index.intersects(geometry.boundingBox()):
        if geometry.intersects(zone_geometries[zone_id]):
            print(parcel["parcel_id"], "hits zone", zone_id)
            break

Breakdown: QgsSpatialIndex(zones.getFeatures()) builds an R-tree over the bounding boxes in one pass. index.intersects(bbox) returns candidate feature IDs whose boxes overlap — a fast approximation that can include false positives, which is why the exact geometry.intersects() test still runs. Caching the geometries in a dictionary avoids a second provider round-trip per candidate; on very large layers, prefer QgsSpatialIndex(..., flags=QgsSpatialIndex.FlagStoreFeatureGeometries) and index.geometry(id) so the geometries live in the index itself. Build and Use a Spatial Index in PyQGIS covers the memory trade-off and nearest-neighbour queries.

Validity: why operations fail

GEOS, the library behind these operations, requires valid geometry. A self-intersecting polygon, a ring that touches itself, or a polygon whose hole falls outside its shell will make buffer(), intersection(), or difference() raise or return an empty result — often with a terse message that names a coordinate rather than a feature.

from qgis.core import QgsGeometry, QgsProject

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

for feature in layer.getFeatures():
    geometry = feature.geometry()
    errors = geometry.validateGeometry()
    if errors:
        print(feature.id(), errors[0].what())
        repaired = geometry.makeValid()
        print("  repaired:", repaired.isGeosValid())

Breakdown: validateGeometry() returns a list of QgsGeometry.Error objects, each with a message and often a location — far more actionable than a boolean. makeValid() returns a repaired copy using the GEOS validity routines; it can change the geometry type (a self-intersecting polygon may become a multi-polygon), so check the result before writing it back. isGeosValid() is the cheap boolean check when you do not need the details. The full repair workflow, including what to do when makeValid() changes the part count, is in Find and Fix Invalid Geometries in PyQGIS.

Reading and writing geometry as text

Debugging geometry problems is much easier when you can see the coordinates, and interoperating with other tools usually means exchanging a text or binary encoding rather than a QGIS object. QgsGeometry converts in both directions.

from qgis.core import QgsGeometry

wkt = "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))"
geometry = QgsGeometry.fromWkt(wkt)
if geometry.isNull():
    raise ValueError("could not parse that WKT")

print(geometry.asWkt(precision=1))
print(geometry.asJson(precision=6))          # GeoJSON geometry object
print(len(geometry.asWkb()))                 # compact binary form, in bytes

Breakdown: fromWkt() returns a null geometry rather than raising when the text is malformed, so isNull() is the required check — a null geometry passed downstream produces confusing empty results much later. asWkt(precision=...) rounds the output, which is what makes it readable in a console or a test assertion. asJson() emits the GeoJSON geometry fragment, ready to embed in a feature; note it always assumes the coordinates are already in the CRS you intend to publish, so reproject first. asWkb() is the compact binary encoding used for database round-trips.

This is also the most practical way to write a regression test for a geometry routine: build the input from WKT, run the operation, and compare asWkt(precision=3) against an expected string. That avoids depending on a data file and makes failures readable. The same pattern underpins the plugin test suites described in Testing and CI for QGIS Plugins.

One conversion trap is worth naming: WKT carries no CRS. A polygon read from WKT is just numbers, and QGIS will happily treat it as being in whatever CRS you attach it to. When a geometry arrives from an external system, confirm the CRS out of band and reproject explicitly rather than assuming it matches the layer you are about to compare it against.

Key takeaways

  • QgsGeometry is a value-type handle. Copying one is cheap and mutating it never edits the source layer; constGet() reaches the coordinates when you need vertex-level access.
  • Normalise single versus multi-part early. Iterating parts() with a one-element fallback removes an entire class of "only the first polygon was processed" bugs.
  • Buffer and simplify distances are in map units. In a geographic CRS that means degrees, so reproject to a metric CRS before measuring anything.
  • Touching intersects but does not overlap, and every geometry is within itself. Those two facts explain most spatial-join surprises.
  • Use QgsDistanceArea for real-world measurements. Set both the source CRS and the ellipsoid; without the ellipsoid it quietly stays planar.
  • Index before joining. A QgsSpatialIndex turns a quadratic scan into a bounding-box lookup plus a few exact tests, and the exact test is still required.
  • Validate before operating. validateGeometry() reports what is wrong and where; makeValid() repairs it but may change the geometry type.

Frequently Asked Questions

Why does buffer() produce a strange shape or an empty geometry? Almost always invalid input or the wrong CRS. Run geometry.isGeosValid() first and repair with makeValid() if it fails. If the geometry is valid but the buffer looks enormous or microscopic, the layer is probably in a geographic CRS and the distance is being read as degrees.

What is the difference between intersects() and overlaps()?intersects() is true whenever the geometries share any point at all, including a single boundary touch. overlaps() additionally requires that their interiors share area and that neither one contains the other. Two parcels sharing an edge intersect but do not overlap.

Should I use centroid() or pointOnSurface() for labelling?pointOnSurface(), because it is guaranteed to fall inside the geometry. The centroid of a concave or ring-shaped polygon can land outside it entirely, which puts the label in the wrong place or over a neighbour.

How do I get area in hectares instead of square map units? Use QgsDistanceArea with the source CRS and project ellipsoid set, call measureArea() for square metres, then convertAreaMeasurement() with AreaUnit.AreaHectares. Calling geometry.area() directly returns square map units, which are square degrees in EPSG:4326.

Does a spatial index give me the answer directly? No. It returns candidate feature IDs whose bounding boxes intersect your query box, which is an approximation that includes false positives. You still run the exact predicate on each candidate — the index only shortens the list.