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

Anti-Patterns (Raster Analysis)

1. Ignoring NoData in Raster Calculator

# WRONG: Division without NoData guard — produces NaN or infinity at NoData pixels
expression = '"nir@1" / "red@1"'

# CORRECT: Guard against zero and NoData propagation
expression = '("nir@1" - "red@1") / ("nir@1" + "red@1" + 0.0001)'

WHY: NoData values propagate through arithmetic. Dividing by a band that contains zero or NoData pixels produces corrupt output. ALWAYS add a small epsilon or use conditional expressions to handle edge cases.

---

2. Mismatched Entry Reference Names

# WRONG: entry.ref does not match the expression reference
entry = QgsRasterCalculatorEntry()
entry.ref = 'my_dem@1'
entry.raster = rlayer
entry.bandNumber = 1

calc = QgsRasterCalculator(
    '"dem@1" * 2',  # "dem@1" does NOT match "my_dem@1"
    output_path, 'GTiff', extent, width, height, [entry]
)
# Result: all-zero output with NO error message

# CORRECT: entry.ref MUST match the reference in the expression exactly
entry.ref = 'dem@1'
calc = QgsRasterCalculator('"dem@1" * 2', ...)

WHY: QgsRasterCalculator silently produces all-zero output when references do not match. There is no error or warning. ALWAYS verify that entry.ref strings match the double-quoted references in the expression.

---

3. Forgetting triggerRepaint After Setting Renderer

# WRONG: Renderer set but canvas not updated
renderer = QgsSingleBandPseudoColorRenderer(provider, 1, shader)
rlayer.setRenderer(renderer)
# Map canvas still shows old rendering

# CORRECT: ALWAYS call triggerRepaint
rlayer.setRenderer(renderer)
rlayer.triggerRepaint()

WHY: setRenderer() changes the internal renderer object but does NOT trigger a canvas redraw. Without triggerRepaint(), the user sees stale visualization.

---

4. Using Band 0 Instead of Band 1

# WRONG: Band numbering starts at 1, not 0
stats = provider.bandStatistics(0, QgsRasterBandStats.All)  # Crashes or returns garbage

# CORRECT: Bands are 1-indexed
stats = provider.bandStatistics(1, QgsRasterBandStats.All)

WHY: QGIS raster bands are 1-indexed. Band 0 does NOT exist. Using band 0 causes crashes or returns invalid statistics with no clear error message.

---

5. Wrong Data Type in ContrastEnhancement

# WRONG: Hardcoding data type instead of reading from provider
ce = QgsContrastEnhancement(Qgis.Byte)  # Assumes 8-bit, but DEM may be Float32

# CORRECT: Read actual data type from provider
ce = QgsContrastEnhancement(rlayer.dataProvider().dataType(1))

WHY: Mismatched data types cause incorrect contrast stretching. A Float32 DEM treated as Byte clips all values to 0-255, destroying the data range.

---

6. Not Checking processCalculation Return Value

# WRONG: Assuming calculation succeeded
calc.processCalculation()
result_layer = QgsRasterLayer(output_path, "Result")  # May load empty/corrupt file

# CORRECT: Check return code
result_code = calc.processCalculation()
if result_code != 0:
    error_map = {
        1: "Failed to create output file",
        2: "Input layer error",
        3: "Calculation canceled",
        4: "Expression parser error",
        5: "Memory allocation error",
        6: "Band number error"
    }
    raise RuntimeError(f"Raster calc failed: {error_map.get(result_code, 'Unknown')}")

WHY: processCalculation() returns an integer error code, not a boolean. Code 0 means success. Non-zero values indicate specific failures, but the method does NOT raise exceptions.

---

7. Using 'TEMPORARY_OUTPUT' for GDAL Algorithms

# WRONG: TEMPORARY_OUTPUT is for native algorithms, not GDAL
result = processing.run("gdal:hillshade", {
    'INPUT': dem_layer,
    'OUTPUT': 'TEMPORARY_OUTPUT'  # May fail or produce unexpected path
})

# CORRECT: Use explicit path for GDAL algorithms
result = processing.run("gdal:hillshade", {
    'INPUT': dem_layer,
    'OUTPUT': '/tmp/hillshade.tif'
})

