Practical QGIS & PyQGIS Workflow Tutorials - Hydrology & Siting

Tutorial 1: Automated Hydrological Watershed Delineation from DEM

Hydrological catchment modeling requires extracting continuous drainage channels and watershed boundary polygons from high-resolution raster elevation data. In civil drainage design, flood hazard mitigation, and agricultural runoff modeling, manual delineation of ridge lines introduces unacceptable subjectivity. The following multi-step workflow utilizes native QGIS and GRASS algorithms to transform a raw Digital Elevation Model into vector catchment basins.

1

Sink Identification and Surface Pit Filling

Raw DEMs derived from photogrammetry or spaceborne radar (such as Copernicus DEM or SRTM) contain spurious depressions (sinks) where water cannot escape. To establish a continuous hydraulic gradient:

Verify your DEM is projected into a planar coordinate reference system (e.g., UTM Zone 33N, EPSG:32633). Never execute hydrological tools on geographic coordinates (EPSG:4326).
Open Processing Toolbox > SAGA > Terrain Analysis - Hydrology > Fill Sinks (Wang & Liu).
Specify minimum slope threshold: 0.01 degrees. This guarantees a slight downward hydraulic gradient across flat lakebeds and flooded valleys.
Output: Filled DEM.

2

Flow Direction and D8 Flow Accumulation

Water follows the path of steepest downhill descent. The deterministic 8-neighbor (D8) algorithm calculates flow direction by routing runoff from each pixel to one of its eight neighboring cells:

Open Processing Toolbox > GRASS > Raster > r.watershed.
Supply the Filled DEM as input elevation.
Set minimum size of exterior watershed basin: 500 cells (adjust based on cell resolution, e.g., 500 cells of 10m DEM = 50,000 m2 minimum catchment).
Enable calculation of Flow Accumulation (number of upstream contributing cells) and Drainage Direction.

3

Channel Network Extraction and Stream Order

Channels initiate when upslope accumulation exceeds a critical physical threshold:

Open Raster > Raster Calculator.
Enter the boolean channel extraction expression:
"accumulation@1" >= 1000
This assigns value 1 to active streams and 0 to dry hillsides.
Execute Processing Toolbox > GRASS > Raster > r.stream.order to assign Strahler or Shreve hierarchical ordering numbers to each tributary.
Convert raster streams to vector polylines using Processing Toolbox > GRASS > Vector > r.to.vect with feature type set to line.

4

Catchment Basin Vectorization

To generate closed polygon boundaries for municipal watershed management:

Open Processing Toolbox > GRASS > Raster > r.watershed and check the Unique label for each drainage basin output.
Convert the resulting basin raster to vector polygons with Raster > Conversion > Polygonize (Raster to Vector).
Execute Vector Geometry > Fix Geometries to resolve any self-intersecting polygon boundary rings.

Tutorial 2: Multi-Criteria Land Suitability Analysis for Utility Solar

Siting utility-scale ground-mounted photovoltaic power plants requires combining multiple spatial criteria: topography, solar irradiation, distance to infrastructure, and environmental zoning exclusions. Multi-Criteria Evaluation (MCE) via Weighted Linear Combination (WLC) normalizes continuous geographic layers into a unified suitability index between 0.0 (unsuitable) and 1.0 (optimal).

Suitability Criteria Matrix

Evaluation Factor	Optimal Condition	Exclusion Boundary	Normalization Function	Assigned Weight
Terrain Slope	0 to 3 degrees (flat ground)	> 12 degrees (unbuildable)	Linear decay: $(12 - Slope) / 12$	25% (0.25)
Aspect / Orientation	160 to 200 degrees (True South)	315 to 45 degrees (North slopes)	Cosine bell curve centered at 180 deg	20% (0.20)
Grid Distance	< 500 meters to 20 kV line	> 5,000 meters (prohibitive CAPEX)	Inverse distance decay: $(5000 - Dist) / 5000$	30% (0.30)
Road Accessibility	< 300 meters from arterial roads	> 2,500 meters	Linear decay: $(2500 - Dist) / 2500$	15% (0.15)
Environmental Buffers	Non-protected farmland / wasteland	Natura 2000, wetlands, river buffers	Binary mask: 0 (inside) or 1 (outside)	Boolean Exclusion

Raster Calculator Execution Recipe

Once individual factor rasters are generated and resampled to a consistent 10-meter grid resolution, the master Weighted Linear Combination is computed in the QGIS Raster Calculator:

        # Master Solar PV Suitability Expression in QGIS Raster Calculator
(
  (0.25 * ("norm_slope@1")) +
  (0.20 * ("norm_aspect@1")) +
  (0.30 * ("norm_grid_dist@1")) +
  (0.15 * ("norm_road_dist@1"))
) * ("env_exclusion_mask@1")
      

