Qgis Agents Analysis Orchestrator

qgis-agents-analysis-orchestrator

Quick Reference

This skill is an orchestrator — it guides analysis method selection and workflow ordering. It does NOT contain implementation details. For implementation, refer to the specific skill indicated by the decision trees.

When to Use This Skill

---

Critical Warnings

NEVER start coding before completing the workflow plan. ALWAYS determine: (1) analysis type, (2) CRS, (3) input format, (4) algorithm chain, (5) output format.

NEVER mix vector and raster operations without explicit conversion steps. ALWAYS include rasterize/vectorize as a named step.

NEVER chain algorithms with mismatched CRS. ALWAYS reproject to a common CRS as the FIRST step.

NEVER use geographic CRS (EPSG:4326) for distance or area calculations. ALWAYS reproject to a projected CRS first.

NEVER assume third-party provider algorithms (GRASS, SAGA) are available. ALWAYS check availability with QgsApplication.processingRegistry().algorithmById() and fall back to native: algorithms.

NEVER hardcode file paths. ALWAYS use QgsProcessing.TEMPORARY_OUTPUT or 'memory:' for intermediate results.

ALWAYS wrap processing.run() in try/except for QgsProcessingException.

ALWAYS check layer.isValid() after loading any layer.

ALWAYS create spatial indexes (native:createspatialindex) before spatial queries on large datasets.

---

Decision Tree 1: Analysis Type Selection

What is your primary question?
│
├─ "How are features distributed spatially?"
│  └─ VECTOR ANALYSIS → Decision Tree 2
│
├─ "What is the value at this location?" / "How does a surface vary?"
│  └─ RASTER ANALYSIS → Decision Tree 3
│
├─ "What is the shortest/fastest route?" / "What areas are reachable?"
│  └─ NETWORK ANALYSIS → Decision Tree 4
│
├─ "How do two datasets relate spatially?"
│  ├─ Both datasets are vector → VECTOR OVERLAY → Decision Tree 2
│  ├─ Both datasets are raster → RASTER CELL STATISTICS → Decision Tree 3
│  └─ One vector + one raster → HYBRID ANALYSIS:
│     ├─ Extract raster values at vector locations → `native:rastersampling`
│     ├─ Summarize raster within vector zones → `native:zonalstatisticsfb`
│     └─ Convert vector to raster first → `gdal:rasterize` then Decision Tree 3
│
├─ "Create a continuous surface from point samples?"
│  └─ INTERPOLATION:
│     ├─ Sparse, irregular points → `native:tininterpolation`
│     ├─ Dense, regular points → `native:idwinterpolation`
│     └─ Event density visualization → `native:heatmapkerneldensityestimation`
│
└─ "Classify or cluster features?"
   └─ CLUSTERING:
      ├─ Known number of groups → `native:kmeansclustering`
      ├─ Unknown groups, density-based → `native:dbscanclustering`
      └─ Spatiotemporal data → `native:stdbscanclustering`

---

Decision Tree 2: Vector Analysis Method

What vector operation do you need?
│
├─ PROXIMITY ANALYSIS
│  ├─ Buffer around features → `native:buffer`
│  ├─ Distance between all pairs → `native:distancematrix`
│  ├─ Nearest feature connection → `native:distancetonearesthublines`
│  ├─ Shortest line between features → `native:shortestline`
│  └─ Select within distance → `native:extractwithindistance`
│
├─ OVERLAY ANALYSIS
│  ├─ Keep area in BOTH layers → `native:intersection`
│  ├─ Keep area in EITHER layer → `native:union`
│  ├─ Keep area in A but NOT B → `native:difference`
│  ├─ Cut A to shape of B → `native:clip`
│  └─ Keep area in A OR B but NOT both → `native:symmetricaldifference`
│
├─ SPATIAL JOIN
│  ├─ Join by shared location → `native:joinattributesbylocation`
│  ├─ Join by location with stats → `native:joinattributesbylocationsummary`
│  ├─ Join by nearest feature → `native:joinattributesbynearest`
│  └─ Join by matching field value → `native:joinattributesbyfieldvalue`
│
├─ AGGREGATION
│  ├─ Merge geometries by attribute → `native:dissolve`
│  ├─ Aggregate with expressions → `native:aggregate`
│  ├─ Count points in polygons → `native:countpointsinpolygon`
│  └─ Sum line lengths in polygons → `native:sumlinelengths`
│
├─ GEOMETRY TRANSFORMATION
│  ├─ Simplify (reduce vertices) → `native:simplifygeometries`
│  ├─ Convert polygon to line → `native:polygonstolines`
│  ├─ Convert line to polygon → `native:linestopolygons`
│  ├─ Multi-part to single → `native:multiparttosingleparts`
│  ├─ Fix invalid geometries → `native:fixgeometries`
│  ├─ Calculate centroids → `native:centroids`
│  └─ Convex hull → `native:convexhull`
│
└─ SELECTION / EXTRACTION
   ├─ By attribute value → `native:extractbyattribute`
   ├─ By expression → `native:extractbyexpression`
   ├─ By spatial relationship → `native:extractbylocation`
   └─ Random sample → `native:randomextract`

