Raster Hazard Layers & Cloud-Optimized GeoTIFF

A flood model finishes at 03:40 and writes a 4-gigabyte depth grid covering three counties at two-metre resolution. Six forward command posts need it, each working one division, each on an uplink that would take four hours to move the whole file — and the forecast that produced it will be superseded in six. Publishing that raster as a plain GeoTIFF gives every consumer an all-or-nothing choice. Publishing it as a Cloud-Optimized GeoTIFF (COG) lets each post read the division it is working in, at the scale it is displaying, in a few hundred kilobytes.

Problem Framing

Vector layers on this site have a well-developed story for constrained links: clip to the incident footprint, serialise deterministically, ship deltas, and let a seekable format serve a working extent without the archive. Raster hazard products — flood depth grids, fire progression rasters, smoke dispersion surfaces, damage-probability layers — have historically had no equivalent, so they are either shipped whole, downsampled beyond usefulness, or rendered server-side into pictures that discard the values responders need to query.

The failure this causes is specific and quiet. A division supervisor asked “how deep is the water on this block?” needs a value, not a colour, and a pre-rendered picture cannot answer. The workaround — reading the value off a legend — introduces a rounding step nobody records, and it fails completely where the underlying cell was nodata rather than zero. This topic implements the raster half of the core architecture and data standards contract.

Prerequisites

• GDAL 3.1 or newer, which can both write COGs through the COG driver and validate them, along with rasterio built against the same GDAL. The pinned-binary discipline from Dockerized GIS environments applies here as much as to vector work.
• A settled horizontal CRS for the incident, per the coordinate reference system standard. Raster reprojection is lossy in a way vector reprojection is not, so the target CRS should be chosen once and written at production time rather than per consumer.
• An object store or web server that honours HTTP range requests. The entire benefit depends on partial reads; a store that returns whole objects turns a COG back into a plain GeoTIFF with extra steps.
• A declared vertical datum and unit for any depth or elevation product, which is a separate decision from the horizontal CRS and is the one most often left implicit.

What Makes a GeoTIFF Cloud-Optimized

A COG is not a new format. It is a GeoTIFF that satisfies layout constraints a range-reading client can exploit, and every one of those constraints exists to make some read cheap that would otherwise be expensive.

The three properties are independent and all three are needed. Front-loading the header lets a client discover the structure in one small read. Internal tiling means a rectangular window corresponds to contiguous byte ranges rather than to slivers of every scanline it crosses. Internal overviews mean a zoomed-out view reads a small pre-computed image rather than reading full-resolution pixels and throwing most of them away.

Dropping any one of them degrades gracefully in appearance and catastrophically in cost, which is why validation matters more here than the format’s simplicity suggests — a file can be a perfectly valid GeoTIFF, open correctly in every tool, and still force every client to download all of it.

The overview column is the one worth defending in a review. A production pipeline under time pressure will often skip overview generation because it roughly doubles write time and the output looks identical when opened locally. The cost lands entirely on the remote consumer, and it lands hardest on the situational-awareness view — the county-wide picture an operations chief looks at most often becomes the most expensive read in the system.

What Must Survive the Conversion

The habit of treating a raster as an image is where hazard products go wrong, because three of their properties have no equivalent in an image and are silently discarded by tooling that assumes one.

The nodata case deserves the emphasis. A flood model solves over a domain and leaves cells outside it — or cells where the solver failed to converge — as nodata. Convert with a tool that has no concept of nodata, or write a format that cannot express it, and those cells become zero. Zero is a perfectly good depth value meaning dry, so the raster now asserts that unmodelled ground is safe, and it does so in exactly the areas where the model had trouble.

Production Python Implementation

The writer below produces a validated COG from a model output array, preserving nodata, stamping the vertical datum and model run into the file’s own metadata, and building overviews sized to the display scales the field application actually uses.

from __future__ import annotations

import logging
from dataclasses import dataclass, asdict
from pathlib import Path

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

logger = logging.getLogger("incidentgis.hazard_raster")

# Overview factors chosen to match the display scales the field app renders,
# not the generic powers of two. An unused level is pure write cost.
OVERVIEW_FACTORS = (2, 4, 8, 16, 32)


