This guide is part of Automated Vector & Raster Cleaning Workflows.

Raster Mosaicking and Tiling for Pipeline Outputs

The Problem: A Product Made of Scenes Is Not Yet a Product

An imagery pipeline that ends with a directory of processed scenes has produced inputs, not outputs. Scenes overlap, they have different acquisition dates, their edges do not align with anything a consumer cares about, and there is no way to ask for “the value at this coordinate” without knowing which file to open. Turning that into a product means three decisions: how the pieces are combined, what the addressable unit is, and how the result is written so that a reader can fetch part of it cheaply.

Get those wrong and the failures are expensive rather than subtle. A single continent-sized GeoTIFF has to be rewritten in full every time one input scene is reprocessed. A mosaic assembled without an explicit priority rule changes between runs because the file listing order changed. A product written without overviews forces every zoomed-out client to read full-resolution pixels, which turns a map pan into a bandwidth incident.

Prerequisites and Environment

pip install "rasterio>=1.3" "rio-cogeo>=5.3" "shapely>=2.0" "geopandas>=1.0" "numpy>=1.26"

GDAL supplies the VRT and warping machinery; rio-cogeo handles cloud-optimized output and validation. Confirm the GDAL version supports the COG driver directly rather than relying on the older two-step creation:

import rasterio

assert rasterio.__gdal_version__ >= "3.4", "COG driver and its predictors need GDAL 3.4+"
with rasterio.Env() as env:
    assert "COG" in rasterio.drivers.raster_driver_extensions().values() or True

Version and Compatibility Matrix

Pipeline Architecture: Scenes → Grid → Product

The stages are the same regardless of scale, and each one has a natural artefact.

Index. Build a footprint table of every input scene: geometry, acquisition datetime, cloud cover, source path, checksum. This is a GeoDataFrame and it is the only thing later stages need to consult in order to decide what belongs in a tile.

Assign. For each tile in the output grid, query the index for intersecting scenes and order them by the priority rule. The result is a small manifest per tile — a list of source paths in the order they should be painted.

Render. Warp each contributing scene into the tile’s grid and combine according to the seam rule. This is the only stage that touches pixels, and it is embarrassingly parallel across tiles.

Publish. Write each tile as a validated COG, refresh the virtual mosaic and the catalogue, and record the manifest that produced each tile alongside it.

Keeping the index and the manifests as data — rather than as logic inside the render step — is what makes the product reproducible. A tile can be rebuilt months later from its manifest without re-running the assignment, and a change in the priority rule shows up as a diff in the manifests before a single pixel is written.

Deterministic Seams

Where two scenes overlap, something has to decide which one wins, and “whichever GDAL saw first” is not a decision — it is an accident that changes when the file listing changes.

import geopandas as gpd


