Qgis Impl Vector Analysis

qgis-impl-vector-analysis

Quick Reference

Core Classes

Key Processing Algorithm IDs

---

Critical Warnings

NEVER pass layers with different CRS to overlay operations without reprojecting first. Processing algorithms reproject automatically, but PyQGIS geometry methods do NOT. ALWAYS verify CRS match before direct QgsGeometry operations.

NEVER call QgsGeometry methods on NULL geometries. ALWAYS check feature.hasGeometry() before ANY geometry operation. NULL geometries cause silent failures or crashes.

NEVER forget to call del writer after using QgsVectorFileWriter.create(). The file is NOT written to disk until the writer object is destroyed.

NEVER use 'OUTPUT': 'memory:' for large datasets in batch processing. Memory layers consume RAM and are lost when QGIS closes. ALWAYS write to GeoPackage for persistence.

NEVER assume processing.run() modifies the input layer. It ALWAYS creates a new output layer. Access results via result['OUTPUT'].

ALWAYS use native:fixgeometries before overlay operations when input data comes from external sources. Invalid geometries cause overlay algorithms to fail silently or produce incomplete results.

ALWAYS initialize the Processing framework before calling processing.run() in standalone scripts:

from qgis.core import QgsApplication
import processing
from processing.core.Processing import Processing
Processing.initialize()

---

Decision Tree: Which Overlay Operation to Use

What is your goal?
|
+-- Keep ONLY the area where both layers overlap
|   --> native:intersection
|
+-- Combine ALL areas from both layers into one
|   --> native:union
|
+-- Remove areas of layer B from layer A
|   +-- Keep remainder of A only
|   |   --> native:difference
|   +-- Keep remainder of BOTH A and B (exclude overlap)
|       --> native:symmetricaldifference
|
+-- Cut layer A to the boundary of layer B
|   --> native:clip
|
+-- Transfer attributes from one layer to another
    +-- Based on matching field values
    |   --> native:joinattributestable
    +-- Based on spatial relationship
    |   --> native:joinattributesbylocation
    +-- Based on nearest feature
        --> native:joinattributesbynearest

Clip vs Intersection

---

Essential Patterns

Pattern 1: Spatial Query with Feature Request

from qgis.core import QgsFeatureRequest, QgsRectangle

# Filter by bounding rectangle with exact geometry test
area = QgsRectangle(100.0, -1.0, 101.0, 0.0)
request = QgsFeatureRequest().setFilterRect(area)
request.setFlags(QgsFeatureRequest.ExactIntersect)
request.setLimit(100)

for feature in layer.getFeatures(request):
    print(feature.id(), feature.geometry().asWkt())

Pattern 2: Spatial Index for Fast Lookups

from qgis.core import QgsSpatialIndex, QgsPointXY

index = QgsSpatialIndex(layer.getFeatures())

# Nearest neighbor: returns feature IDs
nearest_ids = index.nearestNeighbor(QgsPointXY(15.5, 47.1), 5)

# Bounding box intersection: returns feature IDs
bbox = QgsRectangle(14.0, 46.0, 17.0, 49.0)
intersecting_ids = index.intersects(bbox)

# Retrieve actual features by ID
for fid in nearest_ids:
    feature = layer.getFeature(fid)

Pattern 3: Buffer Analysis via Processing

import processing

result = processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'SEGMENTS': 5,
    'END_CAP_STYLE': 0,   # 0=Round, 1=Flat, 2=Square
    'JOIN_STYLE': 0,       # 0=Round, 1=Miter, 2=Bevel
    'MITER_LIMIT': 2,
    'DISSOLVE': False,
    'OUTPUT': 'memory:'
})
buffered_layer = result['OUTPUT']

Pattern 4: Buffer via QgsGeometry (Single Feature)

geom = feature.geometry()
if not geom.isNull():
    buffered = geom.buffer(100.0, 5)  # distance, segments

Pattern 5: Overlay Operation (Intersection)

import processing

result = processing.run("native:intersection", {
    'INPUT': layer_a,
    'OVERLAY': layer_b,
    'INPUT_FIELDS': [],
    'OVERLAY_FIELDS': [],
    'OVERLAY_FIELDS_PREFIX': '',
    'OUTPUT': 'memory:'
})
intersection_layer = result['OUTPUT']

Pattern 6: Attribute Join by Field Value