@dataclass(frozen=True)
class HazardProvenance:
    """Everything needed to identify which model run produced this raster."""
    model: str                  # e.g. "HEC-RAS 2D"
    run_id: str                 # the modelling system's own run identifier
    forecast_cycle: str         # ISO 8601 UTC of the forecast this run consumed
    vertical_datum: str         # e.g. "NAVD88"
    value_units: str            # e.g. "metres"
    quantity: str               # e.g. "water_surface_depth"


def write_hazard_cog(
    values: np.ndarray,
    *,
    transform,
    crs,
    nodata: float,
    provenance: HazardProvenance,
    destination: Path,
) -> Path:
    """Write a validated Cloud-Optimized GeoTIFF for a hazard surface.

    Raises rather than degrading: a hazard raster that silently lost its nodata
    mask or its vertical datum is more dangerous than a job that failed.
    """
    if not np.issubdtype(values.dtype, np.floating):
        raise TypeError("hazard values must be floating point to carry nodata")
    if np.isnan(nodata):
        # NaN nodata round-trips badly through several readers and cannot be
        # compared with ==; insist on a sentinel the format can store.
        raise ValueError("use an explicit sentinel nodata value, not NaN")
    if np.any(values[values != nodata] < 0.0):
        raise ValueError("negative depth outside the nodata mask — check solver output")

    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "count": 1,
        "height": values.shape[0],
        "width": values.shape[1],
        "transform": transform,
        "crs": crs,
        "nodata": nodata,
        "tiled": True,
        "blockxsize": 512,
        "blockysize": 512,
        "compress": "deflate",
        "predictor": 3,          # floating-point predictor; lossless
    }

    staging = destination.with_suffix(".staging.tif")
    try:
        with rasterio.open(staging, "w", **profile) as dst:
            dst.write(values.astype("float32"), 1)
            # Provenance lives in the file, not in the filename. A raster that
            # is copied, renamed or re-served must still say which forecast
            # cycle produced it.
            dst.update_tags(**{
                "INCIDENTGIS_" + key.upper(): str(val)
                for key, val in asdict(provenance).items()
            })
            dst.build_overviews(OVERVIEW_FACTORS, Resampling.average)
            dst.update_tags(ns="rio_overview", resampling="average")

        # The COG driver rewrites the staged file into the required layout:
        # header first, tiles contiguous, overviews internal.
        rio_copy(
            staging, destination, driver="COG",
            compress="deflate", predictor="YES",
            overview_resampling="average", blocksize=512,
        )
    except (RasterioIOError, ValueError) as exc:
        logger.error("hazard_cog_write_failed", exc_info=exc,
                     extra={"destination": str(destination)})
        raise
    finally:
        staging.unlink(missing_ok=True)

    _assert_cog_contract(destination, nodata=nodata)
    logger.info("hazard_cog_written", extra={
        "destination": str(destination),
        "run_id": provenance.run_id,
        "forecast_cycle": provenance.forecast_cycle,
    })
    return destination


def _assert_cog_contract(path: Path, *, nodata: float) -> None:
    """Fail closed if the written file lost a property consumers depend on."""
    with rasterio.open(path) as src:
        if src.nodata is None or src.nodata != nodata:
            raise ValueError("nodata not preserved: " + repr(src.nodata))
        if not src.profile.get("tiled", False):
            raise ValueError("not internally tiled — range reads will be useless")
        if not src.overviews(1):
            raise ValueError("no internal overviews — every scale costs full resolution")
        tags = src.tags()
        for field in ("VERTICAL_DATUM", "VALUE_UNITS", "RUN_ID"):
            if not tags.get("INCIDENTGIS_" + field):
                raise ValueError("missing provenance tag " + field)
    logger.info("hazard_cog_contract_ok", extra={"path": str(path)})

The _assert_cog_contract call is the load-bearing part. Every check in it corresponds to a way the file can be valid, open cleanly, and be useless or dangerous to a downstream consumer — and none of them is visible by looking at the raster.

Configuration Reference

Verification and Smoke Test