# OR use memory: for native algorithms
result = processing.run("native:hillshade", {
    'INPUT': dem_layer,
    'OUTPUT': 'memory:'
})

WHY: GDAL algorithms write to disk via GDAL drivers. They do NOT support in-memory outputs the same way native algorithms do. ALWAYS provide a file path for GDAL algorithms, or use the equivalent native: algorithm with 'memory:'.

---

8. Not Validating Raster Layer Before Use

# WRONG: Using layer without validity check
rlayer = QgsRasterLayer("/nonexistent/path.tif", "DEM")
stats = rlayer.dataProvider().bandStatistics(1)  # Crashes: provider is None

# CORRECT: ALWAYS validate before use
rlayer = QgsRasterLayer("/data/dem.tif", "DEM")
if not rlayer.isValid():
    raise ValueError(f"Failed to load: {rlayer.error().message()}")

WHY: An invalid QgsRasterLayer has a None data provider. Any method call on the provider causes an AttributeError or segfault with no useful error context.

---

9. Resampling Method Mismatch

# WRONG: Using nearest neighbor for continuous DEM data
result = processing.run("gdal:warpreproject", {
    'INPUT': dem_layer,
    'TARGET_CRS': 'EPSG:28992',
    'RESAMPLING': 0,  # Nearest neighbor — creates staircase artifacts in DEMs
    'OUTPUT': '/output/dem_rd.tif'
})

# CORRECT: Use bilinear or cubic for continuous data
result = processing.run("gdal:warpreproject", {
    'INPUT': dem_layer,
    'TARGET_CRS': 'EPSG:28992',
    'RESAMPLING': 1,  # Bilinear — smooth interpolation for continuous data
    'OUTPUT': '/output/dem_rd.tif'
})

WHY: Nearest neighbor resampling preserves exact pixel values (correct for categorical data like land use), but creates visible staircase artifacts in continuous data like DEMs and temperature grids. ALWAYS use bilinear (1) or cubic (2) resampling for continuous rasters, and nearest neighbor (0) for categorical/classified rasters.

---

10. Loading Result Without Checking Output Path

# WRONG: Assuming result dict key matches expected output
result = processing.run("gdal:hillshade", {
    'INPUT': '/data/dem.tif',
    'OUTPUT': '/output/hillshade.tif'
})
layer = QgsRasterLayer('/output/hillshade.tif', 'Hillshade')  # Fragile

# CORRECT: Use the result dictionary to get actual output path
result = processing.run("gdal:hillshade", {
    'INPUT': '/data/dem.tif',
    'OUTPUT': '/output/hillshade.tif'
})
layer = QgsRasterLayer(result['OUTPUT'], 'Hillshade')
if not layer.isValid():
    raise ValueError("Hillshade output failed to load")

WHY: The actual output path may differ from the requested path (e.g., GDAL may append extensions or modify the path). ALWAYS use result['OUTPUT'] to get the confirmed output path.

Working Code Examples (Raster Analysis)

Example 1: Load DEM and Compute Statistics

from qgis.core import QgsRasterLayer, QgsRasterBandStats, QgsProject

# Load raster
rlayer = QgsRasterLayer("/data/dem.tif", "DEM")
if not rlayer.isValid():
    raise ValueError("Raster layer failed to load")

QgsProject.instance().addMapLayer(rlayer)

# Get statistics
provider = rlayer.dataProvider()
stats = provider.bandStatistics(1, QgsRasterBandStats.All, rlayer.extent(), 0)

print(f"Dimensions: {rlayer.width()} x {rlayer.height()} pixels")
print(f"Bands: {rlayer.bandCount()}")
print(f"Min elevation: {stats.minimumValue}")
print(f"Max elevation: {stats.maximumValue}")
print(f"Mean elevation: {stats.mean}")
print(f"Std deviation: {stats.stdDev}")

---

Example 2: NDVI Calculation with Raster Calculator

from qgis.core import QgsRasterLayer
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

# Load multispectral image
rlayer = QgsRasterLayer("/data/sentinel2.tif", "Sentinel2")
if not rlayer.isValid():
    raise ValueError("Layer failed to load")

