Resampling and Overviews When Writing Cloud-Optimized GeoTIFFs

A Cloud-Optimized GeoTIFF is an ordinary GeoTIFF whose bytes are arranged so a client can fetch a useful piece of it without downloading the rest, and three creation choices decide whether that arrangement actually holds: the internal block size, the overview factor set, and the resampling kernel used to build those overviews. This guide is for anyone publishing analysis-ready rasters with Raster Data Handling with Rasterio in Mastering Core Geospatial Python Libraries; the read side of the same contract is covered in windowed reads from a Cloud Optimized GeoTIFF.

Why This Approach / What Goes Wrong

The format specification is short: the image must be internally tiled, it must carry a pyramid of reduced-resolution copies, and the directory structures that index them must sit at the front of the file, ahead of the pixel data they point at. A reader opens the file with one small ranged request, learns where every tile of every level lives, and then asks for exactly the tiles it needs. Break any of the three and the file still opens perfectly in QGIS or rasterio.open() — it simply stops being cheap to read remotely, silently, with no error anywhere in your pipeline.

Three distinct failures produce that outcome, and they need different fixes. The first is a file that is tiled but has no pyramid: every zoomed-out request decodes full-resolution tiles and downsamples them on the client, so a continental preview costs gigabytes of transfer. The second is a pyramid written in the wrong place. Opening a plain GTiff for writing, calling build_overviews(), and closing appends the overview directories and their pixels after the main image, because GDAL has already streamed the full-resolution tiles to disk by then. The bytes are all there, the pyramid works locally, and rio cogeo validate rejects the file — correctly, because a remote client would have to seek to the tail before it could plan a single read.

The third failure is structural perfection with wrong numbers in it. Overviews are new pixels, computed by a kernel you chose, and the kernel has to match what the values mean. Averaging four land-cover codes produces a code that names no class; averaging four elevations produces exactly the elevation the level should hold. Nothing validates this — the file is a valid COG either way, and the damage only shows up as a zoomed-out map full of classes nobody defined.

Byte layout of a valid COG compared with a GeoTIFF that had overviews appended Two file ribbons drawn as byte sequences from left to right. The valid COG, written by the COG driver, starts with the header and all image file directories, followed by the overview tiles from coarsest to finest, and ends with the full-resolution tiles; a bracket under the header segment notes that one sixteen kilobyte ranged request retrieves every tile offset in the file. The second ribbon is a plain GeoTIFF written first and given overviews afterwards: the full-resolution tiles occupy the front of the file and the overview directories and overview tiles are appended at the end, so a client must seek to the tail before it can discover the pyramid, and rio cogeo validate rejects the file. Both files contain identical pixels and identical overviews. Byte order decides whether a client can plan its reads Valid COG · directories first, pyramid before full resolution header + IFDs overview tiles · coarsest first full-resolution tiles one ~16 KB ranged request → every tile offset in the file the client then fetches only the tiles it needs, at the level it needs Plain GTiff written, then build_overviews() — same pixels, wrong order full-resolution tiles overview IFDs overview tiles the pyramid is only discoverable after a seek to the end of the file — rio cogeo validate rejects it Identical pixels, identical overviews, identical checksums per tile — only one of the two is a COG.
The pyramid being present is not enough; it has to be indexed and stored ahead of the full-resolution tiles, which is what the COG driver and COPY_SRC_OVERVIEWS guarantee.

Prerequisites

conda install -c conda-forge "rasterio=1.3.*" "rio-cogeo=5.0.*" "numpy=1.26.*"

Install the whole stack from conda-forge in one solve — Rasterio links GDAL, which links PROJ, and the compression codecs are compiled into that same GDAL. A pip wheel mixed with a system GDAL is the usual reason COMPRESS=ZSTD raises at write time on one machine and works on another.

Step-by-Step Implementation

1. Inspect the source and confirm it is georeferenced. A COG without a CRS is a thumbnail, not a dataset, so fail loudly here rather than shipping one.

