Validating Raster Alignment and Resolution with Rasterio

Alignment defects are the quietest failures in a raster pipeline. A tile with a pixel size of 0.9999998 metres, or an origin at 412,000.5 instead of 412,000, opens fine, renders fine and produces a mosaic that is imperceptibly soft. Every downstream product built from it has been resampled once more than it should have been. This guide asserts the grid properties directly from the affine transform with rasterio — pixel size, rotation, origin snapping, tile adjacency and seam continuity — all of it from headers, which makes it cheap enough to run on every tile of every delivery. It implements the alignment rules introduced in Raster and Elevation Data Quality Checks.

Prerequisites

  • rasterio 1.3+ and Python 3.10+. The affine package ships with rasterio and provides the transform object used throughout.
  • A product grid definition: the pixel size in CRS units, and the grid interval that origins must snap to. These are usually the same number for a simple product (a 1 m raster on a 1 m grid) and different for a tiled scheme (a 0.5 m raster on a 1,000 m tile grid).
  • A tile index or naming convention if you validate adjacency. Deriving neighbours from bounds alone works but is O(n²); a tile identifier that encodes row and column makes it linear.
  • Read access to headers only. These checks never read pixel data except the optional seam comparison in Step 5, so they run happily against /vsicurl/ sources without transferring the raster.

Step-by-Step Procedure

Step 1 — Read the transform and turn it into checkable facts

# align/step1_transform.py
import rasterio


def grid_facts(path: str) -> dict:
    """The six affine terms, expressed as the properties a specification talks about."""
    with rasterio.open(path) as src:
        t = src.transform
        return {
            "path": path,
            "pixel_width": t.a,
            "row_rotation": t.b,
            "origin_x": t.c,
            "col_rotation": t.d,
            "pixel_height": t.e,      # negative for a north-up raster
            "origin_y": t.f,
            "width": src.width,
            "height": src.height,
            "bounds": tuple(src.bounds),
            "epsg": src.crs.to_epsg() if src.crs else None,
        }

Verification: on a standard north-up product, row_rotation and col_rotation are exactly 0.0 and pixel_height is the negative of pixel_width. Print the facts for one known-good tile and keep it as the reference in your test fixtures.

The five checks, in the order they can failChain of five checks running from the affine transform: resolution, orientation, origin snapping, tile adjacency and finally seam continuity, which is the only one that reads pixels.Resolutionpixel sizevs the specOrientationrotation = 0north-upSnappingorigin onthe gridAdjacencybounds abutexactlySeampixels agreereads data
Four of the five checks are pure header arithmetic; only the last one costs a read.

Step 2 — Assert resolution and orientation

# align/step2_resolution.py
TOL = 1e-6


def check_resolution(facts: dict, target_px: float) -> list[str]:
    problems = []
    if abs(abs(facts["pixel_width"]) - target_px) > TOL:
        problems.append(
            f"pixel width {facts['pixel_width']!r} != target {target_px}")
    if abs(abs(facts["pixel_height"]) - target_px) > TOL:
        problems.append(
            f"pixel height {facts['pixel_height']!r} != target {target_px}")
    if facts["pixel_height"] > 0:
        problems.append("raster is south-up (positive pixel height)")
    if facts["row_rotation"] != 0.0 or facts["col_rotation"] != 0.0:
        problems.append(
            f"raster is rotated ({facts['row_rotation']}, {facts['col_rotation']})")
    return problems

Verification: run against a tile you deliberately warp with gdalwarp -tr 0.999 0.999. The resolution check must fire. Note that the tolerance is deliberately tight: a pixel size that is nearly the target is a resampled raster, not a rounding artefact.

Step 3 — Test origin snapping against the product grid

This is the check most pipelines lack, and the one that prevents silent resampling downstream.

# align/step3_snapping.py
TOL = 1e-6


def snap_offset(value: float, interval: float) -> float:
    """Distance from `value` to the nearest multiple of `interval`."""
    r = value % interval
    return min(r, interval - r)


def check_snapping(facts: dict, grid_interval: float) -> list[str]:
    problems = []
    dx = snap_offset(facts["origin_x"], grid_interval)
    dy = snap_offset(facts["origin_y"], grid_interval)
    if dx > TOL:
        problems.append(
            f"origin_x {facts['origin_x']} is {dx:.6f} off the {grid_interval} grid")
    if dy > TOL:
        problems.append(
            f"origin_y {facts['origin_y']} is {dy:.6f} off the {grid_interval} grid")
    return problems