result = processing.run("native:joinattributestable", {
    'INPUT': layer_a,
    'FIELD': 'id',
    'INPUT_2': layer_b,
    'FIELD_2': 'foreign_id',
    'FIELDS_TO_COPY': [],
    'METHOD': 0,               # 0=one-to-many, 1=first match only
    'DISCARD_NONMATCHING': False,
    'PREFIX': '',
    'OUTPUT': 'memory:'
})

Pattern 7: Spatial Join by Location

result = processing.run("native:joinattributesbylocation", {
    'INPUT': target_layer,
    'PREDICATE': [0],          # 0=intersects, 1=contains, 2=equals,
                               # 3=touches, 4=overlaps, 5=within, 6=crosses
    'JOIN': join_layer,
    'JOIN_FIELDS': [],
    'METHOD': 0,               # 0=one-to-many, 1=one-to-first, 2=largest overlap
    'DISCARD_NONMATCHING': False,
    'PREFIX': '',
    'OUTPUT': 'memory:'
})

Pattern 8: Dissolve by Attribute

result = processing.run("native:dissolve", {
    'INPUT': layer,
    'FIELD': ['province'],     # Dissolve field(s); empty list = dissolve all
    'SEPARATE_DISJOINT': False,
    'OUTPUT': 'memory:'
})

---

Common Operations

Field Calculation (Edit Buffer)

with edit(layer):
    field_idx = layer.fields().indexOf('area_m2')
    for feature in layer.getFeatures():
        if feature.hasGeometry():
            area = feature.geometry().area()
            layer.changeAttributeValue(feature.id(), field_idx, area)

Field Calculation (Processing)

result = processing.run("native:fieldcalculator", {
    'INPUT': layer,
    'FIELD_NAME': 'area_m2',
    'FIELD_TYPE': 0,           # 0=Float, 1=Integer, 2=String, 3=Date
    'FIELD_LENGTH': 10,
    'FIELD_PRECISION': 3,
    'FORMULA': '$area',
    'OUTPUT': 'memory:'
})

Feature Selection

# Select by expression
layer.selectByExpression('"type" = \'highway\'')

# Iterate selected features
for feature in layer.selectedFeatures():
    print(feature.id())

# Clear selection
layer.removeSelection()

Extract Features by Expression

result = processing.run("native:extractbyexpression", {
    'INPUT': layer,
    'EXPRESSION': '"population" > 50000',
    'OUTPUT': 'memory:'
})

Extract Features by Location

result = processing.run("native:extractbylocation", {
    'INPUT': layer,
    'PREDICATE': [0],          # 0=intersects
    'INTERSECT': reference_layer,
    'OUTPUT': 'memory:'
})

Write Output to GeoPackage

from qgis.core import QgsVectorFileWriter, QgsCoordinateTransformContext

save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.driverName = "GPKG"
save_options.fileEncoding = "UTF-8"

error = QgsVectorFileWriter.writeAsVectorFormatV3(
    layer,
    "/path/to/output.gpkg",
    QgsCoordinateTransformContext(),
    save_options
)

Write Features Individually

from qgis.core import QgsVectorFileWriter, QgsWkbTypes, QgsCoordinateTransformContext

save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.driverName = "GPKG"

writer = QgsVectorFileWriter.create(
    "/path/to/output.gpkg",
    layer.fields(),
    QgsWkbTypes.Polygon,
    layer.crs(),
    QgsCoordinateTransformContext(),
    save_options
)
for feature in layer.getFeatures():
    writer.addFeature(feature)
del writer  # ALWAYS delete to flush and close the file

Aggregate with Statistics

result = processing.run("native:aggregate", {
    'INPUT': layer,
    'GROUP_BY': '"province"',
    'AGGREGATES': [
        {'aggregate': 'sum', 'delimiter': ',', 'input': '"population"',
         'length': 10, 'name': 'total_pop', 'precision': 0, 'type': 6}
    ],
    'OUTPUT': 'memory:'
})

---

Reference Links

• references/methods.md -- API signatures for QgsGeometry, QgsVectorFileWriter, QgsFeatureRequest, QgsSpatialIndex, and processing algorithm parameters
• references/examples.md -- Complete vector analysis workflows with realistic data
• references/anti-patterns.md -- Common vector analysis mistakes and how to avoid them

Official Sources