import rasterio

src_path = "dem_utm33n.tif"          # 20 000 x 20 000 float32 elevation model

with rasterio.open(src_path) as src:
    if src.crs is None:
        raise ValueError(f"{src_path} has no CRS — assign one before writing a COG")
    print(src.crs, src.dtypes[0], src.nodata)   # EPSG:32633 float32 -9999.0
    print(src.width, src.height)                # 20000 20000
    print(src.block_shapes[0])                  # (1, 20000) — striped source
    print(src.overviews(1))                     # [] — no pyramid yet

Keep the source CRS. A COG destined for a web map does not need to be in Web Mercator: rio-tiler-style tile servers reproject on the fly, and re-gridding a metric UTM archive into EPSG:3857 resamples every pixel and destroys the equal-ground-spacing property that elevation and reflectance analysis depends on. Reproject a display copy if you need one; never the archive.

2. Derive the overview factor set from the raster's size. The rule is mechanical: keep halving until the coarsest level fits inside a single internal block, so that a whole-scene preview is one tile fetch. Stopping early is the most common under-build — a pyramid ending at 4× on a 20 000-pixel raster still forces a 5 000-pixel decode for a thumbnail.

def overview_factors(width: int, height: int, blocksize: int = 512) -> list[int]:
    """Halve until the coarsest level fits inside one internal block."""
    factors, level = [], 1
    while max(width, height) // 2 ** (level - 1) > blocksize:
        factors.append(2 ** level)
        level += 1
    return factors


print(overview_factors(20000, 20000))   # [2, 4, 8, 16, 32, 64]
print(overview_factors(3000, 1800))     # [2, 4, 8]

Use 512 × 512 blocks as the default for imagery and elevation. 256 makes each request smaller but multiplies the number of tile offsets in the header — a 20 000-pixel scene jumps from 1 600 to 6 400 entries per level — while 1 024 wastes bandwidth on any client that wants a small window. The block size and the pyramid depth are the same decision: the pyramid stops when a level fits one block. On multi-band files keep the default pixel interleave as well, so one tile fetch returns every band for that footprint; band interleave forces a separate range request per band and turns a four-band patch read into four.

Overview levels of a 20 000 pixel scene and the share of bytes each holds Seven horizontal bars, one per stored level of a twenty thousand by twenty thousand pixel raster tiled at 512 pixels. Full resolution holds 1600 tiles and 75 percent of the stored bytes. The 2x level holds 400 tiles and 18.8 percent, the 4x level 100 tiles and 4.7 percent, the 8x level 25 tiles and 1.2 percent, the 16x level 9 tiles and 0.3 percent, the 32x level 4 tiles and 0.07 percent, and the 64x level a single 312 pixel tile and 0.02 percent. The 16x level is highlighted as the one a 1024 pixel viewport of the whole scene resolves to. Together the six overviews add roughly a third to the file size. One pyramid: what each level costs and which one answers a read level stored raster at that level share of stored bytes 1:1 full res 20 000 × 20 000 px · 1 600 tiles 75.0% 2× 10 000 × 10 000 px · 400 tiles 18.8% 4× 5 000 × 5 000 px · 100 tiles 4.7% 8× 2 500 × 2 500 px · 25 tiles 1.2% 16× 1 250 × 1 250 px · 9 tiles ← a 1 024 px viewport of the whole scene lands here 0.3% 32× 625 × 625 px · 4 tiles 0.07% 64× 312 × 312 px · 1 tile — fits one block 0.02% Each level is a quarter of the one above, so the six overviews together add about 33% to the file — the cheapest third you will ever spend.
The pyramid is geometric, so depth is nearly free: every level after the first two rounds to noise in the file size while removing a full decode from the client.

3. Choose the resampling kernel from what the band means. This is the only step in the recipe that no validator can check for you, and the only one that changes the data rather than its packaging.

from rasterio.enums import Resampling

