PyQGIS Engineering & Automation

Reading Raster Pixel Values with QgsRasterDataProvider

Sample raster values in PyQGIS: QgsRasterDataProvider.identify at a point, block reads with QgsRasterBlock, band selection, nodata handling, and map-to-pixel…

TL;DR: For a single point use provider.identify(QgsPointXY(x, y), QgsRaster.IdentifyFormatValue).results() and read the band-keyed dictionary; for a region use provider.block(band, extent, cols, rows) to get a QgsRasterBlock and read block.value(row, col), always testing block.isNoData(row, col) before trusting a cell.

This page is part of the Vector and Raster Data Access Patterns in PyQGIS guide, which sits inside the broader PyQGIS Core Architecture & Data Handling reference. Raster sampling looks trivial until you meet nodata sentinels, multi-band datasets, and the moment your sample point and your raster disagree about which coordinate reference system they live in. This deep-dive gives you two reliable primitives — point identify and windowed block reads — and the conversion maths that ties map coordinates to pixel indices.

Complete Runnable Sampling Utility

Drop this module into any headless script or the QGIS Python console. It loads a raster, samples a single value through identify(), then reads a small window through block() with full nodata filtering. Replace the path and coordinate to match your data.

python
"""
raster_sampler.py

Read raster pixel values in PyQGIS via QgsRasterDataProvider.
Covers single-point identify(), windowed block() reads, nodata
handling, and map-coordinate to pixel-index conversion.

Compatible with QGIS 3.28+ LTR.
"""
from __future__ import annotations

from typing import Optional

from qgis.core import (
    QgsPointXY,
    QgsRaster,
    QgsRasterBlock,
    QgsRasterDataProvider,
    QgsRasterLayer,
    QgsRectangle,
)


def load_provider(path: str, name: str = "raster") -> QgsRasterDataProvider:
    """
    Load a raster layer and return its data provider.

    Raises:
        ValueError: if the raster fails to load or has no provider.
    """
    layer = QgsRasterLayer(path, name)
    if not layer.isValid():
        raise ValueError(f"Raster failed to load: {path}")
    provider = layer.dataProvider()
    if provider is None:
        raise ValueError(f"Raster has no data provider: {path}")
    return provider


def sample_point(
    provider: QgsRasterDataProvider,
    x: float,
    y: float,
    band: int = 1,
) -> Optional[float]:
    """
    Sample a single raster value at a map coordinate.

    The coordinate MUST already be in the raster's CRS. identify()
    returns a results dict keyed by 1-based band number; a value of
    None means the point fell on nodata or outside the extent.

    Returns:
        The band value, or None for nodata / out-of-extent.
    """
    point = QgsPointXY(x, y)
    result = provider.identify(point, QgsRaster.IdentifyFormatValue)
    if not result.isValid():
        return None
    # results() keys are band numbers (1-based); values may be None.
    return result.results().get(band)


def map_to_pixel(
    provider: QgsRasterDataProvider,
    x: float,
    y: float,
) -> tuple[int, int]:
    """
    Convert a map coordinate to a (row, col) pixel index.

    Uses the provider extent and raster dimensions to derive the
    per-pixel resolution, then floors the offset from the top-left
    origin. Y increases downward in raster space, hence (top - y).
    """
    extent: QgsRectangle = provider.extent()
    width: int = provider.xSize()
    height: int = provider.ySize()

    x_res: float = extent.width() / width
    y_res: float = extent.height() / height

    col = int((x - extent.xMinimum()) / x_res)
    row = int((extent.yMaximum() - y) / y_res)
    return row, col


def read_window(
    provider: QgsRasterDataProvider,
    extent: QgsRectangle,
    cols: int,
    rows: int,
    band: int = 1,
) -> list[list[Optional[float]]]:
    """
    Read a rectangular window of pixels into a nested list.

    Requests a QgsRasterBlock resampled to the given cols x rows grid
    over `extent`, then reads each cell with isNoData() filtering.
    Nodata cells are returned as None.
    """
    block: QgsRasterBlock = provider.block(band, extent, cols, rows)
    grid: list[list[Optional[float]]] = []
    for r in range(rows):
        line: list[Optional[float]] = []
        for c in range(cols):
            if block.isNoData(r, c):
                line.append(None)
            else:
                line.append(block.value(r, c))
        grid.append(line)
    return grid


