Use Rasterio and NumPy with QGIS Rasters

QGIS's raster calculator and Processing algorithms cover a great deal, but some raster work is simply clearer as array arithmetic: a vegetation index with conditional masking, a reclassification driven by a lookup table, a moving-window statistic, a model trained in scikit-learn applied pixel by pixel. Rasterio gives a Pythonic interface to GDAL — windows, masks, profiles — and NumPy does the arithmetic. QGIS remains the place to load, inspect and style the inputs and outputs.

This recipe belongs to PyQGIS and the Python Data Stack. It opens the file behind a QGIS raster layer with rasterio, reads masked arrays and windows, computes an index and a reclassification, writes the result with a correct profile, and loads and styles it in QGIS.

QGIS layer to array and backA QGIS raster layer points to a file. Rasterio opens that file and reads bands as masked NumPy arrays, respecting NoData. NumPy computes the result. Rasterio writes it with a profile copied from the input, updated for the output data type and NoData value. QGIS loads the new file as a layer and applies a style.File in, arrays through, file outQGIS layersource() → pathrasterio readmasked arraysNumPyindex, reclasswrite + loadprofile, NoDatageoreferencing travels in the profile — copy it, then change only what differs

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series, with rasterio installed into QGIS's Python. On conda-based and Linux installations rasterio shares QGIS's GDAL; on Windows, install it with the OSGeo4W shell or pip as described in installing Python packages into QGIS.
  • A raster layer backed by a file — GeoTIFF, VRT, COG. Rasters from WMS or other services have no file for rasterio to open.

Open the file behind a layer

A QGIS raster layer's source is usually a file path, possibly with provider options appended. Rasterio needs just the path.

import rasterio
from qgis.core import QgsProject, QgsProviderRegistry

layer = QgsProject.instance().mapLayersByName("sentinel2_2026_06")[0]
parts = QgsProviderRegistry.instance().decodeUri("gdal", layer.source())
path = parts["path"]

with rasterio.open(path) as src:
    print(src.count, "bands", src.width, "x", src.height, src.dtypes[0])
    print(src.crs, src.res, src.nodata)
    print(src.descriptions)          # band names, if the file has them

Breakdown: decodeUri with the gdal provider key separates the file path from any layer options, which is safer than using layer.source() directly. The dataset object reports band count, size, data type, CRS, pixel size and NoData — compare them with QGIS's Layer Properties to confirm both see the same file. Using with closes the file when done, which matters on Windows where an open handle blocks QGIS from overwriting or deleting the file.

Read masked arrays

NoData pixels must not take part in arithmetic. Reading with masked=True returns NumPy masked arrays in which NoData is masked out automatically.

import numpy as np

with rasterio.open(path) as src:
    red = src.read(4, masked=True).astype("float32")      # Sentinel-2 B04
    nir = src.read(8, masked=True).astype("float32")      # Sentinel-2 B08
    profile = src.profile

print(red.mask.sum(), "masked pixels;", float(red.mean()), "mean red")

Breakdown: Masked arrays carry their mask through arithmetic, so a NoData pixel in either band stays masked in the result. Converting to float before computing an index avoids integer overflow and integer division on reflectance stored as uint16. Keeping the source profile — CRS, transform, size, dtype, compression — is how the output inherits the input's georeferencing. Band numbers in rasterio start at 1, like in QGIS and GDAL.

Compute an index

A normalised difference index is the canonical example: two bands, one formula, a few lines, with masking and division by zero handled.

NDVI from two bandsThe NIR and red bands are read as masked float arrays. NDVI is computed as NIR minus red divided by NIR plus red. Pixels where the denominator is zero are masked to avoid division by zero. The result ranges from minus one to one: water and bare ground near or below zero, sparse vegetation around 0.2 to 0.4, dense vegetation above 0.6.(NIR − red) / (NIR + red)NIR bandB08 · maskedred bandB04 · maskedNDVImask whereNIR + red = 0≤ 0 water0.2–0.4 sparse> 0.6 dense−1+1

denominator = nir + red
ndvi = np.ma.where(denominator == 0, np.ma.masked, (nir - red) / denominator)
ndvi = np.ma.clip(ndvi, -1, 1)
print("NDVI range", float(ndvi.min()), float(ndvi.max()),
      "| vegetated share", float((ndvi > 0.4).sum() / ndvi.count()))

Breakdown: Masking where the denominator is zero avoids division warnings and infinite values over pixels that are dark in both bands. Clipping guards against values just outside the valid range caused by sensor noise or scaling. ndvi.count() counts unmasked pixels only, so the vegetated share is a proportion of valid land, not of the whole rectangle. The same expression in QGIS's raster calculator would work too, as shown in the raster calculator recipe, but would not mask zero denominators as cleanly.

Reclassify with a lookup table

Reclassification maps ranges or codes to new values. With NumPy, digitize handles ranges and a lookup array handles codes, both in one vectorised step.

breaks = np.array([0.0, 0.2, 0.4, 0.6])           # class edges
classes = np.ma.masked_array(np.digitize(ndvi.filled(-2), breaks), mask=ndvi.mask)
# 0: < 0 (water/bare), 1: 0–0.2, 2: 0.2–0.4, 3: 0.4–0.6, 4: > 0.6
counts = {int(c): int((classes == c).sum()) for c in range(5)}
print(counts)

with rasterio.open("/data/landcover/corine_2018.tif") as lc:
    codes = lc.read(1)
lookup = np.zeros(1000, dtype="uint8")             # code → simplified class
lookup[[111, 112]] = 1                             # urban
lookup[[311, 312, 313]] = 2                        # forest
simplified = lookup[codes]