OVERVIEW_RESAMPLING = {
    "elevation":   Resampling.average,   # continuous float — mean is meaningful
    "reflectance": Resampling.average,
    "landcover":   Resampling.mode,      # class codes — majority, never a mean
    "cloud_mask":  Resampling.nearest,   # binary 0/1 — keep it binary
    "rgb_basemap": Resampling.average,   # appearance only
}

kernel = OVERVIEW_RESAMPLING["elevation"]
print(kernel.name)                       # average

One dependency turns average from correct into corrupting: nodata. GDAL excludes nodata pixels from an averaged overview only when the band declares a nodata value or carries an internal mask. On a scene where fill is -9999 but src.nodata is None, every coarse level near the edge averages real elevations with fill and produces a dark fringe that spreads one pixel wider at each level. Set nodata before you build, or write an internal mask.

Which overview resampling kernel matches which band semantics A four-row matrix mapping band semantics to a resampling kernel. Continuous float bands such as elevation, reflectance and NDVI take Resampling.average, because the mean of four neighbours is a value the sensor could have measured. Classified integer bands such as land cover and cloud classes take Resampling.mode, keeping the majority class, with nearest as a faster fallback. Binary masks take Resampling.nearest, because any interpolation turns a hard edge into fractional values that no longer round-trip to zero or one. Visual RGB byte imagery takes Resampling.average, or RMS for amplitude data, since only appearance matters. A warning panel below shows two class codes, three and five, averaged into a class code four that appears in no legend. Pick the kernel from what the numbers mean, not from what looks smooth band semantics overview resampling why Continuous float DEM, reflectance, NDVI Resampling.average bilinear also defensible The mean of four neighbours is a value the sensor could have measured; the level mean stays honest. Classified integer land cover, cloud classes Resampling.mode nearest when speed wins Codes are labels, not magnitudes. Mode keeps the majority; a mean invents codes that do not exist. Binary mask 0 / 1 validity layer Resampling.nearest mode is equivalent here Interpolation turns a hard edge into fractions that no longer round-trip to 0 or 1 after casting. Visual RGB byte basemap, orthophoto Resampling.average rms for SAR amplitude Only appearance matters, and average removes the aliasing that nearest shows on zoomed-out roofs. 3 5 4 average over a land-cover pyramid produced class 4 … which is in no legend. The file still validates as a COG; only the map is wrong.
Structure can be validated by a tool; kernel choice cannot, which is why a classified raster with an averaged pyramid passes every check and still ships invented classes.

4. Choose compression and predictor together. The predictor is a reversible pre-transform that makes neighbouring values compress better, and it is typed: predictor=2 (horizontal differencing) for integers, predictor=3 (floating-point differencing) for float32/float64, predictor=1 for byte imagery and anything already compressed. Setting 2 on floats gains nothing and may raise; setting 3 on integers is invalid.

DEFLATE with the right predictor is the safe default — every GDAL, browser-side reader and cloud tiler in existence decodes it. ZSTD compresses slightly better and writes two to three times faster, which matters when you are re-encoding thousands of scenes. LERC is the specialist: it is designed for continuous rasters and takes a max_z_error bound, so you can trade a guaranteed maximum per-pixel error for a large size reduction on elevation or reflectance. max_z_error=0 is lossless.

Compression level is worth far less tuning than it looks. Pushing ZSTD_LEVEL from 9 to 22 typically buys a few percent for several times the write time, and ZLEVEL above 9 on DEFLATE behaves the same way; the predictor and the block size move the number much further. The overviews inherit the main image's codec, so whatever you choose here is also what every zoomed-out read has to decode — one more reason to keep the archive on something universally supported and reserve LERC for the products whose error budget you actually control.

