Aligning Rasterized Labels to an Existing Feature Grid

Diagnose and fix half-pixel offsets between a label raster and a feature stack: compare affine transforms, detect sub-pixel shifts, and rebuild labels on the reference grid.

Compare src.transform for the feature raster and the label raster term by term. If the c (x origin) or f (y origin) terms differ by anything other than an exact multiple of the pixel size, the two grids are offset and every boundary pixel in your training set is mislabelled. The fix is not to shift the array — it is to re-burn the labels from the original vector geometries onto the feature raster’s transform.

This page covers detection and repair. The burn itself is in burning polygon labels into a training raster, and the surrounding design decisions are in Rasterizing Vector Layers for Model Inputs.

A half-pixel offset makes every label straddle four feature pixels A feature raster grid is drawn with solid lines and a label raster grid is drawn over it with dashed lines, offset by half a cell in both directions. Each label cell therefore overlaps four feature cells, so no sample is trained against exactly the pixel it describes. Same extent, same resolution, wrong origin label origin feature origin solid = features dashed = labels offset = 30 m on a 60 m grid = exactly half a pixel Every label cell now covers parts of four different feature cells.

Why This Fails in Geospatial ML Pipelines

Offsets are born in small conveniences. Someone computes the label grid from gdf.total_bounds instead of reading it from the features, and the origin lands wherever the northernmost vertex happens to be. Someone windows a raster with rasterio.windows.Window and reuses the parent transform instead of rasterio.windows.transform(window, src.transform). Someone reprojects the labels separately and lets GDAL pick the output grid. All three produce files with the right CRS, the right resolution and the right approximate extent — which is exactly what a quick visual check confirms.

The damage is concentrated at class boundaries. Interior pixels of a large field are still labelled correctly, so overall accuracy hardly moves. What degrades is the model’s ability to place edges: it is trained on boundary pixels whose spectra belong half to one class and half to another, so it learns a smeared decision surface. In segmentation this shows up as rounded, bloated shapes; in per-pixel classification it shows up as a fringe of misclassification around every polygon that no amount of hyperparameter tuning removes. If you have ever tuned a model for a week and never beaten a boundary artefact, this is the first thing to check.

Sub-pixel offsets are also invisible to most tooling. gdalinfo reports the transform but nobody compares two of them digit by digit; QGIS overlays the layers and a half-pixel shift is a fraction of a screen pixel at any sensible zoom. Only an explicit numerical assertion catches it, which is why that assertion belongs in the pipeline rather than in a checklist.

Core Principles

  • Grids are equal or they are not. Compare all six affine terms, plus width, height and CRS. “Close enough” is not a category.
  • Re-burn, do not shift. Rasterizing again from the source geometries onto the reference transform is exact; translating an array is a guess.
  • Resample categorical data only with nearest neighbour. Bilinear on class codes invents classes.
  • Derive windowed transforms with windows.transform. Never reuse a parent transform for a subset.
  • Assert alignment at the training-data boundary. The check costs microseconds and prevents a week of chasing edge artefacts.
  • Record the reference grid identity. A hash of the transform plus shape, written into both files’ tags, makes mismatches obvious in a listing.

Production-Ready Code

from __future__ import annotations

import hashlib
import logging
import numpy as np
import rasterio
from rasterio.affine import Affine
from rasterio.enums import Resampling
from rasterio.warp import reproject

logger = logging.getLogger(__name__)


def grid_fingerprint(path: str) -> str:
    """Stable short hash of (crs, transform, width, height) — two files match or they don't."""
    with rasterio.open(path) as src:
        payload = "|".join([
            src.crs.to_string() if src.crs else "NOCRS",
            ",".join(f"{v:.9g}" for v in src.transform[:6]),
            f"{src.width}x{src.height}",
        ])
    return hashlib.sha256(payload.encode()).hexdigest()[:12]


