Burning Polygon Labels into a Training Raster

Turn a labelled polygon layer into a per-pixel training raster with rasterio.features.rasterize: class codes, ignore values, burn order and an area check that proves it worked.

Call rasterio.features.rasterize(shapes, out_shape=..., transform=..., fill=IGNORE, dtype="uint8") with the transform and shape taken from the feature raster, and with the shape list sorted so priority classes burn last. That is the whole operation. Everything that goes wrong goes wrong in the arguments: a transform invented rather than inherited, a fill value that collides with a class code, or an unordered shape list that makes overlaps non-deterministic.

This page is the minimal end-to-end burn. For the wider decisions — burn rules, class priority design, performance at national scale — see Rasterizing Vector Layers for Model Inputs.

The grid comes from the features; only the values come from the vectors A polygon layer supplies geometries and class codes while a feature raster supplies the transform, width, height and coordinate reference system. Both feed a rasterize call whose output is a label raster with exactly the same grid as the feature raster. Never construct the grid — inherit it polygon layer geometry + class_name reprojected to the grid CRS feature raster transform, width, height crs, nodata rasterize() shapes sorted by priority fill = 255 (ignore) all_touched = False label raster same shape, same transform uint8, nodata = 255 one band, deflate + tiled A half-pixel offset here trains every sample against its neighbour’s spectra — and never raises.

Why This Fails in Geospatial ML Pipelines

The failure that costs the most time is an empty output. rasterize does not reproject, does not warn, and does not check that the geometries fall anywhere near the transform’s extent. Hand it a layer in EPSG:4326 and a transform in UTM metres and it returns an array full of the fill value — a perfectly valid raster that trains a model to predict one class everywhere.

The second failure is the fill collision. fill=0 is the default, and 0 is the most natural code for a background class. When they coincide, “nobody labelled this” and “this is background” become the same integer, and the loss function happily learns background across every unlabelled region. Since unlabelled areas are usually the majority of a scene, the model’s apparent accuracy goes up while its usefulness goes down.

Third, burn order is invisible. Overlapping polygons are resolved by list position — last write wins — and if the list order comes from a GeoDataFrame’s row order, an unrelated edit that reorders rows silently changes the labels. Two runs of “the same” pipeline then produce different training sets and different scores, and there is nothing in the diff to explain it. Sorting explicitly is what makes the label raster a function of the data rather than of the file.

Core Principles

  • Take transform, width, height and crs from the feature raster. Never build them from geometry bounds.
  • Reproject the vectors, not the raster. Vector reprojection is exact; raster reprojection resamples.
  • Reserve a dedicated ignore value. 255 in uint8, set as the raster’s nodata, distinct from every class code.
  • Sort the shape list before burning. Broad classes first, priority classes last, from an explicit configuration.
  • Reconcile areas after burning. Rasterized class area should match vector class area within a few percent.
  • Burn once and version the artifact. The label raster is data, not a derived value to recompute per run.

Production-Ready Code

from __future__ import annotations

import logging
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import rasterize
from shapely import make_valid

logger = logging.getLogger(__name__)

CLASS_CODES: dict[str, int] = {
    "background": 0, "cropland": 1, "forest": 2, "water": 3, "built": 4,
}
BURN_PRIORITY = ["background", "cropland", "forest", "built", "water"]  # last wins
IGNORE = 255