# Band 4 = Red, Band 8 = NIR (Sentinel-2)
entries = []

red_entry = QgsRasterCalculatorEntry()
red_entry.ref = 'sentinel@4'
red_entry.raster = rlayer
red_entry.bandNumber = 4
entries.append(red_entry)

nir_entry = QgsRasterCalculatorEntry()
nir_entry.ref = 'sentinel@8'
nir_entry.raster = rlayer
nir_entry.bandNumber = 8
entries.append(nir_entry)

# NDVI = (NIR - Red) / (NIR + Red)
expression = '("sentinel@8" - "sentinel@4") / ("sentinel@8" + "sentinel@4")'

calc = QgsRasterCalculator(
    expression,
    '/output/ndvi.tif',
    'GTiff',
    rlayer.extent(),
    rlayer.width(),
    rlayer.height(),
    entries
)

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

# Load result
ndvi_layer = QgsRasterLayer('/output/ndvi.tif', 'NDVI')
QgsProject.instance().addMapLayer(ndvi_layer)

---

Example 3: Complete Terrain Analysis Pipeline

import processing
from qgis.core import QgsRasterLayer, QgsProject

dem_layer = QgsRasterLayer("/data/dem.tif", "DEM")
if not dem_layer.isValid():
    raise ValueError("DEM failed to load")

QgsProject.instance().addMapLayer(dem_layer)

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

# Slope
slope_result = processing.run("native:slope", {
    'INPUT': dem_layer,
    'Z_FACTOR': 1.0,
    'OUTPUT': 'memory:'
})
QgsProject.instance().addMapLayer(slope_result['OUTPUT'])
slope_result['OUTPUT'].setName('Slope')

# Aspect
aspect_result = processing.run("native:aspect", {
    'INPUT': dem_layer,
    'Z_FACTOR': 1.0,
    'OUTPUT': 'memory:'
})
QgsProject.instance().addMapLayer(aspect_result['OUTPUT'])
aspect_result['OUTPUT'].setName('Aspect')

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

---

Example 4: DEM Visualization with Pseudocolor Renderer

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

rlayer = QgsRasterLayer("/data/dem.tif", "DEM")
if not rlayer.isValid():
    raise ValueError("DEM failed to load")

QgsProject.instance().addMapLayer(rlayer)

# Get actual min/max from data
stats = rlayer.dataProvider().bandStatistics(1, QgsRasterBandStats.All, rlayer.extent(), 0)

# Build color ramp shader
fcn = QgsColorRampShader()
fcn.setColorRampType(QgsColorRampShader.Interpolated)
fcn.setColorRampItemList([
    QgsColorRampShader.ColorRampItem(stats.minimumValue, QColor(0, 100, 0), "Valley"),
    QgsColorRampShader.ColorRampItem(stats.mean, QColor(255, 255, 100), "Mid"),
    QgsColorRampShader.ColorRampItem(stats.maximumValue, QColor(139, 69, 19), "Peak"),
])

shader = QgsRasterShader()
shader.setRasterShaderFunction(fcn)

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

---

Example 5: Satellite Image RGB Composite

from qgis.core import (QgsRasterLayer, QgsMultiBandColorRenderer,
                        QgsContrastEnhancement, QgsProject)

rlayer = QgsRasterLayer("/data/satellite.tif", "Satellite")
if not rlayer.isValid():
    raise ValueError("Layer failed to load")

QgsProject.instance().addMapLayer(rlayer)

# True color: R=Band3, G=Band2, B=Band1 (common for Landsat)
renderer = QgsMultiBandColorRenderer(rlayer.dataProvider(), 3, 2, 1)

# Apply contrast enhancement per band
for band_no in [3, 2, 1]:
    ce = QgsContrastEnhancement(rlayer.dataProvider().dataType(band_no))
    ce.setContrastEnhancementAlgorithm(QgsContrastEnhancement.StretchToMinimumMaximum)
    stats = rlayer.dataProvider().bandStatistics(band_no)
    ce.setMinimumValue(stats.minimumValue)
    ce.setMaximumValue(stats.maximumValue)
    if band_no == 3:
        renderer.setRedContrastEnhancement(ce)
    elif band_no == 2:
        renderer.setGreenContrastEnhancement(ce)
    else:
        renderer.setBlueContrastEnhancement(ce)