Validate the output rather than the process. gdalinfo reports the layout, and a deliberate partial read proves the property the format exists for:

# Layout: expects "Block=512x512" and an Overviews line on band 1.
gdalinfo hazard_depth.tif | grep -E 'Block=|Overviews|NoData'

# Prove the range-read path: read one window over HTTP and count the ranges.
CPL_CURL_VERBOSE=YES CPL_VSIL_CURL_USE_HEAD=NO \
  gdal_translate -srcwin 4000 4000 512 512 \
  /vsicurl/https://example.invalid/hazard_depth.tif /tmp/window.tif 2>&1 \
  | grep -c 'Range: bytes'

A COG that needs more than a handful of range requests for a single window is tiled wrongly or has its overviews in a sidecar. Both open fine locally and both defeat the purpose.

Integration With Adjacent Workflows

The published COG is consumed the same way the vector cache is. A forward node reads windows over its uplink while connected and falls back to a clipped local copy when it is not, which is the offline caching pattern applied to raster. The provenance tags feed the same metadata governance gate vector layers pass, and the model run identity is what lets a depth value quoted in an ICS-209 be traced to a forecast cycle.

Troubleshooting

Symptom: every window read downloads the whole file. The file is not internally tiled, or overviews live in a .ovr sidecar the client never requests. Check gdalinfo for Block= with equal dimensions and an Overviews: line under band 1, not a separate file on disk.

Symptom: unmodelled areas render as zero depth. The nodata value was lost in a conversion step, most often by a tool that read the array and wrote a new file without carrying the mask. Assert src.nodata after every write, not only at the end of the pipeline.

Symptom: depths disagree with a partner agency’s product by roughly a metre. Almost always a vertical datum difference rather than a modelling difference. Compare the INCIDENTGIS_VERTICAL_DATUM tags before investigating the models.

Symptom: overviews look blocky and wrong on a categorical hazard layer. average resampling on class codes produces values that are not classes. Rebuild with nearest for any categorical raster.

Symptom: file size roughly triples after conversion. The predictor is wrong for the data type — 3 on integer data or 2 on floats both defeat compression. Match the predictor to the dtype.

Frequently Asked Questions

Why not just render hazard rasters to image tiles on the server? Because a responder needs a value, not a colour. Asked how deep the water is on a given block, a pre-rendered tile can only be read against a legend, which introduces a rounding step nobody records and fails entirely where the underlying cell was nodata rather than zero. A Cloud-Optimized GeoTIFF lets the client fetch the actual depth values for a small window, so the number quoted on the radio is the number the model produced. Server-side rendering remains useful for a basemap layer, but it is a display convenience rather than a substitute for the data.

What actually makes a GeoTIFF cloud-optimized? Three layout guarantees, all of which serve partial reads. The header and image file directory sit at the front of the file so a client learns the structure in one small read. Pixel data is internally tiled rather than stored as scanlines, so a rectangular window corresponds to contiguous byte ranges instead of fragments of every row it crosses. Overviews are stored inside the same file rather than in a sidecar, so a zoomed-out view reads a small pre-computed image. A file can be a perfectly valid GeoTIFF, open in every tool, and satisfy none of these, in which case every client downloads all of it.

What is the most common way a hazard raster is silently corrupted? Losing the nodata mask during a conversion. A model leaves cells outside its solved domain, or where the solver failed to converge, as nodata. A tool that reads the array and writes a new file without carrying the mask turns those cells into zero, and zero is a valid depth meaning dry. The raster then asserts that unmodelled ground is safe, and it does so precisely where the model had difficulty. Asserting that nodata survived every write is a two-line check that prevents an entire class of dangerous output.

Related

• Offline GIS Data Caching Strategies — the clipped local copy a forward node falls back to when its uplink drops.
• FlatGeobuf vs GeoPackage for Offline Caching — the same partial-read argument, made on the vector side.
• Coordinate Reference Systems for Disaster Zones — the horizontal CRS contract a raster inherits, and why vertical datum is a separate decision.
• Emergency Metadata Standards — the lineage gate the provenance tags written here are designed to satisfy.

Up: Core Emergency GIS Architecture & Data Standards