---

Decision Tree 3: Raster Analysis Method

What raster operation do you need?
│
├─ TERRAIN ANALYSIS (from DEM)
│  ├─ Slope angle → `native:slope`
│  ├─ Aspect direction → `native:aspect`
│  ├─ Hillshade visualization → `native:hillshade`
│  ├─ Terrain ruggedness → `native:ruggednessindex`
│  ├─ Fill sinks (hydrology) → `native:fillsinks`
│  └─ Relief rendering → `native:relief`
│
├─ RASTER CALCULATION
│  ├─ Band math / map algebra → `native:rastercalculator`
│  ├─ Reclassify values → `native:reclassifybytable`
│  ├─ Rescale values → `native:rescaleraster`
│  └─ Fill NoData cells → `native:fillnodata`
│
├─ RASTER STATISTICS
│  ├─ Global statistics → `native:rasterlayerstatistics`
│  ├─ Zonal stats per polygon → `native:zonalstatisticsfb`
│  ├─ Zonal histogram → `native:zonalhistogram`
│  ├─ Sample at point locations → `native:rastersampling`
│  └─ Cell statistics across stack → `native:cellstatistics`
│
└─ RASTER COMPARISON
   ├─ Boolean AND across layers → `native:rasterbooleanand`
   ├─ Boolean OR across layers → `native:rasterbooleanor`
   └─ Frequency analysis → `native:equaltofrequency`

---

Decision Tree 4: Network Analysis Method

What network question do you have?
│
├─ "Shortest path between two points"
│  └─ `native:shortestpathpointtopoint`
│
├─ "Shortest path from point to multiple destinations"
│  └─ `native:shortestpathpointtolayer`
│
├─ "Shortest path from multiple origins to one point"
│  └─ `native:shortestpathlayertopoint`
│
├─ "What area is reachable within X minutes/meters?"
│  ├─ From a single point → `native:serviceareafrompoint`
│  └─ From multiple points → `native:serviceareafromlayer`
│
└─ PREREQUISITES (ALWAYS complete these first):
   ├─ Network layer MUST be a line layer
   ├─ Configure direction field for one-way streets
   ├─ Configure speed/cost field for weighted analysis
   └─ NEVER assume bidirectional — check direction attributes

---

CRS Selection Guide

CRS Selection Decision

Where is your study area?
│
├─ Single country/region
│  └─ Use the national projected CRS (e.g., NL=EPSG:28992, UK=EPSG:27700)
│
├─ Crosses UTM zones (narrow east-west)
│  └─ Use the UTM zone covering the majority of the area
│
├─ Continental or global
│  ├─ Area calculations → Equal-area (e.g., EPSG:6933 World Equal Area)
│  ├─ Distance calculations → Equidistant projection
│  └─ Display only → EPSG:4326 or EPSG:3857
│
└─ ALWAYS verify with:
   crs = QgsCoordinateReferenceSystem("EPSG:28992")
   assert crs.isValid(), "CRS is invalid"

---

Output Format Selection Guide

Format Decision

