PyQGIS Engineering & Automation

Buffering and Spatial Predicates with QgsGeometry

Compute buffers and evaluate spatial predicates in PyQGIS — QgsGeometry.buffer, intersects, contains, within, and disjoint, with correct CRS units and the…

TL;DR: geom.buffer(distance, segments) returns a new QgsGeometry in the layer’s CRS units, so buffer in a projected CRS if you want the distance to mean metres; predicates such as a.intersects(b), a.contains(b), and a.within(b) return a plain bool computed by GEOS against the DE-9IM intersection matrix.

This page is part of the Geometry Operations and Validation in PyQGIS guide, which sits inside the broader PyQGIS Core Architecture & Data Handling reference. Buffering and spatial predicates are the two operations you reach for most often when a task is phrased as “everything within X metres of Y” — and they are the two operations most easily broken by an unprojected CRS or an invalid input geometry.

Complete Runnable Example

The script below buffers every feature in a point layer by a fixed distance in metres, then selects the polygons that intersect any buffer. It is a drop-in template: point it at two loaded layers by name and it will populate the polygon layer’s selection. It also demonstrates the full predicate family and the segments parameter.

python
"""
buffer_and_select.py

Buffer a point layer by N metres and select intersecting polygons.
Demonstrates QgsGeometry.buffer and the spatial-predicate family.
Compatible with QGIS 3.28+ (LTR).
"""
from __future__ import annotations

from typing import Iterator

from qgis.core import (
    Qgis,
    QgsFeature,
    QgsGeometry,
    QgsProject,
    QgsSpatialIndex,
    QgsVectorLayer,
    QgsMessageLog,
    QgsWkbTypes,
)


def require_projected(layer: QgsVectorLayer) -> None:
    """
    Guard against buffering in a geographic CRS.

    Raises:
        ValueError: if the layer's CRS is geographic (degrees), where a
            'metre' buffer distance would be silently misinterpreted.
    """
    if layer.crs().isGeographic():
        raise ValueError(
            f"Layer '{layer.name()}' uses a geographic CRS "
            f"({layer.crs().authid()}). Reproject to a projected CRS so "
            "buffer distances are metres, not degrees."
        )


def buffer_geometry(geom: QgsGeometry, distance: float, segments: int = 8) -> QgsGeometry:
    """
    Return a new geometry offset outward by `distance` CRS units.

    Args:
        geom:     Source geometry (must be valid).
        distance: Offset in the layer's CRS units (metres for a metric CRS).
        segments: Line segments used to approximate each quarter circle;
            higher is smoother and slower. 8 is a common default.

    Returns:
        A new QgsGeometry; the input is left unmodified.
    """
    if not geom.isGeosValid():
        # Buffering an invalid geometry can raise or produce garbage.
        geom = geom.makeValid()
    return geom.buffer(distance, segments)


def select_polygons_near_points(
    points_layer: QgsVectorLayer,
    polygons_layer: QgsVectorLayer,
    distance_m: float,
    segments: int = 8,
) -> list[int]:
    """
    Select polygons that intersect a buffer around any point.

    Uses a spatial index over the polygon layer to shortlist candidates
    by bounding box, then applies the exact `intersects` predicate.

    Returns:
        The list of selected polygon feature IDs.
    """
    require_projected(points_layer)
    require_projected(polygons_layer)

    # Build a bounding-box index over the polygons once, up front.
    index = QgsSpatialIndex(polygons_layer.getFeatures())
    poly_geoms: dict[int, QgsGeometry] = {
        f.id(): QgsGeometry(f.geometry()) for f in polygons_layer.getFeatures()
    }

    selected: set[int] = set()
    for point in points_layer.getFeatures():
        buffer_geom = buffer_geometry(point.geometry(), distance_m, segments)
        # Shortlist by bounding box — cheap integer/float comparisons.
        for fid in index.intersects(buffer_geom.boundingBox()):
            candidate = poly_geoms[fid]
            # Exact test via GEOS / DE-9IM.
            if buffer_geom.intersects(candidate):
                selected.add(fid)

    polygons_layer.selectByIds(list(selected))
    QgsMessageLog.logMessage(
        f"Selected {len(selected)} polygon(s) within {distance_m} m of a point.",
        "buffer_and_select",
        Qgis.Info,
    )
    return list(selected)