• https://docs.qgis.org/3.34/en/docs/pyqgis_developer_cookbook/vector.html
• https://docs.qgis.org/3.34/en/docs/pyqgis_developer_cookbook/geometry.html
• https://qgis.org/pyqgis/3.34/core/QgsGeometry.html
• https://qgis.org/pyqgis/3.34/core/QgsVectorFileWriter.html
• https://qgis.org/pyqgis/3.34/core/QgsSpatialIndex.html

qgis-impl-vector-analysis — Anti-Patterns

AP-01: Operating on NULL Geometries

Wrong

for feature in layer.getFeatures():
    area = feature.geometry().area()  # CRASHES if geometry is NULL
    buffered = feature.geometry().buffer(100, 5)  # Returns empty geometry silently

Right

for feature in layer.getFeatures():
    if not feature.hasGeometry():
        continue
    geom = feature.geometry()
    if geom.isNull() or geom.isEmpty():
        continue
    area = geom.area()
    buffered = geom.buffer(100, 5)

Why

Vector layers from external sources frequently contain features with NULL or empty geometries. Calling geometry methods on NULL geometries causes crashes or produces silently wrong results. ALWAYS check hasGeometry() before ANY geometry operation.

---

AP-02: Mixing CRS in Direct Geometry Operations

Wrong

# layer_a is EPSG:4326 (geographic), layer_b is EPSG:28992 (projected)
for feat_a in layer_a.getFeatures():
    for feat_b in layer_b.getFeatures():
        if feat_a.geometry().intersects(feat_b.geometry()):  # WRONG CRS comparison
            print("Intersects!")

Right

from qgis.core import QgsCoordinateTransform, QgsCoordinateReferenceSystem, QgsProject

transform = QgsCoordinateTransform(
    layer_a.crs(),
    layer_b.crs(),
    QgsProject.instance()
)

for feat_a in layer_a.getFeatures():
    geom_a = QgsGeometry(feat_a.geometry())
    geom_a.transform(transform)  # Transform to layer_b's CRS
    for feat_b in layer_b.getFeatures():
        if geom_a.intersects(feat_b.geometry()):
            print("Intersects!")

Why

QgsGeometry methods perform pure coordinate comparison with NO automatic CRS transformation. Comparing geometries in different CRS produces completely wrong results. Processing algorithms handle CRS automatically, but direct QgsGeometry operations do NOT. ALWAYS reproject to a common CRS before direct geometry operations.

---

AP-03: Forgetting to Delete QgsVectorFileWriter

Wrong

writer = QgsVectorFileWriter.create(
    "/path/to/output.gpkg", fields, QgsWkbTypes.Polygon,
    crs, QgsCoordinateTransformContext(), save_options
)
for feature in layer.getFeatures():
    writer.addFeature(feature)
# File may be incomplete — writer is still open

Right

writer = QgsVectorFileWriter.create(
    "/path/to/output.gpkg", fields, QgsWkbTypes.Polygon,
    crs, QgsCoordinateTransformContext(), save_options
)
for feature in layer.getFeatures():
    writer.addFeature(feature)
del writer  # ALWAYS delete to flush buffers and close file

Why

QgsVectorFileWriter buffers data in memory. The file is NOT fully written until the writer object is destroyed. Without del writer, the output file may be incomplete, truncated, or locked by the process.

---

AP-04: Using Memory Layers for Persistent Results

Wrong

result = processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': 'memory:'  # Lost when QGIS closes
})

Right (for persistent results)

result = processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': '/path/to/output/buffers.gpkg'
})

When Memory IS Appropriate

Use 'memory:' ONLY for intermediate results in a processing chain where the final output is written to disk. NEVER use 'memory:' as the final output of an analysis that must persist.

---

AP-05: Not Fixing Geometries Before Overlay Operations

Wrong

# External data may contain invalid geometries
result = processing.run("native:intersection", {
    'INPUT': external_layer,  # May have self-intersections, duplicate vertices
    'OVERLAY': boundary,
    'OUTPUT': 'memory:'
})
# Algorithm may fail or produce incomplete results

Right

# ALWAYS fix geometries from external sources first
fixed = processing.run("native:fixgeometries", {
    'INPUT': external_layer,
    'OUTPUT': 'memory:'
})['OUTPUT']

result = processing.run("native:intersection", {
    'INPUT': fixed,
    'OVERLAY': boundary,
    'OUTPUT': 'memory:'
})

Why