What is the output purpose?
│
├─ Intermediate / temporary result
│  └─ Use 'memory:' or QgsProcessing.TEMPORARY_OUTPUT
│
├─ Final file output
│  ├─ Single dataset → GeoPackage (ALWAYS default choice)
│  ├─ Multiple related layers → GeoPackage with layername parameter
│  ├─ Web API / JavaScript → GeoJSON
│  └─ Legacy system requirement → Shapefile (truncate field names to 10 chars)
│
├─ Multi-user / enterprise
│  └─ PostGIS (see qgis-impl-postgis skill)
│
└─ NEVER use Shapefile as default — ALWAYS prefer GeoPackage

---

Workflow Patterns

Pattern 1: Standard Analysis Chain

1. LOAD      → Load input layers, check isValid()
2. VALIDATE  → Fix geometries (native:fixgeometries)
3. REPROJECT → Reproject all to common CRS (native:reprojectlayer)
4. INDEX     → Create spatial index (native:createspatialindex)
5. ANALYZE   → Run analysis algorithm(s)
6. EXPORT    → Save to target format

Pattern 2: Multi-Layer Overlay

1. LOAD      → Load all input layers
2. VALIDATE  → Fix geometries on ALL layers
3. REPROJECT → Reproject ALL to common projected CRS
4. INDEX     → Create spatial index on ALL layers
5. OVERLAY   → Run overlay operation (intersection/union/difference)
6. CLEAN     → Remove slivers, fix topology
7. EXPORT    → Save result

Pattern 3: Raster-Vector Hybrid

1. LOAD      → Load raster + vector layers
2. REPROJECT → Reproject vector to match raster CRS
3. EXTRACT   → Sample raster at vector locations (native:rastersampling)
              OR Zonal statistics (native:zonalstatisticsfb)
4. ANALYZE   → Further vector analysis on enriched data
5. EXPORT    → Save result

Pattern 4: Iterative Processing (Batch)

# ALWAYS use processing.runAndLoadResults() for final output only
# ALWAYS use 'memory:' for intermediate results
layers = QgsProject.instance().mapLayersByName("input")
for layer in layers:
    if not layer.isValid():
        continue
    try:
        buffered = processing.run("native:buffer", {
            'INPUT': layer,
            'DISTANCE': 100,
            'OUTPUT': 'memory:'
        })['OUTPUT']
        # Chain next operation...
    except QgsProcessingException as e:
        QgsMessageLog.logMessage(f"Failed: {e}", "Analysis")

---

Performance Optimization Checklist

---

MCP Server Integration

Two existing MCP servers enable Claude to control a running QGIS instance:

jjsantos01/qgis_mcp (855+ stars)

• Purpose: Direct QGIS control from Claude Desktop
• Capabilities: Load layers, run processing algorithms, manage projects, create layouts
• Setup: Install via QGIS Plugin Manager + configure Claude Desktop MCP settings
• URL: https://github.com/jjsantos01/qgis_mcp

nkarasiak/qgis-mcp (51+ tools)

• Purpose: Comprehensive QGIS MCP with 51 tools
• Capabilities: Layer management, styling, processing, spatial queries, project management
• URL: https://github.com/nkarasiak/qgis-mcp

When to Use MCP vs PyQGIS Scripts

Is QGIS currently running and accessible?
│
├─ YES → Use MCP server for interactive control
│  ├─ Layer manipulation → MCP tools
│  ├─ Visual feedback needed → MCP (sees map canvas)
│  └─ Complex multi-step analysis → MCP + this orchestrator skill
│
└─ NO → Generate standalone PyQGIS scripts
   ├─ Headless processing → Use qgis.core initialization
   ├─ Batch processing → Standalone script with QgsApplication
   └─ Plugin development → Follow plugin skill patterns

---

Code Validation Checklist

ALWAYS verify generated PyQGIS code against this checklist before presenting it:

• [ ] Layer validity: Every QgsVectorLayer / QgsRasterLayer creation is followed by isValid() check
• [ ] CRS consistency: All layers in the same operation share a CRS, or explicit reprojection is included
• [ ] CRS for measurements: Distance/area calculations use a projected CRS, NEVER geographic
• [ ] Transform context: Every QgsCoordinateTransform includes QgsProject.instance().transformContext()
• [ ] Edit sessions: Feature modifications use with edit(layer): context manager
• [ ] Error handling: processing.run() is wrapped in try/except
• [ ] Algorithm existence: Third-party algorithms are checked with algorithmById() before use
• [ ] No GUI in threads: processAlgorithm() uses feedback object, NEVER QMessageBox or iface
• [ ] Temp outputs: Intermediate results use 'memory:' or QgsProcessing.TEMPORARY_OUTPUT
• [ ] No hardcoded paths: Paths use os.path.join() or pathlib.Path, NEVER backslashes
• [ ] Spatial index: Large datasets have native:createspatialindex before spatial operations
• [ ] Feature request flags: NoGeometry flag used when geometry is not needed
• [ ] NULL geometry check: geom.isNull() checked before geometry operations
• [ ] Expression validation: hasParserError() checked after QgsExpression creation
• [ ] Memory layer updates: updateExtents() called after adding features to memory layers
• [ ] Field updates: updateFields() called after modifying field structure

---

Reference Links

• references/methods.md — Algorithm selection methods and workflow construction patterns
• references/examples.md — Complete workflow examples for common analysis scenarios
• references/anti-patterns.md — Workflow-level anti-patterns and planning mistakes

Official Sources

• https://docs.qgis.org/latest/en/docs/pyqgis_developer_cookbook/
• https://docs.qgis.org/latest/en/docs/user_manual/processing/
• https://qgis.org/pyqgis/master/
• https://github.com/jjsantos01/qgis_mcp
• https://github.com/nkarasiak/qgis-mcp

anti-patterns.md — Analysis Orchestrator

Workflow Planning Anti-Patterns

AP-1: Skipping CRS Alignment

WRONG: Running overlay operations on layers with different CRS.

# WRONG — layers have different CRS, result is geometrically incorrect
result = processing.run("native:intersection", {
    'INPUT': layer_epsg4326,
    'OVERLAY': layer_epsg28992,
    'OUTPUT': 'memory:'
})

CORRECT: ALWAYS reproject all layers to a common CRS before ANY multi-layer operation.

# CORRECT — reproject first
reprojected = processing.run("native:reprojectlayer", {
    'INPUT': layer_epsg4326,
    'TARGET_CRS': 'EPSG:28992',
    'OUTPUT': 'memory:'
})['OUTPUT']

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

---

AP-2: Using Geographic CRS for Distance/Area

WRONG: Calculating buffer distance in degrees (EPSG:4326).

# WRONG — distance 0.01 is in degrees, not meters
result = processing.run("native:buffer", {
    'INPUT': points_epsg4326,
    'DISTANCE': 0.01,  # What does 0.01 degrees mean in meters? It varies!
    'OUTPUT': 'memory:'
})

CORRECT: ALWAYS reproject to a projected CRS for distance/area operations.

# CORRECT — reproject, then buffer in meters
reproj = processing.run("native:reprojectlayer", {
    'INPUT': points_epsg4326,
    'TARGET_CRS': 'EPSG:28992',
    'OUTPUT': 'memory:'
})['OUTPUT']

result = processing.run("native:buffer", {
    'INPUT': reproj,
    'DISTANCE': 1000,  # 1000 meters — unambiguous
    'OUTPUT': 'memory:'
})

---

AP-3: No Error Handling in Algorithm Chains

WRONG: Chaining algorithms without error handling — one failure crashes the entire chain with an unhelpful traceback.

# WRONG — no error handling
buffered = processing.run("native:buffer", params1)['OUTPUT']
clipped = processing.run("native:clip", {'INPUT': buffered, ...})['OUTPUT']
dissolved = processing.run("native:dissolve", {'INPUT': clipped, ...})['OUTPUT']

CORRECT: ALWAYS wrap each step in try/except with meaningful error messages.

# CORRECT — each step has error handling
try:
    buffered = processing.run("native:buffer", params1)['OUTPUT']
except QgsProcessingException as e:
    raise RuntimeError(f"Step 1 (buffer) failed: {e}")

try:
    clipped = processing.run("native:clip", {
        'INPUT': buffered, 'OVERLAY': mask, 'OUTPUT': 'memory:'
    })['OUTPUT']
except QgsProcessingException as e:
    raise RuntimeError(f"Step 2 (clip) failed: {e}")

---

AP-4: Missing Layer Validity Checks

