Rasterize a Vector Layer in PyQGIS

Some analyses work only on grids. A cost surface for least-cost paths needs roads, rivers and land cover as raster values; a suitability model multiplies raster layers together; a mask for clipping imagery must match the imagery's pixels; a habitat model takes protected areas as a 0/1 grid. Rasterizing burns vector features into pixels, and the result is only as good as three decisions: which grid the pixels follow, which value each pixel gets, and which pixels count as "covered" by a feature.

This recipe belongs to Raster Analysis Workflows. It rasterizes polygons with values from a field, aligns the output to an existing raster, chooses between all-touched and centre-pixel rules, sets NoData and data types, and burns several layers into one grid.

Three decisions behind every rasterizationRasterizing needs a grid, defined by extent, pixel size and CRS, ideally copied from an existing raster; a value per feature, from a field or a constant; and a coverage rule deciding which pixels a feature burns, either only pixels whose centre lies inside or every pixel the feature touches. Pixels no feature covers get the NoData value.Grid, value, coverage rulegridextentpixel sizeCRScopy from a rastervaluefrom a fieldor a constantdata type to fitcoveragepixel centre insideor all touchedthin features!uncovered pixels get NoData

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series, with the GDAL Processing provider.
  • A vector layer in a projected CRS, and ideally a reference raster whose grid the output should match.

Rasterize with values from a field

The simplest call burns a numeric field into a new raster with a chosen pixel size.

import processing
from qgis.core import QgsProject, QgsRasterLayer

landuse = QgsProject.instance().mapLayersByName("landuse")[0]
out = processing.run("gdal:rasterize", {
    "INPUT": landuse,
    "FIELD": "cost",                # friction per land-use class
    "BURN": 0,
    "USE_Z": False,
    "UNITS": 1,                     # 1 = georeferenced units (pixel size)
    "WIDTH": 10, "HEIGHT": 10,      # 10 m pixels
    "EXTENT": landuse.extent(),
    "NODATA": -9999,
    "DATA_TYPE": 5,                 # Float32
    "INIT": None,
    "INVERT": False,
    "OPTIONS": "COMPRESS=DEFLATE|TILED=YES",
    "OUTPUT": "/data/model/landuse_cost.tif",
})["OUTPUT"]

cost = QgsRasterLayer(out, "land-use cost")
print(cost.width(), "x", cost.height())

Breakdown: With UNITS set to georeferenced units, WIDTH and HEIGHT are pixel sizes in map units rather than pixel counts. FIELD gives each feature's value; leave it empty and set BURN to burn a constant instead. NODATA marks pixels no feature covers. The data type must hold the values: Float32 for fractional costs, a small integer type for class codes, Byte for 0/1 masks. The extent here is the layer's own, which is fine for a stand-alone raster but rarely what an analysis needs — the next section aligns it.

Align to an existing grid

Rasters combined in an analysis must share a grid: same CRS, same pixel size, same origin, same extent. Rasterizing onto the grid of a reference raster avoids resampling later.

Same grid or notLeft: a rasterized layer with its own extent and origin is offset from the DEM grid by a fraction of a pixel, so pixel-by-pixel arithmetic needs resampling and edges blur. Right: rasterized with the DEM's extent and pixel size, every pixel lines up with a DEM pixel and raster calculator expressions work directly.Copy the reference grid exactlyown extent: offsetreference grid: aligned

dem = QgsRasterLayer("/data/dem/dem_10m.tif", "dem")
ext = dem.extent()
aligned = processing.run("gdal:rasterize", {
    "INPUT": landuse, "FIELD": "cost", "UNITS": 1,
    "WIDTH": dem.rasterUnitsPerPixelX(), "HEIGHT": dem.rasterUnitsPerPixelY(),
    "EXTENT": f"{ext.xMinimum()},{ext.xMaximum()},{ext.yMinimum()},{ext.yMaximum()} [{dem.crs().authid()}]",
    "NODATA": -9999, "DATA_TYPE": 5, "OPTIONS": "COMPRESS=DEFLATE|TILED=YES",
    "OUTPUT": "/data/model/landuse_cost_aligned.tif"})["OUTPUT"]

