Watershed Modeling Python for hydrology automation

Writing Cloud-Optimized GeoTIFFs for Hydrology Pipelines

A Cloud-Optimized GeoTIFF is an ordinary GeoTIFF with its bytes arranged so a client can fetch the part it needs without reading the rest. That property is what makes chunked and tiled hydrologic processing viable, and it is entirely a matter of how the file was written. This guide covers the write options that matter for elevation and hydrology rasters, as part of the tiled and parallel watershed processing topic within watershed delineation and catchment synchronization.

Prerequisites

  • rasterio built against GDAL 3.1 or later, which includes the COG driver.
  • rio-cogeo is optional but its validator is worth having.
  • Knowledge of the raster’s dtype and semantics — the compression choice differs sharply between a float elevation surface and an integer direction grid.

Core Technique: Internal Tiling and Byte Order

Three properties make a GeoTIFF cloud-optimized, and only the first two matter for local pipeline performance.

What a Windowed Read Costs in Each Layout A stripe-oriented file stores each full-width row group contiguously, so reading a small square window transfers every row it touches across the whole raster width. A tiled file stores 512 by 512 blocks, so the same window transfers only the four blocks it overlaps. stripe-oriented the window is 52 × 32 the transfer is two full-width strips ≈ 20× more bytes than needed tiled, 512 × 512 blocks the same window transfers four blocks ≈ 3× more bytes than needed

Internal tiling stores the raster as independently compressed blocks rather than full-width strips. This is the property that makes a windowed read cheap, and it matters just as much on local disk as on object storage.

Header-first layout puts the image file directory and the block offsets at the front, so a client can learn where every block lives in one small read. This matters over HTTP and is irrelevant on a local filesystem.

Overviews hold decimated copies for fast display. They cost roughly a third of the file size and are ignored entirely by full-resolution processing.

Annotated Code Example

python 
import logging

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

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

# Predictor 3 is the floating-point predictor: it differences adjacent values
# before compressing, which is extremely effective on elevation surfaces where
# neighbours differ by centimetres. Predictor 2 is the integer equivalent.
PREDICTOR_BY_KIND = {
    "elevation": 3,     # float32/float64 continuous surface
    "direction": 2,     # uint8 D8 codes
    "accumulation": 3,  # float32/float64, huge dynamic range
    "labels": 2,        # int32 basin ids
    "mask": 2,          # uint8 boolean
}


def write_hydrology_cog(
    array: np.ndarray,
    profile: dict,
    out_path: str,
    kind: str = "elevation",
    blocksize: int = 512,
    compress: str = "DEFLATE",
    build_overviews: bool = False,
) -> dict:
    """
    Write a hydrology raster as a Cloud-Optimized GeoTIFF.

    Parameters
    ----------
    array           : The raster data, 2-D.
    profile         : A rasterio profile carrying crs, transform and nodata.
    out_path        : Destination path.
    kind            : One of PREDICTOR_BY_KIND; selects the compression predictor.
    blocksize       : Internal block size; 512 is the right default for DEMs.
    compress        : DEFLATE, ZSTD or LZW. Never a lossy codec.
    build_overviews : Only worth it if the file will be viewed, not just processed.
    """
    if compress.upper() in {"JPEG", "WEBP"}:
        raise ValueError(f"{compress} is lossy — a lossy DEM produces a different "
                         "watershed, silently")
    if kind not in PREDICTOR_BY_KIND:
        raise ValueError(f"unknown raster kind {kind!r}")

    predictor = PREDICTOR_BY_KIND[kind]
    prof = profile.copy()
    prof.update(
        driver="GTiff",
        count=1,
        dtype=array.dtype,
        tiled=True,
        blockxsize=blocksize,
        blockysize=blocksize,
        compress=compress,
        predictor=predictor,
        # BIGTIFF is needed above 4 GB; IF_SAFER lets GDAL decide, which avoids
        # a failed write at the very end of a long job.
        BIGTIFF="IF_SAFER",
        num_threads="ALL_CPUS",
    )

    tmp_path = out_path + ".tmp.tif"
    with rasterio.open(tmp_path, "w", **prof) as dst:
        dst.write(array, 1)
        if build_overviews:
            # Powers of two down to roughly a 256-pixel thumbnail. Average is
            # right for continuous surfaces; use nearest for categorical data.
            factors = []
            size = max(array.shape)
            f = 2
            while size / f > 256:
                factors.append(f)
                f *= 2
            resampling = (Resampling.nearest if kind in {"direction", "labels", "mask"}
                          else Resampling.average)
            dst.build_overviews(factors, resampling)
            dst.update_tags(ns="rio_overview", resampling=resampling.name)
            log.info("Built %d overview level(s) with %s resampling",
                     len(factors), resampling.name)

    # --- The COG driver rewrites the file with the header first, which is what
    # makes a range-read client efficient. Copying is cheaper than writing
    # through the COG driver directly for a large array. ---
    rio_copy(tmp_path, out_path, driver="COG", compress=compress,
             predictor=predictor, blocksize=blocksize, overview_resampling="average")

    import os
    size_mb = os.path.getsize(out_path) / 1e6
    raw_mb = array.nbytes / 1e6
    os.remove(tmp_path)

    log.info("Wrote %s: %.1f MB from %.1f MB raw (%.1f× compression, "
             "%s predictor %d, %d px blocks)",
             out_path, size_mb, raw_mb, raw_mb / size_mb, compress, predictor, blocksize)
    return {"path": out_path, "size_mb": size_mb, "ratio": raw_mb / size_mb}