WRONG: Using a layer without checking if it loaded correctly.

# WRONG — layer may be invalid (wrong path, missing provider)
layer = QgsVectorLayer("/data/missing_file.gpkg", "data", "ogr")
result = processing.run("native:buffer", {'INPUT': layer, ...})
# Crashes with cryptic error

CORRECT: ALWAYS check isValid() immediately after layer creation.

# CORRECT
layer = QgsVectorLayer("/data/missing_file.gpkg", "data", "ogr")
if not layer.isValid():
    raise RuntimeError(f"Failed to load layer from: /data/missing_file.gpkg")

---

AP-5: Hardcoded Temporary Paths

WRONG: Using platform-specific temporary paths.

# WRONG — /tmp does not exist on Windows
result = processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': '/tmp/buffer_result.gpkg'
})

CORRECT: ALWAYS use 'memory:' for intermediate results or QgsProcessing.TEMPORARY_OUTPUT for disk-based temp files.

# CORRECT — cross-platform
result = processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': 'memory:'  # Or QgsProcessing.TEMPORARY_OUTPUT
})

---

AP-6: Assuming Third-Party Algorithms Exist

WRONG: Using GRASS or SAGA algorithms without checking provider availability.

# WRONG — GRASS may not be installed
result = processing.run("grass7:v.clean", {
    'input': layer,
    'type': [0, 1, 2, 3, 4, 5, 6],
    'tool': [0],
    'output': 'memory:'
})

CORRECT: ALWAYS check algorithm availability and provide a native: fallback.

# CORRECT — check first, fallback to native
registry = QgsApplication.processingRegistry()
if registry.algorithmById('grass7:v.clean') is not None:
    result = processing.run("grass7:v.clean", params)
else:
    # Fallback to native equivalent
    result = processing.run("native:fixgeometries", {
        'INPUT': layer,
        'OUTPUT': 'memory:'
    })

---

AP-7: Missing Spatial Index on Large Datasets

WRONG: Running spatial operations on large datasets without an index.

# WRONG — spatial join on 500,000 features without index takes 100x longer
result = processing.run("native:joinattributesbylocation", {
    'INPUT': large_layer,    # 500,000 features, no spatial index
    'JOIN': another_layer,
    'OUTPUT': 'memory:'
})

CORRECT: ALWAYS create a spatial index before spatial operations on datasets with more than a few thousand features.

# CORRECT — index first
processing.run("native:createspatialindex", {'INPUT': large_layer})
processing.run("native:createspatialindex", {'INPUT': another_layer})

result = processing.run("native:joinattributesbylocation", {
    'INPUT': large_layer,
    'JOIN': another_layer,
    'OUTPUT': 'memory:'
})

---

AP-8: Wrong Analysis Type Selection

WRONG: Using vector overlay when raster analysis is more appropriate.

Scenario: "What percentage of each municipality is forested?"

WRONG approach: Vectorize forest raster → intersect with municipalities → calculate areas
(Vectorizing a high-resolution raster creates millions of polygons, takes hours)

CORRECT approach: Use native:zonalstatisticsfb directly on the raster with municipality polygons
(Completes in seconds, no intermediate conversion needed)

Rule: If one input is already raster and the question is "summarize raster within vector zones", ALWAYS use native:zonalstatisticsfb. NEVER vectorize the raster first.

---

AP-9: Ignoring Geometry Validity

WRONG: Running overlay operations on layers with invalid geometries.

# WRONG — invalid geometries cause silent failures or wrong results
result = processing.run("native:intersection", {
    'INPUT': layer_with_self_intersections,
    'OVERLAY': another_layer,
    'OUTPUT': 'memory:'
})
# Result may be missing features or have corrupt geometries

CORRECT: ALWAYS fix geometries before overlay operations.

# CORRECT — fix first
fixed = processing.run("native:fixgeometries", {
    'INPUT': layer_with_self_intersections,
    'OUTPUT': 'memory:'
})['OUTPUT']

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

---

AP-10: Using Shapefile as Default Output

WRONG: Defaulting to Shapefile format for output.

# WRONG — Shapefile truncates field names to 10 chars, 2GB limit, no multi-layer
processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': '/data/output/result.shp'  # Field names will be truncated
})