check = QgsRasterLayer(aligned, "check")
assert (check.width(), check.height()) == (dem.width(), dem.height())
assert check.extent() == dem.extent()
print("aligned to DEM grid")

Breakdown: Passing the DEM's pixel size and its exact extent — as a string with the CRS — makes the output grid identical to the DEM's, so the two can be combined in the raster calculator without any resampling. The assertions confirm it. The vector layer must be in the same CRS as the reference raster, or features are burned in the wrong place; reproject first if needed. When the reference grid matters across many rasters, resampling and aligning rasters covers bringing existing rasters onto it too.

All touched or pixel centre

By default a pixel is burned only if its centre lies inside a polygon (or on a line). Narrow features — roads, rivers, hedges — can then vanish or break into dashes where they pass between pixel centres. The all-touched rule burns every pixel a feature touches.

Centre rule versus all-touchedA thin diagonal road crosses a grid of pixels. With the default centre rule, only pixels whose centres the line passes close enough to are burned, producing a broken line of isolated pixels. With all-touched, every pixel the line crosses is burned, producing a continuous but thicker line. Polygons grow by up to one pixel at their edges with all-touched.Thin features need all-touchedcentre rule: brokenall touched: continuous

roads = QgsProject.instance().mapLayersByName("roads")[0]
road_mask = processing.run("gdal:rasterize", {
    "INPUT": roads, "FIELD": None, "BURN": 1, "UNITS": 1,
    "WIDTH": dem.rasterUnitsPerPixelX(), "HEIGHT": dem.rasterUnitsPerPixelY(),
    "EXTENT": dem.extent(), "NODATA": 0, "DATA_TYPE": 0,        # Byte
    "INIT": 0, "EXTRA": "-at",
    "OUTPUT": "/data/model/roads_mask.tif"})["OUTPUT"]

Breakdown: The -at flag in EXTRA switches GDAL to all-touched mode. Lines become continuous chains of pixels, which matters for connectivity in cost-distance and network-on-raster analyses. The trade-off is that polygons grow by up to one pixel along their edges, so areas computed from an all-touched raster overestimate; use the centre rule for area-faithful polygon rasters and all-touched for thin features and masks that must not miss anything. INIT fills the raster with 0 before burning, so uncovered pixels are an explicit "no road" rather than NoData — useful for masks.

Burn several layers into one grid

A cost surface usually combines several sources — land use as the base, roads as cheap corridors on top, rivers as barriers. Rasterizing into an existing raster burns each layer over the previous ones.

import shutil

shutil.copy(aligned, "/data/model/cost_surface.tif")
for layer_name, value in (("roads", 1.0), ("rivers", 9999.0)):
    lyr = QgsProject.instance().mapLayersByName(layer_name)[0]
    processing.run("gdal:rasterize_over_fixed_value", {
        "INPUT": lyr, "INPUT_RASTER": "/data/model/cost_surface.tif",
        "BURN": value, "ADD": False, "EXTRA": "-at"})
print("cost surface: land use base, roads and rivers burned over it")

Breakdown: gdal:rasterize_over_fixed_value writes into an existing raster in place, so the land-use costs remain wherever no road or river is burned. Order matters: later layers overwrite earlier ones, so rivers, burned last, act as barriers even where a road crosses them — burn bridges after rivers if they should be passable. Copying the aligned base first keeps the original intact. The result feeds least-cost path tools in GRASS or SAGA, run from Python as in running GRASS and SAGA algorithms.

Combine rasterized layers in a model

Rasterizing is usually a step towards a model in which several layers are combined pixel by pixel. Because the rasters share one grid, the raster calculator can combine them directly — here a simple suitability score that excludes protected areas, rewards proximity to roads and penalises steep slopes.