def validate_cog(path: str) -> dict:
    """Check the properties a tiled pipeline actually depends on."""
    with rasterio.open(path) as src:
        blocks = src.block_shapes[0]
        tiled = src.profile.get("tiled", False)
        overviews = src.overviews(1)
        report = {
            "tiled": bool(tiled),
            "block_shape": blocks,
            "compression": str(src.compression),
            "overview_levels": len(overviews),
            "dtype": str(src.dtypes[0]),
            "nodata": src.nodata,
        }
    if not report["tiled"]:
        log.error("%s is NOT internally tiled — every windowed read will "
                  "transfer full-width strips", path)
    if report["nodata"] is None:
        log.warning("%s has no nodata value; edge cells will be treated as "
                    "terrain by the conditioning step", path)
    log.info("COG report for %s: %s", path, report)
    return report


# --- Example usage ---
# with rasterio.open("dem_filled.tif") as src:
#     write_hydrology_cog(src.read(1), src.profile, "dem_filled_cog.tif",
#                         kind="elevation", compress="ZSTD")
# validate_cog("dem_filled_cog.tif")

Parameter Reference

Raster dtype Compression Predictor Typical ratio
Elevation (DEM) float32 ZSTD or DEFLATE 3 2.5–4×
Flow direction uint8 DEFLATE 2 6–12×
Flow accumulation float32 ZSTD 3 2–3×
Basin labels int32 DEFLATE 2 10–40×
Stream mask uint8 DEFLATE 2 30–100×

The predictor is where most of the benefit sits and it is the option most often left at its default of 1. On a float32 DEM, predictor 3 typically halves the file relative to the same compressor with no predictor, because neighbouring elevations differ by a few centimetres and differencing them leaves mostly zeros.

The Predictor Does Most of the Work File sizes for one 100 million cell float32 DEM. Uncompressed is 400 megabytes, LZW without a predictor 340, DEFLATE without a predictor 312, LZW with predictor 3 is 190, DEFLATE with predictor 3 is 148 and ZSTD with predictor 3 is 132. uncompressed 400 LZW, no predictor 340 DEFLATE, no predictor 312 LZW, predictor 3 190 DEFLATE, predictor 3 148 ZSTD, predictor 3 132 file size (MB) — one 100 M cell float32 DEM the predictor matters more than the compressor

Compression is CPU work, and on a large write it is most of the wall-clock time. Letting GDAL use every core usually costs nothing and halves the write.

Compression Threads Halve the Write Writing a 100 million cell float32 raster with ZSTD takes 212 seconds on one thread, 118 on two, 71 on four and 54 on eight, after which the write becomes I/O bound and further threads do not help. 212 s 118 s 71 s 54 s 1 thread 2 threads 4 threads 8 threads Beyond eight the write is I/O bound; NUM_THREADS=ALL_CPUS is a safe default. write time

Gotchas and Edge Cases

  • Lossy compression on elevation. JPEG or WEBP on a DEM changes elevations by fractions of a metre, which changes flow directions on flat ground and therefore changes the watershed. There is no acceptable use of a lossy codec here.
  • Predictor 3 on integer data. GDAL will accept it and the result is often larger than no predictor at all. Match the predictor to the dtype.
  • Missing nodata. The conditioning step treats undefined edge cells as terrain and fills against them, producing a rim of spurious drainage around the raster.
  • Block size above 1024. Windowed reads decompress far more than they use. The gain in per-block overhead is negligible by comparison.
  • Overviews on intermediate products. A third more bytes written and read on files no human will ever look at.
  • num_threads left at 1. Compression is the slow part of writing a large raster, and GDAL parallelises it well. ALL_CPUS often halves the write time.
  • Writing then converting in place. Some workflows write a plain GeoTIFF and translate to COG afterwards, doubling disk traffic. Where the array fits in memory, going straight through the COG driver is cheaper.