Sample Raster Values at Points in PyQGIS

"What is the elevation at each of my survey stations?" is one of the most common questions in applied GIS, and it has two very different answers in PyQGIS depending on how many points there are. For a whole layer, one Processing algorithm does the job and writes the values into new fields. For a handful of coordinates, a direct provider call returns a number in microseconds. Using the wrong one turns a two-second job into a twenty-minute one.

This recipe belongs to Raster Analysis Workflows in PyQGIS. It covers both routes, the CRS mismatch that silently returns nothing, nodata handling, and when a point sample should really have been a zonal statistic.

What point sampling actually readsPoints are drawn over a raster grid. Each point takes the value of the single cell it falls inside, regardless of how close it is to a boundary. One point falls on a nodata cell and returns no value. One point lies outside the raster extent and returns an invalid result rather than zero.One point, one cell — proximity to the edge changes nothingnodata112118124131109127135104108121130inside a valid cell — returns that cell's value112, even though 118 is two pixels awayinside a nodata cell — returns no valuethe field is NULL, not zerooutside the raster extent — invalidcheck the validity flag, never the value

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • A point layer and a raster. They do not have to share a CRS for the Processing route, but they do for the direct API.
  • Knowledge of the raster's nodata value — provider.sourceNoDataValue(1) reports it.

Sample a whole layer

import processing

result = processing.run("native:rastersampling", {
    "INPUT": "/data/stations.gpkg|layername=stations",
    "RASTERCOPY": "/data/dem.tif",
    "COLUMN_PREFIX": "elev_",
    "OUTPUT": "/data/output/stations_elev.gpkg",
})

Breakdown: The algorithm adds one field per raster band, named with the prefix and the band number — elev_1, elev_2 and so on — so a single-band DEM produces elev_1. Points outside the raster or on nodata cells get a NULL rather than a fabricated value, which is the behaviour you want and the reason to prefer this over a hand-written loop. Reprojection between the point layer's CRS and the raster's happens automatically, using the project's transform context.

For several rasters — elevation, slope, rainfall — run it once per raster, feeding each output into the next. The fields accumulate, and one pass is added per raster rather than one pass per point.

Sample a single coordinate

from qgis.core import QgsPointXY, QgsRasterLayer

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

value, ok = provider.sample(QgsPointXY(433250, 189400), 1)
if ok:
    print(f"elevation {value:.1f} m")
else:
    print("no data at that location")

Breakdown: sample() takes a point in the raster's own CRS and a one-based band number, and returns a (value, valid) tuple. The validity flag is the part that matters: a point outside the extent or on a nodata cell returns (nan, False), and code that reads only the first element will happily record a NaN as an elevation. This call is a direct read with no algorithm overhead, which makes it the right tool inside an interactive map tool — see Identify the Feature Under the Cursor in PyQGIS for the canvas side of that pattern.

Transform coordinates before sampling

Unlike the algorithm, sample() does no reprojection. Give it a point in the wrong CRS and it will return values from wherever those numbers happen to land — usually outside the extent, occasionally inside it, which is the worse outcome.

from qgis.core import QgsCoordinateTransform, QgsCoordinateReferenceSystem, QgsProject

transform = QgsCoordinateTransform(
    QgsCoordinateReferenceSystem("EPSG:4326"),
    raster.crs(),
    QgsProject.instance(),
)

point = transform.transform(QgsPointXY(-1.4701, 53.3811))
value, ok = provider.sample(point, 1)

Breakdown: Building the transform once outside the loop matters — constructing it per point is one of the classic profile-topping mistakes described in Profile Slow PyQGIS Code. Passing the project as the third argument gives the transform access to the configured datum transformations, so the result matches what the canvas shows. The mechanics of transformation are covered in Transform Point Coordinates in PyQGIS.

Sampling many points efficiently

When the points are already a layer, the algorithm wins. When they are generated — a grid, a route, model output — sample directly, and keep the loop lean.

from qgis.core import QgsFeatureRequest

request = QgsFeatureRequest().setSubsetOfAttributes([])
samples = []

for feature in points.getFeatures(request):
    point = feature.geometry().asPoint()
    value, ok = provider.sample(point, 1)
    samples.append(value if ok else None)

print(f"{sum(1 for s in samples if s is None)} points had no value")

Breakdown: Requesting no attributes halves the iteration cost when only geometry is needed, as explained in Speed Up Feature Iteration with QgsFeatureRequest. Storing None rather than a sentinel keeps the missing values honestly missing when the list is later written to a field. Counting the failures and reporting them is the habit that catches a CRS mismatch immediately: if every point failed, the coordinates are in the wrong system.

Point sample, interpolate, or summarise an areaSampling one cell suits a value that genuinely belongs to a location, such as a land cover class. Interpolating between neighbouring cells suits a continuous surface where the point sits near a cell boundary. Zonal statistics suit a value that belongs to an area rather than a point, such as mean rainfall over a catchment.Is the answer really a single cell?point samplenative:rastersamplingcategorical rastersland cover, soil classaveraging classes is meaninglessinterpolatewarp to a finer grid firstcontinuous surfaceselevation, temperaturesmooths the cell-edge stepzonal statisticsnative:zonalstatisticsfbthe value belongs to an areamean rainfall per catchmenta centroid sample is not a mean

