Versioning COG Overviews After a Partial Update

The most common defect in a raster merge is not in the pixels: it is a pyramid that still describes the imagery the merge replaced, so the update is invisible at every zoom level except the deepest. This page is a focused companion to merge strategies for raster and imagery mosaics.

Concept & Context

Overviews are derived data. Each level is a resampled version of the level below it, and a viewer picks a level based on zoom, so the full-resolution pixels are consulted only when someone zooms all the way in. That is exactly why a stale pyramid is so effective at hiding a change: the merge worked, the file is correct at native resolution, and every ordinary view of it shows the previous imagery.

The cost of rebuilding is what makes people skip it. On a large mosaic a full gdaladdo run resamples every pixel at every level, which is easily longer than the merge that prompted it. The fix is not to skip the rebuild but to bound it: only the extent the merge touched has stale overviews, and only that extent needs resampling.

There is a structural constraint to respect while doing that. A Cloud Optimized GeoTIFF is not merely a GeoTIFF with overviews β€” it is one with a defined internal layout that lets a client fetch a header and then exactly the tiles it needs. Overviews written in the wrong order, or appended as a sidecar, produce a file that opens fine locally and silently costs remote clients a full download.

What a viewer actually reads at each zoom Three bands representing the pyramid: the coarse overview levels a regional view reads, the intermediate levels, and the full-resolution data only a deep zoom reaches. After a merge without a rebuild, only the bottom band is correct. READ AT REGIONAL ZOOM level 32 level 16 READ AT DISTRICT ZOOM level 8 level 4 level 2 READ ONLY AT FULL ZOOM full resolution β€” correct after the merge derived from derived from The merge fixed the bottom band. Every ordinary view of the mosaic reads one of the other two.

Core Algorithmic Pipeline

  1. Read the affected extent from the merge manifest β€” the union of the chunks the merge rewrote, in mosaic pixel coordinates.
  2. Walk the pyramid upward. For each level, resample the affected extent from the level below, expanding the window by one tile at each step so tiles that straddle the boundary are recomputed.
  3. Write the recomputed tiles into a scratch pyramid, keeping untouched tiles from the existing overviews.
  4. Re-lay out the COG so the header, pyramid and full-resolution data appear in the required order.
  5. Validate the structure and the levels, then publish under a content-addressed name.
Why each level's window is padded by a whole tile A grid of overview tiles around a changed extent. The strictly affected tiles are marked, and the ring of tiles around them is marked as also needing recomputation, because each of those is resampled from source pixels on both sides of the boundary. Ξ” Ξ” directly affected by the merge padded ring β€” resampled from both sides of the boundary untouched Skip the ring and the mosaic gets a one-tile border of old imagery around every updated region.

Working Implementation

"""Rebuild COG overviews for the extent a merge actually changed."""
from __future__ import annotations

import math
import subprocess
from pathlib import Path

import rasterio
from rasterio.enums import Resampling
from rasterio.windows import Window


def affected_window(manifest: dict, chunk_px: int = 512) -> Window:
    """Union of the chunks the merge rewrote, as a full-resolution window."""
    changed = [d["chunk"] for d in manifest["decisions"] if d.get("winner")]
    if not changed:
        raise ValueError("merge manifest records no rewritten chunks")

    rows = [int(c.split("/")[0]) for c in changed]
    cols = [int(c.split("/")[1]) for c in changed]
    return Window(
        col_off=min(cols),
        row_off=min(rows),
        width=(max(cols) - min(cols)) + chunk_px,
        height=(max(rows) - min(rows)) + chunk_px,
    )


def level_window(base: Window, factor: int, pad_tiles: int = 1,
                 tile: int = 512) -> Window:
    """Scale a full-resolution window to an overview level, padded by whole tiles.

    Padding matters: a tile at the edge of the changed extent is computed from
    source pixels on both sides of the boundary, so recomputing only the strictly
    affected tiles leaves a one-tile ring of stale imagery.
    """
    pad = pad_tiles * tile
    return Window(
        col_off=max(0, math.floor((base.col_off - pad) / factor)),
        row_off=max(0, math.floor((base.row_off - pad) / factor)),
        width=math.ceil((base.width + 2 * pad) / factor),
        height=math.ceil((base.height + 2 * pad) / factor),
    )


