Write a NumPy Array to a Raster in PyQGIS

Reading cells into NumPy is half of a custom raster workflow. The other half is getting the result back out as a file that lines up with everything else: the same cell size, the same origin, the same CRS, a sensible data type, and NoData where the calculation had nothing to say. Get the georeferencing half a cell wrong and the output looks fine on its own and shifts by a metre against the input — the kind of error that survives review.

This recipe belongs to Raster Analysis Workflows in PyQGIS. It writes arrays through QgsRasterFileWriter so no second library is needed in a plugin, matches the output grid to a source layer, handles data types and NoData, writes large results tile by tile, and sets GeoTIFF creation options.

Values are the easy partA NumPy array of computed values sits on the left. Around it, four other things the file needs: the grid, taken from the source layer as extent, width, height and CRS; the data type, such as Float32; the NoData value, such as minus 9999; and creation options such as deflate compression and tiling. All five go to QgsRasterFileWriter, which creates a provider, receives one or more blocks, and closes to produce a GeoTIFF.Five things the output file needsvalues: NumPy arraygrid: extent, size, CRSdata type: Float32NoData: -9999QgsRasterFileWriter+ creation optionsCOMPRESS=DEFLATEresult.tiflines up with source

Prerequisites

  • QGIS 3.40 LTR or newer, or the QGIS 4 series.
  • A source raster whose grid the output should match, and a NumPy array computed from it — for example with the block reading shown in reading raster pixels with QgsRasterBlock.
  • Write access to the output folder.

Write a single-band raster that matches its source

The writer creates an empty raster with a given data type, size, extent and CRS and returns a data provider for it. You then put the array into a QgsRasterBlock and write the block.

import numpy as np
from qgis.core import (
    Qgis, QgsProject, QgsRasterFileWriter, QgsRasterBlock,
)

source = QgsProject.instance().mapLayersByName("dem_2m")[0]
width, height = source.width(), source.height()

elevation = read_whole_band(source)            # float64 array with NaN for NoData
slope_proxy = np.hypot(*np.gradient(elevation, 2.0))

NODATA = -9999.0
out = np.where(np.isnan(slope_proxy), NODATA, slope_proxy).astype(np.float32)

writer = QgsRasterFileWriter("/data/work/gradient.tif")
writer.setOutputFormat("GTiff")
provider = writer.createOneBandRaster(
    Qgis.DataType.Float32, width, height, source.extent(), source.crs())
if provider is None or not provider.isValid():
    raise RuntimeError(f"could not create output: {writer.error()}")

provider.setNoDataValue(1, NODATA)
provider.setEditable(True)
block = QgsRasterBlock(Qgis.DataType.Float32, width, height)
block.setData(out.tobytes())
if not provider.writeBlock(block, 1, 0, 0):
    raise RuntimeError("writeBlock failed")
provider.setEditable(False)
del provider

Breakdown: Passing the source layer's extent(), width(), height() and crs() reproduces its grid exactly, which is the single most important step — the cell size is implied by extent divided by size, so the same four values guarantee alignment. The array is cast to float32 to match the declared data type; tobytes() then produces exactly the byte layout the block expects, row by row from the top. NoData must be replaced with a real number before writing because a NaN in an integer file is impossible and in a float file is often mishandled by other software. setEditable(False) flushes and closes the dataset, and deleting the provider releases the file handle so the output can be opened straight away.

read_whole_band stands for whichever reading approach fits the raster's size — a single block for small rasters, the tiled reader for large ones.

Pick the smallest type that holds the answerA table of five data types with value range, typical use and uncompressed size for a ten thousand by ten thousand cell raster. Byte, 0 to 255, class maps, 95 megabytes. Int16, minus 32,768 to 32,767, elevations in decimetres, 191 megabytes. UInt16, 0 to 65,535, scaled reflectance, 191 megabytes. Float32, about seven significant digits, continuous values, 381 megabytes. Float64, fifteen digits, rarely needed, 763 megabytes.Data type sets range, precision and sizetyperangetypical use10k × 10kByte0 – 255classes, masks95 MBInt16±32,767heights in decimetres191 MBUInt160 – 65,535scaled reflectance191 MBFloat32~7 significant digitscontinuous measurements381 MBFloat64~15 significant digitsrarely justified763 MB