CORRECT: ALWAYS use GeoPackage as the default output format.

# CORRECT — GeoPackage has no field name limits, no size limit, supports multi-layer
processing.run("native:buffer", {
    'INPUT': layer,
    'DISTANCE': 100,
    'OUTPUT': '/data/output/result.gpkg'
})

---

AP-11: No Cancellation Check in Processing Loops

WRONG: Processing features in a loop without checking for cancellation.

# WRONG — user cannot cancel, UI appears frozen
def processAlgorithm(self, parameters, context, feedback):
    for feature in source.getFeatures():
        # Heavy processing...
        pass

CORRECT: ALWAYS check feedback.isCanceled() at the top of every loop iteration.

# CORRECT
def processAlgorithm(self, parameters, context, feedback):
    total = 100.0 / source.featureCount() if source.featureCount() else 0
    for current, feature in enumerate(source.getFeatures()):
        if feedback.isCanceled():
            break
        # Heavy processing...
        feedback.setProgress(int(current * total))

---

AP-12: Missing Transform Context in CRS Operations

WRONG: Creating coordinate transforms without project transform context.

# WRONG — ignores datum transformation preferences, silent precision loss
xform = QgsCoordinateTransform(crs_source, crs_target)

CORRECT: ALWAYS include the transform context.

# CORRECT
xform = QgsCoordinateTransform(
    crs_source, crs_target,
    QgsProject.instance().transformContext()
)

examples.md — Analysis Orchestrator

Example 1: Site Suitability Analysis (Vector Overlay Chain)

Question: Find areas within 500m of a train station, within a residential zone, and NOT in a flood risk area.

Workflow Plan: 1. Load: stations (point), zoning (polygon), flood_risk (polygon) 2. Validate + Reproject to EPSG:28992 3. Buffer stations by 500m 4. Intersect buffer with residential zones 5. Difference with flood risk areas 6. Export to GeoPackage

import processing
from qgis.core import (
    QgsVectorLayer, QgsProject, QgsProcessingException, QgsMessageLog
)

# Step 1: Load
stations = QgsVectorLayer("/data/stations.gpkg", "stations", "ogr")
zoning = QgsVectorLayer("/data/zoning.gpkg|layername=residential", "zoning", "ogr")
flood = QgsVectorLayer("/data/flood_risk.gpkg", "flood", "ogr")

for lyr in [stations, zoning, flood]:
    assert lyr.isValid(), f"Layer {lyr.name()} failed to load"

# Step 2: Reproject all to EPSG:28992
target_crs = "EPSG:28992"
reprojected_layers = {}
for name, layer in [("stations", stations), ("zoning", zoning), ("flood", flood)]:
    try:
        result = processing.run("native:reprojectlayer", {
            'INPUT': layer,
            'TARGET_CRS': target_crs,
            'OUTPUT': 'memory:'
        })
        reprojected_layers[name] = result['OUTPUT']
    except QgsProcessingException as e:
        QgsMessageLog.logMessage(f"Reproject failed for {name}: {e}", "Analysis")
        raise

# Step 3: Buffer stations
try:
    buffered = processing.run("native:buffer", {
        'INPUT': reprojected_layers['stations'],
        'DISTANCE': 500,
        'SEGMENTS': 16,
        'DISSOLVE': True,
        'OUTPUT': 'memory:'
    })['OUTPUT']
except QgsProcessingException as e:
    raise RuntimeError(f"Buffer failed: {e}")

# Step 4: Intersect with residential zones
try:
    suitable = processing.run("native:intersection", {
        'INPUT': buffered,
        'OVERLAY': reprojected_layers['zoning'],
        'OUTPUT': 'memory:'
    })['OUTPUT']
except QgsProcessingException as e:
    raise RuntimeError(f"Intersection failed: {e}")

# Step 5: Remove flood risk areas
try:
    final = processing.run("native:difference", {
        'INPUT': suitable,
        'OVERLAY': reprojected_layers['flood'],
        'OUTPUT': '/data/output/suitable_sites.gpkg'
    })['OUTPUT']
except QgsProcessingException as e:
    raise RuntimeError(f"Difference failed: {e}")

