Spatial Epi Public Health GIS

Writing Cloud-Optimized GeoTIFFs for Exposure Rasters

This guide, part of Cloud-Native Geospatial Storage for Surveillance, covers writing the raster products from the exposure section so that the thousands of small reads an assignment job performs are cheap rather than catastrophic.

Problem Context & Constraints

An exposure surface has an unusual access profile. It is written once and read enormously: a cohort assignment samples it at a million points, a zonal statistic reads a few thousand windows, and a map viewer requests overview tiles. Every one of those is a small, scattered read.

A striped GeoTIFF — the historical default — stores data in full-width rows. Reading a 3-by-3 window at a point therefore requires reading three complete rows of the state, which at 30 m resolution is several megabytes for nine pixels. Over a million points that is the difference between minutes and days, and over network storage it is the difference between feasible and not.

Three properties make a raster cloud-optimized in practice. Internal tiling, so a window read fetches a block rather than a row. Internal overviews, so a viewer drawing a thumbnail reads a downsampled copy rather than the full resolution. And a header layout that places all metadata at the front, so a reader learns the file’s structure in one request rather than seeking through it.

The fourth constraint is specific to exposure work: the prediction and its variance belong in the same file as separate bands, because separating them guarantees that somewhere downstream only one of the two gets carried.

What a Point Sample Costs in Each Layout Two representations of the same raster. In the striped layout, data is stored as full-width rows, so sampling a small window near the centre requires reading three complete rows spanning the entire state, shown as three wide bands. In the tiled layout, the same sample reads a single square block containing the point, shown as one small square. The caption gives the byte cost of each for a thirty-metre state-wide raster, several megabytes against a few kilobytes, and multiplies both by a million sample points. The same sample, two layouts, a thousandfold difference in bytes Striped 3 full rows read: about 4 MB for 9 pixels of data Tiled 1 block read: about 4 kB for the same 9 pixels

Prerequisites

  • python 3.11
  • rasterio 1.3.10, numpy 1.26.4
  • The exposure surface and its prediction variance as arrays or an intermediate raster
  • A declared projected coordinate reference system and a nodata value
  • Object storage that supports ranged reads, which is where the layout pays off

Step-by-Step Solution

Write the file tiled with a sensible block size, build overviews internally, declare nodata and band descriptions, and then validate that the result is genuinely cloud-optimized rather than merely tiled.

# Write and validate a cloud-optimized exposure raster with prediction and variance.
# Pinned: rasterio==1.3.10, numpy==1.26.4
import logging

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.shutil import copy as rio_copy

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("storage.cog")

BLOCK = 512
NODATA = -9999.0


def write_exposure_cog(prediction: np.ndarray, variance: np.ndarray,
                       transform, crs, dst_path: str,
                       clear_count: np.ndarray | None = None,
                       tags: dict | None = None) -> None:
    """Write prediction, variance and an optional support band as one COG."""
    if prediction.shape != variance.shape:
        raise ValueError("prediction and variance grids differ")
    bands = [("prediction", prediction), ("prediction_variance", variance)]
    if clear_count is not None:
        bands.append(("support_count", clear_count.astype("float32")))

    profile = {
        "driver": "GTiff", "height": prediction.shape[0], "width": prediction.shape[1],
        "count": len(bands), "dtype": "float32", "crs": crs, "transform": transform,
        "nodata": NODATA, "tiled": True, "blockxsize": BLOCK, "blockysize": BLOCK,
        "compress": "deflate", "predictor": 2, "BIGTIFF": "IF_SAFER",
    }
    tmp = dst_path + ".tmp.tif"
    with rasterio.open(tmp, "w", **profile) as dst:
        for i, (name, arr) in enumerate(bands, start=1):
            a = np.where(np.isfinite(arr), arr, NODATA).astype("float32")
            dst.write(a, i)
            dst.set_band_description(i, name)
        dst.build_overviews([2, 4, 8, 16, 32], Resampling.average)
        dst.update_tags(ns="rio_overview", resampling="average")
        if tags:
            dst.update_tags(**{k: str(v) for k, v in tags.items()})

    # Rewrite through the COG driver so the header and overviews land in the
    # required order; a tiled file with overviews is NOT automatically a COG.
    rio_copy(tmp, dst_path, driver="COG", compress="DEFLATE", predictor="YES",
             blocksize=BLOCK, overview_resampling="average", BIGTIFF="IF_SAFER")
    log.info("wrote COG %s bands=%s block=%d", dst_path,
             [b[0] for b in bands], BLOCK)