Data types and NoData choices

The output type should be the smallest that holds every value the calculation can produce, with room for a NoData value outside that range. Too small and values wrap or clip silently; too large and files are two to four times bigger than they need to be.

classes = np.digitize(slope_proxy, bins=[2, 5, 15, 30]).astype(np.uint8)
classes[np.isnan(slope_proxy)] = 255

writer = QgsRasterFileWriter("/data/work/slope_classes.tif")
provider = writer.createOneBandRaster(
    Qgis.DataType.Byte, width, height, source.extent(), source.crs())
provider.setNoDataValue(1, 255)
provider.setEditable(True)
block = QgsRasterBlock(Qgis.DataType.Byte, width, height)
block.setData(classes.tobytes())
provider.writeBlock(block, 1, 0, 0)
provider.setEditable(False)
del provider

Breakdown: A classified output with five classes needs one byte per cell, a quarter of Float32, and it compresses far better. Using 255 as NoData keeps 0 available as a real class. np.digitize returns int64 by default, so the explicit astype(np.uint8) is what makes the bytes match the declared type — writing int64 bytes into a Byte block produces a striped, eight-times-too-wide mess. For continuous results Float32 is almost always enough: its seven significant digits exceed the precision of any elevation or reflectance measurement.

Write large outputs in tiles

When the result is too large to hold as one array, write it the way you read it: tile by tile. The provider accepts blocks at any row and column offset, so each tile is computed, written and discarded.

Read, compute, write, repeatA cycle of three steps repeated for every tile: read a block from the source at row0 and col0, compute the result array for that tile, and write the result block to the output provider at the same row0 and col0. A grid on the right shows tiles already written in green, the current tile in orange, and remaining tiles in grey.One tile in memory, however big the outputread tilerow0, col0computeNumPy on the tilewriteBlocksame row0, col0next tilewritten · current · to do

TILE = 2048
writer = QgsRasterFileWriter("/data/work/gradient_national.tif")
writer.setCreateOptions(["COMPRESS=DEFLATE", "PREDICTOR=3", "TILED=YES",
                         "BLOCKXSIZE=512", "BLOCKYSIZE=512", "BIGTIFF=IF_SAFER"])
out_provider = writer.createOneBandRaster(
    Qgis.DataType.Float32, source.width(), source.height(),
    source.extent(), source.crs())
out_provider.setNoDataValue(1, NODATA)
out_provider.setEditable(True)

for row0, col0, src_block in iter_tiles(source, tile=TILE):
    arr = block_to_array(src_block)
    result = np.hypot(*np.gradient(arr, 2.0))
    result = np.where(np.isnan(result), NODATA, result).astype(np.float32)
    h, w = result.shape
    blk = QgsRasterBlock(Qgis.DataType.Float32, w, h)
    blk.setData(result.tobytes())
    out_provider.writeBlock(blk, 1, col0, row0)

out_provider.setEditable(False)
del out_provider

Breakdown: writeBlock(block, band, xOffset, yOffset) takes the column offset before the row offset — the reverse of the row0, col0 order the reader yields, and an easy swap to make. iter_tiles and block_to_array are the helpers from the reading guide. The creation options produce an internally tiled, deflate-compressed GeoTIFF; PREDICTOR=3 is the floating-point predictor that makes continuous surfaces compress well, and BIGTIFF=IF_SAFER avoids the 4 GB classic TIFF limit on large outputs. One caveat is specific to this example: a gradient computed per tile has edge effects at tile boundaries, because each tile cannot see its neighbours' cells. Read tiles with a one-cell overlap and trim the result when a calculation depends on neighbours.

Verify that the output lines up

A written raster that opens and looks right can still be misaligned by a fraction of a cell, and that is invisible at normal zoom. Three checks catch it in a second, and they are worth running at the end of every script that writes rasters — especially one that other people's analysis will depend on.

from qgis.core import QgsRasterLayer