def burn_labels(vector_path: str, feature_raster: str, out_path: str,
                class_col: str = "class_name", all_touched: bool = False) -> np.ndarray:
    """Rasterize a labelled polygon layer onto the feature raster's exact grid.

    Args:
        vector_path: Polygon layer carrying a class column.
        feature_raster: The raster whose grid the labels must match.
        out_path: Destination GeoTIFF for the label raster.
        class_col: Column holding class names present in CLASS_CODES.
        all_touched: True only for line-like or sub-pixel geometries.

    Returns:
        The label array that was written.

    Raises:
        ValueError: on unmapped classes, a missing CRS, or an all-fill result.
    """
    with rasterio.open(feature_raster) as src:
        transform, width, height, crs = src.transform, src.width, src.height, src.crs
        profile = src.profile.copy()

    gdf = gpd.read_file(vector_path)
    if gdf.crs is None:
        raise ValueError(f"{vector_path}: no CRS — cannot align to the feature grid")
    if gdf.crs != crs:
        logger.info("reprojecting labels %s -> %s", gdf.crs.to_string(), crs.to_string())
        gdf = gdf.to_crs(crs)

    # Repair, then keep only polygonal parts (make_valid can emit collections).
    invalid = ~gdf.geometry.is_valid
    if invalid.any():
        logger.info("repairing %d invalid geometries", int(invalid.sum()))
        gdf.loc[invalid, "geometry"] = gdf.loc[invalid, "geometry"].apply(make_valid)
    gdf = gdf[gdf.geometry.geom_type.isin(["Polygon", "MultiPolygon"])]
    gdf = gdf[~gdf.geometry.is_empty]

    unknown = set(gdf[class_col]) - set(CLASS_CODES)
    if unknown:
        raise ValueError(f"unmapped class labels: {sorted(unknown)}")

    rank = {name: i for i, name in enumerate(BURN_PRIORITY)}
    ordered = gdf.assign(_rank=gdf[class_col].map(rank)).sort_values("_rank", kind="stable")
    shapes = [(geom, CLASS_CODES[name])
              for geom, name in zip(ordered.geometry, ordered[class_col])]

    arr = rasterize(shapes, out_shape=(height, width), transform=transform,
                    fill=IGNORE, all_touched=all_touched, dtype="uint8")

    if (arr == IGNORE).all():
        raise ValueError(
            "every pixel is the fill value — the geometries do not intersect the grid. "
            f"labels bounds {tuple(gdf.total_bounds)} vs raster bounds "
            f"{tuple(rasterio.open(feature_raster).bounds)}"
        )

    profile.update(count=1, dtype="uint8", nodata=IGNORE, compress="deflate",
                   predictor=2, tiled=True, blockxsize=512, blockysize=512)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(arr, 1)
        dst.set_band_description(1, "label")
        dst.update_tags(classes=",".join(f"{k}={v}" for k, v in CLASS_CODES.items()),
                        ignore_value=str(IGNORE), all_touched=str(all_touched))
    return arr

Step-by-Step Walkthrough

Step 1 — Burn and read the class histogram

labels = burn_labels("field_labels.gpkg", "features.tif", "labels.tif")

values, counts = np.unique(labels, return_counts=True)
inv = {v: k for k, v in CLASS_CODES.items()}
for v, c in zip(values, counts):
    name = "IGNORE" if v == IGNORE else inv.get(int(v), f"?{v}")
    print(f"{name:<12} {c:>12,}  {c / labels.size:6.2%}")

The histogram is the first diagnostic. A run where IGNORE is 99% means the geometries barely intersect the grid; a run where one class is 0% means it lost every overlap or its polygons are sub-pixel.

Step 2 — Confirm the grids are identical

def assert_grid_identical(a_path: str, b_path: str) -> None:
    with rasterio.open(a_path) as a, rasterio.open(b_path) as b:
        assert a.crs == b.crs, f"CRS differs: {a.crs} vs {b.crs}"
        assert (a.width, a.height) == (b.width, b.height), "shape differs"
        assert a.transform.almost_equals(b.transform, precision=1e-6), "transform differs"


assert_grid_identical("features.tif", "labels.tif")

Step 3 — Reconcile area, class by class

gdf = gpd.read_file("field_labels.gpkg").to_crs(rasterio.open("features.tif").crs)
px_area = abs(rasterio.open("features.tif").transform.a *
              rasterio.open("features.tif").transform.e)

for name, code in CLASS_CODES.items():
    if code == 0:
        continue
    vec = gdf.loc[gdf["class_name"] == name].geometry.area.sum()
    ras = float((labels == code).sum()) * px_area
    if vec:
        print(f"{name:<10} vector {vec:>12,.0f} m2  raster {ras:>12,.0f} m2  "
              f"{(ras - vec) / vec:+6.1%}")

Uniform over-estimate across all classes means all_touched=True; a deficit concentrated in one class means it is losing overlaps and needs to move later in BURN_PRIORITY.

