Resample and Align Rasters in PyQGIS

Raster arithmetic assumes that pixel (row, column) in one raster covers the same ground as pixel (row, column) in another. That holds only when both share a CRS, pixel size, origin and extent — a common grid. Real inputs rarely do: a 25 m DEM, a 10 m land-cover map, a 30 m satellite product and a 1 km climate grid, each with its own extent. Combining them correctly means choosing one target grid and bringing every input onto it, with a resampling method that respects what each raster's values mean.

This recipe belongs to Raster Analysis Workflows. It explains what a grid is, resamples rasters to a new resolution, aligns them to a reference raster, picks the right resampling method for continuous and categorical data, and verifies alignment before arithmetic.

Four properties define a gridA raster grid is defined by its CRS, its pixel size in x and y, its origin which is the coordinate of the top-left corner, and its extent or number of rows and columns. Two rasters are aligned only when all four match. A half-pixel shift in origin, with everything else equal, already misaligns every pixel.Same CRS, size, origin, extentCRSsame systemsame unitspixel sizex and ye.g. 10 morigintop-left corneron the same latticeextentrows × columnssame windowa half-pixel shift in origin misaligns every pixel

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series, with the GDAL provider.
  • A reference raster that defines the target grid — usually the most important or finest input of the analysis.

Inspect the grids you have

Before resampling anything, compare the grids. A small function reporting CRS, pixel size, origin and size makes the differences obvious.

from qgis.core import QgsRasterLayer

def grid(path):
    r = QgsRasterLayer(path, "g")
    e = r.extent()
    return {"crs": r.crs().authid(),
            "px": (round(r.rasterUnitsPerPixelX(), 6), round(r.rasterUnitsPerPixelY(), 6)),
            "origin": (round(e.xMinimum(), 3), round(e.yMaximum(), 3)),
            "size": (r.width(), r.height())}

inputs = {"dem": "/data/model/dem_25m.tif",
          "landcover": "/data/model/landcover_10m.tif",
          "rain": "/data/model/rain_1km.tif"}
for name, path in inputs.items():
    print(f"{name:<10}", grid(path))

Breakdown: Pixel size, origin and size together pin down every pixel's position. The origin uses the top-left corner (minimum x, maximum y), which is how GDAL stores the geotransform. Rounding hides floating-point noise in the printout without changing the comparison. Two rasters with the same CRS and pixel size but origins that differ by a non-multiple of the pixel size are misaligned even though they look identical at a glance.

Choose the target grid

The target grid is an analysis decision. The usual choices are the grid of the most important input, or a grid at the coarsest resolution among the inputs, which avoids inventing detail that coarse data does not contain.

Finest or coarsestResampling coarse data to a fine grid multiplies pixels without adding information and inflates file sizes, while preserving the fine inputs. Resampling fine data to a coarse grid loses detail but keeps every input honest about its resolution. A common compromise is the resolution of the key input, with coarser inputs smoothed by bilinear interpolation.Detail you have, or detail you invent?finest grid (10 m)keeps fine inputs1 km rain → 10,000× pixelsdetail is inventedcoarsest grid (1 km)honest about resolutionfine inputs aggregateddetail is lost

There is no universally right answer: a land-cover analysis at parcel level needs the 10 m grid and must accept that rainfall is constant over each square kilometre; a regional climate model needs the 1 km grid and aggregates land cover into proportions. Decide based on the question and the resolution of the data that matters most, then bring everything onto that grid explicitly.

Resample to the reference grid

gdal:warpreproject reprojects and resamples in one step. Supplying the reference raster's CRS, pixel size and extent produces an output on exactly its grid.

import processing

ref = QgsRasterLayer(inputs["dem"], "dem")
e = ref.extent()
target_extent = f"{e.xMinimum()},{e.xMaximum()},{e.yMinimum()},{e.yMaximum()} [{ref.crs().authid()}]"