rlayer.setRenderer(renderer)
rlayer.triggerRepaint()

---

Example 6: Land Use Classification Visualization

from qgis.core import (QgsRasterLayer, QgsPalettedRasterRenderer, QgsProject)
from qgis.PyQt.QtGui import QColor

rlayer = QgsRasterLayer("/data/landuse.tif", "Land Use")
if not rlayer.isValid():
    raise ValueError("Layer failed to load")

QgsProject.instance().addMapLayer(rlayer)

classes = [
    QgsPalettedRasterRenderer.Class(1, QColor(255, 0, 0), "Urban"),
    QgsPalettedRasterRenderer.Class(2, QColor(0, 128, 0), "Forest"),
    QgsPalettedRasterRenderer.Class(3, QColor(0, 0, 255), "Water"),
    QgsPalettedRasterRenderer.Class(4, QColor(255, 255, 0), "Agriculture"),
    QgsPalettedRasterRenderer.Class(5, QColor(128, 128, 128), "Bare Soil"),
]

renderer = QgsPalettedRasterRenderer(rlayer.dataProvider(), 1, classes)
rlayer.setRenderer(renderer)
rlayer.triggerRepaint()

---

Example 7: Raster Reprojection and Compression

import processing
from qgis.core import QgsRasterLayer, QgsProject

# Reproject from WGS84 to Dutch RD New with compression
result = processing.run("gdal:warpreproject", {
    'INPUT': '/data/dem_wgs84.tif',
    'SOURCE_CRS': 'EPSG:4326',
    'TARGET_CRS': 'EPSG:28992',
    'RESAMPLING': 1,               # Bilinear (good for continuous data like DEMs)
    'NODATA': -9999,
    'TARGET_RESOLUTION': 25,
    'OPTIONS': 'COMPRESS=DEFLATE|PREDICTOR=2|ZLEVEL=9',
    'OUTPUT': '/output/dem_rd.tif'
})

reprojected = QgsRasterLayer(result['OUTPUT'], 'DEM_RD')
QgsProject.instance().addMapLayer(reprojected)

---

Example 8: Raster to Vector and Back

import processing

# Polygonize a classified raster
poly_result = processing.run("gdal:polygonize", {
    'INPUT': '/data/classified.tif',
    'BAND': 1,
    'FIELD': 'class_id',
    'EIGHT_CONNECTEDNESS': False,
    'OUTPUT': '/output/classes.gpkg'
})

# Rasterize a vector layer
raster_result = processing.run("gdal:rasterize", {
    'INPUT': '/data/buildings.gpkg',
    'FIELD': 'height',
    'UNITS': 1,          # Georeferenced units
    'WIDTH': 0.5,        # 0.5m pixel width
    'HEIGHT': 0.5,       # 0.5m pixel height
    'NODATA': 0,
    'OUTPUT': '/output/building_heights.tif'
})

---

Example 9: Pixel-Level Analysis with QgsRasterBlock

from qgis.core import QgsRasterLayer, QgsRasterBandStats

rlayer = QgsRasterLayer("/data/dem.tif", "DEM")
provider = rlayer.dataProvider()

# Read entire raster as block
block = provider.block(1, rlayer.extent(), rlayer.width(), rlayer.height())

# Count pixels above threshold
threshold = 500.0
above_count = 0
total_valid = 0

for row in range(block.height()):
    for col in range(block.width()):
        if not block.isNoData(row, col):
            total_valid += 1
            if block.value(row, col) > threshold:
                above_count += 1

if total_valid > 0:
    pct = (above_count / total_valid) * 100
    print(f"Pixels above {threshold}m: {above_count} ({pct:.1f}%)")

---

Example 10: Merge Raster Tiles with NoData Handling

import processing

# Build virtual raster first (lightweight check)
vrt_result = processing.run("gdal:buildvirtualraster", {
    'INPUT': ['/data/tile_01.tif', '/data/tile_02.tif', '/data/tile_03.tif'],
    'RESOLUTION': 0,       # Average resolution
    'INPUT_NODATA_VALUE': -9999,
    'OUTPUT': '/output/mosaic.vrt'
})