The resulting raster layer ranks the territory on a continuous scale from 0.0 to 0.90. Pixels scoring above 0.70 represent prime development parcels suitable for ground-mounted solar plant feasibility surveys.

Tutorial 3: Headless PyQGIS Batch Automation Script

Geospatial production pipelines often require processing hundreds of satellite scenes or regional cadastral tiles without opening the desktop GUI. The PyQGIS library allows Python scripts to execute the full suite of QGIS C++ processing algorithms from the terminal.

The following standalone script initializes a headless QGIS instance, iterates through an input directory of GeoTIFF rasters, reprojects them to UTM Zone 33N (EPSG:32633) with high-order cubic resampling, clips each raster against a project boundary shapefile, and generates Cloud-Optimized GeoTIFF (COG) outputs:

headless_pyqgis_pipeline.py Python 3.10+

import os
import glob
import sys

# 1. Configure QGIS Installation Environment Paths
# On Linux: /usr / On Windows: C:/OSGeo4W/apps/qgis-ltr
QGIS_PREFIX = os.environ.get("QGIS_PREFIX_PATH", "/usr")

from qgis.core import (
    QgsApplication,
    QgsCoordinateReferenceSystem,
    QgsVectorLayer,
    QgsRasterLayer
)

# Initialize QgsApplication without GUI enabled (first argument False)
qgs = QgsApplication([], False)
qgs.setPrefixPath(QGIS_PREFIX, True)
qgs.initQgis()

# Import processing module after initializing QgsApplication
sys.path.append(os.path.join(QGIS_PREFIX, "share", "qgis", "python", "plugins"))
import processing
from processing.core.Processing import Processing
Processing.initialize()

def batch_process_rasters(input_dir, mask_shapefile, output_dir):
    """
    Automated batch reprojection and masking pipeline utilizing PyQGIS.
    """
    os.makedirs(output_dir, exist_ok=True)
    target_crs = QgsCoordinateReferenceSystem("EPSG:32633")
    
    raster_files = glob.glob(os.path.join(input_dir, "*.tif"))
    print(f"Discovered {len(raster_files)} raster files to process.")

    for idx, r_path in enumerate(raster_files, start=1):
        filename = os.path.basename(r_path)
        warped_temp = os.path.join(output_dir, f"temp_warped_{filename}")
        final_cog = os.path.join(output_dir, f"proc_{filename}")
        
        print(f"[{idx}/{len(raster_files)}] Processing {filename}...")

        # Step A: Reproject raster to target projected CRS
        warp_params = {
            "INPUT": r_path,
            "TARGET_CRS": target_crs,
            "RESAMPLING": 2, # Cubic convolution
            "NODATA": -9999,
            "TARGET_RESOLUTION": 10.0,
            "OPTIONS": "COMPRESS=DEFLATE|TILED=YES",
            "OUTPUT": warped_temp
        }
        processing.run("gdal:warpreproject", warp_params)

        # Step B: Clip raster by vector polygon boundary mask
        clip_params = {
            "INPUT": warped_temp,
            "MASK": mask_shapefile,
            "SOURCE_CRS": target_crs,
            "TARGET_CRS": target_crs,
            "NODATA": -9999,
            "ALPHA_BAND": False,
            "CROP_TO_CUTLINE": True,
            "KEEP_RESOLUTION": True,
            "OPTIONS": "COMPRESS=DEFLATE|TILED=YES|PREDICTOR=2",
            "OUTPUT": final_cog
        }
        processing.run("gdal:cliprasterbymasklayer", clip_params)

        # Clean intermediate scratch file
        if os.path.exists(warped_temp):
            os.remove(warped_temp)

    print("Batch processing pipeline completed successfully.")

if __name__ == "__main__":
    try:
        batch_process_rasters(
            input_dir="/data/raw_satellite_tiles",
            mask_shapefile="/data/municipal_boundary.shp",
            output_dir="/data/output_cogs"
        )
    finally:
        # Exit QGIS runtime cleanly to flush file handles and memory locks
        qgs.exitQgis()

Tutorial 4: Building Reusable Graphical Modeler Pipelines

The QGIS Graphical Modeler provides a visual Directed Acyclic Graph (DAG) interface for chaining multiple processing algorithms into a single custom tool. Once constructed, a model can be executed from the Processing Toolbox like any native tool, run in batch mode over multiple input files, or exported to a Python script.

Step-by-Step Model Construction Workflow