if __name__ == "__main__":
    prov = load_provider("/data/dem/elevation.tif", "dem")

    # 1. Single-point sample (coordinate in the raster's CRS).
    value = sample_point(prov, 13.4050, 52.5200, band=1)
    print(f"point value: {value}")

    # 2. Windowed read: a 5x5 grid over a small extent.
    ext = prov.extent()
    window = QgsRectangle(
        ext.xMinimum(),
        ext.yMaximum() - ext.height() * 0.1,
        ext.xMinimum() + ext.width() * 0.1,
        ext.yMaximum(),
    )
    values = read_window(prov, window, cols=5, rows=5, band=1)
    for line in values:
        print(line)

Map Coordinate to Pixel Index

The one conversion that trips people up is turning a map coordinate into a (row, col) pixel index. The provider extent gives you the geographic bounds; the raster dimensions give you the pixel grid; dividing yields per-pixel resolution. The Y axis is inverted — row 0 sits at the top (yMaximum), so you subtract the coordinate from the top edge rather than the bottom.

Map Coordinate to Pixel Index Conversion A raster grid with its top-left origin at xMinimum, yMaximum. A sample point maps to a column via (x minus xMinimum) divided by x resolution, and to a row via (yMaximum minus y) divided by y resolution. origin (xMinimum, yMaximum) (x, y) col = (x − xMinimum) / xRes row = (yMaximum − y) / yRes xRes = extent.width() / xSize() yRes = extent.height() / ySize() Row 0 is the TOP row: raster Y increases downward, map Y increases upward. int() floors to the enclosing cell index.

You rarely need map_to_pixel() when you use identify() or block() directly — both accept map coordinates and a map extent respectively. Reach for it when you already hold a QgsRasterBlock for a known window and want to address a specific cell inside it without a second provider round-trip.

Architecture Breakdown

identify() — the point-sampling contract

QgsRasterDataProvider.identify(point, format) is the fastest way to pull one value. The format argument is a QgsRaster.IdentifyFormat enum — use QgsRaster.IdentifyFormatValue for numeric band values. Other formats exist (IdentifyFormatHtml for the desktop identify panel, IdentifyFormatFeature for WMS getFeatureInfo), but only IdentifyFormatValue returns clean numbers.

The call returns a QgsRasterIdentifyResult. Two members matter: isValid() tells you whether the query succeeded at all, and results() returns a Python dict keyed by 1-based band number whose values are floats — or None when the point lands on nodata or falls outside the raster extent. That None is the reason you cannot skip the validity check: a syntactically successful identify can still yield {1: None}.

python
from qgis.core import QgsPointXY, QgsRaster

result = provider.identify(QgsPointXY(x, y), QgsRaster.IdentifyFormatValue)
if result.isValid():
    band_values: dict = result.results()  # {1: 42.0, 2: 17.5, ...}
    band_1 = band_values.get(1)

The point coordinate must be expressed in the raster’s own CRS. identify() performs no reprojection — a point in EPSG:4326 sampled against a raster in EPSG:25832 simply misses, returning None or an invalid result with no error. Confirm alignment before you sample; the Coordinate Transformations and CRS Handling in PyQGIS guide covers building a QgsCoordinateTransform to bring the point into the raster CRS first.

block() and QgsRasterBlock — reading regions

For anything larger than a handful of points, identify() in a loop is wasteful — every call re-enters the provider. provider.block(band, extent, cols, rows) reads a whole rectangular window in one request and returns a QgsRasterBlock. You pass a QgsRectangle extent in the raster CRS plus the target grid size (cols × rows); the provider resamples the source pixels onto that grid, so you can read data at native resolution (match cols/rows to the extent’s pixel count) or downsample deliberately.