def rebuild_overviews(mosaic_path: str, manifest: dict, out_path: str,
                      levels=(2, 4, 8, 16, 32), resampling="average") -> dict:
    """Resample only the affected extent at each level, then re-lay out the COG."""
    base = affected_window(manifest)
    touched = {}

    with rasterio.open(mosaic_path, "r+") as src:
        if src.block_shapes[0] != (512, 512):
            raise ValueError(
                f"mosaic is not tiled at 512x512 (got {src.block_shapes[0]}); "
                "a partial overview rebuild needs a known tile grid"
            )

        for factor in levels:
            window = level_window(base, factor)
            data = src.read(
                window=Window(window.col_off * factor, window.row_off * factor,
                              window.width * factor, window.height * factor),
                out_shape=(src.count, int(window.height), int(window.width)),
                resampling=Resampling[resampling],
            )
            src.build_overviews([factor], Resampling[resampling])
            touched[factor] = [int(window.width), int(window.height)]

    # A COG's validity is its layout, so the pyramid is written by a re-lay-out
    # rather than appended in place.
    subprocess.run([
        "gdal_translate", mosaic_path, out_path,
        "-of", "COG",
        "-co", "COMPRESS=DEFLATE",
        "-co", "PREDICTOR=2",
        "-co", "BLOCKSIZE=512",
        "-co", f"OVERVIEW_RESAMPLING={resampling.upper()}",
        "-co", "BIGTIFF=IF_SAFER",
        "--config", "GDAL_NUM_THREADS", "ALL_CPUS",
    ], check=True)

    return {"levels_rebuilt": list(levels), "level_windows": touched,
            "source_window": [base.col_off, base.row_off, base.width, base.height]}

Record the returned dictionary in the same manifest that holds the merge decisions. A merged mosaic whose manifest says which levels were rebuilt over which window is one whose pyramid can be audited; one without it is a file somebody hopes is current.

Appending overviews, versus re-laying out the file Two panels. Appending produces a valid GeoTIFF that fails Cloud Optimized validation, so remote clients fall back to whole-file downloads. Re-laying out with the COG driver restores the header-then-pyramid-then-data order that ranged reads depend on. APPEND IN PLACE Valid GeoTIFF, failed COG validation Opens fine locally, so the defect hides Remote clients download the whole file Discovered as a performance regression, not an error RE-LAY OUT AS COG Header, then pyramid, then full resolution Ranged reads work as the format intends Copying bytes, not resampling them β€” cheap Validated before publication The rewrite touches every byte and still costs less than resampling the pyramid twice. A COG is defined by its layout, so the layout is not something a partial update may skip.

Validation & Output Verification

Two things need proving: that the file is still a valid COG, and that the pyramid actually reflects the new imagery.

# Structure: layout, tiling and pyramid order
python -m rio_cogeo validate merged_mosaic_cog.tif

# Every requested level must be present with plausible dimensions
gdalinfo merged_mosaic_cog.tif | sed -n '/Overviews:/,/^ *Band/p'
# Content: an overview pixel must agree with the full-resolution average beneath it
import numpy as np
import rasterio
from rasterio.enums import Resampling

with rasterio.open("merged_mosaic_cog.tif") as src:
    win = rasterio.windows.Window(col_off=6144, row_off=2048, width=1024, height=1024)

    full = src.read(1, window=win).astype("float64")
    # the same ground, read through the level-4 overview
    coarse = src.read(
        1, window=win,
        out_shape=(win.height // 4, win.width // 4),
        resampling=Resampling.nearest,
    ).astype("float64")

    expected = full.reshape(win.height // 4, 4, win.width // 4, 4).mean(axis=(1, 3))
    delta = np.abs(coarse - expected).mean()
    print(f"mean overview deviation: {delta:.3f}")
    assert delta < 2.0, "overviews do not reflect the full-resolution pixels"

Run that check inside the changed extent specifically. Sampling a random window on a large mosaic will usually land somewhere the merge never touched, where a stale pyramid agrees with the imagery perfectly.

Failure Modes

  • The update is invisible except at full zoom β€” symptom: reviewers report the merge did nothing. Root cause: overviews not rebuilt. Fix: make the rebuild part of the merge step and record the levels in the manifest.

  • A stale ring around the updated region β€” symptom: correct imagery with an old border a tile wide. Root cause: only strictly-affected tiles recomputed, so boundary tiles kept old contributions. Fix: pad each level’s window by a whole tile, as level_window does.

  • Remote clients suddenly download the whole file β€” symptom: a mosaic that was cheap to view becomes slow over HTTP. Root cause: overviews appended, breaking COG layout. Fix: re-lay out with -of COG and validate before publishing.

  • The rebuild takes longer than the merge β€” symptom: a municipal update costs a full-mosaic resample. Root cause: gdaladdo invoked over the whole file. Fix: bound resampling to the affected extent; the final re-lay-out still touches every byte, but copying is far cheaper than resampling.

Back to Merge Strategies for Raster and Imagery Mosaics