Vector-Raster Hybrid Processing

A model cannot read a raster. What it can read is a statement derived from one — “seventy per cent of this parcel is impermeable surface, from thirty-metre land-cover data captured in 2023” — and the whole of this topic is about producing statements like that correctly. The pixels have to be clipped to the right shape, sampled at a defensible resolution, summarised without inventing precision, and labelled with what they came from.

This topic belongs to spatial LLM architecture and core concepts and pairs with context-window optimization for maps, which decides how much room the resulting statements get. Which raster to use in the first place is a catalog question, answered by spatial metadata and catalog indexing.

Foundational Principles

Alignment precedes everything. A raster and a vector mask in different frames, or on grids offset by half a cell, produce a plausible number that is systematically wrong. Reproject the mask to the raster’s frame rather than the reverse, because resampling pixels degrades them and reprojecting a polygon does not.

Clip to the shape, never to its extent. A parcel’s bounding box can be twice its area for an elongated or L-shaped feature, and the pixels in the difference belong to somewhere else. The bounding box is a pre-filter for reading, not a substitute for the mask.

Report the resolution with the number. “Sixty-eight per cent impermeable” from thirty-metre data means something different from the same figure at one metre, and only one of the two supports a statement about a single building. The resolution travels with the statistic or the statistic is not usable.

Step-by-Step Implementation Pipeline

1. Read the raster’s grid definition before anything else

Every later decision depends on the transform, the frame and the nodata value. Read them once, and fail early if any is missing rather than discovering it as a wrong number later.

import logging
from dataclasses import dataclass
from typing import Optional

log = logging.getLogger("vector_raster")


@dataclass(frozen=True)
class Grid:
    epsg: int
    pixel_size_m: float
    nodata: Optional[float]
    captured_year: Optional[int]


class RasterUnusable(ValueError):
    """The raster cannot be used for a defensible statistic."""


def read_grid(dataset) -> Grid:
    """Extract the grid definition, refusing anything that cannot be interpreted."""
    epsg = dataset.crs.to_epsg() if dataset.crs else None
    if epsg is None:
        raise RasterUnusable("raster has no declared reference frame")
    if dataset.res[0] <= 0 or dataset.res[1] <= 0:
        raise RasterUnusable(f"implausible pixel size {dataset.res}")
    if abs(dataset.res[0] - dataset.res[1]) > 1e-6 * dataset.res[0]:
        log.info("non-square pixels %s — area statistics will use cell area", dataset.res)
    return Grid(epsg, float(dataset.res[0]), dataset.nodata,
                dataset.tags().get("captured_year"))

2. Reproject the mask into the raster’s frame

Moving the vector is cheap and lossless; moving the raster resamples every pixel and changes the values being counted. The direction of this transformation is the single most consequential decision in the whole pipeline.

from shapely.ops import transform as shapely_transform
from pyproj import CRS, Transformer


def mask_in_raster_frame(mask, mask_epsg: int, grid: Grid):
    """Bring the vector to the pixels, never the pixels to the vector."""
    if mask_epsg == grid.epsg:
        return mask
    tf = Transformer.from_crs(CRS.from_epsg(mask_epsg), CRS.from_epsg(grid.epsg),
                              always_xy=True)
    try:
        return shapely_transform(tf.transform, mask)
    except Exception as exc:
        raise RasterUnusable(f"mask could not be reprojected to EPSG:{grid.epsg}: {exc}") from exc

3. Check that the shape is big enough for the grid

A thirty-metre grid cannot say anything about a shape smaller than a few cells. Computing a percentage from four pixels is arithmetic, not measurement, and the honest response is to decline rather than to round.

MIN_CELLS_FOR_PROPORTION = 12


def cells_covered(mask, grid: Grid) -> float:
    """Approximate cell count for the mask at this resolution."""
    return mask.area / (grid.pixel_size_m ** 2)


def resolution_adequate(mask, grid: Grid) -> tuple[bool, str]:
    n = cells_covered(mask, grid)
    if n < MIN_CELLS_FOR_PROPORTION:
        return False, (f"the shape covers about {n:.0f} cells at "
                       f"{grid.pixel_size_m:g} m; too few for a proportion")
    return True, ""