Relative file size, write cost and reader reach for four compression settings Four horizontal bars for the same float32 elevation model. DEFLATE with predictor 3 is the baseline at 100 percent size, 1.0 times write cost, decodable by any GDAL. ZSTD level 9 with predictor 3 lands at 92 percent and roughly three times faster writes but needs GDAL 2.3 or newer. LERC with max_z_error zero, still lossless, lands at 65 percent at roughly twice the write speed and needs GDAL 3.4 or newer. LERC with max_z_error of one centimetre lands at 30 percent, lossy within that hard bound. A note reminds that predictor 2 suits integers, 3 suits floats and 1 suits byte imagery or already-compressed data. Same float32 DEM, four codec settings relative file size write cost reader support DEFLATE · predictor=3 100% 1.0× baseline any GDAL ZSTD level 9 · predictor=3 92% ≈3× faster GDAL ≥ 2.3 LERC · max_z_error=0 65% · still lossless ≈2× faster GDAL ≥ 3.4 LERC · max_z_error=0.01 30% · lossy within ±1 cm ≈2× faster GDAL ≥ 3.4 predictor=2 for integer bands · predictor=3 for float bands · predictor=1 for byte imagery and already-compressed data Indicative on one 20 000 × 20 000 float32 elevation model. LERC’s bound is per pixel, so ±1 cm is a guarantee, not an average.
DEFLATE buys reach, ZSTD buys throughput, LERC buys size on continuous data — and only LERC lets you name the error you are willing to accept.

5. Write the COG in one pass with the COG driver. GDAL 3.1+ ships a driver that handles tiling, pyramid construction and byte ordering itself. Drive it through rasterio.shutil.copy, which issues a straight CreateCopy from the source dataset.

from rasterio.shutil import copy as rio_copy

rio_copy(
    src_path,
    "dem_cog.tif",
    driver="COG",
    blocksize=512,
    compress="DEFLATE",
    predictor="YES",              # COG driver: YES/NO/STANDARD/FLOATING_POINT
    overview_resampling="AVERAGE",
    num_threads="ALL_CPUS",
    bigtiff="IF_SAFER",
)

Two details save an afternoon here. PREDICTOR on the COG driver takes a named value, not the numeric GTiff one — "YES" resolves to horizontal differencing for integers and floating-point differencing for floats, so it is the portable spelling. And prefer rio_copy over rasterio.open("dem_cog.tif", "w", driver="COG", ...): the COG driver is copy-only, so Rasterio silently falls back to a buffered writer that holds the entire dataset in memory until close, which a 20 000² raster will not survive.

6. When you need an explicit factor set, build the pyramid yourself and hand it over. Use this when the array is produced by your own pipeline, when the overviews are expensive enough to build exactly once, or when the kernel differs from what you would want the driver to guess — a classified raster resampled with mode being the standard case.

import rasterio
from rasterio.enums import Resampling
from rasterio.io import MemoryFile
from rasterio.shutil import copy as rio_copy

with rasterio.open("landcover_utm33n.tif") as src:
    classes = src.read(1)                       # uint8 class codes
    profile = src.profile.copy()

profile.update(
    driver="GTiff", tiled=True, blockxsize=512, blockysize=512,
    compress="DEFLATE", predictor=2,            # numeric here — GTiff, integers
)
factors = overview_factors(profile["width"], profile["height"])

with MemoryFile() as memfile:
    with memfile.open(**profile) as mem:
        mem.write(classes, 1)
        mem.build_overviews(factors, Resampling.mode)     # majority, never a mean
        mem.update_tags(ns="rio_overview", resampling="mode")
        rio_copy(
            mem,
            "landcover_cog.tif",
            driver="COG",
            overviews="FORCE_USE_EXISTING",     # reuse the pyramid, do not rebuild
            blocksize=512,
            compress="DEFLATE",
            predictor="YES",
            num_threads="ALL_CPUS",
        )

MemoryFile keeps the intermediate in RAM, so use it only for rasters that fit; for anything larger write the tiled GTiff to a scratch path and copy from there. The band-indexing rules for multi-band sources are covered in reading multi-band TIFFs with Rasterio — read and write band by band rather than materialising a full (bands, rows, cols) stack.