def describe_misalignment(feature_path: str, label_path: str) -> dict:
    """Quantify how two grids differ. Returns a report; never raises."""
    with rasterio.open(feature_path) as f, rasterio.open(label_path) as l:
        ft, lt = f.transform, l.transform
        report = {
            "crs_match": f.crs == l.crs,
            "shape_match": (f.width, f.height) == (l.width, l.height),
            "pixel_size_match": (np.isclose(ft.a, lt.a) and np.isclose(ft.e, lt.e)),
            "dx_m": lt.c - ft.c,
            "dy_m": lt.f - ft.f,
            "feature_shape": (f.height, f.width),
            "label_shape": (l.height, l.width),
        }
    px, py = abs(ft.a), abs(ft.e)
    report["dx_px"] = report["dx_m"] / px if px else float("nan")
    report["dy_px"] = report["dy_m"] / py if py else float("nan")
    report["aligned"] = (
        report["crs_match"] and report["shape_match"] and report["pixel_size_match"]
        and abs(report["dx_px"] - round(report["dx_px"])) < 1e-6
        and abs(report["dy_px"] - round(report["dy_px"])) < 1e-6
    )
    return report


def assert_aligned(feature_path: str, label_path: str) -> None:
    """Fail loudly, with numbers, when the two grids are not identical."""
    r = describe_misalignment(feature_path, label_path)
    if r["aligned"] and r["dx_px"] == 0 and r["dy_px"] == 0:
        return
    raise ValueError(
        "label grid does not match the feature grid: "
        f"crs_match={r['crs_match']} shape_match={r['shape_match']} "
        f"pixel_size_match={r['pixel_size_match']} "
        f"offset=({r['dx_px']:.3f}, {r['dy_px']:.3f}) px "
        f"shapes {r['feature_shape']} vs {r['label_shape']}"
    )


def snap_categorical(label_path: str, feature_path: str, out_path: str) -> None:
    """Last-resort repair when the source vectors are gone: nearest-neighbour snap.

    Prefer re-burning from geometry. This preserves class codes exactly but cannot
    recover boundary detail that the original burn placed on the wrong grid.
    """
    with rasterio.open(feature_path) as ref, rasterio.open(label_path) as src:
        profile = src.profile.copy()
        profile.update(transform=ref.transform, width=ref.width, height=ref.height,
                       crs=ref.crs)
        dst_arr = np.full((ref.height, ref.width), src.nodata or 255, dtype=src.dtypes[0])
        reproject(
            source=rasterio.band(src, 1),
            destination=dst_arr,
            src_transform=src.transform, src_crs=src.crs,
            dst_transform=ref.transform, dst_crs=ref.crs,
            resampling=Resampling.nearest,          # never anything else for classes
            src_nodata=src.nodata, dst_nodata=src.nodata or 255,
        )
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(dst_arr, 1)
    logger.info("snapped %s onto the grid of %s", label_path, feature_path)

Step-by-Step Walkthrough

Step 1 — Fingerprint both files

print("features:", grid_fingerprint("features.tif"))
print("labels:  ", grid_fingerprint("labels.tif"))

Identical fingerprints mean you are done. Different ones tell you nothing about how they differ, which is what the next step is for.

Step 2 — Quantify the difference

r = describe_misalignment("features.tif", "labels.tif")
for k, v in r.items():
    print(f"{k:>18}: {v}")

A typical bad report reads dx_px: 0.5, dy_px: -0.5 with everything else True — the classic half-pixel case, usually from a grid built out of geometry bounds. A report with pixel_size_match: False is a different problem: the label raster was written at a different resolution and needs rebuilding, not snapping.

Step 3 — Re-burn from the source geometries

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

This is the correct repair whenever the vector source still exists, which it usually does. It is exact: the geometries are evaluated against the reference transform directly, so boundary pixels land where they belong.

Step 4 — If the vectors are gone, snap and record the compromise

snap_categorical("labels.tif", "features.tif", "labels_snapped.tif")
assert_aligned("features.tif", "labels_snapped.tif")

with rasterio.open("labels_snapped.tif", "r+") as dst:
    dst.update_tags(alignment="nearest-neighbour snap from labels.tif",
                    caveat="boundary pixels inherited from a misaligned burn")

Nearest-neighbour snapping moves the labels onto the right grid but cannot recover the boundary detail lost in the original burn. Recording that in the file tags stops a future reader from assuming the labels are pristine.