QgsMessageLog.logMessage("Site suitability analysis complete", "Analysis")

---

Example 2: Raster-Vector Hybrid (Zonal Statistics)

Question: Calculate average elevation per municipality from a DEM.

Workflow Plan: 1. Load: DEM (raster), municipalities (polygon) 2. Reproject municipalities to match DEM CRS 3. Run zonal statistics 4. Export enriched vector layer

import processing
from qgis.core import QgsVectorLayer, QgsRasterLayer, QgsProcessingException

# Step 1: Load
dem = QgsRasterLayer("/data/dem_25m.tif", "dem")
municipalities = QgsVectorLayer("/data/gemeenten.gpkg", "gemeenten", "ogr")

assert dem.isValid(), "DEM failed to load"
assert municipalities.isValid(), "Municipalities failed to load"

# Step 2: Reproject vector to match raster CRS
try:
    reproj_muni = processing.run("native:reprojectlayer", {
        'INPUT': municipalities,
        'TARGET_CRS': dem.crs(),
        'OUTPUT': 'memory:'
    })['OUTPUT']
except QgsProcessingException as e:
    raise RuntimeError(f"Reproject failed: {e}")

# Step 3: Zonal statistics
try:
    result = processing.run("native:zonalstatisticsfb", {
        'INPUT': reproj_muni,
        'INPUT_RASTER': dem,
        'RASTER_BAND': 1,
        'COLUMN_PREFIX': 'elev_',
        'STATISTICS': [0, 1, 2],  # Count, Sum, Mean
        'OUTPUT': '/data/output/municipalities_elevation.gpkg'
    })
except QgsProcessingException as e:
    raise RuntimeError(f"Zonal statistics failed: {e}")

---

Example 3: Network Analysis (Service Area)

Question: Find all areas reachable within 10 minutes driving from a fire station.

Workflow Plan: 1. Load: road network (line), fire station (point) 2. Validate road network has speed attribute 3. Run service area analysis 4. Export result

import processing
from qgis.core import QgsVectorLayer, QgsProcessingException

# Step 1: Load
roads = QgsVectorLayer("/data/roads.gpkg", "roads", "ogr")
station = QgsVectorLayer("/data/fire_station.gpkg", "station", "ogr")

assert roads.isValid(), "Roads failed to load"
assert station.isValid(), "Station failed to load"

# Step 2: Verify speed field exists
field_names = [f.name() for f in roads.fields()]
assert 'speed_kmh' in field_names, "Road layer MUST have 'speed_kmh' field"

# Step 3: Service area (10 minutes = 600 seconds)
try:
    result = processing.run("native:serviceareafrompoint", {
        'INPUT': roads,
        'START_POINT': f"{station.getFeature(1).geometry().asPoint().x()},"
                       f"{station.getFeature(1).geometry().asPoint().y()}"
                       f" [{roads.crs().authid()}]",
        'TRAVEL_COST': 600,
        'STRATEGY': 1,  # 0=Shortest, 1=Fastest
        'DEFAULT_SPEED': 50,
        'SPEED_FIELD': 'speed_kmh',
        'DEFAULT_DIRECTION': 2,  # Both directions
        'OUTPUT': '/data/output/fire_service_area.gpkg'
    })
except QgsProcessingException as e:
    raise RuntimeError(f"Service area analysis failed: {e}")

---

Example 4: Clustering Analysis

Question: Find clusters of crime incidents using DBSCAN.

Workflow Plan: 1. Load incident points 2. Reproject to projected CRS for accurate distance 3. Create spatial index 4. Run DBSCAN clustering 5. Export with cluster IDs

import processing
from qgis.core import QgsVectorLayer, QgsProcessingException

# Step 1: Load
incidents = QgsVectorLayer("/data/crime_incidents.gpkg", "incidents", "ogr")
assert incidents.isValid(), "Incidents layer failed to load"

# Step 2: Reproject to projected CRS
try:
    reproj = processing.run("native:reprojectlayer", {
        'INPUT': incidents,
        'TARGET_CRS': 'EPSG:28992',
        'OUTPUT': 'memory:'
    })['OUTPUT']
except QgsProcessingException as e:
    raise RuntimeError(f"Reproject failed: {e}")