# Physical merge with compression
merge_result = processing.run("gdal:merge", {
    'INPUT': ['/data/tile_01.tif', '/data/tile_02.tif', '/data/tile_03.tif'],
    'NODATA_INPUT': -9999,
    'NODATA_OUTPUT': -9999,
    'OPTIONS': 'COMPRESS=DEFLATE',
    'OUTPUT': '/output/mosaic.tif'
})

API Signatures Reference (Raster Analysis)

QgsRasterLayer

class QgsRasterLayer(QgsMapLayer):
    def __init__(self, path: str, baseName: str = '', providerType: str = 'gdal') -> None
    def isValid(self) -> bool
    def width(self) -> int                    # Pixel columns
    def height(self) -> int                   # Pixel rows
    def extent(self) -> QgsRectangle
    def bandCount(self) -> int
    def bandName(self, bandNo: int) -> str    # 1-indexed
    def rasterType(self) -> int               # 0=GrayOrUndefined, 1=Palette, 2=Multiband
    def crs(self) -> QgsCoordinateReferenceSystem
    def dataProvider(self) -> QgsRasterDataProvider
    def renderer(self) -> QgsRasterRenderer
    def setRenderer(self, renderer: QgsRasterRenderer) -> None
    def triggerRepaint(self) -> None
    def setContrastEnhancement(self, algorithm: QgsContrastEnhancement.ContrastEnhancementAlgorithm,
                                limits: QgsRasterMinMaxOrigin.Limits = ...,
                                extent: QgsRectangle = ...,
                                sampleSize: int = ...,
                                generateLookupTableFlag: bool = True) -> None

---

QgsRasterDataProvider

class QgsRasterDataProvider(QgsDataProvider, QgsRasterInterface):
    def sample(self, point: QgsPointXY, band: int) -> Tuple[float, bool]
    def identify(self, point: QgsPointXY, format: QgsRaster.IdentifyFormat) -> QgsRasterIdentifyResult
    def block(self, bandNo: int, extent: QgsRectangle, width: int, height: int) -> QgsRasterBlock
    def bandStatistics(self, bandNo: int, stats: int = QgsRasterBandStats.All,
                       extent: QgsRectangle = ..., sampleSize: int = 0) -> QgsRasterBandStats
    def histogram(self, bandNo: int, binCount: int = 0) -> QgsRasterHistogram
    def dataType(self, bandNo: int) -> Qgis.DataType
    def sourceNoDataValue(self, bandNo: int) -> float
    def sourceHasNoDataValue(self, bandNo: int) -> bool
    def setEditable(self, enabled: bool) -> bool
    def writeBlock(self, block: QgsRasterBlock, band: int, xOffset: int = 0, yOffset: int = 0) -> bool

---

QgsRasterBlock

class QgsRasterBlock:
    def __init__(self, dataType: Qgis.DataType = ..., width: int = 0, height: int = 0) -> None
    def value(self, row: int, column: int) -> float
    def setValue(self, row: int, column: int, value: float) -> bool
    def isNoData(self, row: int, column: int) -> bool
    def setIsNoData(self, row: int, column: int) -> bool
    def setNoDataValue(self, noDataValue: float) -> None
    def width(self) -> int
    def height(self) -> int
    def setData(self, data: bytes) -> None

---

QgsRasterCalculator

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

class QgsRasterCalculatorEntry:
    ref: str                     # Reference name used in expression, e.g. 'dem@1'
    raster: QgsRasterLayer       # Source raster layer
    bandNumber: int              # Band number (1-indexed)

class QgsRasterCalculator:
    def __init__(self,
                 formulaString: str,           # Expression with refs in double quotes
                 outputFile: str,              # Output file path
                 outputFormat: str,            # e.g. 'GTiff'
                 outputExtent: QgsRectangle,   # Output extent
                 nOutputColumns: int,          # Output width in pixels
                 nOutputRows: int,             # Output height in pixels
                 rasterEntries: List[QgsRasterCalculatorEntry]) -> None

    def processCalculation(self, feedback: QgsRasterCalcFeedback = None) -> int
    # Returns: 0=Success, 1=CreateOutputError, 2=InputLayerError,
    #          3=Canceled, 4=ParserError, 5=MemoryError, 6=BandError