Overlay operations (intersection, union, difference) rely on valid topology. Invalid geometries (self-intersections, duplicate vertices, unclosed rings) cause algorithms to fail silently, produce incomplete output, or raise errors. ALWAYS run native:fixgeometries on external data before overlay operations.

---

AP-06: Assuming processing.run() Modifies the Input Layer

Wrong

processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': 'memory:'
})
# Expecting 'layer' to now contain buffers — it does NOT
for feature in layer.getFeatures():
    print(feature.geometry().area())  # Still the original geometries

Right

result = processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': 'memory:'
})
buffered_layer = result['OUTPUT']  # This is the new layer with buffers
for feature in buffered_layer.getFeatures():
    print(feature.geometry().area())

Why

processing.run() NEVER modifies the input layer. It ALWAYS creates a new output. The result dictionary contains the output layer under the 'OUTPUT' key. Ignoring the return value means your analysis results are discarded.

---

AP-07: Using Nested Feature Loops Without Spatial Index

Wrong

# O(n*m) complexity — extremely slow for large datasets
for feat_a in layer_a.getFeatures():
    for feat_b in layer_b.getFeatures():
        if feat_a.geometry().intersects(feat_b.geometry()):
            process(feat_a, feat_b)

Right

from qgis.core import QgsSpatialIndex

# Build index — O(m log m)
index_b = QgsSpatialIndex(layer_b.getFeatures())

# Query index — O(n log m) total
for feat_a in layer_a.getFeatures():
    if not feat_a.hasGeometry():
        continue
    bbox = feat_a.geometry().boundingBox()
    candidate_ids = index_b.intersects(bbox)
    for fid in candidate_ids:
        feat_b = layer_b.getFeature(fid)
        if feat_a.geometry().intersects(feat_b.geometry()):
            process(feat_a, feat_b)

Why

Nested feature loops without a spatial index have O(n*m) complexity. For two layers of 10,000 features each, this means 100 million geometry comparisons. A spatial index reduces this to approximately O(n log m), making the operation orders of magnitude faster. ALWAYS use QgsSpatialIndex for spatial queries across layers.

---

AP-08: Using Shapefile for New Projects

Wrong

save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.driverName = "ESRI Shapefile"  # Legacy format with many limitations

Right

save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.driverName = "GPKG"  # Modern, no limitations

Why

Shapefile has severe limitations: field names truncated to 10 characters, single geometry type per file, no NULL value support, 2GB file size limit, requires multiple sidecar files (.shx, .dbf, .prj). GeoPackage has NONE of these limitations. ALWAYS use GeoPackage for new projects. Only use Shapefile when required by external systems.

---

AP-09: Not Initializing Processing Framework in Standalone Scripts

Wrong

# In a standalone Python script (NOT the QGIS Python console)
import processing
result = processing.run("native:buffer", {...})  # ERROR: Processing not initialized

Right

from qgis.core import QgsApplication
import sys

qgs = QgsApplication([], False)
qgs.initQgis()

# Initialize Processing AFTER QgsApplication
from processing.core.Processing import Processing
Processing.initialize()

import processing
result = processing.run("native:buffer", {...})

qgs.exitQgis()

Why

The Processing framework is NOT automatically available in standalone scripts. It requires explicit initialization after QgsApplication.initQgis(). In the QGIS Python console, Processing is already initialized. This distinction catches many developers who test code in the console and then deploy it as a standalone script.

---

AP-10: Buffer Distance in Wrong Units

Wrong

# Layer is in EPSG:4326 (degrees)
result = processing.run("native:buffer", {
    'INPUT': layer_4326,
    'DISTANCE': 100,       # This is 100 DEGREES, not 100 meters!
    'OUTPUT': 'memory:'
})

Right

# Option 1: Reproject to a projected CRS first
reprojected = processing.run("native:reprojectlayer", {
    'INPUT': layer_4326,
    'TARGET_CRS': 'EPSG:32632',  # UTM zone 32N (meters)
    'OUTPUT': 'memory:'
})['OUTPUT']

result = processing.run("native:buffer", {
    'INPUT': reprojected,
    'DISTANCE': 100,       # Now correctly 100 meters
    'OUTPUT': 'memory:'
})

# Option 2: Use QgsDistanceArea for unit conversion
from qgis.core import QgsDistanceArea, QgsCoordinateReferenceSystem
d = QgsDistanceArea()
d.setSourceCrs(layer_4326.crs(), QgsProject.instance().transformContext())
d.setEllipsoid('WGS84')
# Then use d.measureLine() for accurate distances

