Rasterizing Vector Polygons onto a Raster Grid

To burn polygons into an array that lines up exactly with an existing raster, take the grid from that raster and pass geometry–value pairs to features.rasterize:

import geopandas as gpd
import rasterio
from rasterio.features import rasterize

with rasterio.open("ndvi_20230615.tif") as src:
    zones = gpd.read_file("fields.gpkg").to_crs(src.crs)
    labels = rasterize(
        ((geom, int(fid)) for geom, fid in zip(zones.geometry, zones["field_id_int"])),
        out_shape=(src.height, src.width),
        transform=src.transform,      # the grid comes from the raster, never from the vectors
        fill=0,
        dtype="int32",
        all_touched=False,
    )

This is the mask-building step behind Zonal Statistics and Vector–Raster Integration, used here to produce a persistent layer rather than a throwaway mask.

Why This Arises in Remote Sensing Workflows

Rasterising is how vector knowledge enters an array pipeline. A field register becomes a label raster so that per-parcel statistics can be computed with array operations instead of geometry loops. A cloud polygon becomes a mask. A training set becomes a target layer for a classifier. In every case the requirement is the same and is easy to state: the output must sit on exactly the same grid as the imagery it will be used with.

“Exactly” is the operative word. A layer that is half a pixel out looks correct at any zoom level a human would inspect, and quietly assigns edge pixels to the wrong parcel. The failure surfaces later as parcels whose statistics include a strip of the neighbouring field, which is indistinguishable from a genuine agronomic signal unless someone checks the geometry.

The defence is to copy the grid rather than to derive one, which is the same principle applied to warping in Reprojecting a Raster from UTM to WGS84.

Environment & Setup

pip install "rasterio>=1.3.0" "geopandas>=0.13" "shapely>=2.0"

Complete Working Example

This function burns an integer identifier per polygon onto the grid of a reference raster and writes the result as a tiled GeoTIFF, so the label layer is itself efficient to read later.

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import rasterize


def rasterize_to_reference(
    vector_path: str,
    reference_raster: str,
    dst_path: str,
    *,
    value_field: str = "field_id_int",
    fill: int = 0,
    all_touched: bool = False,
    dtype: str = "int32",
) -> None:
    """Burn a vector layer's values onto the reference raster's exact grid."""
    with rasterio.open(reference_raster) as ref:
        gdf = gpd.read_file(vector_path)
        if gdf.crs != ref.crs:
            gdf = gdf.to_crs(ref.crs)         # geometries move; the grid never does

        if gdf[value_field].isna().any():
            raise ValueError(f"{value_field} contains nulls; every geometry needs a burn value")

        shapes = ((geom, int(val)) for geom, val in zip(gdf.geometry, gdf[value_field]))

        labels = rasterize(
            shapes,
            out_shape=(ref.height, ref.width),
            transform=ref.transform,
            fill=fill,
            all_touched=all_touched,
            dtype=dtype,
        )

        profile = ref.profile | {
            "dtype": dtype,
            "count": 1,
            "nodata": fill,
            "compress": "deflate",
            "tiled": True,
            "blockxsize": 512,
            "blockysize": 512,
        }

    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(labels, 1)
        dst.build_overviews([2, 4, 8, 16], rasterio.enums.Resampling.nearest)
        dst.update_tags(ns="rio_overview", resampling="nearest")
        dst.update_tags(source_vector=vector_path, all_touched=str(all_touched))


if __name__ == "__main__":
    rasterize_to_reference("fields.gpkg", "ndvi_20230615.tif", "field_labels.tif")
    with rasterio.open("field_labels.tif") as src:
        arr = src.read(1)
        print("distinct labels:", len(np.unique(arr)) - 1, "unlabelled pixels:", int((arr == 0).sum()))

Three choices in that profile are worth naming. nodata=fill means unlabelled pixels are honestly declared rather than silently zero-valued. nearest overviews are mandatory because these are class codes, for the reasons in Adding Internal Overviews with the Right Resampling. And recording all_touched in the tags means the convention travels with the file.

Variant Patterns

1. Overlaps and priority

rasterize burns geometries in sequence, so the last one written wins wherever polygons overlap. That is a usable rule if you control the order.

# Burn in ascending priority so the most important class ends up on top
priority = {"water": 1, "forest": 2, "urban": 3, "fire_scar": 4}
gdf = gdf.assign(prio=gdf["class"].map(priority)).sort_values("prio")