Open Processing > Graphical Modeler.
Assign a model name (e.g., Cadastral_Parcel_Zoning_Validator) and group name (e.g., Planning_Automation).
From the Inputs tab on the left panel, drag a Vector Layer input named Input Parcels (Geometry type: Polygon).
Drag a second Vector Layer input named Zoning Restrictions (Geometry type: Polygon).
Drag a Number input named Buffer Distance (Meters) with minimum 0 and default 50.
Switch to the Algorithms tab. Search for Buffer and drag it onto the canvas. Set input layer to Zoning Restrictions and distance parameter to Buffer Distance (Meters). Name the output Buffered Restrictions.
Search for Intersection and drag it onto the canvas. Set Input layer to Input Parcels and Overlay layer to Algorithm output: Buffered Restrictions.
Label the final output Non_Compliant_Parcels.
Click the Validate Model button. If green, click Save Model to Project. The custom tool immediately appears in your QGIS Processing Toolbox under the Project models category.

Common Pitfalls and Field Troubleshooting in QGIS

Desktop GIS workflows frequently stall due to obscure geometric, coordinate, or memory errors:

1. On-The-Fly Reprojection Distortion Trap

When you add a raster in geographic coordinates (EPSG:4326) and vector data in UTM (EPSG:32633), QGIS projects them together seamlessly on screen. However, spatial algorithms execute on the layer's native coordinate reference system, not the canvas projection. Running Slope on an unprojected raster treats one degree of longitude ($~111\text{ km}$ at the equator) as equivalent to 1 vertical meter, producing slopes that erroneously read $89.9^\circ$ across completely flat farmland. Always re-project layers to a projected metric CRS before running spatial tools.

2. Self-Intersecting Polygons and Bowtie Errors

Vector digitizing errors frequently create self-intersecting boundary rings (figure-eight or bowtie shapes). These geometries cause tools like Intersection, Union, or Dissolve to abort with an GEOS Geoprocessing Error: TopologyException. Resolve this instantly by inserting the Fix Geometries algorithm into your workflow ahead of spatial overlays.

3. Nodata Value Bleeding and Black Borders

Reprojecting a rectangular GeoTIFF into a rotated UTM zone creates black empty wedges along the image edges. If the nodata value is unassigned, QGIS interprets these black pixels as valid zero-reflectance data, contaminating raster statistics and index histograms. In Warp (Reproject), explicitly define -dstnodata -9999 or -dstnodata 0 to guarantee complete transparency for non-data boundaries.

Frequently Asked Questions on QGIS & PyQGIS Automation

Why does QGIS on-the-fly CRS transformation cause measurement errors in raster spatial analysis?

On-the-fly CRS transformation dynamically re-projects vector and raster layers in the map canvas for visual display without altering the underlying raw file coordinates. When executing raster analysis tools (such as slope, aspect, or buffer calculations) on a layer in geographic coordinates (EPSG:4326), horizontal units are measured in angular degrees while elevation is in meters. This mismatch distorts slope calculations by orders of magnitude unless layers are permanently re-projected into an appropriate conformal projected coordinate system (such as UTM).

How does pit filling prevent disconnected flow networks in DEM hydrological modeling?

Raw digital elevation models frequently contain spurious single-pixel depressions or data artifacts caused by sensor noise, interpolation errors, or dense vegetation canopy occlusions. If uncorrected, simulated surface water flows into these artificial sinks and cannot exit, terminating the flow path prematurely. Pit filling algorithms raise the elevation of these depression pixels to the level of their lowest pour point, ensuring continuous hydraulic flow across the digital surface.

How can self-intersecting polygon geometry errors be resolved automatically in QGIS?

Self-intersecting polygons or bowtie errors violate OGC Simple Features specifications and cause spatial overlay algorithms to fail. In QGIS, navigate to Processing Toolbox > Vector Geometry > Fix Geometries. This tool applies modern topological repair algorithms that split self-intersecting rings into valid multi-part geometries without modifying vertex coordinates.

What is the primary architectural difference between running PyQGIS inside QGIS and running headless PyQGIS scripts?

Running scripts in the QGIS Python Console operates within an already initialized GUI application context with an active canvas, project instance, and loaded layer tree. Headless PyQGIS scripts run from a standard operating system shell without launching the GUI. Headless execution requires explicitly importing QgsApplication, setting the QGIS installation prefix path, and invoking qgs.initQgis() before processing, followed by qgs.exitQgis() on completion.

How do virtual rasters (VRT) in QGIS optimize working with massive multi-gigabyte aerial orthomosaic datasets?

A Virtual Raster (VRT) is an XML file that catalogs and references hundreds of individual GeoTIFF tiles as a single continuous mosaic layer without copying or duplicating pixel data. QGIS reads the VRT file and dynamically streams only the tiles and pyramid overviews within the active viewport, drastically reducing memory consumption and project loading latency.