Most geospatial training data arrives in two incompatible shapes. The predictors are rasters — reflectance bands, spectral indices, elevation derivatives — laid out on a regular grid. The labels are vectors: field boundaries digitised by an agronomist, building footprints from a cadastre, flood extents traced from an aerial survey. Before a per-pixel model can be trained, the vectors have to be burned onto exactly the same grid the predictors live on. That is rasterization, and getting it wrong is one of the quietest ways to poison a training set.
This topic is part of Spatial Feature Engineering for Machine Learning. Rasterization is the mirror image of zonal statistics for polygon aggregation: zonal statistics move information from a grid up into polygons for tabular models, while rasterization pushes information from polygons down onto a grid for per-pixel models such as the ones in deep learning segmentation for satellite imagery. Both directions depend on the same prerequisite — a single agreed grid — and both fail silently when that agreement breaks.
Problem Framing
Rasterization has exactly one hard requirement and several soft ones, and teams routinely get the hard one wrong. The hard requirement is that the label raster and the feature raster share a grid: identical CRS, identical affine transform, identical width and height. Not “roughly the same extent” — identical. If the label grid is offset by half a pixel, every sample is trained against the spectral signature of its neighbour, and the model still converges, just to a worse optimum than it should. Nothing in the stack raises an error.
The soft requirements are all judgement calls that change what the model learns:
- Burn rule. The default rule fills a pixel only when the pixel centre falls inside the geometry. This preserves area faithfully but erases anything narrower than one pixel — every road, every stream, every hedgerow.
all_touched=Truefills any pixel the geometry intersects, which rescues linear features at the cost of inflating polygon areas by roughly half a pixel around each perimeter. - Overlap priority.
rasterio.features.rasterizewrites shapes in sequence, so the last one wins. When a water polygon and a wetland polygon overlap, which class survives is decided purely by list order — an implicit decision that should be made explicit. - Fill value. Pixels touched by nothing take the
fillvalue. Iffill=0and one of your classes is also coded0, background and that class become indistinguishable, and the confusion only shows up as an oddly high recall for class zero. - Dtype.
uint8is the natural choice for a label raster, but it caps you at 255 classes and silently wraps if a class id exceeds that. For instance ids — one integer per building — you needint32.
A related decision is what to do with the pixels no polygon covers. In a segmentation problem, “unlabelled” and “background” are different: unlabelled pixels should be excluded from the loss, background pixels should be learned. That distinction has to survive into the label raster, usually as a dedicated ignore value.
Prerequisites & Environment Setup
# Pinned requirements for vector rasterization
rasterio==1.3.10
geopandas==0.14.4
shapely==2.0.5
numpy==1.26.4
pyproj==3.6.1
Install with:
pip install "rasterio==1.3.10" "geopandas==0.14.4" "shapely==2.0.5" \
"numpy==1.26.4" "pyproj==3.6.1"GDAL/PROJ system dependencies (Ubuntu/Debian):
sudo apt-get install -y gdal-bin libgdal-dev libproj-devShapely 2.x matters here. Its make_valid is a first-class function rather than the old buffer(0) trick, and buffer(0) on a self-intersecting polygon can silently delete the smaller ring — a whole labelled field disappearing from the training set with no warning. Geometry repair is covered end-to-end in fixing invalid geometries before buffering.
Step-by-Step Implementation
Step 1 — Adopt the reference grid from the feature raster
Never construct a grid for the labels. Read it from the raster the model will actually consume, so alignment is true by construction rather than by coincidence.
import rasterio
import geopandas as gpd
def load_reference_grid(raster_path: str) -> dict:
"""Return the grid definition the label raster must match exactly."""
with rasterio.open(raster_path) as src:
grid = {
"transform": src.transform,
"width": src.width,
"height": src.height,
"crs": src.crs,
"bounds": src.bounds,
}
if grid["crs"] is None:
raise ValueError("Reference raster has no CRS")
return grid
grid = load_reference_grid("terrain_features.tif")
labels_gdf = gpd.read_file("field_labels.gpkg")
if labels_gdf.crs != grid["crs"]:
labels_gdf = labels_gdf.to_crs(grid["crs"])
assert labels_gdf.crs == grid["crs"], "label CRS still does not match the grid"Step 2 — Repair geometries and drop the unusable ones
rasterize raises on some invalid geometries and silently mis-renders others. Repair first, and log what you dropped rather than dropping it quietly.
from shapely import make_valid
def clean_geometries(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
"""Repair invalid geometries; drop empty and null ones with a count."""
before = len(gdf)
gdf = gdf[~gdf.geometry.isna()].copy()
gdf = gdf[~gdf.geometry.is_empty].copy()
invalid = ~gdf.geometry.is_valid
if invalid.any():
gdf.loc[invalid, "geometry"] = gdf.loc[invalid, "geometry"].apply(make_valid)
gdf = gdf[gdf.geometry.geom_type.isin(["Polygon", "MultiPolygon"])].copy()
print(f"cleaned {before} -> {len(gdf)} geometries ({int(invalid.sum())} repaired)")
assert gdf.geometry.is_valid.all(), "geometries still invalid after repair"
return gdf
labels_gdf = clean_geometries(labels_gdf)make_valid can turn a self-intersecting polygon into a GeometryCollection containing stray lines and points. Filtering to polygon types afterwards is what keeps those fragments out of the burn.
Step 3 — Encode classes and set the burn order
Class priority is a modelling decision, so state it as data. Rarer and more specific classes are burned last so they survive overlaps with broad background classes.
CLASS_CODES = {"background": 0, "cropland": 1, "forest": 2, "water": 3, "built": 4}
IGNORE = 255 # excluded from the loss, never predicted
BURN_PRIORITY = ["cropland", "forest", "built", "water"] # last wins
def build_shapes(gdf: gpd.GeoDataFrame, class_col: str = "class_name"):
"""Yield (geometry, burn value) pairs ordered so priority classes are burned last."""
unknown = set(gdf[class_col]) - set(CLASS_CODES)
if unknown:
raise ValueError(f"Unmapped class labels: {sorted(unknown)}")
order = {name: i for i, name in enumerate(BURN_PRIORITY)}
ordered = gdf.assign(_rank=gdf[class_col].map(order)).sort_values("_rank")
return [(geom, CLASS_CODES[name])
for geom, name in zip(ordered.geometry, ordered[class_col])]
shapes = build_shapes(labels_gdf)
assert len(shapes) == len(labels_gdf)Step 4 — Burn the shapes
Every argument here is deliberate. Passing the reference transform and shape is what guarantees alignment; fill=IGNORE keeps unlabelled ground distinguishable from the background class.
import numpy as np
from rasterio.features import rasterize
def rasterize_labels(shapes, grid: dict, fill: int = IGNORE,
all_touched: bool = False) -> np.ndarray:
"""Burn (geometry, value) pairs onto the reference grid."""
arr = rasterize(
shapes=shapes,
out_shape=(grid["height"], grid["width"]),
transform=grid["transform"],
fill=fill,
all_touched=all_touched,
dtype="uint8",
)
if (arr == fill).all():
raise ValueError(
"Every pixel is fill — the geometries do not intersect the grid. "
"Check the CRS and the raster bounds."
)
return arr
label_arr = rasterize_labels(shapes, grid)
print({int(v): int(c) for v, c in zip(*np.unique(label_arr, return_counts=True))})Linear features need a second pass with all_touched=True, burned on top of the polygon pass so a road crossing a field is not erased by it.
roads = gpd.read_file("roads.gpkg").to_crs(grid["crs"])
road_arr = rasterize_labels([(g, CLASS_CODES["built"]) for g in roads.geometry],
grid, fill=0, all_touched=True)
label_arr = np.where(road_arr > 0, road_arr, label_arr)Step 5 — Write the label raster with the reference profile
def write_labels(path: str, arr: np.ndarray, grid: dict, nodata: int = IGNORE) -> None:
profile = {
"driver": "GTiff", "dtype": "uint8", "count": 1, "nodata": nodata,
"width": grid["width"], "height": grid["height"],
"transform": grid["transform"], "crs": grid["crs"],
"compress": "deflate", "tiled": True, "blockxsize": 512, "blockysize": 512,
}
with rasterio.open(path, "w", **profile) as dst:
dst.write(arr, 1)
dst.set_band_description(1, "land_cover_class")
write_labels("labels.tif", label_arr, grid)Verification & Testing
Two checks catch nearly every rasterization defect: a grid identity check and an area reconciliation.
import numpy as np
import rasterio
def assert_grids_identical(feature_path: str, label_path: str) -> None:
"""Fail loudly if the label raster is not pixel-identical to the feature grid."""
with rasterio.open(feature_path) as f, rasterio.open(label_path) as l:
assert f.crs == l.crs, f"CRS differs: {f.crs} vs {l.crs}"
assert (f.width, f.height) == (l.width, l.height), "shape differs"
assert f.transform.almost_equals(l.transform, precision=1e-6), (
f"transform differs:\n{f.transform}\n{l.transform}"
)
def reconcile_class_areas(gdf, label_arr, grid, class_col="class_name", tol=0.05):
"""Compare rasterized class area against the vector area, class by class."""
px_area = abs(grid["transform"].a * grid["transform"].e)
for name, code in CLASS_CODES.items():
if code == 0:
continue
vector_area = gdf.loc[gdf[class_col] == name].geometry.area.sum()
raster_area = float((label_arr == code).sum()) * px_area
if vector_area == 0:
continue
err = abs(raster_area - vector_area) / vector_area
print(f"{name:<10} vector {vector_area:,.0f} m2 raster {raster_area:,.0f} m2 "
f"({err:.1%})")
assert err < tol, f"{name}: rasterized area is {err:.1%} off the vector area"
assert_grids_identical("terrain_features.tif", "labels.tif")
reconcile_class_areas(labels_gdf, label_arr, grid)An area error consistently in the same direction points at the burn rule: a uniform over-estimate is the signature of all_touched=True, a uniform under-estimate means small polygons are being dropped. An error concentrated in one class points at burn order.
Troubleshooting & Common Errors
ValueError: Invalid geometry object — a null, empty or non-polygonal geometry reached rasterize. Run the cleaning step and filter GeometryCollection outputs from make_valid down to polygon members.
The label raster is entirely the fill value — the geometries do not intersect the transform’s extent. Ninety-five percent of the time this is a CRS mismatch; the rest of the time the reference raster covers a different tile than the labels. Print labels_gdf.total_bounds next to grid["bounds"] and compare.
rasterio.errors.WindowError or a shape mismatch when stacking with features — out_shape was derived from the geometry bounds rather than from the reference raster. Always take (height, width) from the feature raster.
Class counts change between runs on unchanged data — the shape list order is coming from an unordered source such as a set or a dict iteration, so overlap resolution is non-deterministic. Sort explicitly, as in Step 3.
Small polygons disappear entirely — expected under the default centre-in-polygon rule when a polygon is smaller than a pixel. Either raise the grid resolution, switch that layer to all_touched=True, or model those features as points sampled from a coarser grid.
Every class is off by one row or column — the transform came from a raster that was later windowed or resampled without updating the transform. Recompute it with rasterio.windows.transform(window, src.transform) rather than reusing the parent transform.
Performance Optimisation
rasterize is a C loop over geometry vertices, so its cost tracks total vertex count far more than polygon count. Three levers matter.
Simplify before burning, not after. A cadastre digitised at 1:1000 carries vertices spaced under a metre; on a 10 m grid they cannot possibly change which pixels are filled. gdf.geometry.simplify(tolerance=cell_size / 4, preserve_topology=True) typically removes 60–80% of vertices with no change in the output raster — verify that claim with the area reconciliation above before trusting it on a new dataset.
Rasterize per tile, not per country. Allocating a national uint8 array is cheap, but the intermediate float work is not. Loop over the feature raster’s block_windows(), subset the GeoDataFrame with its spatial index (gdf.sindex.query(box(*window_bounds))), and burn only the geometries that intersect that window. Indexing technique is covered in speeding up nearest neighbour joins with a spatial index.
Write tiled and compressed. A label raster is mostly runs of identical integers, so compress="deflate" with predictor=2 routinely cuts it to a few percent of its raw size, and tiled layout means the training loader can read a patch without decompressing a full strip.
Finally, rasterize once. The label raster is a data artifact, not a derived value to recompute each run. Store it, hash it, and version it alongside the feature stack so an edit to the source vector layer shows up as a deliberate new dataset version rather than as unexplained variance between two supposedly identical experiments.
FAQ
Why is my rasterized label layer empty?
Either the geometries fall outside the grid — almost always a CRS mismatch — or the burn value equals the fill value so the output looks untouched. Reproject the vector layer to the reference raster’s CRS first, then compare gdf.total_bounds against the raster bounds, and pick a fill value no class uses.
When should all_touched be True?
For linear features (roads, rivers, pipelines) and for polygons smaller than a pixel, where the default centre rule would delete them. For area classes leave it False, because all_touched adds roughly half a pixel around every perimeter and systematically inflates class area — which shows up downstream as a shifted class balance.
How do overlapping polygons get resolved?
By list order alone: rasterize burns sequentially and the last shape covering a pixel wins. Make that explicit by sorting the (geometry, value) pairs so the class that should survive is last. See handling class imbalance in land cover classification for what these ordering choices do to the class distribution the model sees.
Should labels be rasterized once or per training run?
Once. Re-rasterizing per run couples every experiment to the live state of the vector source, so one edited polygon silently changes the training set and two runs stop being comparable. Burn once, version the artifact, and treat a new burn as a new dataset version.
Related
- Zonal Statistics for Polygon Aggregation — the reverse operation, grid up into polygons
- Burning Polygon Labels into a Training Raster — the minimal end-to-end label burn
- Aligning Rasterized Labels to an Existing Feature Grid — diagnosing and fixing half-pixel offsets
- Fixing Invalid Geometries Before Buffering — geometry repair that rasterization also depends on
- Deep Learning Segmentation for Satellite Imagery — the models that consume these label rasters