IDW Interpolation from Points in PyQGIS

Inverse distance weighting fills the space between measurements by averaging nearby samples, weighted so that closer ones count for more. It makes very few assumptions, which is its strength, and it will happily produce a confident-looking value fifty kilometres from the nearest observation, which is its weakness. Getting a useful surface out of it is mostly about constraining where it is allowed to speak.

This recipe belongs to Terrain & Interpolation Analysis in PyQGIS. It covers the awkward layer-specification string the algorithm expects, choosing the distance power and pixel size, masking the output to where the data supports it, and cross-validating so the error has a number attached.

How one interpolated cell gets its valueAn output cell is surrounded by sample points at different distances. Each sample contributes a weight proportional to one over its distance raised to the chosen power. A power of one spreads influence broadly; a power of four concentrates it on the nearest sample, producing plateaus around each observation.The power decides how quickly influence dies awayone output cellnear — heavyfar — lightw = 1 / dᵖthe same data, three powersp = 1 · broadp = 2 · defaultp = 4 · plateausIDW never produces a value outside the range of the observations

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • A point layer with a numeric attribute to interpolate, in a projected CRS. Distances in degrees weight north–south and east–west differently and skew the surface.
  • Enough samples to be worth interpolating. Below roughly thirty points, a surface implies far more than the data can support.

Build the interpolation data string

The algorithm takes its input as an encoded string rather than as parameters, which is the only genuinely unpleasant part of the process.

from qgis.core import QgsProject

layer = QgsProject.instance().mapLayersByName("rain_gauges")[0]
attribute_index = layer.fields().indexOf("annual_mm")
if attribute_index < 0:
    raise LookupError("field 'annual_mm' not found on the layer")

spec = f"{layer.source()}::~::0::~::{attribute_index}::~::0"