Step 4 — Add linear features in a second pass

roads = gpd.read_file("roads.gpkg").to_crs(rasterio.open("features.tif").crs)
with rasterio.open("features.tif") as src:
    road_arr = rasterize([(g, CLASS_CODES["built"]) for g in roads.geometry],
                         out_shape=(src.height, src.width), transform=src.transform,
                         fill=0, all_touched=True, dtype="uint8")
labels = np.where(road_arr > 0, road_arr, labels)

A one-pixel-wide road would vanish under the default centre rule, so it needs all_touched=True — but only for that layer, burned on top of the polygon pass.

fill=0 merges “unlabelled” with “background” Two label rasters are compared. In the first, unlabelled pixels take the value zero, which is also the background class code, so the model is trained to predict background across the entire undigitised area. In the second, unlabelled pixels take the dedicated ignore value 255, which the loss function skips, so only genuinely labelled pixels contribute gradients. fill = 0 fill = 255 (ignore) 0 0 1 1 0 2 1 0 0 0 0 0 Digitised background and never-digitised ground are the same integer. Loss learns “background” everywhere. 255 255 1 1 255 2 1 0 255 255 255 0 Only the digitised zero on the right is a real background label. ignore_index=255 skips the rest. Set the ignore value as the raster’s nodata so every reader inherits the distinction.

Verification

import numpy as np
import rasterio
from rasterio.transform import from_origin
from rasterio.features import rasterize
from shapely.geometry import box


def test_single_square_burns_exact_pixel_count():
    """A 30 m square on a 10 m grid must claim exactly 9 pixels."""
    transform = from_origin(0, 100, 10, 10)      # 10 m pixels, origin top-left
    square = box(20, 40, 50, 70)                 # 30 m x 30 m
    arr = rasterize([(square, 1)], out_shape=(10, 10), transform=transform,
                    fill=255, dtype="uint8")
    assert int((arr == 1).sum()) == 9


def test_last_shape_wins_the_overlap():
    transform = from_origin(0, 100, 10, 10)
    a, b = box(0, 60, 40, 100), box(20, 60, 60, 100)
    arr = rasterize([(a, 1), (b, 2)], out_shape=(10, 10), transform=transform,
                    fill=255, dtype="uint8")
    assert arr[0, 3] == 2, "the shape burned last must own the shared pixels"


def test_empty_result_raises(tmp_path):
    """A CRS mismatch must fail loudly rather than return a raster of fill."""
    try:
        burn_labels("labels_wgs84.gpkg", "features_utm.tif", str(tmp_path / "l.tif"))
    except ValueError as exc:
        assert "do not intersect" in str(exc)
    else:
        raise AssertionError("expected a ValueError on an all-fill result")

The first test is the one to keep: it pins the burn rule, the transform convention and the dtype in six lines, and it fails the moment somebody flips all_touched in a shared configuration.

Area reconciliation tells you which failure you have For four land cover classes the vector area and the rasterized area are drawn as paired bars. Cropland and forest agree closely. Water's rasterized bar is visibly shorter than its vector bar, indicating lost overlaps. Built's rasterized bar is visibly longer, indicating the all_touched rule. Rasterized area vs vector area, by class area cropland −1% forest −2% water −53% lost overlaps built +71% all_touched vector rasterized

FAQ

What value should unlabelled pixels get?

A dedicated ignore value that no class uses — 255 in a uint8 raster — set as the file’s nodata. Unlabelled and background are different things: background is a class to learn, ignore pixels must be excluded from the loss entirely.

Why do class pixel counts not match polygon areas?

Three causes: all_touched=True inflates every polygon by about half a pixel around its perimeter; overlapping polygons hand shared pixels to whichever class burned last; and sub-pixel polygons vanish under the default rule. The reconciliation in Step 3 distinguishes them by the sign and the concentration of the error.

Should the label raster be uint8 or int32?

uint8 for semantic classes — 255 values, excellent compression. int32 for instance labels, where each field or building carries its own id; uint8 will wrap silently past 255 and merge unrelated instances.

Part of: Rasterizing Vector Layers for Model Inputs Part of: Spatial Feature Engineering for Machine Learning