Step 5 — Wire the assertion into the data loader

class LabelledScene:
    def __init__(self, feature_path: str, label_path: str):
        assert_aligned(feature_path, label_path)      # cheap, and it never lies
        self.feature_path, self.label_path = feature_path, label_path
Three routes to the same broken label grid Three causes are shown feeding a single misaligned label raster: deriving the transform from vector bounds rather than the feature raster, reusing a parent transform after taking a window, and reprojecting the labels independently so GDAL chooses its own output grid. Each is paired with its one-line fix. Every offset comes from inventing a grid instead of inheriting one from_bounds(*gdf.total_bounds) origin lands on a vertex, not on the feature grid reusing a parent transform after Window(col_off, row_off) shifts by the offset separate reprojection GDAL picks its own grid offset labels opens fine, trains worse the one fix read transform, width, height and crs from the FEATURE raster then re-burn from geometry and assert before training

Verification

import numpy as np
import rasterio
from rasterio.transform import from_origin


def test_detects_half_pixel_offset(tmp_path):
    """A 5 m shift on a 10 m grid must be reported as 0.5 px, not as aligned."""
    ref = tmp_path / "features.tif"
    off = tmp_path / "labels.tif"
    base = {"driver": "GTiff", "count": 1, "dtype": "uint8", "width": 10, "height": 10,
            "crs": "EPSG:32633"}
    with rasterio.open(ref, "w", transform=from_origin(0, 100, 10, 10), **base) as d:
        d.write(np.zeros((10, 10), "uint8"), 1)
    with rasterio.open(off, "w", transform=from_origin(5, 95, 10, 10), **base) as d:
        d.write(np.zeros((10, 10), "uint8"), 1)

    r = describe_misalignment(str(ref), str(off))
    assert np.isclose(r["dx_px"], 0.5) and np.isclose(r["dy_px"], -0.5)
    assert not r["aligned"]


def test_whole_pixel_shift_is_not_a_subpixel_offset(tmp_path):
    """A clean 1-pixel shift is a windowing bug, not a sampling bug — report it as such."""
    r = describe_misalignment(str(ref_1px_shifted), str(labels))
    assert r["aligned"] is True          # grids are commensurate
    assert abs(r["dx_px"]) == 1.0        # but the extents differ


def test_reburn_reproduces_the_reference_fingerprint(tmp_path):
    out = tmp_path / "labels_aligned.tif"
    burn_labels("field_labels.gpkg", "features.tif", str(out))
    assert grid_fingerprint("features.tif") == grid_fingerprint(str(out))

A useful field check after repair: compute the class-boundary agreement between the old and new label rasters. If interior pixels agree above 99% while boundary pixels agree around 50%, you have just confirmed the offset was real and that the re-burn fixed exactly the pixels it should have.

Overall accuracy hides the offset; boundary accuracy exposes it Two pairs of bars compare a model trained on offset labels against the same model trained on re-burned aligned labels. Overall accuracy is nearly identical between the two. Boundary pixel accuracy is much lower for the offset labels and much higher after realignment. Why the offset survives review: the headline metric barely moves 0.91 0.93 overall accuracy 0.38 0.79 boundary pixels only offset labels re-burned Always report a boundary-restricted metric alongside the overall one.

FAQ

How do I tell whether my labels are offset from my features?

Compare the affine transforms numerically. If the c or f terms differ by anything that is not an exact multiple of the pixel size, the grids are offset. A visual overlay will not show a half-pixel shift at any usable zoom, which is precisely the magnitude that hurts training most.

Can I fix an offset by resampling the label raster?

Only with nearest neighbour, and only as a fallback. Any interpolating resampler invents class codes that never existed. Re-burning from the source geometry onto the reference transform is exact and is almost always available.

Does a half-pixel offset really matter?

At class boundaries, badly. Every edge pixel is trained against its neighbour’s spectra, so the model learns a blurred decision surface — bloated shapes in segmentation, a persistent misclassification fringe in per-pixel models. Interior pixels are unaffected, which is why overall accuracy barely moves and the problem survives review.

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