Verification: take a compliant tile, shift it by half a pixel with gdal_translate -a_ullr, and confirm both offsets are reported. Report the offset value, not just a boolean — a consistent half-pixel offset across a delivery points at a different production convention (pixel-centre versus pixel-corner registration) rather than at random error.

Step 4 — Verify tile adjacency from bounds

# align/step4_adjacency.py
from itertools import combinations

TOL = 1e-6


def adjacency_report(all_facts: list[dict]) -> list[dict]:
    """Flag tile pairs that nearly abut but do not abut exactly."""
    issues = []
    for a, b in combinations(all_facts, 2):
        ax0, ay0, ax1, ay1 = a["bounds"]
        bx0, by0, bx1, by1 = b["bounds"]

        vertical_overlap = min(ay1, by1) - max(ay0, by0) > TOL
        horizontal_overlap = min(ax1, bx1) - max(ax0, bx0) > TOL

        if vertical_overlap:
            gap = bx0 - ax1
            if 0 < abs(gap) < a["pixel_width"] * 2 and abs(gap) > TOL:
                issues.append({"a": a["path"], "b": b["path"],
                               "kind": "horizontal gap/overlap", "delta": gap})
        if horizontal_overlap:
            gap = ay0 - by1
            if 0 < abs(gap) < abs(a["pixel_height"]) * 2 and abs(gap) > TOL:
                issues.append({"a": a["path"], "b": b["path"],
                               "kind": "vertical gap/overlap", "delta": gap})
    return issues

Verification: on a clean delivery this returns an empty list. A delta equal to one pixel size means the tiling scheme double-counts or skips an edge row — common when a producer switches between inclusive and exclusive tile bounds.

Step 5 — Measure seam continuity where it matters

Bounds arithmetic proves the tiles touch; only pixels prove they agree.

# align/step5_seam.py
import numpy as np
import rasterio


def seam_difference(left: str, right: str) -> dict:
    """Mean and max absolute difference along the shared column of two tiles."""
    with rasterio.open(left) as l, rasterio.open(right) as r:
        if abs(l.bounds.right - r.bounds.left) > 1e-6:
            return {"comparable": False, "reason": "tiles do not abut"}
        rows = min(l.height, r.height)
        a = l.read(1, window=((0, rows), (l.width - 1, l.width))).astype("float64")
        b = r.read(1, window=((0, rows), (0, 1))).astype("float64")
        nodata = l.nodata

    if nodata is not None:
        valid = (a != nodata) & (b != nodata)
        a, b = a[valid], b[valid]
    if a.size == 0:
        return {"comparable": False, "reason": "no overlapping valid data"}

    diff = np.abs(a - b)
    return {
        "comparable": True,
        "compared_cells": int(a.size),
        "mean_abs_diff": float(diff.mean()),
        "max_abs_diff": float(diff.max()),
        "median_signed_diff": float(np.median(b - a)),
    }

Verification: the median_signed_diff is the diagnostic worth watching. Random noise gives a median near zero with a non-trivial mean absolute difference; a constant vertical offset between production batches gives a median equal to that offset, which no amount of edge blending will fix.

Interpreting Results

Result What it means Action
Pixel size off by < 1e-6 Floating-point representation of the correct value Pass
Pixel size off by a visible fraction The raster has been resampled since production Reject; source the original
Origin offset = half a pixel Pixel-centre versus pixel-corner registration mismatch Fix the convention at source; do not warp
Origin offset arbitrary Grid built from an ad-hoc extent, not the product grid Regenerate with a snapped transform
Non-zero rotation Raw sensor geometry, or a warp without -tap Re-warp with target-aligned pixels
Adjacency delta = one pixel Inclusive/exclusive bounds confusion in the tiler Fix the tiler; re-cut the affected row
Seam median offset constant Calibration or datum difference between batches Escalate to the producer with both tile identifiers

The distinction that matters most is between defects you may repair locally and defects that require the source. Anything that would need a resample to fix — resolution drift, arbitrary origins, rotation — must go back to production, because repairing it locally means resampling, which is the very thing the check exists to prevent.

