Writing and Validating Cloud-Optimized GeoTIFFs

Reading a Cloud-Optimized GeoTIFF efficiently is a property of the file, not of the reader. Every technique in Understanding Cloud-Optimized GeoTIFF Structure — ranged reads, overview-backed previews, header-only inventories — depends on choices made at write time, and a file written carelessly cannot be read cleverly. This topic covers those choices and the validation pass that proves they were made, within the broader context of Core Raster Fundamentals & STAC Mapping.

The problem is specific: a pipeline produces derived rasters — composites, indices, masks, mosaics — and those outputs become somebody else’s inputs. If they are written as untiled, overview-free GeoTIFFs, every downstream windowed read degrades into a whole-file read, every preview decimates full-resolution data, and the egress bill quietly triples. The write step is where that is decided.

Prerequisites

pip install "rasterio>=1.3.0" "rio-cogeo>=5.0" "numpy>=1.23"

Conceptually, you need the byte layout from Understanding Cloud-Optimized GeoTIFF Structure, the codec trade-off from Choosing COG Compression: ZSTD vs DEFLATE, and the resampling rules from Handling Pixel Resolution and Scaling, because overview construction is a resampling operation like any other.

What “cloud-optimized” adds to a plain GeoTIFF

A GeoTIFF becomes cloud-optimized through three properties, and all three are independent decisions that a naive write gets wrong.

Tiling makes pixels independently addressable. Overviews make zoomed-out reads proportional to what is displayed rather than to what is stored. Layout ordering makes the metadata reachable in one request instead of several. A file can have any subset of the three, which is why “is it a COG?” is a validation question rather than a naming convention.

Step-by-step workflow

1. Decide the profile from the read pattern

Start from how the file will be consumed. Analysis pipelines read moderate windows and benefit from 512-pixel tiles; tile servers read 256-pixel map tiles and benefit from matching that. Archives read by unknown consumers should use DEFLATE; internal intermediates can use ZSTD and save write time.

import rasterio
from rasterio.enums import Resampling

COG_PROFILE = {
    "driver": "GTiff",
    "tiled": True,
    "blockxsize": 512,        # square tiles; must match blockysize
    "blockysize": 512,
    "compress": "deflate",    # widest reader support
    "predictor": 2,           # horizontal differencing: integer data only
    "interleave": "pixel",
    "BIGTIFF": "IF_SAFER",    # avoids the 4 GB limit without forcing BigTIFF
}

The predictor is worth a sentence because it is free compression that is easy to get wrong: predictor=2 (horizontal differencing) helps integer data substantially, predictor=3 is the floating-point variant, and applying the integer predictor to float data hurts both size and speed.

2. Write tiled pixel data

import numpy as np
import rasterio


def write_tiled(dst_path: str, array: np.ndarray, profile: dict) -> None:
    """Write a single-band array as a tiled GeoTIFF using a COG-ready profile."""
    out = profile.copy()
    out.update(
        height=array.shape[-2],
        width=array.shape[-1],
        count=1,
        dtype=array.dtype,
    )
    with rasterio.open(dst_path, "w", **out) as dst:
        dst.write(array, 1)

Any profile that omits tiled=True produces strips regardless of the block sizes given, and strips are what force a windowed read to pull entire rows — the effect measured in Optimizing rasterio Window Reads for Memory Efficiency.

3. Build overviews with the right resampling

import rasterio
from rasterio.enums import Resampling


def add_overviews(path: str, *, categorical: bool = False) -> None:
    """Append decimated overview levels sized so the coarsest fits one screen."""
    with rasterio.open(path, "r+") as dst:
        # Levels chosen so the smallest overview is roughly 300-500 px on a side
        factors = [2, 4, 8, 16, 32]
        method = Resampling.nearest if categorical else Resampling.average
        dst.build_overviews(factors, method)
        dst.update_tags(ns="rio_overview", resampling=method.name)

average is correct for reflectance and elevation; nearest is the only defensible choice for class maps, because averaging class codes produces classes that do not exist. The same rule governs every other resampling decision in the pipeline, as set out in Choosing the Right Resampling Method for Sentinel-2.

4. Reorder into COG layout

Writing tiles and appending overviews leaves the overviews at the end of the file. The COG layout requires them at the front, which means a second pass that copies the file in the right order. GDAL’s COG driver does this in one step:

import rasterio
from rasterio.shutil import copy as rio_copy


def to_cog(src_path: str, dst_path: str, *, overview_resampling: str = "average") -> None:
    """Rewrite a tiled GeoTIFF into true COG layout, header and overviews first."""
    rio_copy(
        src_path,
        dst_path,
        driver="COG",
        compress="DEFLATE",
        predictor="YES",              # COG driver picks 2 or 3 by dtype
        blocksize=512,
        overview_resampling=overview_resampling,
        BIGTIFF="IF_SAFER",
    )

The COG driver builds overviews itself when the source has none, so for a fresh output you can skip step 3 entirely and let the driver do both. Keep step 3 when you need per-level control the driver does not expose.

5. Validate before publishing

import rasterio
from rio_cogeo.cogeo import cog_validate


def assert_cog(path: str, *, min_levels: int = 3) -> None:
    """Fail loudly if a file is not a usable COG."""
    is_valid, errors, warnings = cog_validate(path)
    if not is_valid:
        raise ValueError(f"{path} is not a valid COG: {errors}")

    with rasterio.open(path) as src:
        block = src.block_shapes[0]
        if block[0] != block[1] or block[0] < 128:
            raise ValueError(f"{path} has unusable block shape {block}")
        levels = src.overviews(1)
        if len(levels) < min_levels:
            raise ValueError(f"{path} has only {len(levels)} overview level(s)")
    print(f"{path}: valid COG, blocks {block}, overviews {levels}, warnings={len(warnings)}")

Validation belongs in the pipeline, not in a person’s habits. A single assertion at the write boundary is what prevents an entire archive from being republished six months later.

Parameter reference

The cost of each write choice

Write-time settings trade three quantities against each other: file size, write time and read latency. The table above says what to set; the chart says what it costs.

The asymmetry in that last line is the whole argument for spending write time. A pipeline that writes ten thousand outputs pays the compression cost ten thousand times; a downstream that reads them a hundred times each pays the read cost a million times.

Verification and testing

Three assertions cover almost every real regression, and all three are cheap enough to run on every output:

import rasterio


def audit(path: str) -> dict:
    """Return the structural facts that decide whether a file reads efficiently."""
    with rasterio.open(path) as src:
        return {
            "tiled": src.profile.get("tiled", False),
            "block": src.block_shapes[0],
            "overviews": src.overviews(1),
            "compress": src.profile.get("compress"),
            "dtype": src.dtypes[0],
            "nodata": src.nodata,
            "crs": str(src.crs),
        }


facts = audit("composite_2023.tif")
assert facts["tiled"], "not tiled — every window read becomes a strip read"
assert len(facts["overviews"]) >= 3, "too few overview levels for wide zoom-outs"
assert facts["nodata"] is not None, "no declared fill value"

The fourth check needs a network: open the published object over HTTPS and confirm that reading its profile costs the number of requests you expect. A file that validates locally but triggers a dozen requests on open has a layout problem the local validator will not see, and the settings that control that are covered in How to Read COG Headers Without Downloading Full Files.

Troubleshooting

blockxsize is set but the file is still striped

tiled=True was omitted. GDAL accepts the block-size keywords silently and writes strips. Check src.profile["tiled"] on the output rather than trusting the write call.

Overviews exist but a preview is still slow

The overview levels stop too early — often at 1/4 — so a full-extent preview falls back to decimating a much larger level. Extend the factor list until the coarsest level is a few hundred pixels on a side.

predictor=2 made the file larger

The data is floating point. Use predictor=3 for float32 and float64, or drop the predictor entirely; the integer predictor de-correlates the wrong thing and defeats the compressor.

The COG driver refuses a float32 array with predictor=YES

Older GDAL builds map YES to the integer predictor. Set the predictor explicitly by dtype rather than relying on the automatic mode.

A valid COG reads slowly from object storage

Structure is not the problem: check the region. Bytes crossing a region boundary are billed and slow regardless of layout, as quantified in Reducing S3 Egress Costs in Raster Pipelines.

Writing COGs from a labelled-array pipeline

Pipelines that compute in xarray hit one extra wrinkle: rio.to_raster writes a plain GeoTIFF unless it is told otherwise, and the array’s chunk layout leaks into the file’s block layout.