---

QgsRasterBandStats

class QgsRasterBandStats:
    All = 63                    # Flag to compute all statistics
    minimumValue: float
    maximumValue: float
    mean: float
    stdDev: float
    sum: float
    range: float
    elementCount: int
    sumOfSquares: float

---

QgsRasterIdentifyResult

class QgsRasterIdentifyResult:
    def isValid(self) -> bool
    def results(self) -> Dict[int, Any]    # {bandNo: value}
    def error(self) -> QgsError

---

Raster Renderers

QgsSingleBandGrayRenderer

class QgsSingleBandGrayRenderer(QgsRasterRenderer):
    def __init__(self, input: QgsRasterInterface, grayBand: int) -> None
    def setContrastEnhancement(self, ce: QgsContrastEnhancement) -> None
    def setGradient(self, gradient: int) -> None  # 0=BlackToWhite, 1=WhiteToBlack

QgsSingleBandPseudoColorRenderer

class QgsSingleBandPseudoColorRenderer(QgsRasterRenderer):
    def __init__(self, input: QgsRasterInterface, band: int,
                 shader: QgsRasterShader) -> None
    def shader(self) -> QgsRasterShader
    def setBand(self, bandNo: int) -> None

QgsMultiBandColorRenderer

class QgsMultiBandColorRenderer(QgsRasterRenderer):
    def __init__(self, input: QgsRasterInterface,
                 redBand: int, greenBand: int, blueBand: int) -> None
    def setRedBand(self, band: int) -> None
    def setGreenBand(self, band: int) -> None
    def setBlueBand(self, band: int) -> None
    def setRedContrastEnhancement(self, ce: QgsContrastEnhancement) -> None
    def setGreenContrastEnhancement(self, ce: QgsContrastEnhancement) -> None
    def setBlueContrastEnhancement(self, ce: QgsContrastEnhancement) -> None

QgsPalettedRasterRenderer

class QgsPalettedRasterRenderer(QgsRasterRenderer):
    def __init__(self, input: QgsRasterInterface, bandNumber: int,
                 classes: List[QgsPalettedRasterRenderer.Class]) -> None

    class Class:
        def __init__(self, value: float, color: QColor, label: str = '') -> None

QgsHillshadeRenderer

class QgsHillshadeRenderer(QgsRasterRenderer):
    def __init__(self, input: QgsRasterInterface, band: int,
                 lightAzimuth: float, lightAngle: float) -> None
    def setAzimuth(self, azimuth: float) -> None
    def setAltitude(self, altitude: float) -> None
    def setZFactor(self, factor: float) -> None
    def setMultiDirectional(self, isMultiDirectional: bool) -> None

---

QgsRasterShader / QgsColorRampShader

class QgsRasterShader:
    def setRasterShaderFunction(self, function: QgsRasterShaderFunction) -> None

class QgsColorRampShader(QgsRasterShaderFunction):
    Interpolated = 0
    Discrete = 1
    Exact = 2

    def setColorRampType(self, type: int) -> None
    def setColorRampItemList(self, items: List[QgsColorRampShader.ColorRampItem]) -> None

    class ColorRampItem:
        def __init__(self, value: float, color: QColor, label: str = '') -> None

---

QgsContrastEnhancement

class QgsContrastEnhancement:
    NoEnhancement = 0
    StretchToMinimumMaximum = 1
    StretchAndClipToMinimumMaximum = 2
    ClipToMinimumMaximum = 3

    def __init__(self, dataType: Qgis.DataType) -> None
    def setContrastEnhancementAlgorithm(self, algorithm: int) -> None
    def setMinimumValue(self, value: float) -> None
    def setMaximumValue(self, value: float) -> None

---

Processing Algorithm Parameter Reference

gdal:warpreproject

gdal:hillshade

gdal:slope

gdal:contour

gdal:polygonize

gdal:rasterize

gdal:translate