Why

Buffer distances are specified in the CRS units of the input layer. For geographic CRS (EPSG:4326), units are degrees. A buffer of 100 in EPSG:4326 means 100 degrees — roughly 11,000 km. ALWAYS reproject to a projected CRS (meters) before buffer operations, or use QgsDistanceArea for ellipsoidal calculations.

qgis-impl-vector-analysis — Examples

Example 1: Complete Buffer and Clip Workflow

Buffer a point layer by 500m and clip to a study area polygon.

import processing
from qgis.core import QgsProject

# Load layers
points = QgsProject.instance().mapLayersByName('sample_points')[0]
study_area = QgsProject.instance().mapLayersByName('study_boundary')[0]

# Step 1: Fix geometries (ALWAYS do this for external data)
fixed = processing.run("native:fixgeometries", {
    'INPUT': points,
    'OUTPUT': 'memory:'
})['OUTPUT']

# Step 2: Buffer points by 500 meters
buffered = processing.run("native:buffer", {
    'INPUT': fixed,
    'DISTANCE': 500,
    'SEGMENTS': 5,
    'END_CAP_STYLE': 0,
    'JOIN_STYLE': 0,
    'MITER_LIMIT': 2,
    'DISSOLVE': False,
    'OUTPUT': 'memory:'
})['OUTPUT']

# Step 3: Clip buffers to study area
clipped = processing.run("native:clip", {
    'INPUT': buffered,
    'OVERLAY': study_area,
    'OUTPUT': '/path/to/output/buffered_clipped.gpkg'
})['OUTPUT']

# Step 4: Add to map
QgsProject.instance().addMapLayer(clipped)

---

Example 2: Spatial Join — Transfer Zone Data to Points

Transfer land use zone attributes to point features based on spatial location.

import processing

points = QgsProject.instance().mapLayersByName('buildings')[0]
zones = QgsProject.instance().mapLayersByName('land_use_zones')[0]

result = processing.run("native:joinattributesbylocation", {
    'INPUT': points,
    'PREDICATE': [5],              # 5 = within
    'JOIN': zones,
    'JOIN_FIELDS': ['zone_type', 'max_height', 'fsi'],
    'METHOD': 1,                   # First match only
    'DISCARD_NONMATCHING': False,
    'PREFIX': 'zone_',
    'OUTPUT': 'memory:'
})

joined_layer = result['OUTPUT']
joined_layer.setName('buildings_with_zones')
QgsProject.instance().addMapLayer(joined_layer)

# Verify join results
total = joined_layer.featureCount()
matched = 0
for f in joined_layer.getFeatures():
    if f['zone_zone_type'] is not None:
        matched += 1
print(f"Matched {matched} of {total} features")

---

Example 3: Dissolve and Aggregate Statistics

Dissolve parcels by municipality and calculate total area and count.

import processing

parcels = QgsProject.instance().mapLayersByName('parcels')[0]

# Simple dissolve by field
dissolved = processing.run("native:dissolve", {
    'INPUT': parcels,
    'FIELD': ['municipality'],
    'SEPARATE_DISJOINT': False,
    'OUTPUT': 'memory:'
})['OUTPUT']

# Aggregate with statistics
aggregated = processing.run("native:aggregate", {
    'INPUT': parcels,
    'GROUP_BY': '"municipality"',
    'AGGREGATES': [
        {'aggregate': 'first_value', 'delimiter': ',', 'input': '"municipality"',
         'length': 100, 'name': 'municipality', 'precision': 0, 'type': 10},
        {'aggregate': 'sum', 'delimiter': ',', 'input': '"area_ha"',
         'length': 10, 'name': 'total_area_ha', 'precision': 2, 'type': 6},
        {'aggregate': 'count', 'delimiter': ',', 'input': '"id"',
         'length': 10, 'name': 'parcel_count', 'precision': 0, 'type': 2}
    ],
    'OUTPUT': 'memory:'
})['OUTPUT']

aggregated.setName('municipality_summary')
QgsProject.instance().addMapLayer(aggregated)

---

Example 4: Multi-Layer Intersection Analysis

Find areas where flood zones overlap with residential land use.

import processing
from qgis.core import QgsProject, QgsVectorFileWriter, QgsCoordinateTransformContext