7. Or drive the whole thing from rio-cogeo for batch work. The CLI is the shortest correct path when you are re-encoding a directory, and its Python API takes a profile dict you can version-control.

rio cogeo create dem_utm33n.tif dem_cog.tif \
    --cog-profile deflate \
    --overview-resampling average \
    --blocksize 512
from rio_cogeo.cogeo import cog_translate
from rio_cogeo.profiles import cog_profiles

dst_profile = cog_profiles.get("zstd")          # deflate | zstd | lerc | lzw | webp
dst_profile.update(blockxsize=512, blockysize=512, predictor=3)   # float32

cog_translate(
    "dem_utm33n.tif",
    "dem_cog_zstd.tif",
    dst_profile,
    overview_resampling="average",
    quiet=False,
)

Verification

Two things need proving: that the file satisfies the layout rules, and that the pyramid you asked for is the pyramid that landed. The validator answers the first, Rasterio answers the second.

rio cogeo validate dem_cog.tif
# dem_cog.tif is a valid cloud optimized GeoTIFF

rio cogeo info dem_cog.tif
# IFD
#     Id      Size            BlockSize     Decimation
#     0       20000x20000     512x512       0
#     1       10000x10000     512x512       2
#     ...
import numpy as np
import rasterio
from rasterio.windows import Window
from rio_cogeo.cogeo import cog_validate

is_valid, errors, warnings = cog_validate("dem_cog.tif")
assert is_valid, errors
print("warnings:", warnings)                    # warnings: []

with rasterio.open("dem_cog.tif") as cog:
    assert cog.block_shapes[0] == (512, 512), cog.block_shapes[0]

    factors = cog.overviews(1)
    print("overview factors:", factors)         # overview factors: [2, 4, 8, 16, 32, 64]
    assert factors, "no pyramid — every zoomed-out read decodes full resolution"

    coarsest_px = max(cog.width, cog.height) / factors[-1]
    assert coarsest_px <= 512, f"coarsest level is {coarsest_px:.0f} px — add a level"

    # Lossless codec: the full-resolution pixels must survive the re-encode
    win = Window(10000, 10000, 512, 512)
    with rasterio.open(src_path) as src:
        assert np.array_equal(src.read(1, window=win), cog.read(1, window=win))
print("dem_cog.tif: valid COG, 6 overview levels, pixels unchanged")

With a lossy max_z_error, replace the equality check with the bound you paid for — assert np.abs(src.read(1, window=win) - cog.read(1, window=win)).max() <= 0.01 — which is the whole point of LERC over a generic lossy codec.

Edge Cases & Debugging

Frequently Asked Questions

How much does the overview pyramid add to the file? About a third. Each level holds a quarter of the pixels of the one above, so the series converges to 1/3 of the full-resolution payload regardless of how deep you go. That makes pyramid depth effectively free: stopping at 4× to "save space" saves under 5% of the file while forcing every client to decode a 5 000-pixel image for a thumbnail.

Can I add overviews to an existing COG without rewriting it? No — not while keeping it a COG. Overview directories have to sit ahead of the full-resolution pixel data, and there is no room to insert them into a finished file, so GDAL appends them and the layout stops validating. Re-copy the file through the COG driver instead; on a lossless codec the operation is a decode-and-re-encode with no data change, which the pixel-equality assertion in the verification step confirms.

Should the overviews live in a sidecar .ovr file? Not for a COG. External .ovr pyramids are a fine local optimisation, but a remote client that only has a URL to the .tif will never find them, and rio cogeo validate will report the file as having no overviews at all. Keep everything internal — that is the format's entire premise.

Does a COG have to be in Web Mercator to serve a web map? No. Tile servers reproject on read, so a COG in its native UTM zone serves web tiles perfectly well, and keeping the archive in a metric projected CRS preserves the equal-ground-spacing that analysis needs. Use rio cogeo create --web-optimized only for a display copy whose blocks must align exactly with the EPSG:3857 tile grid, and keep it alongside the UTM original rather than replacing it.

Related