Qgis Impl Raster Analysis

qgis-impl-raster-analysis

Quick Reference

Raster Layer Properties

Renderer Types

Terrain Analysis Algorithms

Key GDAL Processing Algorithms

Critical Warnings

NEVER ignore NoData values in raster calculations. NoData pixels propagate through arithmetic -- a single NoData input produces NoData output. ALWAYS check and set NoData handling explicitly.

NEVER use QgsRasterCalculator without verifying that entry.ref matches the reference string in the expression (e.g., 'dem@1'). Mismatched references cause silent failures with all-zero output.

NEVER forget rlayer.triggerRepaint() after calling rlayer.setRenderer(). The map canvas does NOT update automatically.

ALWAYS check rlayer.isValid() after loading a raster layer. Invalid layers produce cryptic downstream errors.

ALWAYS verify processCalculation() returns 0 (success). Non-zero return values indicate errors but provide no error message.

ALWAYS match the data type in QgsContrastEnhancement to the raster band's actual data type via provider.dataType(bandNo).

NEVER assume band numbering starts at 0. QGIS bands are 1-indexed. Band 0 does NOT exist.

ALWAYS use 'memory:' (with colon) for in-memory processing outputs, NOT 'memory' or 'TEMPORARY_OUTPUT' for GDAL algorithms.

Decision Tree

Which Raster Analysis Approach?

Need raster math (add, multiply, threshold)?
├── Simple expression → QgsRasterCalculator (native PyQGIS)
├── Complex multi-band → gdal:rastercalculator (processing.run)
└── Pixel-by-pixel custom logic → QgsRasterBlock read/write loop

Need terrain derivatives?
├── Hillshade → processing.run("native:hillshade", ...) or "gdal:hillshade"
├── Slope → processing.run("native:slope", ...) or "gdal:slope"
├── Aspect → processing.run("native:aspect", ...) or "gdal:aspect"
├── Contours → processing.run("gdal:contour", ...)
└── Ruggedness → processing.run("native:ruggednessindex", ...)

Need format conversion?
├── Reproject → processing.run("gdal:warpreproject", ...)
├── Change format → processing.run("gdal:translate", ...)
├── Raster to vector → processing.run("gdal:polygonize", ...)
└── Vector to raster → processing.run("gdal:rasterize", ...)

Need visualization?
├── Single-band grayscale → QgsSingleBandGrayRenderer
├── Color ramp (elevation/temperature) → QgsSingleBandPseudoColorRenderer
├── RGB satellite → QgsMultiBandColorRenderer
├── Classified categories → QgsPalettedRasterRenderer
└── Live hillshade → QgsHillshadeRenderer

Native vs GDAL Algorithms?

Use native: algorithms when:
├── Simpler parameter set is sufficient
├── In-memory output needed (memory:)
└── Fewer dependencies preferred

Use gdal: algorithms when:
├── Need advanced options (SCALE, AS_PERCENT, creation options)
├── Need specific GDAL creation options (COMPRESS=DEFLATE)
├── Working with large rasters (GDAL is optimized for I/O)
└── Need format-specific features

Essential Patterns

Load and Validate Raster Layer

from qgis.core import QgsRasterLayer, QgsProject

rlayer = QgsRasterLayer("/path/to/dem.tif", "DEM")
if not rlayer.isValid():
    raise ValueError(f"Raster layer failed to load: {rlayer.error().message()}")

QgsProject.instance().addMapLayer(rlayer)

Query Pixel Values

# Single band, single point
provider = rlayer.dataProvider()
val, result = provider.sample(QgsPointXY(20.50, -34.0), 1)  # point, band
if result:
    print(f"Value: {val}")

# All bands at a point
from qgis.core import QgsRaster
ident = provider.identify(QgsPointXY(20.5, -34.0), QgsRaster.IdentifyFormatValue)
if ident.isValid():
    print(ident.results())  # {1: 323.0, 2: 127.0, ...}

Raster Calculator (Map Algebra)

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

entries = []
entry = QgsRasterCalculatorEntry()
entry.ref = 'dem@1'
entry.raster = rlayer
entry.bandNumber = 1
entries.append(entry)

calc = QgsRasterCalculator(
    '"dem@1" * 2 + 100',         # Expression (refs in double quotes)
    '/path/to/output.tif',       # Output path
    'GTiff',                     # Output format
    rlayer.extent(),             # Output extent
    rlayer.width(),              # Output columns
    rlayer.height(),             # Output rows
    entries                      # List of entries
)

result = calc.processCalculation()
if result != 0:
    raise RuntimeError(f"Raster calculation failed with code {result}")

Band Statistics

from qgis.core import QgsRasterBandStats

stats = rlayer.dataProvider().bandStatistics(
    1,                           # Band number (1-indexed)
    QgsRasterBandStats.All,      # Compute all statistics
    rlayer.extent(),             # Extent to analyze
    0                            # Sample size (0 = all pixels)
)
print(f"Min: {stats.minimumValue}, Max: {stats.maximumValue}")
print(f"Mean: {stats.mean}, StdDev: {stats.stdDev}")
print(f"Sum: {stats.sum}, Range: {stats.range}")

Raster Block Access (Pixel-Level)

provider = rlayer.dataProvider()
block = provider.block(1, rlayer.extent(), rlayer.width(), rlayer.height())

for row in range(block.height()):
    for col in range(block.width()):
        if not block.isNoData(row, col):
            value = block.value(row, col)
            # Process pixel value