def align_to_reference(path, method, out):
    return processing.run("gdal:warpreproject", {
        "INPUT": path,
        "SOURCE_CRS": None,
        "TARGET_CRS": ref.crs(),
        "RESAMPLING": method,
        "NODATA": None,
        "TARGET_RESOLUTION": ref.rasterUnitsPerPixelX(),
        "TARGET_EXTENT": target_extent,
        "TARGET_EXTENT_CRS": ref.crs(),
        "MULTITHREADING": True,
        "OPTIONS": "COMPRESS=DEFLATE|TILED=YES",
        "OUTPUT": out,
    })["OUTPUT"]

rain_25 = align_to_reference(inputs["rain"], 1, "/data/model/rain_25m.tif")          # bilinear
lc_25 = align_to_reference(inputs["landcover"], 6, "/data/model/landcover_25m.tif")   # mode

Breakdown: Target CRS, resolution and extent together define the output grid; when the extent is the reference's own extent and the resolution its pixel size, the output's origin and size match the reference exactly. The resampling enum follows GDAL's order in the Processing dialog — 0 nearest, 1 bilinear, 2 cubic, 5 average, 6 mode on current releases; confirm with processing.algorithmHelp("gdal:warpreproject"). Rain is continuous, so bilinear interpolation gives smooth values; land cover is categorical, so mode picks the most common class among source pixels. Multithreading speeds up large warps.

Pick the resampling method by data type

The single most common resampling error is using an interpolating method on categorical data, which produces class codes that do not exist.

Method by data meaningCategorical data such as land-cover codes must use nearest neighbour when the grid is similar or finer, and mode when aggregating to a coarser grid. Continuous data such as elevation or temperature uses bilinear or cubic when refining and average when aggregating. Counts and totals such as population per pixel must be summed when aggregating, not averaged.What do the values mean?categoricalland cover, soil classfiner: nearestcoarser: modecontinuouselevation, rainfallfiner: bilinear/cubiccoarser: averagecountspopulation per pixelcoarser: sumnever average

from qgis.core import QgsApplication

alg = QgsApplication.processingRegistry().algorithmById("gdal:warpreproject")
options = alg.parameterDefinition("RESAMPLING").options()
METHOD = {}
for i, label in enumerate(options):
    METHOD.setdefault(label.split()[0].lower(), i)      # "Bilinear (2x2 kernel)" → "bilinear"
print(METHOD)

def method_for(kind, refining):
    if kind == "categorical":
        return METHOD["nearest"] if refining else METHOD["mode"]
    if kind == "count":
        return METHOD["sum"]
    return METHOD["bilinear"] if refining else METHOD["average"]

for name, kind in (("rain", "continuous"), ("landcover", "categorical")):
    src = QgsRasterLayer(inputs[name], name)
    refining = src.rasterUnitsPerPixelX() > ref.rasterUnitsPerPixelX()
    print(name, "→", method_for(kind, refining))

Breakdown: Reading the method list from the algorithm's own parameter definition maps names to positions on whatever GDAL version is installed, so the script never depends on a hard-coded enum order. Encoding the rule in a function makes the choice explicit and reviewable rather than a magic number per call. "Refining" means going to smaller pixels. Counts are a special case: population per 1 km pixel aggregated to 5 km must be summed, or totals shrink by a factor of 25; refining counts properly needs redistribution, not resampling. The sum method requires GDAL 3.1 or newer; if it is missing from the printed list, the installed GDAL is too old for count aggregation in one step.

Verify alignment before arithmetic

A check that every aligned raster matches the reference grid exactly turns silent misalignment into an immediate error.

def assert_aligned(paths, reference):
    g0 = grid(reference)
    for p in paths:
        g = grid(p)
        if g != g0:
            raise ValueError(f"{p} not aligned:\n  {g}\n  vs reference {g0}")
    print(len(paths), "rasters aligned with", reference)

assert_aligned([rain_25, lc_25], inputs["dem"])

Breakdown: Comparing the full grid description — CRS, pixel size, origin, size — catches every kind of mismatch, including the half-pixel origin shift that is invisible on the map. Run it at the start of any script that combines rasters; the raster calculator will otherwise happily resample on the fly or produce results offset by half a pixel.

Snap an extent to the grid

