Geospatial Orchestration Hub

Collecting Mapped Results into a Single Mosaic

The reduce step of a tile fan-out should never hold more than one block of pixels at a time. Build a VRT over the successful tile paths, translate it once into the output format, then add overviews — and account explicitly for what is missing, because a mosaic assembled from 398 of 412 tiles is a valid raster with holes and nothing about it says so. Memory stays flat, the merge is one sequential pass, and the coverage figure travels with the product.

When to Use This Pattern

  • A fan-out produced many tiles that must become one addressable raster or one published layer.
  • The output is larger than memory, which for any national mosaic it is.
  • Partial failure is possible, so the reduce step has to decide what an incomplete set means.
  • Consumers read windows rather than the whole file, which makes tiling and overviews part of the merge rather than an afterthought.

Complete Working Example

The reduce step accounts first and merges second, so an unacceptable coverage never becomes an output file at all.

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Optional, Sequence


@dataclass(frozen=True)
class MosaicResult:
    path: Optional[Path]
    expected: int
    present: int
    missing_keys: list[str]

    @property
    def coverage(self) -> float:
        return self.present / self.expected if self.expected else 0.0


def collect_mosaic(
    manifest: Sequence["TileItem"],
    tile_paths: Sequence[Optional[str]],
    dest: Path,
    min_coverage: float = 0.995,
    nodata: float = -9999.0,
) -> MosaicResult:
    """Merge what arrived, or refuse — but always say which."""
    present = [p for p in tile_paths if p and Path(p).exists()]
    missing = [item.work_key for item, p in zip(manifest, tile_paths)
               if not p or not Path(p).exists()]

    result = MosaicResult(None, len(manifest), len(present), missing)
    if result.coverage < min_coverage:
        # Refusing to publish is a decision, and it is usually the right one:
        # a mosaic with holes is indistinguishable from a complete one downstream.
        return result

    # A file list keeps the command line short — a VRT over 400 paths would
    # otherwise exceed the argument limit on some systems.
    listing = dest.with_suffix(".files.txt")
    listing.write_text("\n".join(present))

    vrt = dest.with_suffix(".vrt")
    run_gdal([
        "gdalbuildvrt",
        "-input_file_list", str(listing),
        "-srcnodata", str(nodata), "-vrtnodata", str(nodata),
        # Later files win where tiles overlap; ordering the list deterministically
        # is what makes the seam reproducible between runs.
        "-resolution", "highest",
        str(vrt),
    ])

    tmp = dest.with_suffix(".part.tif")
    run_gdal([
        "gdal_translate",
        "-co", "COMPRESS=DEFLATE", "-co", "PREDICTOR=2",
        "-co", "TILED=YES", "-co", "BLOCKXSIZE=512", "-co", "BLOCKYSIZE=512",
        "-co", "BIGTIFF=IF_SAFER",       # over 4 GB a plain GeoTIFF silently fails
        "-a_nodata", str(nodata),
        str(vrt), str(tmp),
    ])

    # Overviews once, on the finished file, with the same resampling the tiles used.
    run_gdal(["gdaladdo", "-r", "average", "--config", "COMPRESS_OVERVIEW", "DEFLATE",
              str(tmp), "2", "4", "8", "16", "32"])

    tmp.replace(dest)                     # the commit point
    return MosaicResult(dest, len(manifest), len(present), missing)

Recording coverage with the artefact is what stops an incomplete mosaic from being indistinguishable from a complete one three systems downstream:

def publish(s3, bucket: str, key: str, result: MosaicResult) -> None:
    s3.upload_file(str(result.path), bucket, key, ExtraArgs={"Metadata": {
        "tile-coverage": f"{result.coverage:.4f}",
        "tiles-expected": str(result.expected),
        "tiles-present": str(result.present),
    }})
Two ways to merge, two memory profilesReading all tiles into arrays and concatenating peaks at the full mosaic size. A VRT followed by a translate holds one block at a time and stays flat.26 GB13 GB0progress through the mergeconcatenate arrays — peaks at the whole mosaicVRT + translate — flat at 340 MBSame inputs, same output bytes. The upper approach also needs a worker sized for a number that grows with the data.
The flat line is what lets the reduce step run on the same worker class as the tiles, rather than needing a special large machine for the last step alone.

The ordering of the three GDAL steps is not interchangeable, and each one is doing something the others cannot. gdalbuildvrt decides the geometry — the union extent, the resolution, and which tile wins in an overlap — without touching a pixel, which is why it runs first and costs nothing. gdal_translate is the only step that reads and writes the actual data, so it is where compression, tiling and block size are chosen, and it is the step whose cost scales with the mosaic. gdaladdo then reads the finished file to build the pyramid, which requires the data to already be in its final form. Reordering them either loses the creation options or builds overviews of something that is about to change.

What each merge step decidesgdalbuildvrt decides extent, resolution and overlap order at negligible cost. gdal_translate decides compression, tiling and block size and costs one full pass. gdaladdo builds the pyramid from the finished file.gdalbuildvrtextent, resolution,overlap orderreads no pixelsgdal_translatecompression, tiling,block size, BIGTIFFone full sequential passgdaladdothe pyramid, fromthe finished filere-reads the outputSwap the last two and the overviews describe a file that no longer exists. Skip the first and the merge has no geometry to work from.
Only the middle step costs time proportional to the data, which is why the creation options belong there and nowhere else.

Parameter & Option Reference