def predicate_report(a: QgsGeometry, b: QgsGeometry) -> dict[str, bool]:
    """
    Evaluate the full topological predicate family for two geometries.

    Every method returns a plain bool computed by GEOS against the
    DE-9IM intersection matrix of `a` and `b`.
    """
    return {
        "intersects": a.intersects(b),   # share at least one point
        "disjoint": a.disjoint(b),       # share no points (== not intersects)
        "contains": a.contains(b),       # b lies entirely inside a
        "within": a.within(b),           # a lies entirely inside b
        "touches": a.touches(b),         # boundaries meet, interiors do not
        "overlaps": a.overlaps(b),       # partial same-dimension overlap
        "crosses": a.crosses(b),         # interiors cross, lower-dim result
        "equals": a.equals(b),           # geometrically identical
    }


if __name__ == "__main__":
    project = QgsProject.instance()
    pts = project.mapLayersByName("survey_points")[0]
    polys = project.mapLayersByName("parcels")[0]
    fids = select_polygons_near_points(pts, polys, distance_m=250.0, segments=12)
    print(f"{len(fids)} parcels selected.")

The predicate family returned by predicate_report is the practical toolkit: intersects answers “do these touch at all?”, contains and within are exact inverses of each other, and touches, overlaps, and crosses distinguish how two geometries relate once you know they are not disjoint.

The Predicate Model

Every predicate below is derived from one 3×3 matrix — the DE-9IM (Dimensionally Extended 9-Intersection Model) — that records the dimension of the intersection between the interior, boundary, and exterior of two geometries. A named predicate is just a shorthand for a particular pattern in that matrix. The diagram shows which predicate is true for a common geometric arrangement.

Spatial Predicates for Two Shapes Four small panels each showing two shapes A and B and naming the predicate that is true: A contains B when B sits fully inside A; A intersects B when they partly overlap; A touches B when only their edges meet; A disjoint B when they are fully separate. A B A.contains(B) B.within(A) A B A.intersects(B) A.overlaps(B) A B A.touches(B) edges meet only A B A.disjoint(B) share no points buffer distance = CRS units source geometry (valid input) buffer(dist, seg) new geometry predicate(candidate) bool via DE-9IM

Architecture Breakdown

QgsGeometry.buffer() — the contract and CRS units

The signature is buffer(distance: float, segments: int) -> QgsGeometry. It returns a new geometry — the original is never mutated — representing every point within distance of the source. A positive distance grows the geometry outward; a negative distance on a polygon shrinks it inward (an erosion), which can collapse thin slivers to empty geometry.

The single most important fact about buffer() is that the distance is expressed in the geometry’s CRS units, not metres. A QgsGeometry carries no CRS of its own — it inherits the CRS of the layer or feature it came from. If that CRS is geographic (EPSG:4326, units of degrees), then buffer(250, 8) builds a 250-degree buffer, which is nonsensical and wraps the globe. The fix is to work in a projected CRS whose map units are metres. Reproject the layer, or transform the individual geometry, before buffering. The mechanics of building and applying those transforms belong to the Coordinate Transformations and CRS Handling in PyQGIS guide; the rule here is simply: confirm layer.crs().isGeographic() is False before you trust a metric distance.

The segments parameter controls how many straight line segments approximate each quarter of a circular arc at convex corners. A value of 8 (the common default) gives 32 segments around a full circle — smooth enough for display and most analysis. Raising it improves roundness at the cost of vertex count and downstream predicate speed; lowering it to 1 or 2 produces a chamfered, blocky buffer that is cheaper to test but visibly faceted.

The predicate family and DE-9IM

The predicate methods — intersects, disjoint, contains, within, touches, overlaps, crosses, and equals — each return a plain Python bool. Internally GEOS computes the DE-9IM intersection matrix once and each named predicate checks whether that matrix matches a specific pattern. Two identities are worth memorising: a.disjoint(b) is always the logical negation of a.intersects(b), and a.contains(b) is exactly b.within(a). Choosing the more natural direction keeps calling code readable.

The distinctions between the non-disjoint predicates are precise. touches is true only when the geometries share boundary points but their interiors do not meet — two parcels sharing a fence line. overlaps requires the geometries to be of the same dimension and to share interior area while neither contains the other. crosses applies when interiors intersect in a result of lower dimension than at least one input — a road line crossing a district polygon. When you only need “are these related at all?”, intersects is both the cheapest to reason about and the safest default, because it is true for every arrangement except full separation.

