Use the Raster Calculator in PyQGIS

Band maths is where raster analysis starts: an index from two bands, a mask from a threshold, a difference between two dates. The calculation itself is arithmetic; what makes it go wrong is everything around it — nodata treated as a number, two rasters that do not line up pixel for pixel, and an output type that quietly truncates a ratio to zero.

This recipe belongs to Raster Analysis Workflows in PyQGIS. It covers the two calculator APIs, expression syntax, protecting nodata, aligning misaligned inputs, and choosing an output data type that can hold the answer.

How a band-maths expression is evaluatedTwo input bands are read pixel by pixel. The expression combines them into a result. A nodata mask taken from the inputs is applied so that any pixel missing in either input becomes nodata in the output rather than an arithmetic result computed from a fill value.A missing pixel must stay missing, not become −9999NIR bandFloat32, nodata −9999red bandFloat32, nodata −9999expression(nir − red) / (nir + red)evaluated per pixelNDVI, Float32values from −1 to 1nodata carried throughgaps stay gapsAn Int16 output here would round every NDVI value to −1, 0 or 1

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • One or more rasters that share a CRS. If they do not, reproject first — see Batch Reprojecting Raster Datasets.
  • Knowledge of each input's nodata value, which layer.dataProvider().sourceNoDataValue(1) reports.

The Processing route

For most work native:rastercalc is the right call: it takes a plain expression, handles the block iteration and writes the output.

import processing

processing.run("native:rastercalc", {
    "LAYERS": ["/data/sentinel.tif"],
    "EXPRESSION": '("sentinel@8" - "sentinel@4") / ("sentinel@8" + "sentinel@4")',
    "EXTENT": None,
    "CELL_SIZE": None,
    "CRS": None,
    "OUTPUT": "/data/output/ndvi.tif",
})

Breakdown: Band references take the form "layername@band", where the layer name is the file's base name for a path input, and the band number is one-based. Leaving EXTENT, CELL_SIZE and CRS as None inherits them from the first input, which is what you want when every input already matches. The output is Float32 by default — appropriate for a ratio, and the reason an index calculated with the GDAL calculator's default Int16 comes out as three distinct values.

For a threshold mask, comparison operators return 1 or 0:

processing.run("native:rastercalc", {
    "LAYERS": ["/data/dem.tif"],
    "EXPRESSION": '"dem@1" > 250',
    "OUTPUT": "/data/output/above_250.tif",
})

Breakdown: A boolean result is a raster of ones and zeros — useful directly as a mask for clipping another raster, and convertible to polygons with gdal:polygonize when the areas are what you need. Because the values are small integers, this is one of the few cases where forcing an Int16 or Byte output is worth doing to save space.

Protect nodata explicitly

The calculator propagates nodata from its inputs, but only where the input's nodata value is actually declared. A raster whose gaps are stored as -9999 with no nodata flag set will happily contribute -9999 to the arithmetic.

from qgis.core import QgsRasterLayer

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

if not provider.sourceHasNoDataValue(1):
    provider.setNoDataValue(1, -9999)
    print("declared nodata for band 1")

Breakdown: sourceHasNoDataValue() asks whether the file itself declares a nodata value; setNoDataValue() sets it on the layer for this session, which is enough for the calculator to honour it. Making this check part of the loading step is the cheapest way to avoid a class of results that look plausible and are wrong — a mean elevation of −4 000 metres is obvious, but a five percent bias from a handful of fill pixels is not.

Where a condition must exclude values without relying on the nodata flag, write it into the expression:

expression = '("dem@1" > 250) * ("dem@1" != -9999)'

Breakdown: Multiplying two boolean sub-expressions is the raster equivalent of an AND, and it keeps the arithmetic explicit. This idiom is worth knowing because the calculator has no if — every conditional is expressed as multiplication by a boolean mask.

Align inputs before combining them

Two rasters from different sources rarely share an origin and pixel size. The calculator resamples to the first layer's grid, which silently changes values by nearest-neighbour sampling.