Breakdown: np.digitize returns, for each value, the index of the bin it falls in; filling masked pixels with an out-of-range value first and re-applying the mask afterwards keeps NoData separate. For categorical rasters, indexing a lookup array with the code array reclassifies every pixel in one step, far faster than a loop or chained where calls; the array just needs to be at least as long as the largest code. QGIS's own approach is covered in reclassifying raster values.

Write the result correctly

The output needs the input's georeferencing but its own data type, band count and NoData value. Copying the profile and updating only what differs is the reliable pattern.

Copy the profile, change what differsThe input profile holds driver, width, height, CRS, transform, data type, band count, NoData and compression. For the NDVI output, width, height, CRS and transform are kept; data type changes to float32, count to one, NoData to minus 9999, and compression and tiling are set for a cloud-friendly GeoTIFF.Keep the grid, change the payloadkept from inputwidth · heightcrs · transformsame grid, same placeupdated for outputdtype float32 · count 1nodata −9999compress, tiled

out_profile = profile.copy()
out_profile.update(dtype="float32", count=1, nodata=-9999.0,
                   compress="deflate", predictor=3, tiled=True,
                   blockxsize=512, blockysize=512, driver="GTiff")

out_path = "/data/results/ndvi_2026_06.tif"
with rasterio.open(out_path, "w", **out_profile) as dst:
    dst.write(ndvi.filled(-9999.0).astype("float32"), 1)
    dst.set_band_description(1, "NDVI")
    dst.update_tags(source=path, formula="(B08-B04)/(B08+B04)")

Breakdown: Keeping crs, transform, width and height from the input places the output exactly on the input grid. filled(-9999) replaces masked pixels with the declared NoData value before writing. Deflate compression with predictor 3 suits floating-point data; tiled 512-pixel blocks make the file efficient for QGIS to render at any zoom and are a step towards a Cloud Optimized GeoTIFF, covered in exporting a COG. The band description and tags record what the file contains and how it was made.

Load and style the result in QGIS

The new file is an ordinary GeoTIFF; QGIS loads it like any other and a pseudocolour renderer makes it readable.

from qgis.core import (QgsRasterLayer, QgsColorRampShader, QgsRasterShader,
                       QgsSingleBandPseudoColorRenderer, QgsStyle)

ndvi_layer = QgsRasterLayer(out_path, "NDVI June 2026")
shader_fn = QgsColorRampShader(-0.2, 0.9, QgsStyle.defaultStyle().colorRamp("RdYlGn"))
shader_fn.classifyColorRamp(9)
shader = QgsRasterShader()
shader.setRasterShaderFunction(shader_fn)
ndvi_layer.setRenderer(QgsSingleBandPseudoColorRenderer(ndvi_layer.dataProvider(), 1, shader))
QgsProject.instance().addMapLayer(ndvi_layer)

Breakdown: A fixed range of −0.2 to 0.9 rather than the layer's own minimum and maximum keeps colours comparable between dates — the same NDVI value gets the same colour in June and August. The red–yellow–green ramp is the convention for vegetation indices. For a reusable look, save the style as QML once and apply it to every new date, as in applying a colour ramp to a raster.

Process large rasters in windows

A national mosaic does not fit in memory. Rasterio's windows let you read, compute and write block by block with the same code.

with rasterio.open(path) as src, rasterio.open(out_path, "w", **out_profile) as dst:
    for _, window in src.block_windows(1):
        r = src.read(4, window=window, masked=True).astype("float32")
        n = src.read(8, window=window, masked=True).astype("float32")
        d = n + r
        v = np.ma.where(d == 0, np.ma.masked, (n - r) / d)
        dst.write(np.ma.clip(v, -1, 1).filled(-9999.0).astype("float32"), 1, window=window)

Breakdown: block_windows iterates the file's internal blocks, so each read is aligned with how the data is stored and costs one disk read. Memory use is bounded by the block size regardless of raster size. The output must use the same grid; with matching block sizes in the output profile, writes are aligned too. For per-pixel operations like this the result is identical to the whole-array version.

QGIS version compatibility

The rasterio and NumPy code is independent of QGIS version; the QGIS parts work on 3.34 LTR, 3.40 LTR and QGIS 4. Rasterio must be built against a GDAL compatible with the one QGIS uses — mixing them in one process can crash on import. On conda, install both from conda-forge into the same environment; with OSGeo4W, use its packages.

Troubleshooting

  • ImportError or a crash when importing rasterio. It was built against a different GDAL; install it from the same distribution as QGIS.
  • The output is shifted or stretched. The transform or size in the profile was changed; copy them unchanged from the input.
  • NoData shows as black. The NoData value was not set in the profile or the masked array was not filled with it.
  • QGIS cannot overwrite the output. The file is still open in a rasterio handle or loaded in the project; close or remove it first.

Conclusion

Get the file path from the layer with decodeUri, read masked float arrays, compute with NumPy while keeping masks, copy the input profile and change only data type, band count, NoData and compression, write with a description and tags, load and style with a fixed range, and use block windows for rasters that do not fit in memory.

Frequently Asked Questions

Can I skip rasterio and use only PyQGIS? Yes. QgsRasterBlock reads arrays and QgsRasterFileWriter writes them, as in writing a NumPy array to a raster. Rasterio is more concise for heavy array work.

Does rasterio read COGs over HTTP? Yes, with /vsicurl/ or an https:// path, using GDAL's range requests.

How do I apply a scikit-learn model to every pixel? Reshape bands to a 2-D array of pixels by features, predict, and reshape back to the grid.

Is xarray an alternative? Yes; rioxarray adds rasterio-backed I/O to xarray and suits time series of rasters.