Aligning Raster Tiles with Vector Masks for LLM Context

A shape rarely sits inside one tile, and the tiles it spans do not always share a grid. Getting the alignment right is what separates a statistic that describes the shape from one that describes a rectangle near it, and the failure produces no error at any point. This guide handles the alignment and the reading, as the mechanical core of vector-raster hybrid processing.

When to Use This Approach

Use it whenever a mask spans more than one file, or whenever the mask and the raster arrive in different frames — which together cover almost every real request.

The third row is where judgement is required. Unaligned grids mean the pixels are not comparable, and pretending otherwise by resampling one to the other is a decision with consequences that should be recorded rather than absorbed.

Implementation

The reader checks alignment, rasterizes the mask per tile, and accumulates values with their coverage.

import logging
from dataclasses import dataclass
from typing import Iterable, Sequence

import numpy as np

log = logging.getLogger("tile_alignment")


@dataclass(frozen=True)
class TileGrid:
    epsg: int
    pixel_size_m: float
    origin: tuple[float, float]
    nodata: float | None


class Unalignable(ValueError):
    """The tiles cannot be combined without a decision the caller must make."""


def grids_align(a: TileGrid, b: TileGrid, tolerance: float = 1e-3) -> bool:
    """Same frame, same pixel size, origins offset by a whole number of cells."""
    if a.epsg != b.epsg:
        return False
    if abs(a.pixel_size_m - b.pixel_size_m) > 1e-6 * a.pixel_size_m:
        return False
    dx = (a.origin[0] - b.origin[0]) / a.pixel_size_m
    dy = (a.origin[1] - b.origin[1]) / a.pixel_size_m
    return abs(dx - round(dx)) < tolerance and abs(dy - round(dy)) < tolerance


def read_masked(datasets: Sequence, mask, mask_epsg: int,
                rasterize) -> tuple[np.ndarray, float]:
    """Read every pixel of the mask across the tiles. Raises if grids disagree."""
    if not datasets:
        raise Unalignable("no tiles cover the mask")

    grids = [TileGrid(d.crs.to_epsg(), float(d.res[0]), (d.transform.c, d.transform.f),
                      d.nodata) for d in datasets]
    reference = grids[0]
    for other in grids[1:]:
        if not grids_align(reference, other):
            raise Unalignable(
                f"tiles do not share a grid ({reference.epsg}/{reference.pixel_size_m:g} m "
                f"against {other.epsg}/{other.pixel_size_m:g} m) — resample deliberately "
                "or report per tile")

    projected = mask if mask_epsg == reference.epsg else _reproject(mask, mask_epsg,
                                                                   reference.epsg)
    collected, inside_total = [], 0
    for dataset, grid in zip(datasets, grids):
        try:
            window = dataset.window(*projected.bounds)
            block = dataset.read(1, window=window, masked=True)
            inside = rasterize(projected, dataset, window)      # True where the mask covers
        except Exception as exc:
            log.warning("tile read failed, skipping it: %s", exc)
            continue
        values = block[inside]
        inside_total += int(inside.sum())
        usable = values.compressed() if hasattr(values, "compressed") else values
        if grid.nodata is not None:
            usable = usable[usable != grid.nodata]
        collected.append(usable)

    if not collected or inside_total == 0:
        raise Unalignable("the mask selected no pixels in any tile")
    stacked = np.concatenate(collected)
    return stacked, float(stacked.size) / float(inside_total)

Raising rather than resampling when grids disagree is the decision this whole guide turns on. Resampling is sometimes the right answer and is never one a reader should discover by accident, so the function refuses and names the alternatives; the caller decides, and records the decision.

Skipping a tile that fails to read, rather than aborting, is the opposite instinct and is also right: a mask spanning six tiles where one is corrupt still has five tiles of real data, and the coverage fraction reports the loss. What must not happen is the loss being invisible, which is why coverage is a return value rather than a log line.

Validation & Testing

def test_offset_grids_are_refused():
    a = TileGrid(27700, 10.0, (0.0, 0.0), None)
    b = TileGrid(27700, 10.0, (5.0, 0.0), None)      # half a cell
    assert not grids_align(a, b)


def test_whole_cell_offsets_are_accepted():
    a = TileGrid(27700, 10.0, (0.0, 0.0), None)
    b = TileGrid(27700, 10.0, (30.0, -20.0), None)
    assert grids_align(a, b)


def test_mixed_frames_raise_rather_than_reprojecting_the_raster():
    try:
        read_masked([TILE_27700, TILE_4326], MASK, 27700, rasterize)
    except Unalignable as exc:
        assert "do not share a grid" in str(exc)
        return
    raise AssertionError("mixed frames must be refused, not silently resampled")


def test_unreadable_tile_reduces_coverage_without_aborting():
    values, coverage = read_masked([GOOD, GOOD, BROKEN], MASK, 27700, rasterize)
    assert values.size > 0 and coverage < 1.0


def test_empty_selection_raises():
    try:
        read_masked([TILE_ELSEWHERE], MASK, 27700, rasterize)
    except Unalignable:
        return
    raise AssertionError("a mask selecting no pixels must not return an empty statistic")

The third test is the one that will be argued about, and it is worth keeping. Resampling on the caller’s behalf is convenient, it works, and it silently changes categorical values into interpolated ones that correspond to no class — which is a corruption rather than a loss.

Build the fixtures from real tile metadata rather than from constructed grids. Origin values in production data carry the floating-point residue of whatever produced them, and a tolerance that looks generous against clean numbers can be too tight against real ones — which presents as a refusal to combine tiles that every other tool combines happily.

Gotchas & Edge Cases

Reprojecting the raster to match the mask. Resamples every pixel and, for categorical data, invents values. Move the vector; it is cheap and lossless.

Window computed from the mask’s bounds in the wrong frame. Produces a window somewhere else entirely and usually reads nothing, which presents as “no data for this shape”. Reproject the mask before computing the window, not after.

Coverage computed against the window rather than the mask. The window is a rectangle and the mask is not, so dividing by the window’s pixel count understates coverage badly for irregular shapes. Divide by the pixels inside the mask.

Tiles with different nodata values. Concatenating values from tiles whose nodata conventions differ leaves one tile’s nodata counted as data. Read each tile’s own nodata from its metadata rather than assuming a shared constant.

Overlapping tiles double-counting. Many tile schemes overlap by a few pixels at the edges, so a mask spanning a boundary counts those pixels twice. Deduplicate by position, or clip each tile’s window to its exclusive extent.

A mask spanning tiles from two product versions. They can share a grid and disagree about what a class code means, so concatenation produces proportions that are arithmetically clean and semantically mixed. Check the product version alongside the grid, and treat a mismatch exactly as you would treat a grid mismatch — as a refusal rather than as something to average over.

Alignment tolerance too tight. Origins stored as floating-point values accumulate tiny errors, and an exact comparison rejects grids that are aligned in every practical sense. A thousandth of a cell is a workable tolerance and a whole cell is not.

Frequently Asked Questions

Related

• Up to the parent topic: Vector-Raster Hybrid Processing
• Summarising Raster Statistics for Model Prompts
• Related technique: Indexing Catalog Collections for Agent Retrieval
• Related topic: Coordinate Reference System Normalization