def validate_cog(path: str) -> dict:
    """Assert the properties that actually make windowed reads cheap."""
    problems = []
    with rasterio.open(path) as src:
        prof = src.profile
        if not prof.get("tiled", False):
            problems.append("not tiled: every window read fetches full rows")
        bs = (prof.get("blockxsize"), prof.get("blockysize"))
        if bs[0] != bs[1] or not bs[0] or bs[0] < 128:
            problems.append(f"unsuitable block size {bs}")
        if not src.overviews(1):
            problems.append("no internal overviews: a thumbnail reads full resolution")
        if src.nodata is None:
            problems.append("nodata not declared")
        if src.crs is None:
            problems.append("no CRS")
        missing_desc = [i for i in range(1, src.count + 1)
                        if not src.descriptions[i - 1]]
        if missing_desc:
            problems.append(f"bands without descriptions: {missing_desc}")
        info = {"tiled": prof.get("tiled"), "block": bs,
                "overviews": src.overviews(1), "bands": src.count,
                "descriptions": list(src.descriptions), "nodata": src.nodata,
                "crs": src.crs.to_string() if src.crs else None}
    if problems:
        log.error("COG validation failed for %s: %s", path, "; ".join(problems))
    else:
        log.info("COG validated: %s", info)
    return {"ok": not problems, "problems": problems, "info": info}

The rewrite through the dedicated driver is the step most implementations skip. A file that is tiled and has overviews still is not cloud-optimized if its header and overview blocks are scattered through the file, because a reader then needs several round trips to discover the structure before reading any data. The driver lays them out so one ranged request at the start tells the reader everything it needs.

Carrying prediction, variance and a support band together is what keeps the downstream analysis honest. The assignment code in Zonal Statistics for Raster Exposure at Address Points reads band descriptions rather than positions, so a file without them is a file whose bands will eventually be swapped by somebody.

Header Placement Decides How Many Round Trips a Read Costs Two file layouts drawn as horizontal bars. In the cloud-optimized layout, the header and the internal overviews sit at the front of the file followed by full-resolution tiles, so a reader learns the structure and can fetch a thumbnail in one or two ranged requests. In the scattered layout, the header is at the front but overviews are appended at the end and tile offsets are interleaved, so a reader makes several round trips before any pixel arrives. The note records that both files are tiled and have overviews, which is why tiling alone is not the test. Tiled with overviews is necessary, not sufficient Cloud-optimized layout header overviews full-resolution tiles thumbnail: 2 ranged requests Scattered layout header tiles overviews more tiles thumbnail: 5 or more round trips before the first pixel

Validation & Edge Cases

Validation that only checks tiling. A file can be tiled, have overviews, and still require many round trips. Validate the layout with a dedicated checker, and confirm by timing a windowed read over the network rather than locally.

Nodata confused with a real value. A sentinel of −9999 in a float band is fine when declared and dangerous when not. Declare it, and additionally range-check the band on write so a genuine value cannot collide with the sentinel.

Compression that costs more than it saves. A high compression level slows every read for a modest size reduction. Deflate with a horizontal predictor suits continuous surfaces; measure decompression time before choosing something more aggressive.

Overviews built with the wrong resampling. Nearest-neighbour overviews of a continuous surface produce visibly wrong thumbnails; average resampling is right for exposure values and wrong for a categorical band. Where a file mixes both, build separate products.

Very large rasters and the classic format limit. A state-wide 10 m multi-band float raster can exceed the four-gigabyte offset limit. Setting the big-file flag conditionally avoids both the failure and the incompatibility of forcing it unnecessarily.

Band order assumed downstream. Always set band descriptions, and read by description rather than index in every consumer.

INFO storage.cog wrote COG /archive/exposure/pm25_2026.tif bands=['prediction', 'prediction_variance', 'support_count'] block=512
ERROR storage.cog COG validation failed for /archive/exposure/ndvi_2026.tif: no internal overviews: a thumbnail reads full resolution; bands without descriptions: [2]

The overview resampling choice in that list is visible in a thumbnail, which is exactly where users first judge whether a surface looks right:

Overview Resampling Changes What Users See First Two thumbnails of the same continuous exposure surface at a coarse overview level. Built with nearest-neighbour resampling, the thumbnail picks one source pixel per block and shows speckle and false hot spots where an isolated high cell happened to be chosen. Built with average resampling, the thumbnail shows the smooth gradient the full-resolution surface actually has. Average suits continuous exposure bands; nearest neighbour suits only categorical bands such as land cover. Average for continuous bands; nearest only for categories Nearest neighbour speckle and false hot spots Average the gradient the surface actually has same surface, same overview level, different resampling mixing continuous and categorical bands needs separate products

Compliance Notes

Write provenance into the file’s own tags: the run identifier, the configuration digest, the source dataset vintages, the model parameters and the validation statistics. A raster is the artifact most likely to be copied out of a system and sent to a partner, and tags travel with it while a sidecar metadata file does not.

Exposure rasters themselves contain no protected data and should be published where policy allows, which makes them the easiest part of the system to share and the most useful to other agencies. The one caution is resolution: a very fine surface published alongside an areal release can help an adversary narrow a residence to a contour band, so the publication decision belongs with the disclosure review for the products built on it rather than being automatic.