flood_zones = QgsProject.instance().mapLayersByName('flood_zones')[0]
land_use = QgsProject.instance().mapLayersByName('land_use')[0]

# Step 1: Extract residential areas
residential = processing.run("native:extractbyexpression", {
    'INPUT': land_use,
    'EXPRESSION': '"category" = \'residential\'',
    'OUTPUT': 'memory:'
})['OUTPUT']

# Step 2: Intersect flood zones with residential areas
at_risk = processing.run("native:intersection", {
    'INPUT': flood_zones,
    'OVERLAY': residential,
    'INPUT_FIELDS': ['flood_level', 'return_period'],
    'OVERLAY_FIELDS': ['neighborhood', 'dwelling_count'],
    'OVERLAY_FIELDS_PREFIX': 'res_',
    'OUTPUT': 'memory:'
})['OUTPUT']

# Step 3: Calculate affected area
with_area = processing.run("native:fieldcalculator", {
    'INPUT': at_risk,
    'FIELD_NAME': 'affected_area_m2',
    'FIELD_TYPE': 0,
    'FIELD_LENGTH': 15,
    'FIELD_PRECISION': 2,
    'FORMULA': '$area',
    'OUTPUT': 'memory:'
})['OUTPUT']

# Step 4: Write to GeoPackage
save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.driverName = "GPKG"
save_options.fileEncoding = "UTF-8"

QgsVectorFileWriter.writeAsVectorFormatV3(
    with_area,
    '/path/to/flood_risk_residential.gpkg',
    QgsCoordinateTransformContext(),
    save_options
)

---

Example 5: Attribute Join from CSV

Join census data from a CSV table to a polygon layer.

import processing
from qgis.core import QgsProject, QgsVectorLayer

# Load CSV as table (no geometry)
csv_uri = "file:///path/to/census_data.csv?delimiter=,&encoding=UTF-8"
census_table = QgsVectorLayer(csv_uri, "census_data", "delimitedtext")

neighborhoods = QgsProject.instance().mapLayersByName('neighborhoods')[0]

# Join by matching field
result = processing.run("native:joinattributestable", {
    'INPUT': neighborhoods,
    'FIELD': 'neighborhood_code',
    'INPUT_2': census_table,
    'FIELD_2': 'code',
    'FIELDS_TO_COPY': ['population', 'median_income', 'avg_household_size'],
    'METHOD': 1,                   # First match only
    'DISCARD_NONMATCHING': False,
    'PREFIX': 'census_',
    'OUTPUT': 'memory:'
})

joined = result['OUTPUT']
joined.setName('neighborhoods_with_census')
QgsProject.instance().addMapLayer(joined)

---

Example 6: Spatial Query with Index for Performance

Find all buildings within 200m of a river using spatial index for fast lookups.

from qgis.core import (
    QgsSpatialIndex, QgsFeatureRequest, QgsGeometry, QgsProject
)

buildings = QgsProject.instance().mapLayersByName('buildings')[0]
rivers = QgsProject.instance().mapLayersByName('rivers')[0]

# Build spatial index on buildings
building_index = QgsSpatialIndex(buildings.getFeatures())

# Buffer each river segment and find nearby buildings
nearby_ids = set()
for river_feat in rivers.getFeatures():
    if not river_feat.hasGeometry():
        continue
    river_geom = river_feat.geometry()
    buffer_geom = river_geom.buffer(200.0, 5)
    bbox = buffer_geom.boundingBox()

    # Fast bounding box check via index
    candidate_ids = building_index.intersects(bbox)

    # Exact geometry check
    for fid in candidate_ids:
        building = buildings.getFeature(fid)
        if building.hasGeometry() and buffer_geom.intersects(building.geometry()):
            nearby_ids.add(fid)

# Select the results
buildings.selectByIds(list(nearby_ids))
print(f"Found {len(nearby_ids)} buildings within 200m of rivers")

---

Example 7: Field Calculation via Edit Buffer

Calculate population density on an existing field using the edit buffer.

from qgis.core import QgsProject

layer = QgsProject.instance().mapLayersByName('districts')[0]
density_idx = layer.fields().indexOf('pop_density')
pop_idx = layer.fields().indexOf('population')