labels = rasterize(
    ((g, int(p)) for g, p in zip(gdf.geometry, gdf["prio"])),
    out_shape=(ref.height, ref.width), transform=ref.transform,
    fill=0, dtype="uint8", all_touched=False,
)

When several classes genuinely coexist in a pixel, a single label layer cannot express it, and the honest representation is one binary layer per class — or a fractional-coverage layer per class, which is the weighting idea developed in the parent topic.

2. Tiled rasterisation for very large extents

A national parcel layer against a continental grid will not fit in memory as one array. Iterate the reference raster’s own blocks and burn only the geometries that intersect each one.

import rasterio
from rasterio.features import rasterize


def rasterize_tiled(gdf, reference_raster: str, dst_path: str, value_field: str) -> None:
    """Burn geometries block by block, so peak memory is one block, not one continent."""
    sindex = gdf.sindex                       # build the spatial index once
    with rasterio.open(reference_raster) as ref:
        profile = ref.profile | {"dtype": "int32", "count": 1, "nodata": 0,
                                 "compress": "deflate", "tiled": True}
        with rasterio.open(dst_path, "w", **profile) as dst:
            for _, window in ref.block_windows(1):
                bounds = rasterio.windows.bounds(window, ref.transform)
                hits = list(sindex.intersection(bounds))
                if not hits:
                    continue                  # nothing here: leave the fill value
                sub = gdf.iloc[hits]
                block = rasterize(
                    ((g, int(v)) for g, v in zip(sub.geometry, sub[value_field])),
                    out_shape=(int(window.height), int(window.width)),
                    transform=ref.window_transform(window),
                    fill=0, dtype="int32", all_touched=False,
                )
                dst.write(block, 1, window=window)

The spatial index is what makes this viable: without it, every block would test every geometry, and the loop would be quadratic.

3. Rasterising lines and points

Lines burn the pixels they cross, and with all_touched=False a diagonal line can leave gaps because consecutive pixels touch only at corners. For connectivity-sensitive work — stream networks, roads — use all_touched=True and accept the slightly thicker result.

Verifying Alignment

Two assertions are enough, and both are one line each:

import rasterio

with rasterio.open("field_labels.tif") as lab, rasterio.open("ndvi_20230615.tif") as ref:
    assert lab.transform == ref.transform, "grids differ — the labels will not line up"
    assert (lab.width, lab.height) == (ref.width, ref.height), "shapes differ"
    assert lab.crs == ref.crs

Then spot-check semantics: pick a coordinate you can identify on a map, read both rasters there with src.index(x, y), and confirm the label is the parcel you expect. Transform equality proves the grids match; the spot check proves the geometries were not shifted before rasterisation, which equality cannot see.

Common Errors

ValueError: Invalid geometry object

The sequence contains a null or empty geometry — common after a spatial join. Drop them with gdf[~gdf.geometry.is_empty & gdf.geometry.notna()] before burning.

Every pixel is the fill value

The geometries are in a different CRS, so they fall outside the grid. Compare gdf.total_bounds with ref.bounds before rasterising.

The output is float64 and enormous

dtype was not passed, so rasterio inferred it from the burn values. Always set the smallest dtype that holds your codes — uint8 for a handful of classes, int32 for identifiers.

Frequently Asked Questions

Q: What happens where two polygons overlap? The later geometry in the sequence wins, because rasterize burns in order. If priority matters, sort the geometries so the most important class is burned last, or rasterize each class into its own layer.

Q: Why is my rasterized layer offset from the imagery? The transform came from the geometries’ bounds rather than from the reference raster. Always copy transform, width and height from the raster you need to align with; deriving a grid from bounds almost never reproduces it exactly.

Q: How do I rasterize a large national dataset without running out of memory? Rasterize per tile. Iterate the reference raster’s block windows, select the geometries intersecting each window with a spatial index, and burn only those into the window-sized array.

Related

• Zonal Statistics and Vector–Raster Integration — the parent topic, where these masks are consumed.
• Computing Zonal Statistics with rasterstats — the library that builds these masks internally.
• How to Clip Rasters to Irregular Polygon Boundaries — the same rasterisation used to cut imagery rather than to label it.
• Writing and Validating Cloud-Optimized GeoTIFFs — writing the label layer so it is cheap to read back.