def assert_same_grid(a, b, tolerance=1e-6):
    problems = []
    if a.crs() != b.crs():
        problems.append(f"CRS {a.crs().authid()} != {b.crs().authid()}")
    if (a.width(), a.height()) != (b.width(), b.height()):
        problems.append(f"size {a.width()}x{a.height()} != {b.width()}x{b.height()}")
    for name in ("xMinimum", "yMinimum", "xMaximum", "yMaximum"):
        va, vb = getattr(a.extent(), name)(), getattr(b.extent(), name)()
        if abs(va - vb) > tolerance:
            problems.append(f"{name} {va} != {vb}")
    if problems:
        raise AssertionError("; ".join(problems))

check = QgsRasterLayer("/data/work/gradient.tif", "check")
assert_same_grid(source, check)

probe = source.extent().center()
src_val, _ = source.dataProvider().sample(probe, 1)
out_val, _ = check.dataProvider().sample(probe, 1)
print(f"centre cell: source {src_val:.2f}, output {out_val:.4f}")
print("output NoData declared:", check.dataProvider().sourceHasNoDataValue(1))

Breakdown: Comparing CRS, size and all four extent edges proves the two grids are identical, because cell size and origin follow from those values. The tolerance absorbs floating-point noise in coordinates without letting a half-cell shift through. Sampling the same point in both layers is a sanity check on the values rather than the geometry — for a derived product like a gradient the numbers differ, but a NoData output where the source has data, or a value wildly out of range, shows up immediately. Confirming that NoData is declared in the written file guards against the most common downstream complaint, a black border where there should be transparency. Wrapped in a test, the same function belongs in a plugin's test suite, as in testing Processing algorithms with pytest-qgis.

Add the result to the project with a style

A written raster is just a file until someone opens it. Loading it and applying a renderer immediately makes the output reviewable in the same run.

from qgis.core import QgsRasterLayer, QgsSingleBandPseudoColorRenderer, QgsStyle

result_layer = QgsRasterLayer("/data/work/gradient.tif", "gradient")
stats = result_layer.dataProvider().bandStatistics(1)
ramp = QgsStyle.defaultStyle().colorRamp("Viridis")

renderer = QgsSingleBandPseudoColorRenderer(result_layer.dataProvider(), 1)
renderer.setClassificationMin(stats.minimumValue)
renderer.setClassificationMax(stats.maximumValue)
renderer.createShader(ramp)
result_layer.setRenderer(renderer)
QgsProject.instance().addMapLayer(result_layer)

Breakdown: Band statistics come from the file and exclude NoData, so the ramp spans only real values — a check in itself, since a minimum of -9999 would mean NoData was not declared. createShader builds the colour ramp shader between the classification bounds. Styling options beyond this are covered in applying a colour ramp to a raster.

QGIS version compatibility

QgsRasterFileWriter.createOneBandRaster and createMultiBandRaster have been available since QGIS 3.0, and writeBlock with offsets since 2.x. setCreateOptions accepts GDAL creation options as a list of strings on all 3.x releases; on the QGIS 4 series check the method name against the API documentation for your build, as some option setters have gained typed variants. Qgis.DataType is required in its scoped form on QGIS 4.

Troubleshooting

  • The output is shifted against the input. The extent or size did not come from the source layer, so the implied cell size differs.
  • Stripes or noise instead of values. The NumPy dtype does not match the block's data type.
  • The file is locked or empty after the script. The provider was not closed; call setEditable(False) and delete it.
  • NoData shows as black. setNoDataValue was not called, so the value is treated as data.
  • Seams at tile edges. A neighbourhood calculation was run per tile without overlap.

Conclusion

Create the output from the source layer's extent, size and CRS so the grids match, pick the smallest data type with room for NoData, cast arrays to exactly that type before tobytes(), and close the provider when done. For large rasters write tiles at their offsets with compression options set, and add overlap for neighbourhood calculations.

Frequently Asked Questions

How do I write several bands? Use createMultiBandRaster with a band count, then call writeBlock once per band number.

Can I write to formats other than GeoTIFF? Any GDAL driver that supports creation works via setOutputFormat; GeoTIFF is the safest default.

Should I use GDAL's Python bindings instead? They work equally well. The QGIS writer avoids a second API in a plugin and respects QGIS's CRS objects directly.

Why is my output so much larger than the input? Probably Float64 without compression. Use Float32 and COMPRESS=DEFLATE with a predictor.