with edit(layer):
    for feature in layer.getFeatures():
        if not feature.hasGeometry():
            continue
        area_km2 = feature.geometry().area() / 1_000_000  # m2 to km2
        population = feature[pop_idx]
        if area_km2 > 0 and population is not None:
            density = population / area_km2
            layer.changeAttributeValue(feature.id(), density_idx, density)

---

Example 8: Chain Multiple Overlay Operations

Identify parks within 1km of schools that are NOT in flood zones.

import processing

schools = QgsProject.instance().mapLayersByName('schools')[0]
parks = QgsProject.instance().mapLayersByName('parks')[0]
flood_zones = QgsProject.instance().mapLayersByName('flood_zones')[0]

# Step 1: Buffer schools by 1km
school_buffers = processing.run("native:buffer", {
    'INPUT': schools,
    'DISTANCE': 1000,
    'SEGMENTS': 5,
    'END_CAP_STYLE': 0,
    'JOIN_STYLE': 0,
    'DISSOLVE': True,          # Merge overlapping buffers
    'OUTPUT': 'memory:'
})['OUTPUT']

# Step 2: Extract parks within school buffer zones
parks_near_schools = processing.run("native:extractbylocation", {
    'INPUT': parks,
    'PREDICATE': [0],          # intersects
    'INTERSECT': school_buffers,
    'OUTPUT': 'memory:'
})['OUTPUT']

# Step 3: Remove parks in flood zones (difference by location)
safe_parks = processing.run("native:extractbylocation", {
    'INPUT': parks_near_schools,
    'PREDICATE': [2],          # disjoint (no spatial relationship)
    'INTERSECT': flood_zones,
    'OUTPUT': '/path/to/safe_parks_near_schools.gpkg'
})['OUTPUT']

# Note: PREDICATE [2] = disjoint is NOT available in extractbylocation.
# Instead, use a two-step approach:
parks_in_flood = processing.run("native:extractbylocation", {
    'INPUT': parks_near_schools,
    'PREDICATE': [0],          # intersects
    'INTERSECT': flood_zones,
    'OUTPUT': 'memory:'
})['OUTPUT']

# Get IDs of parks in flood zones
flood_park_ids = {f.id() for f in parks_in_flood.getFeatures()}

# Select parks NOT in flood zones
safe_ids = []
for f in parks_near_schools.getFeatures():
    if f.id() not in flood_park_ids:
        safe_ids.append(f.id())

parks_near_schools.selectByIds(safe_ids)

---

Example 9: Write Features to Multiple Formats

Export analysis results to GeoPackage, GeoJSON, and Shapefile.

from qgis.core import (
    QgsVectorFileWriter, QgsCoordinateTransformContext, QgsProject
)

layer = QgsProject.instance().mapLayersByName('analysis_result')[0]
ctx = QgsCoordinateTransformContext()

# GeoPackage (ALWAYS preferred)
opts_gpkg = QgsVectorFileWriter.SaveVectorOptions()
opts_gpkg.driverName = "GPKG"
opts_gpkg.fileEncoding = "UTF-8"
QgsVectorFileWriter.writeAsVectorFormatV3(layer, '/output/result.gpkg', ctx, opts_gpkg)

# GeoJSON (for web applications)
opts_json = QgsVectorFileWriter.SaveVectorOptions()
opts_json.driverName = "GeoJSON"
opts_json.fileEncoding = "UTF-8"
QgsVectorFileWriter.writeAsVectorFormatV3(layer, '/output/result.geojson', ctx, opts_json)

# Shapefile (legacy compatibility only)
opts_shp = QgsVectorFileWriter.SaveVectorOptions()
opts_shp.driverName = "ESRI Shapefile"
opts_shp.fileEncoding = "UTF-8"
QgsVectorFileWriter.writeAsVectorFormatV3(layer, '/output/result.shp', ctx, opts_shp)

qgis-impl-vector-analysis — Methods Reference

QgsGeometry Operations

Spatial Predicates

Geometry Manipulation

Geometry Measurements

Geometry Conversion

Static Constructors

---

QgsFeatureRequest

Flags

---

QgsSpatialIndex

---

QgsVectorFileWriter

Static Methods

SaveVectorOptions

Common Driver Names

---

Processing Algorithm Parameters

native:buffer

native:intersection

native:union

native:difference

native:clip

native:symmetricaldifference

native:dissolve

native:joinattributestable

native:joinattributesbylocation

native:fieldcalculator

native:extractbyexpression

native:extractbylocation