Read individual cells with block.value(row, col), which returns a float. Both indices are 0-based, with row 0 at the top of the extent. Before trusting any cell, call block.isNoData(row, col) — it accounts for the source nodata value, out-of-range results, and unreadable pixels in one test. Reading value() on a nodata cell returns the raw sentinel (often a very large negative number), which will silently poison any statistics you compute.

python
block = provider.block(1, extent, cols, rows)
total, count = 0.0, 0
for r in range(rows):
    for c in range(cols):
        if not block.isNoData(r, c):
            total += block.value(r, c)
            count += 1
mean = total / count if count else float("nan")

Memory is the constraint that decides your window size. A QgsRasterBlock holds the entire requested grid in RAM at the band’s data type width — a 10000 × 10000 Float32 block is roughly 400 MB before you have read a single value. For large datasets, tile the extent and process windows sequentially rather than requesting one enormous block, and drop each block reference as soon as the window is consumed so it can be collected. The tiling and reference-lifetime patterns in Preventing Memory Leaks When Processing Large GeoTIFF Rasters apply directly here.

Band selection and nodata semantics

Band numbers are 1-based everywhere in the provider API — identify() keys, the first argument to block(), and sourceNoDataValue(band). Query provider.bandCount() before assuming a band exists, and provider.dataType(band) to learn the numeric type (for example Qgis.DataType.Float32 or Qgis.DataType.Byte), which governs both value range and the memory footprint of a block.

Nodata handling has two layers. provider.sourceNoDataValue(band) returns the sentinel baked into the file, and provider.sourceHasNoDataValue(band) reports whether one is defined. On top of that, a user or a preceding step can set additional nodata ranges through the provider’s user-nodata settings. block.isNoData() respects both layers, which is exactly why it is more reliable than comparing block.value() against sourceNoDataValue() by hand — a manual comparison misses user-defined nodata and is fragile against floating-point equality on the sentinel.

CRS alignment of the sample point

Neither identify() nor block() reprojects. The sample point and the block extent must be in the raster’s CRS, available from provider.crs() (equivalently layer.crs()). If your input coordinates come from another source — a GPS fix in EPSG:4326, features from a layer in a project CRS — transform them first:

python
from qgis.core import (
    QgsCoordinateTransform,
    QgsCoordinateTransformContext,
    QgsPointXY,
    QgsProject,
)

src_crs = QgsProject.instance().crs()      # e.g. the point's CRS
dst_crs = provider.crs()                    # the raster's CRS
ctx = QgsProject.instance().transformContext()
xform = QgsCoordinateTransform(src_crs, dst_crs, ctx)

raster_point: QgsPointXY = xform.transform(QgsPointXY(lon, lat))
value = provider.identify(raster_point, QgsRaster.IdentifyFormatValue)

A misaligned sample is the single most common cause of “identify keeps returning None”: the point is valid, the raster is valid, but they describe locations on different coordinate systems. Always confirm provider.crs().authid() matches your point CRS, or transform explicitly.

Production Best Practices

  • Confirm layer.isValid() and provider is not None before sampling; a broken path yields an invalid layer whose provider silently returns nothing.
  • Sample in the raster’s CRS. Read provider.crs() and reproject the point with QgsCoordinateTransform if it originates elsewhere — neither identify() nor block() reprojects for you.
  • Check result.isValid() after identify(), then treat a None band value as nodata or out-of-extent rather than an error.
  • Always gate block.value(row, col) behind block.isNoData(row, col). Raw values on nodata cells are sentinels that corrupt sums, means, and min/max.
  • Use 1-based band numbers in identify() results, block(), and sourceNoDataValue(); validate against provider.bandCount().
  • Size blocks against provider.dataType(band). Estimate cols × rows × bytes-per-sample and tile large extents instead of requesting one block that will not fit in RAM.
  • Prefer block() over looping identify() whenever you need more than a few dozen samples — one windowed read beats many single-point round-trips.
  • Release each QgsRasterBlock as soon as its window is processed so large buffers are reclaimed promptly.