aligned = processing.run("gdal:warpreproject", {
    "INPUT": "/data/rainfall.tif",
    "TARGET_CRS": "EPSG:27700",
    "RESAMPLING": 1,                  # bilinear, appropriate for continuous data
    "TARGET_RESOLUTION": 25,
    "TARGET_EXTENT": "432000,436000,187000,191000 [EPSG:27700]",
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

Breakdown: Naming the target extent and resolution explicitly — rather than letting the calculator improvise — makes the alignment a visible, reviewable step. RESAMPLING: 1 is bilinear, correct for continuous surfaces such as elevation or rainfall; use 0 (nearest neighbour) for categorical rasters such as land cover, where averaging class codes produces classes that do not exist.

Why misaligned grids must be warped firstTwo grids are overlaid. The finer grid's origin is offset from the coarser one, so a single coarse cell overlaps four fine cells and no cell centre coincides. After warping to a common origin and cell size, the two grids share every cell boundary and a per-pixel expression is meaningful.A per-pixel expression assumes the pixels are the same pixelsbefore: offset gridsno cell boundary is sharedafter: warped to one gridevery boundary coincides

The direct API, when you need control

QgsRasterCalculator is the class behind the algorithm, and it is worth reaching for when the output geometry has to be set precisely.

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

entry = QgsRasterCalculatorEntry()
entry.ref = "dem@1"
entry.raster = layer
entry.bandNumber = 1

calculator = QgsRasterCalculator(
    '"dem@1" * 0.3048',                    # feet to metres
    "/data/output/dem_metres.tif",
    "GTiff",
    layer.extent(),
    layer.width(),
    layer.height(),
    [entry],
    QgsProject.instance().transformContext(),
)

result = calculator.processCalculation()
if result != QgsRasterCalculator.Success:
    raise RuntimeError(f"calculation failed with code {result}")

Breakdown: Each QgsRasterCalculatorEntry binds a reference name used in the expression to a layer and band, so the names in the expression are yours rather than derived from a filename. Passing the source layer's extent, width and height guarantees the output is pixel-identical to the input — the property that makes results stackable. processCalculation() returns a status code rather than raising, and Success is zero, so the truthiness test people reach for first is inverted.

Avoid writing intermediates you do not need

A chain of three expressions written as three GeoTIFFs reads every pixel three times and leaves two files nobody wanted. A virtual raster describes the calculation instead of materialising it.

import processing

ndvi = processing.run("native:virtualrastercalc", {
    "LAYERS": ["/data/sentinel.tif"],
    "EXPRESSION": '("sentinel@8" - "sentinel@4") / ("sentinel@8" + "sentinel@4")',
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

processing.run("native:rastercalc", {
    "LAYERS": [ndvi],
    "EXPRESSION": '"ndvi@1" > 0.4',
    "OUTPUT": "/data/output/vegetated.tif",
})

Breakdown: native:virtualrastercalc produces a layer that computes its pixels on demand rather than a file on disk, so the intermediate NDVI exists only as a definition. The final step is a real write, and only that pass touches the disk. On a large scene this is the difference between three full read-write cycles and one, and it removes the temporary files that otherwise accumulate in a scratch directory. The trade is that every consumer of the virtual layer recomputes it, so materialise anything read many times.

Where the calculation is a single pass, prefer combining the expressions instead — the calculator is perfectly happy with a compound expression, and one pass beats two:

expression = '(("sentinel@8" - "sentinel@4") / ("sentinel@8" + "sentinel@4")) > 0.4'

Breakdown: Writing the threshold and the index in one statement gives the same result in a single read of the source, with no intermediate at all. Readability is the only cost, and a comment naming what the expression computes usually settles that. As a rule: combine when the expression stays legible, use a virtual raster when it does not, and write a real file only when something downstream needs one.

Two temporary rasters, or noneWriting each stage of a calculation to disk reads and writes every pixel three times and leaves two intermediate files. A virtual raster keeps the intermediate as a definition, so the pixels are read once and only the final result is written.Only the last step has to touch the diskwrittenndvi.tif on diskmasked.tif on diskvegetated.tif — the deliverablethree full read-write cycles, two files to clean upvirtuala definition, computed on demand — nothing on diskvegetated.tif — the deliverableone read of the source, one file written

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.22 LTR3.9qgis:rastercalculator is the algorithm id; native:rastercalc does not exist yet.
3.28 LTR3.9native:rastercalc available and preferred.
3.34 LTR3.12Baseline for this page.
3.40 / 3.443.12Adds native:virtualrastercalc for calculations evaluated lazily without writing a file.

Troubleshooting

  • The output is all zeros or all ones. An integer output type truncated a fractional result. Force Float32.
  • Gaps became large negative numbers. Nodata was not declared on the input. Set it before calculating.
  • "Parser error" on the expression. Band references must be quoted exactly as "name@band", and the layer name must match what QGIS derived from the path.
  • The output extent is not what you expected. It follows the first layer unless EXTENT is given. Set it explicitly when inputs differ.
  • Values shifted slightly compared with a desktop run. The inputs were resampled during alignment. Warp them deliberately with a documented resampling method instead.
  • The calculation is very slow on a large raster. Each expression pass reads every pixel. Combine steps into one expression rather than chaining several passes.

Conclusion

Raster maths in PyQGIS is native:rastercalc with a quoted "layer@band" expression, a Float32 output for anything fractional, and nodata declared on every input before the first calculation. Align inputs deliberately when they come from different sources, and drop to QgsRasterCalculator when the output grid has to match an existing raster exactly.

Frequently Asked Questions

How do I write a conditional? Multiply by a boolean sub-expression: ("dem@1" > 250) * "rainfall@1" yields the rainfall where elevation exceeds 250 and zero elsewhere. There is no if in the calculator's syntax.

Can I use bands from different files in one expression? Yes — list every file in LAYERS and reference each as "basename@band". They must share a grid, or the first layer's grid will be imposed on the rest.

Which output format should I use? GeoTIFF with compression for anything shared, and TEMPORARY_OUTPUT for intermediates. Add COMPRESS=DEFLATE through the creation options when size matters.

Why is my NDVI outside the −1 to 1 range? Either nodata is contributing to the arithmetic, or the bands are the wrong way round. Check both against a known vegetated pixel.

Is the GDAL calculator different?gdal:rastercalculator uses NumPy syntax with letters for inputs and defaults to a narrower output type. It is faster on very large files; the native calculator is easier to read and safer with nodata.