Sample a stack of dated rasters

Monitoring work rarely samples one raster. It samples the same stations against a monthly series, and the result is a table with one column per date.

from pathlib import Path
import processing

rasters = sorted(Path("/data/rainfall").glob("rainfall_*.tif"))
current = "/data/stations.gpkg|layername=stations"

for raster in rasters:
    period = raster.stem.split("_")[-1]          # e.g. 2026_07 -> "07"
    current = processing.run("native:rastersampling", {
        "INPUT": current,
        "RASTERCOPY": str(raster),
        "COLUMN_PREFIX": f"rain_{period}_",
        "OUTPUT": "TEMPORARY_OUTPUT",
    })["OUTPUT"]

processing.run("native:savefeatures", {
    "INPUT": current,
    "OUTPUT": "/data/output/stations_rainfall.gpkg",
})

Breakdown: Each pass takes the previous result as its input, so the columns accumulate on one layer rather than producing a dozen layers to join afterwards. Sorting the paths matters: the column order in the output then follows the chronological order, which is what makes the table readable and any later reshaping straightforward. Keeping every intermediate as TEMPORARY_OUTPUT avoids writing eleven files that exist only to feed the next step, and the single native:savefeatures at the end is the one write to disk.

Two cautions apply to a series. All the rasters must share a grid, or each pass silently samples a slightly different footprint — check with layer.extent() and layer.rasterUnitsPerPixelX() before starting rather than after. And a wide result is often not the shape you want: a table with one row per station per date is far easier to chart or aggregate, which means reshaping the wide output afterwards, most simply by exporting to CSV and pivoting there — see Export an Attribute Table to CSV in PyQGIS.

One layer, one column per dateMonthly rainfall rasters are sampled one after another against the same station layer. Each pass takes the previous result as its input and appends a column named for its month, so the final table carries one row per station and one column per date rather than a dozen separate layers to join.Chain the passes; do not join twelve layers afterwardsrainfall_2026_05.tifrainfall_2026_06.tifrainfall_2026_07.tifrastersamplingonce per rasterstations_rainfall.gpkgstation rain_05_1 rain_06_1 rain_07_1S-014 61.2 48.9 12.4S-021 58.7 51.3 NULLA NULL is a station that fell on nodata that month, not a zero

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.16 LTR3.7native:rastersampling added; sample() present since 3.4.
3.28 LTR3.9Behaviour matches this page.
3.34 LTR3.12Baseline for this page.
3.40 / 3.443.12Unchanged; zonal statistics gained additional statistic types.

Troubleshooting

  • Every sampled value is NULL. The point CRS and raster CRS differ and you used sample(). Transform first, or use the algorithm.
  • Values are wrong but plausible. Same cause, when the coordinates happen to land inside the extent. Check one point against the canvas readout.
  • NaN appears in the output. The validity flag was ignored. Test the second element of the tuple.
  • The band is wrong. sample() is one-based. Band 1 is the first band, not band 0.
  • Sampling a multi-band raster is slow. Each call reads one band; for many bands, native:rastersampling reads them in one pass.
  • A point on a boundary gives an unexpected value. Cell membership follows the raster's geotransform; a point exactly on a boundary belongs to one specific cell. If the boundary matters, interpolate rather than sample.

Conclusion

Sampling a layer is native:rastersampling with a column prefix — it reprojects, honours nodata and writes NULLs where there is no value. Sampling a coordinate is provider.sample(), which is fast, does not reproject, and must have its validity flag checked. If the value logically belongs to an area rather than a location, the question was a zonal statistic all along.

Frequently Asked Questions

Does sampling interpolate between cells? No. Both routes return the value of the cell containing the point. For an interpolated value, warp the raster to a finer grid with bilinear resampling and sample that.

How do I sample several rasters into one layer? Run the algorithm once per raster, chaining each output into the next input. Each pass appends its own prefixed fields.

Can I sample along a line? Convert the line to points first with native:pointsalonglines, then sample those. That is also how a terrain profile is built.

Why is the added field an integer when the raster is float? The output format's field type inference. Write to GeoPackage rather than shapefile, and check the raster's actual data type.

Can I sample a raster that is only available over the network? Yes, if QGIS can open it — a Cloud Optimized GeoTIFF over HTTP works, and only the tiles containing your points are fetched. Expect the first sample to be slow while the header is read, and set a timeout so an unreachable host fails rather than hangs.

How do I record which raster a value came from? Include it in the column prefix, as the dated series above does. A value in a column called elev_1 is anonymous a year later; one in elev_lidar2024_1 explains itself.

What does sampling return on a mosaic with gaps? Nodata, and therefore a NULL. That is correct — a gap is an absence of measurement, not a zero.