The threshold is a judgement rather than a law, and stating it in the refusal is what makes it reviewable. A reader told “this parcel is four pixels across in the available data” understands the limitation immediately; one told “68%” does not.

4. Sample the pixels within the mask

The read uses the bounding box for efficiency and the mask for correctness, with nodata excluded from the denominator rather than counted as a class.

import numpy as np


def sample(dataset, mask, grid: Grid) -> tuple[np.ndarray, float]:
    """Return the values inside the mask and the share of the shape that had data."""
    try:
        window = dataset.window(*mask.bounds)
        block = dataset.read(1, window=window, masked=True)
        shape_mask = rasterize_mask(mask, dataset, window)     # True inside the polygon
    except Exception as exc:
        raise RasterUnusable(f"could not read the raster window: {exc}") from exc

    inside = block[shape_mask]
    if inside.size == 0:
        raise RasterUnusable("the mask selected no pixels")
    valid = inside.compressed() if hasattr(inside, "compressed") else inside
    if grid.nodata is not None:
        valid = valid[valid != grid.nodata]
    coverage = float(valid.size) / float(inside.size)
    return valid, coverage

Returning the coverage fraction alongside the values is what prevents the most common overstatement in this pipeline: a statistic computed over the third of a shape that had data, reported as though it described the whole.

5. Summarise into proportions, at a precision the grid supports

Class proportions are the useful output for categorical rasters and simple statistics for continuous ones. Both should be rounded to a precision the cell count can support — reporting a percentage to one decimal from forty cells implies a resolution that does not exist.

from collections import Counter


def class_proportions(values: np.ndarray, cell_count: int) -> dict[int, float]:
    """Proportions rounded to a precision the sample size supports."""
    if values.size == 0:
        return {}
    digits = 0 if cell_count < 100 else (1 if cell_count < 10_000 else 2)
    counts = Counter(int(v) for v in values)
    total = sum(counts.values())
    return {cls: round(100.0 * n / total, digits) for cls, n in counts.most_common()}

The precision rule is worth stating explicitly because rounding is where invented accuracy usually enters. Forty cells cannot distinguish 67.4% from 67.9%, and printing either one invites a reader to treat the difference as meaningful. The mechanics of turning these into prose are in summarising raster statistics for model prompts.

6. Render the statement with its provenance

The output is a sentence, and the sentence carries the resolution, the capture date, the coverage and the class labels. Omitting any of them produces a statement that reads as more certain than it is.

def to_statement(name: str, proportions: dict[int, float], labels: dict[int, str],
                 grid: Grid, coverage: float) -> str:
    """One sentence a model can quote, with everything a reader needs to judge it."""
    if not proportions:
        return f"{name}: no usable raster values within the shape."
    top = ", ".join(f"{labels.get(cls, f'class {cls}')} {pct:g}%"
                    for cls, pct in list(proportions.items())[:4])
    date = f", captured {grid.captured_year}" if grid.captured_year else ""
    cover = "" if coverage > 0.98 else f", from {coverage * 100:.0f}% of the shape"
    return f"{name}: {top} (from {grid.pixel_size_m:g} m data{date}{cover})."

7. Align tiles when the shape spans more than one

A shape crossing a tile boundary needs values from several files, and those files may not share a grid origin even when they share a frame. Reading them independently and concatenating the values is correct only if the grids align; where they do not, the honest options are to resample to a common grid deliberately or to report per tile. The alignment mechanics are covered in aligning raster tiles with vector masks.

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

8. Refuse rather than resample when the question needs precision

Resampling is sometimes necessary and always a loss, and for categorical data it is worse than a loss — averaging class codes produces values that correspond to no class at all. When a question needs a statistic that the available grid cannot support, the correct output is the refusal, with the reason.