import rioxarray  # noqa: F401  (registers the .rio accessor)
import xarray as xr


def write_cog_from_dataarray(da: xr.DataArray, dst_path: str) -> None:
    """Write a chunked DataArray straight to COG layout."""
    # Chunks become blocks: rechunk to the tile size you actually want on disk
    da = da.chunk({"y": 512, "x": 512})
    da.rio.to_raster(
        dst_path,
        driver="COG",
        compress="DEFLATE",
        blocksize=512,
        overview_resampling="average",
        num_threads="ALL_CPUS",
        BIGTIFF="IF_SAFER",
    )

Two details decide whether the output is usable. The chunk shape must be a multiple of the block size, or GDAL writes partial blocks and the file ends up with ragged tiling that defeats windowed reads. And the array must carry its CRS and nodata on the rio accessor before the write, because anything the accessor does not know about is simply absent from the output — a file with correct pixels and no georeferencing is the most common product of a rushed xarray write. The alignment rules behind that are covered in Band Math Operations with xarray.

For very large outputs, writing per-window rather than materialising the whole array keeps memory flat. Open the destination once, iterate dst.block_windows(), and write each window as it is computed; the pattern is the write-side mirror of the read strategy in Optimizing rasterio Window Reads for Memory Efficiency.

Metadata that belongs in the file, not in a wiki

A structurally perfect COG with no semantic metadata is still a liability, because the next reader has to guess. Four things are cheap to write and expensive to reconstruct later.

The fill value comes first. Declaring nodata costs one profile key and determines whether every downstream statistic is computed over real pixels or over a mixture of pixels and fill. A file without it forces each consumer to invent a convention, and they will not all invent the same one — the failure modes are catalogued in Extracting nodata and dtype from a GeoTIFF.

Scaling comes second. If the array holds scaled integers, write the scale and offset as band attributes rather than relying on a product specification the consumer may not have. If the array holds physical units already, say so in a tag, because “float32 reflectance” and “float32 scaled reflectance” are indistinguishable from the numbers alone when the scene is dark.

Provenance comes third and is the one most often skipped. A handful of tags — the source item identifiers, the processing version, the masking configuration, the resampling method used for overviews — turns an opaque output into one that can be reproduced. dst.update_tags(...) accepts arbitrary key-value pairs and costs nothing in file size.

Band descriptions come fourth. dst.set_band_description(1, "ndvi_median") is what makes a multi-band composite readable by a human six months later, and what stops the second band being misinterpreted as the first in a downstream stack.

None of this is enforced by any validator, which is exactly why it needs to be in the writing function rather than in a checklist. A write_cog helper that always sets nodata, always writes provenance tags and always names its bands makes correctness the default, and every pipeline that imports it inherits the discipline for free.

---

Frequently Asked Questions

Q: Is a tiled GeoTIFF with overviews already a COG? Not necessarily. A COG also requires the layout to be ordered so that the header, the IFDs and the overviews precede the full-resolution tile data. A tiled file with overviews appended at the end reads correctly but forces extra range requests, which is why validation checks ordering rather than only structure.

Q: What tile size should I use? 512 by 512 is a good default for analysis workloads and 256 by 256 for tile-server workloads. Smaller tiles mean more, smaller range requests; larger tiles mean fewer requests that each transfer more than the reader needed. Match the tile size to the typical read window rather than to the file size.

Q: How many overview levels do I need? Enough that the coarsest level fits in roughly one screen — for a 10980 pixel Sentinel-2 band that means levels 2, 4, 8, 16 and 32. Stopping early leaves wide zoom-outs decimating full-resolution data, which is the exact cost overviews exist to remove.

Q: Does adding overviews change the pixel values? It adds new, decimated copies alongside the original data; full-resolution pixels are untouched. The overviews themselves depend entirely on the resampling method, so a class map built with an averaging method will contain class codes that never existed in the source.

---

Related

• Understanding Cloud-Optimized GeoTIFF Structure — the byte layout these write settings produce, and why readers care.
• Choosing COG Compression: ZSTD vs DEFLATE — the codec decision behind the compress argument.
• Converting a GeoTIFF to a COG with rio-cogeo — the one-command version of this workflow.
• Validating COG Structure in CI — turning the assertions above into a build gate.
• Extracting and Parsing Raster Metadata — auditing these properties across an existing archive.