For non-standard relationships GEOS also exposes relate(other, pattern), which tests the DE-9IM matrix against an explicit nine-character pattern string. Reach for it only when the named predicates cannot express the relationship you need — they cover the overwhelming majority of real tasks.

Validity is a precondition

GEOS predicates and the buffer operation both assume topologically valid input. A self-intersecting polygon (a “bowtie”), a ring with duplicate consecutive vertices, or an unclosed exterior ring can make a predicate return a misleading result or make buffer() raise. The example guards this with isGeosValid() and falls back to makeValid() before buffering. In a pipeline processing external data, validate before the first spatial operation rather than debugging a wrong answer later — the failure modes and repair strategies are covered in Fixing Invalid Geometries in PyQGIS. Note that isGeosValid() is the GEOS-backed check used by the predicate engine, which is the relevant one here; a geometry that displays fine in the canvas can still be invalid to GEOS.

Registration and Integration: pairing predicates with a spatial index

A spatial predicate is an exact, per-geometry-pair computation, and it is comparatively expensive. Running buffer.intersects(candidate) across every feature of a large polygon layer is O(n) full geometry comparisons per buffer, which does not scale. The standard remedy is a two-phase filter: use a QgsSpatialIndex to shortlist candidates by bounding box first — an O(log n) R-tree lookup — and only run the exact predicate on the handful of features whose envelopes overlap.

python
from qgis.core import QgsGeometry, QgsSpatialIndex, QgsVectorLayer


def features_intersecting_buffer(
    polygons: QgsVectorLayer,
    buffer_geom: QgsGeometry,
    index: QgsSpatialIndex,
    geom_by_id: dict[int, QgsGeometry],
) -> list[int]:
    """
    Two-phase filter: bounding-box shortlist, then exact predicate.

    The spatial index removes the vast majority of candidates cheaply;
    the exact intersects() test runs only on the survivors.
    """
    hits: list[int] = []
    for fid in index.intersects(buffer_geom.boundingBox()):
        if buffer_geom.intersects(geom_by_id[fid]):
            hits.append(fid)
    return hits

The index returns a superset — every feature whose bounding box overlaps the buffer’s bounding box — so the exact predicate is still required to discard box-overlap-but-geometry-miss cases. The trade-off, index construction cost versus per-query saving, and the choice between an in-memory QgsSpatialIndex and provider-side filtering are the subject of Spatial Indexing and Query Optimization in PyQGIS. For a one-off comparison of two geometries you can skip the index entirely; it earns its keep from roughly the first few hundred candidate features onward.

QgsGeometry.buffer vs the native:buffer algorithm

QgsGeometry.buffer() operates on a single geometry object in memory. The Processing algorithm native:buffer (run via processing.run("native:buffer", {...})) operates on a whole layer, writing a new buffered layer and exposing extra parameters — END_CAP_STYLE, JOIN_STYLE, MITER_LIMIT, and a DISSOLVE flag that merges overlapping buffers into one feature. Use the geometry method when you are already iterating features and need fine-grained control inside a loop; use the algorithm when you want a declarative, one-call layer-to-layer transform that slots into a Processing model or a headless batch. Both ultimately call the same GEOS buffer routine, so results are consistent — the difference is the unit of work and the surrounding parameter surface.

Production Best Practices

  • Buffer in a projected CRS. Confirm layer.crs().isGeographic() is False before treating a buffer distance as metres. For global data, reproject to an appropriate UTM zone or equal-area CRS first.
  • Validate before you buffer or test. Call isGeosValid() and repair with makeValid() up front; GEOS predicates assume valid input and fail quietly on bowties and self-touching rings.
  • Prefer intersects for coarse “are these related” tests, and reserve touches, overlaps, and crosses for cases where the manner of contact matters.
  • Remember contains/within are inverses and disjoint is not intersects — pick the direction that reads naturally and avoid computing both.
  • Pair predicates with a QgsSpatialIndex whenever you test one geometry against many; shortlist by boundingBox() then confirm with the exact predicate.
  • Tune segments to the task, not by habit: high values for smooth cartographic buffers, low values when you need cheap, blocky proximity zones.
  • Reach for native:buffer when you need dissolve, cap, or join styling across a whole layer; use the geometry method for per-feature control inside a loop.