protected = processing.run("gdal:rasterize", {
    "INPUT": QgsProject.instance().mapLayersByName("protected_areas")[0],
    "BURN": 1, "UNITS": 1, "INIT": 0, "NODATA": None, "DATA_TYPE": 0,
    "WIDTH": dem.rasterUnitsPerPixelX(), "HEIGHT": dem.rasterUnitsPerPixelY(),
    "EXTENT": dem.extent(), "OUTPUT": "/data/model/protected.tif"})["OUTPUT"]

score = processing.run("gdal:rastercalculator", {
    "INPUT_A": protected, "BAND_A": 1,
    "INPUT_B": "/data/model/slope_deg.tif", "BAND_B": 1,
    "INPUT_C": road_mask, "BAND_C": 1,
    "FORMULA": "(A == 0) * (B < 15) * (1 + C)",
    "RTYPE": 0, "NO_DATA": 255,
    "OUTPUT": "/data/model/suitability.tif"})["OUTPUT"]

Breakdown: Rasterizing protected areas with INIT 0 and no NoData gives a clean 0/1 grid, so the formula can test it without NoData propagating. The expression is a product of conditions: outside protected land, slope under 15 degrees, with a bonus where a road passes. Because every input is on the DEM's grid — slope derived from the DEM, the masks rasterized onto it — no resampling happens, and each output pixel is computed from exactly corresponding input pixels. More elaborate models follow the same pattern; the raster calculator recipe covers expressions, NoData handling and the QGIS-native calculator.

Check what was burned

Rasterization can silently drop features — too thin for the centre rule, outside the extent, in the wrong CRS. Comparing feature counts with burned values catches it.

import numpy as np

prov = QgsRasterLayer(road_mask, "m").dataProvider()
block = prov.block(1, prov.extent(), prov.xSize(), prov.ySize())
arr = np.frombuffer(bytes(block.data()), dtype=np.uint8).reshape(prov.ySize(), prov.xSize())
burned_km = arr.sum() * dem.rasterUnitsPerPixelX() / 1000
road_km = sum(f.geometry().length() for f in roads.getFeatures()) / 1000
print(f"roads {road_km:.0f} km, burned ≈ {burned_km:.0f} km of pixel length")

Breakdown: For a line mask, the number of burned pixels times the pixel size approximates the burned length; a figure far below the true road length means features fell outside the extent or were lost to the centre rule. For polygon rasters, compare burned area per value with polygon areas, as in the polygonize recipe's area check run in reverse. Reading the whole raster into NumPy is fine for masks of moderate size; tile larger ones.

QGIS version compatibility

gdal:rasterize and gdal:rasterize_over_fixed_value work on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. The EXTRA parameter for additional GDAL switches such as -at has been available since 3.8; on older releases all-touched is not exposed. Data type enums follow GDAL's list and are the same across these releases.

Troubleshooting

  • The raster is empty. The vector and extent are in different CRSs, or UNITS was 0 so width and height were pixel counts.
  • Roads appear as dashes. Use all-touched (-at).
  • Values are truncated. The data type is integer for fractional values; use Float32.
  • Grids do not line up with the DEM. The extent or pixel size was not copied exactly.

Conclusion

Rasterize onto the exact grid of a reference raster, take values from a field with a data type that fits, use the centre rule for area-faithful polygons and all-touched for thin features, initialise masks explicitly, burn layers in a deliberate order, and check burned lengths or areas against the vectors.

Frequently Asked Questions

Can I rasterize point data? Yes — each point burns the pixel it falls in; for densities use kernel density instead.

How do I burn the count of overlapping features? Use ADD in gdal:rasterize_over_fixed_value with a burn value of 1 to accumulate.

Can I rasterize Z values from 3D geometry? Set USE_Z to burn the geometry's Z instead of a field.

What about fractional coverage per pixel? GDAL's rasterize is binary per pixel; for coverage fractions, use zonal tools such as exactextract.