# Step 3: Spatial index
processing.run("native:createspatialindex", {'INPUT': reproj})

# Step 4: DBSCAN clustering
try:
    result = processing.run("native:dbscanclustering", {
        'INPUT': reproj,
        'MIN_SIZE': 5,       # Minimum cluster size
        'EPS': 200,          # 200 meter radius (projected CRS = meters)
        'OUTPUT': '/data/output/crime_clusters.gpkg'
    })
except QgsProcessingException as e:
    raise RuntimeError(f"DBSCAN failed: {e}")

---

Example 5: Batch Processing Multiple Layers

Question: Buffer and dissolve all layers in a GeoPackage.

import processing
from osgeo import ogr
from qgis.core import QgsVectorLayer, QgsProcessingException, QgsMessageLog

gpkg_path = "/data/input_layers.gpkg"

# Discover all layers in GeoPackage
ds = ogr.Open(gpkg_path)
layer_names = [ds.GetLayerByIndex(i).GetName() for i in range(ds.GetLayerCount())]
ds = None

for name in layer_names:
    uri = f"{gpkg_path}|layername={name}"
    layer = QgsVectorLayer(uri, name, "ogr")

    if not layer.isValid():
        QgsMessageLog.logMessage(f"Skipping invalid layer: {name}", "Batch")
        continue

    try:
        buffered = processing.run("native:buffer", {
            'INPUT': layer,
            'DISTANCE': 100,
            'DISSOLVE': True,
            'OUTPUT': 'memory:'
        })['OUTPUT']

        processing.run("native:savefeatures", {
            'INPUT': buffered,
            'OUTPUT': f"/data/output/buffered_{name}.gpkg"
        })
    except QgsProcessingException as e:
        QgsMessageLog.logMessage(f"Failed for {name}: {e}", "Batch")
        continue

QgsMessageLog.logMessage("Batch processing complete", "Batch")

methods.md — Analysis Orchestrator

Workflow Construction Method

Step 1: Classify the Analysis Task

Every spatial analysis task falls into one of these categories. ALWAYS classify FIRST before selecting algorithms.

Step 2: Determine CRS Requirements

def select_crs(study_area_country, analysis_type):
    """
    ALWAYS call this before starting analysis.
    Returns the appropriate EPSG code.
    """
    # National projected CRS lookup (common examples)
    national_crs = {
        'NL': 'EPSG:28992',   # Amersfoort / RD New
        'BE': 'EPSG:31370',   # Belgian Lambert 72
        'DE': 'EPSG:25832',   # ETRS89 / UTM 32N
        'UK': 'EPSG:27700',   # British National Grid
        'FR': 'EPSG:2154',    # RGF93 / Lambert-93
        'US': 'EPSG:5070',    # NAD83 / Conus Albers
        'AU': 'EPSG:3577',    # GDA94 / Australian Albers
    }

    if analysis_type == 'web_display':
        return 'EPSG:3857'
    elif analysis_type == 'global_equal_area':
        return 'EPSG:6933'
    elif study_area_country in national_crs:
        return national_crs[study_area_country]
    else:
        # Fall back to UTM zone
        return determine_utm_zone(study_area_centroid)

Step 3: Build the Algorithm Chain

ALWAYS follow this ordering:

1. Load — Load all input data 2. Validate — Fix geometries, check CRS validity 3. Reproject — Bring all layers to common CRS 4. Index — Create spatial indexes for large datasets 5. Pre-process — Any data cleaning (dissolve, simplify, extract subset) 6. Analyze — Core analysis algorithm(s) 7. Post-process — Join results, calculate additional fields 8. Export — Save to target format

Step 4: Select Algorithms by Task

Proximity Tasks

Overlay Tasks

Aggregation Tasks

Terrain Tasks

Workflow Validation Method

After constructing an algorithm chain, ALWAYS validate:

1. CRS chain: Verify every algorithm receives layers in compatible CRS 2. Field propagation: Verify join/overlay operations propagate needed fields 3. Geometry type chain: Verify output geometry type of step N matches expected input of step N+1 4. Memory management: Verify intermediate results use 'memory:', only final output writes to disk 5. Error propagation: Verify every processing.run() is wrapped in try/except