Breakdown: The four ::~::-separated fields are the layer source, the geometry-type flag (0 for points), the zero-based index of the attribute to interpolate, and a use-Z flag (0 to use the attribute, 1 to use the geometry's Z instead). fields().indexOf() returns -1 for a missing field, and passing -1 into the spec produces an empty raster rather than an error — so the guard above is not optional. Multiple layers can be interpolated together by joining several specs with ;, which is how spot heights and contour vertices are combined into one surface.

Using layer.source() rather than a hard-coded path keeps the script working when the layer came from a GeoPackage with a layer name appended, which a plain path would lose.

Run the interpolation

With the spec built, the call itself is straightforward.

import processing

extent = layer.extent()
extent.grow(500)          # 500 m of margin, in layer units

processing.run("qgis:idwinterpolation", {
    "INTERPOLATION_DATA": spec,
    "DISTANCE_COEFFICIENT": 2.0,
    "EXTENT": extent,
    "PIXEL_SIZE": 100,
    "OUTPUT": "/data/output/rainfall_idw.tif",
})

Breakdown: DISTANCE_COEFFICIENT is the power p. At 2 — the near-universal default — a sample twice as far away counts a quarter as much, which produces a surface that is smooth without ignoring local detail. PIXEL_SIZE is in layer units and should be no finer than roughly a quarter of the mean spacing between samples; going finer produces a beautifully smooth raster whose extra detail is entirely manufactured. Growing the extent slightly gives room for the mask step below without the surface being clipped by its own bounding box.

QGIS's IDW uses every sample for every cell — there is no search radius or neighbour limit in this algorithm. That makes it predictable and slow: cost grows with points multiplied by cells, so ten thousand samples on a fine grid is genuinely expensive. Where a neighbour limit is needed, gdal:gridinversedistancenearestneighbor provides one.

Choosing the power deliberately

The distance power is the one parameter that changes the story the surface tells, and two minutes of comparison settles it better than any rule of thumb.

import processing

for power in (1.0, 2.0, 3.0, 5.0):
    processing.run("qgis:idwinterpolation", {
        "INTERPOLATION_DATA": spec,
        "DISTANCE_COEFFICIENT": power,
        "EXTENT": extent,
        "PIXEL_SIZE": 100,
        "OUTPUT": f"/data/output/idw_p{int(power)}.tif",
    })

Breakdown: Writing four rasters rather than one costs a few minutes of compute and gives something to look at side by side. At a power of 1 the surface tends towards a broad regional average and local measurements barely register. At 5 each sample sits on its own plateau with steep walls between, and the map reads as a set of discs rather than a surface — a pattern usually described as bull's-eyeing, and a reliable sign the power is too high for the sample density.

The choice is not purely aesthetic. A higher power is defensible when the phenomenon genuinely changes sharply over short distances — contamination around point sources, noise around a runway. A lower power suits something that varies gradually and where each sample is a noisy estimate of a smooth field, such as annual rainfall. Where you cannot argue for either, 2 is the convention and requires no defence.

One thing the power cannot do is make the surface exceed the data. IDW is a weighted average, so every output value lies between the minimum and maximum observation. A summit higher than the highest spot height, or a hollow deeper than the deepest borehole, will never appear — which is a genuine limitation when interpolating terrain and a useful safety property when interpolating measurements.

Mask the output to where the data supports it

This is the step that turns a plausible picture into a defensible one, and it is the step almost everyone skips.

buffered = processing.run("native:buffer", {
    "INPUT": layer,
    "DISTANCE": 5000,                  # a defensible influence radius
    "SEGMENTS": 12,
    "DISSOLVE": True,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

processing.run("gdal:cliprasterbymasklayer", {
    "INPUT": "/data/output/rainfall_idw.tif",
    "MASK": buffered,
    "NODATA": -9999,
    "CROP_TO_CUTLINE": True,
    "OUTPUT": "/data/output/rainfall_idw_masked.tif",
})

Breakdown: Buffering every sample by a distance you can justify, dissolving the result, and clipping to it converts the implicit claim "this surface is valid everywhere" into the explicit claim "this surface is valid within 5 km of a gauge". The buffer distance is a judgement about the phenomenon — rainfall varies over kilometres, soil chemistry over metres — and writing it as a named constant with a comment is worth more than any amount of parameter tuning. DISSOLVE: True merges overlapping buffers so the mask is one polygon rather than hundreds.

Masking states what the surface actually coversThe unmasked interpolation fills its entire rectangular extent with values, including large areas with no nearby sample. Clipping to a dissolved buffer around the sample points removes those areas, so the map shows values only where an observation is within a justified distance.The gaps are information, not a defectunmaskedvalues here toono sample within 40 km, but a number anywaymasked to buffered samplesnodatathe map now says where it does not know

Cross-validate, so the error has a number

Holding samples back and testing against them converts an opinion about quality into a statistic.

import processing

holdout = processing.run("native:randomextract", {
    "INPUT": layer,
    "METHOD": 1,                       # percentage
    "NUMBER": 10,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

sampled = processing.run("native:rastersampling", {
    "INPUT": holdout,
    "RASTERCOPY": "/data/output/rainfall_idw.tif",
    "COLUMN_PREFIX": "pred_",
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

errors = [
    abs(f["annual_mm"] - f["pred_1"])
    for f in sampled.getFeatures()
    if f["pred_1"] is not None
]
print(f"mean absolute error: {sum(errors) / len(errors):.1f} mm")

Breakdown: Strictly the held-back points should be excluded from the interpolation before predicting them; interpolating from the full set and then testing on a subset of it flatters the result badly, because IDW passes almost exactly through each observation. Run the interpolation on the complement of the holdout for an honest figure. native:rastersampling writes one column per band with the given prefix, so a single-band raster gives pred_1. Guarding against None handles held-back points that fall outside the mask.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.16 LTR3.7qgis:idwinterpolation present with the ::~:: spec format.
3.22 LTR3.9native:rastersampling available for the validation step.
3.28 LTR3.9native:randomextract percentage method stable.
3.34 LTR3.12Baseline for this page.
3.40+3.12Interpolation algorithms accept layer objects more tolerantly in the spec string.

Troubleshooting

  • The output raster is empty. The attribute index in the spec is wrong, often -1 from a missing field. Print the index before building the string.
  • The surface is a set of bull's-eyes. The power is too high. Drop DISTANCE_COEFFICIENT towards 2.
  • The surface is almost flat. The power is too low, or the samples genuinely have little spatial structure. Compare against a power of 3.
  • It takes forever. Every sample is used for every cell. Coarsen PIXEL_SIZE, or switch to gdal:gridinversedistancenearestneighbor with a neighbour limit.
  • The surface is stretched east–west. The layer is in a geographic CRS. Reproject the points first.
  • Values appear far from any sample. That is IDW working as designed. Mask the output.

Conclusion

Build the spec carefully and check the attribute index, keep the power near 2 unless there is a reason, interpolate no finer than a quarter of the sample spacing, and always clip the result to a buffered mask so the map admits where it has no data. Then hold back a tenth of the samples and publish the error alongside the surface.

Frequently Asked Questions

Can I limit IDW to a search radius? Not in qgis:idwinterpolation, which uses every point. gdal:gridinversedistancenearestneighbor takes a maximum number of neighbours and search radii, at the cost of a slightly different parameter vocabulary.

How do I interpolate the Z coordinate instead of an attribute? Set the last field of the spec to 1 and leave the attribute index at -1. That is the route for a layer of 3D spot heights or contour vertices.

Should I use IDW or TIN? IDW when the phenomenon is smooth and the samples are dense; TIN when the observations are trusted exactly and you want the surface to refuse to extrapolate.

Does IDW handle duplicate points? Two samples at the same coordinate give a zero distance and an infinite weight. QGIS handles it without crashing, but the result is one of the two values arbitrarily. Deduplicate first with count features by attribute or a spatial deduplication step.