Offset patterns and what each one meansGrid of four origin-offset patterns with the cause and the correct response: exactly half a pixel, a constant arbitrary offset, an offset that varies per tile, and no offset at all.CauseResponseExactly half a pixelpixel-centre vs corner registrationfix the convention at sourceConstant, arbitrarywarp without target-aligned pixelsre-warp with -tapVaries per tileeach tile cut from its own extentregenerate from a snapped gridZeroalignedpassNever fix an alignment defect by warping locally — the warp resamples, which is the thing the check exists to prevent.
The offset value is the diagnosis: half a pixel is a convention, an arbitrary number is a missing flag.

Gotchas & Edge Cases

gdalwarp does not snap by default. Warping to a new CRS produces an origin derived from the reprojected corner coordinates, which is almost never on a round grid. Add -tap (target aligned pixels) together with -tr to force the output onto the grid. A pipeline that reprojects without -tap generates unaligned tiles indefinitely.

What a seam comparison actually readsSequence diagram between the seam checker, the left tile, the right tile and the report. The checker confirms the bounds abut, reads one column from each tile, compares the valid overlapping cells and reports the median signed difference.seam checkleft tileright tilereportbounds + pixel sizebounds + pixel sizeread last columnread first columnmean, max and median signed differenceThe median signed difference is the diagnostic: near zero is noise, a constant value is a calibration or datum difference between production batches.
Two one-pixel-wide reads answer a question no header check can: do the tiles actually agree?

Geographic coordinates need a different grid interval. A 1 m grid interval is meaningless in degrees. For an EPSG:4326 product, express the grid in degrees (for example 1/3600 for one arc-second) and keep the tolerance relative — an absolute 1e-6 tolerance in degrees is about 11 cm, which is coarse for a snapping test on fine grids.

Bounds are floating point, so compare with a tolerance, never with equality. left.bounds.right == right.bounds.left looks correct and fails randomly. Every comparison in this guide uses an explicit tolerance for that reason.

A tile can be aligned and still be in the wrong CRS. These checks say nothing about the reference system; run them alongside the CRS assertion from the parent topic. An aligned grid in the wrong projection is a well-organised mistake.

Overviews inherit the transform, so a misaligned raster produces misaligned overviews. Fixing the base raster means rebuilding the overviews; do not assume gdaladdo output survives a transform correction.

When to Escalate

  • Alignment defects across an entire delivery indicate a production configuration problem. Escalate once, with the offending transform values, rather than filing a finding per tile.
  • Seam offsets with a constant median mean two batches were processed against different calibration or datum settings. This needs the producer; edge feathering hides it and preserves the error.
  • Grid interval disagreements between your specification and the delivered data are a contract discussion, not a validation failure. Confirm which grid the product is supposed to use before rejecting thousands of tiles.
  • Rotated rasters arriving in a gridded product feed should be resolved at the ingestion boundary — either the feed is raw imagery and needs a documented warp step with -tap, or the producer is shipping the wrong product.

For the checks that follow alignment — void accounting and structural layout — continue with Detecting Voids and Nodata Gaps in DEMs.

Frequently Asked Questions

Why does a half-pixel offset matter if the raster still renders correctly?

Because every operation that combines the raster with another grid must resample it. A mosaic, a zonal statistic, a raster algebra expression or a tile join all silently interpolate to reconcile the grids, and interpolation changes values. The map still looks right, the numbers move, and nothing in the output records that a resample happened. Asserting alignment at ingestion costs two modulo operations; discovering it after a year of derived products is a re-processing project.

What tolerance should an alignment check use?

Tight — around 1e-6 of a CRS unit for projected coordinates. Alignment is not a measurement with uncertainty; it is an arithmetic property of the transform, and any deviation beyond floating-point representation error means the grid was genuinely built somewhere else. A loose tolerance here hides exactly the defect the check exists to find.

How do I check adjacency without loading pixels?

Bounds arithmetic is enough. Two horizontally adjacent tiles must satisfy left.bounds.right == right.bounds.left within tolerance, and share top and bottom bounds if they are in the same row. Both values come from the header, so an adjacency audit over ten thousand tiles is a header scan rather than a data scan — seconds, not hours.

Is a rotated raster always wrong?

Not intrinsically — a rotated transform is valid, and raw sensor imagery often carries one. It is wrong for a gridded product, because a rotated grid cannot align to a north-up product grid and every consumer will resample it. Treat rotation as acceptable in raw inputs and as a blocker at the point where data enters the product catalogue.

Related

Back to Raster and Elevation Data Quality Checks