Read Raster Pixels with QgsRasterBlock in PyQGIS

Processing algorithms cover most raster work — statistics, reclassification, the raster calculator — and they should be your first choice. But sooner or later a job needs the numbers themselves: a custom index no algorithm computes, a check that every cell in a delivered DEM falls in a sane range, a histogram with bins nobody else uses, a value lookup inside a plugin that has to respond instantly. For that, PyQGIS reads cells through the layer's data provider as a QgsRasterBlock.

This recipe belongs to Raster Analysis Workflows in PyQGIS. It reads blocks for an extent, handles NoData properly, turns blocks into NumPy arrays, and walks a large raster in tiles so memory use stays flat.

A block is a resampled window, not a file readThree inputs feed provider.block: the band number, an extent in the layer's CRS, and a width and height in cells. The provider reads the source cells that fall in the extent and returns a block of exactly the requested size. If width and height match the native resolution of the extent, the values are the original cells. If they are smaller, the provider resamples. The block records its data type and which cells are NoData.You ask for an extent at a size — the provider fills itband1extentQgsRectanglewidth × heightcells in the blockprovider.block()reads + resamplesQgsRasterBlockdataType(), width(), height()value(row, col), isNoData()grey cells: NoData

Prerequisites

  • QGIS 3.40 LTR or newer, or the QGIS 4 series.
  • A raster layer loaded through the gdal provider. The same code works for WMS and WCS providers, but reading cells from a rendered web map rarely gives meaningful numbers.
  • NumPy, which ships with QGIS installers on every platform.

Read a block at native resolution

A block request names a band, an extent and a size in cells. To get the original values rather than a resampled version, compute the size from the extent and the raster's cell size, and align the extent to the cell grid.

import math
from qgis.core import QgsProject, QgsRectangle

dem = QgsProject.instance().mapLayersByName("dem_2m")[0]
provider = dem.dataProvider()

cell_x = dem.rasterUnitsPerPixelX()
cell_y = dem.rasterUnitsPerPixelY()
full = dem.extent()

def aligned_window(xmin, ymin, xmax, ymax):
    col0 = math.floor((xmin - full.xMinimum()) / cell_x)
    col1 = math.ceil((xmax - full.xMinimum()) / cell_x)
    row0 = math.floor((full.yMaximum() - ymax) / cell_y)
    row1 = math.ceil((full.yMaximum() - ymin) / cell_y)
    rect = QgsRectangle(
        full.xMinimum() + col0 * cell_x, full.yMaximum() - row1 * cell_y,
        full.xMinimum() + col1 * cell_x, full.yMaximum() - row0 * cell_y,
    )
    return rect, col1 - col0, row1 - row0

window, width, height = aligned_window(429500, 433000, 430000, 433400)
block = provider.block(1, window, width, height)

print(block.width(), block.height(), block.dataType())
print("top-left cell:", block.value(0, 0), "nodata?", block.isNoData(0, 0))

Breakdown: Rows count down from the top of the raster and columns across from the left, which is why the row numbers are computed from yMaximum. Snapping the window outward to whole cells means every value in the block is a source cell, not an interpolated mix of neighbours; ask for a 250 × 200 block over an extent that is really 251 × 200 cells and the provider silently resamples. value(row, col) returns a float regardless of the underlying type, which is convenient for single lookups and slow for millions of them — the NumPy route below is the one to use for whole blocks.

Snap the window to the gridLeft: a requested window drawn in orange cuts through cells of the raster grid; the provider returns values resampled from partly covered cells. Right: the same area snapped outward to whole cell boundaries, with width and height equal to the cell counts, returns exactly the source values.Half a cell off is a different answerunaligned: resamplededges cut through cellsaligned: source cells5 × 3 cells, requested as 5 × 3

NoData and data types

Every raster band can declare a NoData value, and the block tracks which cells hold it. Ignoring that is the most common source of wrong statistics: an elevation mean dragged down by thousands of -9999 cells looks plausible and is meaningless.

band = 1
print("source NoData:", provider.sourceHasNoDataValue(band),
      provider.sourceNoDataValue(band))
print("data type:", provider.dataType(band))

valid = [block.value(r, c)
         for r in range(block.height())
         for c in range(block.width())
         if not block.isNoData(r, c)]
print(f"{len(valid)} valid of {block.width() * block.height()} cells")
print("min", min(valid), "max", max(valid))

Breakdown: sourceHasNoDataValue reports whether the file declares one; useSourceNoDataValue controls whether QGIS honours it, and user-defined NoData ranges set in the layer properties are applied on top. isNoData accounts for all of those, which is why it is more reliable than comparing values against -9999 yourself. The data type matters beyond correctness: an Int16 DEM with values in decimetres, or a UInt16 satellite band with a scale factor, needs converting before the numbers mean anything — check the band's scale and offset with provider.bandScale(band) and provider.bandOffset(band).

Convert a block to a NumPy array

For anything beyond a few lookups, convert the block's raw bytes to a NumPy array in one step and do the arithmetic there. The byte layout is the provider's native data type, row by row.

import numpy as np
from qgis.core import Qgis

NUMPY_TYPES = {
    Qgis.DataType.Byte: np.uint8,
    Qgis.DataType.Int16: np.int16,
    Qgis.DataType.UInt16: np.uint16,
    Qgis.DataType.Int32: np.int32,
    Qgis.DataType.UInt32: np.uint32,
    Qgis.DataType.Float32: np.float32,
    Qgis.DataType.Float64: np.float64,
}

def block_to_array(block):
    dtype = NUMPY_TYPES[block.dataType()]
    data = np.frombuffer(bytes(block.data()), dtype=dtype)
    array = data.reshape(block.height(), block.width()).astype(np.float64)
    if block.hasNoDataValue():
        array[array == block.noDataValue()] = np.nan
    return array