def statistic_or_refusal(mask, grid: Grid, dataset, name: str, labels) -> str:
    ok, why = resolution_adequate(mask, grid)
    if not ok:
        return f"{name}: cannot be measured from this data — {why}."
    try:
        values, coverage = sample(dataset, mask, grid)
    except RasterUnusable as exc:
        return f"{name}: cannot be measured from this data — {exc}."
    props = class_proportions(values, int(cells_covered(mask, grid)))
    return to_statement(name, props, labels, grid, coverage)

Operating This Stage Over Time

Raster pipelines drift for reasons that have nothing to do with the code. The commonest is a source that changes its class scheme: a land-cover product reissues with an extra category, the numeric codes shift, and every proportion the pipeline produces is now labelled wrongly while remaining arithmetically correct. Pin the class scheme version alongside the data and assert that the codes present in a raster are all in the expected set — an unexpected code is a loud failure and a silent mislabelling otherwise.

The second is resolution creep in the opposite direction from what you would expect. Sources improve, a ten-metre product replaces a thirty-metre one, and the code that refused to measure small parcels starts succeeding. That is good, and it changes the meaning of historical answers: a statement produced last year for the same parcel said something different because it was measured differently. Recording the resolution in the statement, as step 6 does, is what makes those two answers comparable rather than contradictory.

The third is coverage. Cloud gaps, sensor outages and processing failures leave holes that vary over time, so the same query can return a confident number one month and a low-coverage caveat the next. Track the coverage fraction as a monitored value rather than only reporting it: a fall across many queries usually means a data problem upstream rather than a change in the shapes being asked about.

A useful discipline throughout is to treat every number this stage produces as a quotation from a dataset rather than as a measurement of the world. That framing makes the provenance fields feel necessary rather than decorative, and it is the framing a careful reader will apply anyway.

Failure Modes & Root Causes

The bounding-box statistic. A proportion computed over the shape’s extent rather than the shape, inflated by whatever surrounds it. Root cause: using the read window as the analysis unit. Mitigation: rasterize the mask and index by it, as step 4 does.

The resampled category. Class codes averaged during a reprojection, producing values corresponding to no real class. Root cause: reprojecting the raster instead of the mask, with a continuous resampling method. Mitigation: move the vector; where a raster must be resampled, use nearest-neighbour for categorical data and say so.

The confident percentage from four pixels. A statistic that is arithmetically valid and meaningless. Root cause: no minimum sample check. Mitigation: the cell-count gate, with the refusal carrying the cell count.

The vanished nodata. Nodata values counted as a class, producing a large proportion of something that was simply not measured. Root cause: reading without a mask or with the wrong nodata value. Mitigation: exclude explicitly and report the coverage fraction.

Production Validation Protocols

1. Mask-not-box assertion. For an L-shaped fixture, assert the mask statistic differs from the box statistic; if they match, the mask is not being applied.
2. Frame assertion. Assert the mask and raster frames match at the point of sampling, and that the mask was the thing reprojected.
3. Coverage reporting. Assert every emitted statement includes the coverage fraction whenever it falls below the threshold.
4. Class-code invariant. Assert every code encountered is in the expected scheme; an unknown code fails the build rather than being labelled as unknown.
5. Precision rule test. Assert the number of decimal places in an emitted proportion matches the cell count that produced it.
6. Refusal path test. Assert a shape smaller than the minimum cell count produces a refusal naming the cell count, not a number.

Of these, the mask-not-box assertion is the one that catches the most damaging bug for the least effort. It needs one irregular fixture shape and two lines, and it fails loudly the moment somebody optimises the read path by dropping the rasterized mask — an optimisation that looks harmless, speeds the query up noticeably, and changes every number the stage produces.

The class-code invariant deserves the same treatment for a different reason: it is the only check that notices when an upstream product reissues with a changed scheme. Everything else keeps working, the proportions still sum to a hundred, and only the labels are wrong.

Frequently Asked Questions

Related

• Up to the section overview: Spatial LLM Architecture and Core Concepts
• Technique: Aligning Raster Tiles with Vector Masks
• Technique: Summarising Raster Statistics for Model Prompts
• Peer topic: Context-Window Optimization for Maps
• Related topic: Spatial Metadata and Catalog Indexing
• Related topic: Evaluation and Benchmarking for Spatial LLMs