Common Operations

Terrain Analysis via Processing

import processing

# Hillshade
result = processing.run("native:hillshade", {
    'INPUT': dem_layer,
    'Z_FACTOR': 1.0,
    'AZIMUTH': 315,
    'V_ANGLE': 45,
    'OUTPUT': 'memory:'
})
hillshade_layer = result['OUTPUT']

# Slope
result = processing.run("native:slope", {
    'INPUT': dem_layer,
    'Z_FACTOR': 1.0,
    'OUTPUT': 'memory:'
})

# Aspect
result = processing.run("native:aspect", {
    'INPUT': dem_layer,
    'Z_FACTOR': 1.0,
    'OUTPUT': 'memory:'
})

GDAL Terrain Analysis (Advanced Options)

# GDAL hillshade with extra controls
result = processing.run("gdal:hillshade", {
    'INPUT': '/data/dem.tif',
    'BAND': 1,
    'Z_FACTOR': 1,
    'SCALE': 1,
    'AZIMUTH': 315,
    'ALTITUDE': 45,
    'OUTPUT': '/output/hillshade.tif'
})

# GDAL slope with percent option
result = processing.run("gdal:slope", {
    'INPUT': '/data/dem.tif',
    'BAND': 1,
    'SCALE': 1,
    'AS_PERCENT': False,
    'OUTPUT': '/output/slope.tif'
})

Contour Lines from DEM

result = processing.run("gdal:contour", {
    'INPUT': '/data/dem.tif',
    'BAND': 1,
    'INTERVAL': 10,
    'FIELD_NAME': 'ELEV',
    'OUTPUT': '/output/contours.gpkg'
})

Raster Reprojection

result = processing.run("gdal:warpreproject", {
    'INPUT': rlayer,
    'SOURCE_CRS': 'EPSG:4326',
    'TARGET_CRS': 'EPSG:28992',
    'RESAMPLING': 0,        # 0=Nearest, 1=Bilinear, 2=Cubic
    'NODATA': -9999,
    'TARGET_RESOLUTION': 25,
    'OUTPUT': '/output/reprojected.tif'
})

Set Single Band Pseudocolor Renderer

from qgis.core import (QgsSingleBandPseudoColorRenderer,
                        QgsRasterShader, QgsColorRampShader)
from qgis.PyQt.QtGui import QColor

fcn = QgsColorRampShader()
fcn.setColorRampType(QgsColorRampShader.Interpolated)
lst = [
    QgsColorRampShader.ColorRampItem(0, QColor(0, 255, 0), "Low"),
    QgsColorRampShader.ColorRampItem(500, QColor(255, 255, 0), "Medium"),
    QgsColorRampShader.ColorRampItem(1000, QColor(255, 0, 0), "High"),
]
fcn.setColorRampItemList(lst)

shader = QgsRasterShader()
shader.setRasterShaderFunction(fcn)

renderer = QgsSingleBandPseudoColorRenderer(rlayer.dataProvider(), 1, shader)
rlayer.setRenderer(renderer)
rlayer.triggerRepaint()

Set Multiband Color Renderer (RGB)

from qgis.core import QgsMultiBandColorRenderer

renderer = QgsMultiBandColorRenderer(rlayer.dataProvider(), 3, 2, 1)
# Arguments: provider, redBand, greenBand, blueBand
rlayer.setRenderer(renderer)
rlayer.triggerRepaint()

Raster to Vector (Polygonize)

result = processing.run("gdal:polygonize", {
    'INPUT': '/data/classified.tif',
    'BAND': 1,
    'FIELD': 'DN',
    'EIGHT_CONNECTEDNESS': False,
    'OUTPUT': '/output/polygons.gpkg'
})

Vector to Raster (Rasterize)

result = processing.run("gdal:rasterize", {
    'INPUT': '/data/buildings.gpkg',
    'FIELD': 'height',
    'UNITS': 1,       # 1=Georeferenced units
    'WIDTH': 1,       # Pixel width in georef units
    'HEIGHT': 1,      # Pixel height in georef units
    'EXTENT': layer.extent(),
    'OUTPUT': '/output/buildings_raster.tif'
})

Merge and Mosaic Rasters

# Virtual raster (lightweight, no data copy)
result = processing.run("gdal:buildvirtualraster", {
    'INPUT': ['/data/tile1.tif', '/data/tile2.tif'],
    'RESOLUTION': 0,  # 0=Average, 1=Highest, 2=Lowest
    'OUTPUT': '/output/mosaic.vrt'
})

# Physical merge
result = processing.run("gdal:merge", {
    'INPUT': ['/data/tile1.tif', '/data/tile2.tif'],
    'NODATA_INPUT': -9999,
    'NODATA_OUTPUT': -9999,
    'OUTPUT': '/output/merged.tif'
})

Reference Links

• references/methods.md -- API signatures for QgsRasterCalculator, renderers, terrain classes
• references/examples.md -- Complete raster analysis workflows
• references/anti-patterns.md -- Raster analysis pitfalls and fixes

Official Sources

• https://qgis.org/pyqgis/master/core/QgsRasterLayer.html
• https://qgis.org/pyqgis/master/analysis/QgsRasterCalculator.html
• https://qgis.org/pyqgis/master/core/QgsRasterRenderer.html
• https://docs.qgis.org/latest/en/docs/pyqgis_developer_cookbook/raster.html
• https://docs.qgis.org/latest/en/docs/user_manual/processing_algs/gdal/rasteranalysis.html