heights = block_to_array(block)
print("mean", np.nanmean(heights), "p95", np.nanpercentile(heights, 95))
print("cells above 120 m:", int(np.sum(heights > 120)))

Breakdown: block.data() returns a QByteArray of the raw cell bytes; converting it to Python bytes and handing it to np.frombuffer avoids a Python loop entirely. Casting to float64 and replacing NoData with NaN lets the nan-aware NumPy functions ignore missing cells without masks. On a 5,000 × 5,000 block this is thousands of times faster than calling value() per cell. When the source uses a NoData value that also occurs as real data — a DEM with genuine zero heights and NoData of zero — the problem is in the data, not the code, and no conversion can fix it.

Whole raster or tiles?A 40,000 by 40,000 cell Float32 raster needs about 6 GB as one block, more once converted to float64. Reading it in 2,048 by 2,048 tiles needs about 32 MB at a time. Running results such as counts, sums and histograms are accumulated per tile, so the answer is identical while memory stays flat.Same result, very different memoryone block40,000² Float32≈ 6 GB×2 as float64tiles2,048² at a time≈ 32 MBaccumulate per tile

Read a large raster tile by tile

A national DEM at 1 m does not fit in memory as one block. Walk it in fixed-size tiles aligned to the cell grid and accumulate whatever you need — counts, sums, a histogram — as you go.

def iter_tiles(layer, band=1, tile=2048):
    prov = layer.dataProvider()
    ext = layer.extent()
    cols, rows = layer.width(), layer.height()
    cx, cy = layer.rasterUnitsPerPixelX(), layer.rasterUnitsPerPixelY()
    for row0 in range(0, rows, tile):
        for col0 in range(0, cols, tile):
            w = min(tile, cols - col0)
            h = min(tile, rows - row0)
            rect = QgsRectangle(
                ext.xMinimum() + col0 * cx, ext.yMaximum() - (row0 + h) * cy,
                ext.xMinimum() + (col0 + w) * cx, ext.yMaximum() - row0 * cy,
            )
            yield row0, col0, prov.block(band, rect, w, h)

bins = np.arange(0, 1001, 10)
histogram = np.zeros(len(bins) - 1, dtype=np.int64)
count = total = 0

for row0, col0, blk in iter_tiles(dem):
    arr = block_to_array(blk)
    values = arr[~np.isnan(arr)]
    histogram += np.histogram(values, bins=bins)[0]
    count += values.size
    total += values.sum()

print("valid cells:", count, "mean:", total / count)
print("most common 10 m band starts at", bins[histogram.argmax()], "m")

Breakdown: Tiles are computed in cell coordinates and converted to map coordinates, so every tile is exactly aligned and the last row and column are shortened rather than resampled. A tile size that matches the file's internal block size — often 256, 512 or 1024 in a tiled GeoTIFF — lets GDAL read whole internal blocks, which is noticeably faster; gdalinfo reports it. Accumulating a histogram and running sums gives exact results in constant memory; percentiles need either the histogram approach or a sample, because an exact median of forty billion values needs them all. For per-zone results, zonal statistics already does the tiling for you.

Look up single cells quickly

For a handful of points — a plugin that shows the height under the cursor, say — identify or sample avoid building blocks at all.

from qgis.core import QgsPointXY, QgsRaster

point = QgsPointXY(429812.4, 433120.7)
value, ok = provider.sample(point, 1)
print("sample:", value if ok else "outside or NoData")

result = provider.identify(point, QgsRaster.IdentifyFormatValue)
if result.isValid():
    print("all bands:", result.results())

Breakdown: sample returns a single band's value and a success flag, and is the fastest way to read one cell. identify returns every band at once as a dictionary keyed by band number, which suits multispectral imagery. Both expect the point in the layer's CRS; transform it first if it comes from the canvas, as in transforming point coordinates. For many points, sampling raster values at points with the Processing algorithm is simpler than a loop.

QGIS version compatibility

QgsRasterBlock, provider.block and identify are unchanged since QGIS 3.0; provider.sample arrived in 3.4. Qgis.DataType is the scoped spelling used from 3.30 and required on the QGIS 4 series — older code may use Qgis.Float32 directly. QgsRaster.IdentifyFormatValue becomes Qgis.RasterIdentifyFormat.Value on recent releases. The NumPy conversion depends only on the byte layout, which has not changed.

Troubleshooting

  • Values look smoothed. The window was not aligned or the size did not match the cell count, so the provider resampled.
  • The array has the wrong shape or garbage values. The NumPy dtype does not match block.dataType().
  • Statistics are far too low. NoData cells were included; use isNoData or convert to NaN.
  • Reading is slow on a network drive. Tile size does not match the file's internal blocks, or the file is not tiled at all; convert to a Cloud-Optimized GeoTIFF.
  • block.data() is empty. The request fell entirely outside the raster extent.

Conclusion

Request blocks for grid-aligned windows at native size, respect NoData through isNoData or NaN, and convert whole blocks to NumPy rather than reading cells one at a time. Walk large rasters in aligned tiles and accumulate results so memory stays flat, and use sample or identify when you only need a few cells.

Frequently Asked Questions

Can I read a downsampled overview for speed? Yes — request a smaller width and height for the same extent. The provider uses the file's overviews when they exist, which is how renderers read quickly.

Does this work for multiband imagery? Yes. Request one block per band, or use identify for all bands at a point.

Can I write values back through the same block? Not to a read-only layer. Create a writable provider with QgsRasterFileWriter, as in writing a NumPy array to a raster.

Should I use GDAL directly instead? GDAL's ReadAsArray is equally fast for files. The provider route works for any raster QGIS can open, respects layer NoData settings, and needs no second library in a plugin.