def rank_scenes(index: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Order candidate scenes so the winner in any overlap is fully determined."""
    return index.sort_values(
        by=["cloud_cover", "datetime", "scene_id"],
        ascending=[True, False, True],
        kind="mergesort",   # stable, so equal keys keep a defined order
    )

Three keys, applied in order, with the last one guaranteeing that no tie is ever broken by chance. mergesort matters: a stable sort means two scenes identical on every key still order the same way on every machine, which is the difference between a reproducible mosaic and one that differs subtly between the development run and the production one.

The same ordering must be recorded in the tile manifest, because that is what allows the question “why is this pixel from the April scene?” to be answered without re-deriving the whole assignment.

Choosing the Tile Grid

The grid is the product’s addressing scheme, and changing it later is a migration rather than a configuration change. Three properties matter more than the specific choice.

Alignment with the source. A grid aligned to the source scenes’ native footprints — MGRS squares for Sentinel-2, path/row for Landsat — means each tile draws from few scenes and most tiles rebuild after a single scene changes. A grid deliberately misaligned with the source maximises the number of scenes each tile touches.

Alignment with the consumer. A web map wants a quadkey or XYZ scheme so that a tile request maps to a stored object. An analyst wants administrative units. Serving both usually means picking the source-aligned grid for storage and generating the consumer scheme on read, rather than storing two products.

Cell size against feature size. A tile that takes longer than a few minutes to render is a poor unit of work — retries are expensive and progress is invisible. A tile smaller than a few thousand pixels on a side spends more time in metadata than in pixels. Between 2 048 and 8 192 pixels square suits most pipelines, which for 10 m imagery is a 20 to 80 km cell.

The grid definition belongs in configuration, next to the target CRS and the resampling rules, exactly as the fixed grid contract in reprojecting rasters to a common grid with rasterio.warp describes — and it should be treated as equally immutable.

Validation Before Publication

A tile that is written is not yet a tile that is published. Four checks run between the two, and all of them are cheap relative to the render that produced the data.

from rio_cogeo.cogeo import cog_validate
import rasterio


def validate_tile(path: str, expected_transform, expected_shape, nodata) -> None:
    valid, errors, warnings = cog_validate(path)
    if not valid:
        raise ValueError(f"{path} is not a valid COG: {errors}")

    with rasterio.open(path) as src:
        if src.shape != expected_shape:
            raise ValueError(f"{path}: shape {src.shape}, expected {expected_shape}")
        if not src.transform.almost_equals(expected_transform, precision=1e-6):
            raise ValueError(f"{path}: transform does not match the grid contract")
        if src.nodata != nodata:
            raise ValueError(f"{path}: nodata {src.nodata}, expected {nodata}")
        if not src.overviews(1):
            raise ValueError(f"{path}: no overviews — zoomed-out reads will be full resolution")

The COG validation catches structural problems — wrong tiling, missing internal layout — that make a file technically readable and practically useless over HTTP range requests. The grid assertions catch the far more common failure of a tile that rendered correctly onto the wrong grid, which composites into a visible seam that nobody notices until it appears in a published map.

Failure-Mode Reference

Production Integration Notes

Tile rendering is the clearest example in spatial ETL of a task that should be partitioned rather than parallelised ad hoc. One tile is one unit of work: independently retryable, independently addressable, and sized in minutes — the criteria set out in orchestrating spatial ETL pipelines. In Dagster it maps onto a partitioned asset keyed by tile and date; in Airflow onto a mapped task over the tile list; in Prefect onto a bounded .map over the manifests.

Whichever tool runs it, the manifest is what makes a rebuild cheap. When an input scene is reprocessed, query the index for the tiles whose manifests reference it and rebuild exactly those — usually a handful out of hundreds. That query is the entire argument for keeping the manifests as data rather than reconstructing the assignment each run.

Compositing Rules Beyond “Last One Wins”

Painting scenes in priority order and letting the last one win is the simplest seam rule and the right default, but three alternatives matter often enough to be worth naming.

Per-pixel best. Rather than choosing a winner per scene, choose per pixel using a quality band — a cloud mask, a scene classification layer, a per-pixel confidence raster. Each candidate contributes only its good pixels, and the composite has fewer holes than any single scene. The cost is that every candidate must be warped and read in full, rather than only where the previous ones left gaps.

Temporal reduction. Instead of choosing at all, reduce across the stack: median, maximum NDVI, or a percentile. This is standard practice for cloud-affected optical imagery, because the median across a season is far cleaner than any individual acquisition. The output no longer corresponds to a single date, which has to be stated explicitly in the product metadata or downstream users will assume it does.

Feathered blending. Weight contributions near the seam so the transition is gradual rather than a hard edge. It looks better in a rendered map and is usually wrong for analysis, because the blended pixels correspond to no actual observation. Reserve it for visual products and keep an unblended version for anything that will be measured.

The rule belongs in configuration alongside the priority keys, and it belongs in the tile manifest too. A product where some tiles were composited by median and others by last-one-wins — because the rule changed mid-backfill — is nearly impossible to diagnose after the fact and trivial to spot if each tile records what produced it.

Handling Partial Coverage Honestly

Almost every real mosaic has gaps: cloud that no acquisition cleared, an area outside the sensor’s swath, a scene that failed to download and never came back. How those are represented determines whether downstream consumers can reason about them.

Three representations, in increasing order of usefulness. Nodata pixels are the minimum: the gap is visible and arithmetic skips it, provided every reader honours the nodata value. A coverage mask band records per pixel how many observations contributed, which turns “is this a gap?” into “how confident is this pixel?” and costs one extra band. A coverage footprint stored as vector geometry alongside the tile lets a consumer answer the question without opening the raster at all, which is what a catalogue or a validity check actually wants.

The mistake to avoid is filling gaps silently — with a neighbouring value, with a prior date, with zero. Each of those produces a raster that looks complete and contains fabricated observations, and there is no downstream check that can distinguish the fabricated pixels from real ones once the provenance is gone.

Rebuild Cost and Why the Index Earns Its Keep

The economic argument for the index and the manifests is easiest to see in the rebuild case. Suppose one scene out of twelve hundred is reprocessed because its geolocation was corrected.

Without an index, the only safe response is to rebuild everything: twelve hundred scenes warped, four hundred tiles rendered, the whole product rewritten. With the index, the query is one line — which manifests reference this scene id — and the answer is typically three or four tiles. The render cost drops by two orders of magnitude, and more importantly the change is auditable: those four tiles have a new manifest, and every other tile is provably untouched.

The same query answers the reverse question, which is the one auditors ask. Given a suspicious pixel in a published product, the tile’s manifest names the contributing scenes in order, and the index names where each came from and when it was retrieved. That chain — pixel to tile to manifest to scene to source URI — is what makes a raster product defensible rather than merely available, and it is built from two small tables rather than from anything sophisticated.

Naming and Versioning Tiles

A tile’s object key is the closest thing a raster product has to a primary key, and it should carry four things: the grid cell, the temporal period, the product name and a version. products/ndvi_p50/v4/2026-04/31UEQ.tif reads unambiguously and sorts usefully; output/final2/tile_17.tif does not, and every pipeline accumulates a few of the second kind.

The version segment is the one most often omitted and the most valuable. When a compositing rule changes, writing the new tiles under a new version leaves the old ones intact for comparison, lets consumers migrate on their own schedule, and makes the rollback a pointer change rather than a rebuild. The cost is storage for one extra generation, which is trivial next to the cost of discovering a regression with no previous version to compare against.

Keep the period in the key even for products that feel static. A landcover mosaic that is “the current one” today becomes “the 2026 one” the moment a 2027 edition exists, and renaming published objects is the kind of migration that breaks other people’s pipelines.

Related

• Building Virtual Raster Mosaics with GDAL VRT
• Writing Cloud-Optimized GeoTIFFs with rio-cogeo
• Generating XYZ Tiles from a Raster
• Raster Alignment & Resampling Techniques
• Bulk Downloading Satellite Imagery