Parameter Type Default Spatial notes
min_coverage float 0.995 Below this the mosaic is not published. A hole is worse than a late delivery for most products.
-input_file_list — required Four hundred paths on a command line exceeds the argument limit on some systems.
-srcnodata / -vrtnodata -9999 — Both. Omitting either lets nodata read as data at tile edges, producing visible squares.
BIGTIFF=IF_SAFER — on Over 4 GB a plain GeoTIFF fails, sometimes silently. Always set it for a mosaic.
PREDICTOR=2 — on Horizontal differencing. Meaningfully better compression on continuous-tone rasters.
gdaladdo -r average — — Match the resampling used for the tiles; nearest overviews on continuous data look wrong when zoomed out.
overview levels 2–32 — Five levels covers a typical viewer’s zoom range. More costs storage for little benefit.

Verification & Testing

The assertions worth making are about coverage accounting and about the merge not silently accepting less than it was given.

def test_incomplete_set_is_refused(tmp_path, manifest_412) -> None:
    paths = [str(tmp_path / f"{i}.tif") for i in range(398)] + [None] * 14
    result = collect_mosaic(manifest_412, paths, tmp_path / "mosaic.tif")
    assert result.path is None
    assert result.present == 398
    assert len(result.missing_keys) == 14
    assert 0.96 < result.coverage < 0.97


def test_complete_set_merges(tmp_path, manifest_412, rendered_tiles) -> None:
    result = collect_mosaic(manifest_412, rendered_tiles, tmp_path / "mosaic.tif")
    assert result.path and result.path.exists()
    assert result.coverage == 1.0


def test_merge_memory_stays_flat(tmp_path, manifest_412, rendered_tiles) -> None:
    peak = measure_peak_rss(collect_mosaic, manifest_412, rendered_tiles,
                            tmp_path / "mosaic.tif")
    assert peak < 1_500_000_000, "the merge is holding more than one block"

Afterwards, three gdalinfo facts confirm the product is what it claims to be — the extent, the nodata and the presence of overviews. The last is the one most often missing, and its absence is invisible until a viewer tries to render the whole mosaic at once:

gdalinfo mosaic.tif | grep -E "Size is|NoData|Overviews|Block="
# Size is 84213, 62110
# Band 1 Block=512x512 Type=Float32
#   NoData Value=-9999
#   Overviews: 42107x31055, 21054x15528, 10527x7764, 5264x3882, 2632x1941
A hole is invisible in the file and visible in the metadataA mosaic assembled from 398 of 412 tiles has a rectangular gap. The file itself is entirely valid; only the recorded coverage figure distinguishes it from a complete one.a valid GeoTIFF, nodata where the tiles are missingobject metadatatile-coverage: 0.9660tiles-expected: 412tiles-present: 398below min_coverage — do not publishthe only place the incompleteness is stated
Nothing about the raster on the left is invalid. Every downstream consumer will treat its nodata region as genuinely absent data rather than as a pipeline failure.

Common Pitfalls

  • Merging with numpy. Reading every tile into arrays and concatenating rebuilds the memory peak that tiling removed. The VRT exists precisely to avoid this and costs kilobytes.
  • Forgetting BIGTIFF. A mosaic crossing 4 GB fails to write a plain GeoTIFF, sometimes with an error and sometimes with a truncated file. IF_SAFER costs nothing and removes the category.
  • Omitting -vrtnodata. Setting only -srcnodata leaves the VRT without a nodata declaration, so the translate treats the fill value as data and the mosaic acquires visible rectangles at tile boundaries.
  • Publishing without a coverage figure. A mosaic with holes looks exactly like a complete one. Recording expected and present counts in the object’s metadata is three lines and is the only durable evidence.
  • Building overviews per tile instead of on the mosaic. Tile-level overviews are discarded by the merge and cost the fan-out real time. Build them once, on the finished file.

Frequently Asked Questions

Should the VRT be kept, or is it an intermediate?

Keep it. It is a few kilobytes, it documents exactly which tiles composed the mosaic and in what order, and it can be re-translated without re-rendering if the output format or compression needs to change. Deleting it saves nothing and discards the most useful provenance artefact the merge produces.

What if tiles overlap?

gdalbuildvrt lets the later file win in the overlap, so the ordering of the input list decides the seam. Sort the list deterministically — by tile index, or by acquisition date if you want newer imagery on top — and record the ordering rule. An unsorted list means the seam moves between runs, which produces mysterious diffs in otherwise identical outputs.

How do I choose `min_coverage`?

From what the product is for. A visual basemap tolerates a missing tile over open sea and not one over a city, so a flat percentage is a crude proxy; where it matters, weight the check by whether the missing tiles intersect a populated-area mask. For most pipelines, starting at 0.995 and refusing to publish below it is a reasonable default that gets refined the first time it fires.

Should the reduce step retry, and what would it retry?

It should retry, but only the merge — not the tiles. A translate that fails on a full disk or a transient object-store error is worth one more attempt, and because the whole step writes to a .part file and renames, a retry is naturally idempotent. What it must not do is trigger re-rendering: the tile paths it was given either exist or do not, and a missing one is a coverage question rather than something the merge can fix. Keeping that boundary sharp is what stops a reduce failure from cascading into a full re-run of the fan-out.

Can the merge run while tiles are still rendering?

Only if you are willing to re-run it. GDAL will happily build a VRT over whatever exists at the moment you ask, so a merge that starts early produces a mosaic missing the tiles that were still in flight. The dependency in the flow graph exists for exactly this reason, and skipping it to save a minute costs a full re-merge.

← Dynamic Task Mapping for Tile Fan-Out

Related in this section