When clipping or creating new rasters for a sub-area, the extent must also fall on the grid lattice, or the new raster's origin drifts. Snapping the extent outward to whole pixels keeps it aligned.

import math
from qgis.core import QgsRectangle

def snap_extent(rect, reference):
    px = reference.rasterUnitsPerPixelX()
    py = reference.rasterUnitsPerPixelY()
    ox, oy = reference.extent().xMinimum(), reference.extent().yMaximum()
    xmin = ox + math.floor((rect.xMinimum() - ox) / px) * px
    xmax = ox + math.ceil((rect.xMaximum() - ox) / px) * px
    ymax = oy - math.floor((oy - rect.yMaximum()) / py) * py
    ymin = oy - math.ceil((oy - rect.yMinimum()) / py) * py
    return QgsRectangle(xmin, ymin, xmax, ymax)

area = QgsProject.instance().mapLayersByName("catchment")[0].extent()
print(snap_extent(area, ref).toString(3))

Breakdown: Each edge is moved outward to the nearest pixel boundary of the reference lattice, measured from the reference origin, so the snapped rectangle contains the area and starts and ends exactly on pixel edges. Using it as the target extent for clips and warps of a sub-area produces rasters that align with the full reference and with each other. Add from qgis.core import QgsProject if you run this alone.

Aggregate categories into proportions

When a fine categorical raster is aggregated to a coarse grid, mode keeps only the dominant class and throws away the mix. For many models the mix is exactly what matters: the share of forest in each 1 km cell, not whether forest is the most common class. Computing one proportion raster per class preserves it.

forest_mask = processing.run("gdal:rastercalculator", {
    "INPUT_A": inputs["landcover"], "BAND_A": 1,
    "FORMULA": "(A == 311) + (A == 312) + (A == 313)",
    "RTYPE": 5, "NO_DATA": -9999,
    "OUTPUT": "/data/model/forest_10m.tif"})["OUTPUT"]

forest_share = processing.run("gdal:warpreproject", {
    "INPUT": forest_mask, "TARGET_CRS": ref.crs(),
    "RESAMPLING": METHOD["average"], "TARGET_RESOLUTION": 1000,
    "TARGET_EXTENT": target_extent, "TARGET_EXTENT_CRS": ref.crs(),
    "OUTPUT": "/data/model/forest_share_1km.tif"})["OUTPUT"]

Breakdown: The calculator turns the categorical raster into a 0/1 mask for the forest classes, as floating point so averages are fractional. Averaging that mask onto the 1 km grid gives, in each coarse cell, the fraction of fine pixels that were forest — a value between 0 and 1 that carries far more information than a single dominant class. Repeat per class of interest; the shares across all classes sum to one wherever the fine raster has data.

QGIS version compatibility

gdal:warpreproject with target extent, resolution and multithreading works on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. The resampling enum list depends on the bundled GDAL version; current releases include sum, min, max, median and quartile methods. Always verify enum positions on the version that will run the script.

Troubleshooting

  • Land cover has codes that do not exist. Bilinear or cubic was used on categorical data; use nearest or mode.
  • Population totals shrank. Counts were averaged when aggregating; use sum.
  • Results are offset by half a pixel. Origins differ; align with an explicit target extent.
  • The output is much larger than expected. A coarse raster was resampled to a much finer grid; reconsider the target resolution.

Conclusion

Inspect each input's CRS, pixel size, origin and extent, choose a target grid that fits the question, warp every input onto it with explicit resolution and extent, choose nearest or mode for categories, bilinear or average for continuous values and sum for counts, snap sub-area extents to the lattice, and assert alignment before any arithmetic.

Frequently Asked Questions

Does the QGIS raster calculator align automatically? It resamples inputs to the output grid on the fly using nearest neighbour, which is wrong for continuous data. Align explicitly first.

Can I align without reprojecting? Yes — the same call with the input's own CRS only resamples and shifts to the target lattice.

Is cubic better than bilinear? Smoother, but it can overshoot near sharp edges. Bilinear is the safer default for most continuous data.

How do I align many rasters at once? Loop over them with the same reference, as above, or build a VRT on the target grid.