Preserving Nodata When Converting Hazard Rasters

A hydraulic modeller hands over a depth grid at 02:00 and a GIS analyst converts it for the field cache with a one-line rasterio script. The conversion succeeds, the raster opens, the colours look right, and a block where the mesh failed to converge now reads as zero metres of water. At 06:20 a strike team is routed down that block because the depth grid says it is dry.

Root Cause and Operational Impact

A hazard raster carries three distinct cell states — a modelled value, a modelled zero, and no result — and most raster tooling understands two. The nodata mask is what encodes the third, and it is stored as a property of the dataset rather than of the array, so any operation that moves the array without the dataset loses it.

That loss is silent by construction. The sentinel value is an ordinary number in an ordinary floating-point band; once the mask no longer marks it, nothing distinguishes it from data. If the sentinel was -9999 the result is at least conspicuous, but the common defensive move of writing zeros into unmodelled cells “so the raster renders cleanly” produces exactly the failure above: an assertion of dry ground precisely where the model declined to say.

One transect across a flood-depth grid, before and after nodata is lost A cross-section through a flood-depth raster along one street. In the source grid the solver produced depths of 0.4, 1.1 and 1.6 metres over the inundated stretch, zero over ground it solved as dry, and nodata over a block where the mesh failed to converge. After a conversion that dropped the nodata mask, the unconverged block reads as 0.0 metres. The rendered map is identical in both cases over the dry ground and identical in colour over the failed block, because zero and dry share a colour. The result is a map that asserts a block is passable in exactly the place where the model could not say. one transect along a street, through the depth grid source grid 0.00.4 1.11.6 nodata0.0 the mesh did not converge here after a lossy conversion 0.00.4 1.11.6 0.00.0 now indistinguishable from ground the model solved as dry Zero and dry share a colour, so the rendered map is identical — the assertion changed and nothing shows it.

The two transects render identically. Zero and dry share a colour on every depth ramp anyone uses, so no visual review catches it, and the difference only becomes visible when a unit reaches the block. Under NIMS and FEMA traceability expectations, a product that quietly converts “unknown” into “safe” is also unreconstructable after the fact — the audit trail shows a successful conversion.

Tiered Resolution Strategy

Work down these tiers. The first three keep the mask intact; the last two contain the damage when something upstream has already destroyed it.

  1. Never separate the array from its dataset (definitive). Read with masked=True so the mask travels with the values, or copy the source profile wholesale and write through it. Every conversion step must be expressible as “dataset in, dataset out”; a step whose signature takes and returns a bare ndarray is a step where the mask cannot survive.
  2. State the nodata value explicitly at every write. Do not rely on it being inherited. An explicit nodata= on the write profile costs one line and makes the intent auditable in the code rather than implied by the input.
  3. Respect the mask in every resampling operation. Averaging a real depth with a sentinel produces a plausible number derived partly from a cell that had no value. Use mask-aware resampling, or resample the mask separately and reapply it.
  4. Refuse the conversion when the mask cannot be represented (safe default). Some targets — PNG, JPEG, an eight-bit render — have no way to carry nodata. Those are display artefacts, not data products, and the pipeline should mark them as such and refuse to let them feed a routing or assessment step.
  5. Emit an audit record comparing masked-cell counts. Count nodata cells before and after every step and log both. A step that changes the count has either repaired or destroyed something, and either way somebody should know which.
Five conversion paths and whether each preserves the nodata mask Five ways a hazard raster commonly moves between tools, and whether nodata survives. A gdal_translate call with an explicit destination nodata preserves it. A rasterio read-modify-write preserves it only if the profile is copied and the nodata key set explicitly, which is easy to omit. Reading into a plain NumPy array and writing a fresh file loses it, because the array carries no mask. Exporting to PNG or JPEG for a web viewer loses it, since neither format has a nodata concept. Resampling with average over a window containing nodata contaminates real values unless the mask is respected, producing plausible depths derived partly from cells that had none. Only the first path is safe by default; the rest are safe only if someone remembered. where the mask actually gets lost gdal_translate -a_nodata -9999 preserved — the value is stated on the command line safe rasterio read → modify → write(**profile) preserved only if the profile is copied and nodata is set — trivially omitted in a refactor conditional src.read(1) → np.ndarray → new file lost — a bare array carries no mask, so the sentinel becomes an ordinary value lost export to PNG / JPEG for a viewer lost — neither format has a nodata concept, and JPEG also changes the values it keeps lost average resampling across the mask contaminated — a real depth averaged with a sentinel is a plausible number derived from nothing worse

The third row is the one to search the codebase for. src.read(1) returning a plain array is the most natural thing to write in rasterio, appears in every tutorial, and is where the mask is dropped. Its fix is src.read(1, masked=True), which returns a masked array whose fill value can be written back out intact.

The fifth row is subtler and worse. Resampling is not usually thought of as a lossy step for the mask, but an average kernel over a window straddling the mask boundary blends real values with the sentinel, and the output is neither nodata nor a real depth — it is a number with no provenance that will pass every range check you apply to it.

Production Python Implementation

The routine below performs a mask-preserving conversion with an explicit before-and-after audit, and refuses outright when asked to write to a target that cannot carry nodata.

python 
from __future__ import annotations

import logging
from pathlib import Path

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.errors import RasterioIOError

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

# Formats that cannot represent nodata. Writing a hazard product to one of
# these is a rendering step, never a data step.
DISPLAY_ONLY_DRIVERS = frozenset({"PNG", "JPEG", "GIF", "BMP"})


class NodataContractError(RuntimeError):
    """Raised when a conversion would lose or corrupt the nodata mask."""


def convert_preserving_nodata(
    source: Path,
    destination: Path,
    *,
    driver: str = "GTiff",
    scale_factor: float | None = None,
) -> dict[str, int]:
    """Convert a hazard raster, keeping the nodata mask intact and audited.

    Returns the before/after masked-cell counts so the caller can assert on
    them. Raises NodataContractError rather than producing a product whose
    unmodelled cells have become ordinary values.
    """
    if driver.upper() in DISPLAY_ONLY_DRIVERS:
        raise NodataContractError(
            f"{driver} cannot carry nodata; produce it as a render, "
            "not as an input to routing or assessment"
        )

    try:
        with rasterio.open(source) as src:
            if src.nodata is None:
                raise NodataContractError(
                    f"{source} declares no nodata value — refusing to guess one"
                )

            # masked=True keeps the mask attached to the values. A bare
            # src.read(1) here is the single most common way the mask is lost.
            data = src.read(1, masked=True)
            before = int(np.ma.count_masked(data))

            profile = src.profile.copy()
            profile.update(driver=driver, tiled=True,
                           blockxsize=512, blockysize=512)

            if scale_factor is not None:
                # Resample the mask separately and reapply it, so no output
                # cell is an average of real values and the sentinel.
                out_shape = (
                    int(src.height * scale_factor),
                    int(src.width * scale_factor),
                )
                values = src.read(
                    1, out_shape=out_shape, resampling=Resampling.average,
                )
                mask = src.read_masks(1, out_shape=out_shape,
                                      resampling=Resampling.nearest)
                data = np.ma.masked_where(mask == 0, values)
                profile.update(
                    height=out_shape[0], width=out_shape[1],
                    transform=src.transform * src.transform.scale(
                        src.width / out_shape[1], src.height / out_shape[0]
                    ),
                )

            with rasterio.open(destination, "w", **profile) as dst:
                # Writing the filled array plus an explicit nodata keeps the
                # sentinel and the declaration in agreement.
                dst.write(data.filled(src.nodata).astype(profile["dtype"]), 1)
                dst.update_tags(**src.tags())

        with rasterio.open(destination) as check:
            if check.nodata != rasterio.open(source).nodata:
                raise NodataContractError("nodata value changed during conversion")
            after = int(np.ma.count_masked(check.read(1, masked=True)))

    except RasterioIOError as exc:
        logger.error("nodata_conversion_io_failed", exc_info=exc)
        raise

    if scale_factor is None and after != before:
        raise NodataContractError(
            f"masked-cell count changed without resampling: {before} → {after}"
        )

    logger.info("nodata_conversion_ok", extra={
        "source": str(source), "destination": str(destination),
        "masked_before": before, "masked_after": after,
    })
    return {"masked_before": before, "masked_after": after}

The equality assertion on masked-cell counts is what turns a silent class of bug into a loud one. For a straight conversion the count must be identical; if it is not, something read or wrote the sentinel as a value.

What a downstream consumer should do with each of the three cell states A depth raster carries three distinct cell states and each demands different downstream behaviour. A positive depth is a modelled hazard and drives routing closure and evacuation decisions directly. An explicit zero is modelled dry ground and is safe to route across, because the solver considered it and found no water. Nodata is the absence of a result and must propagate as unknown: a route crossing it is not closed and not open, and the correct handling is to flag the segment for reconnaissance rather than to treat it as either. Collapsing nodata into zero destroys the third state, and with it the only signal that would have sent someone to look. three states, three behaviours — the third is the one that gets collapsed depth > 0 depth = 0 nodata modelled hazard modelled dry ground no result at all close the segment drive the evacuation line quote the value on the radio route across it the solver considered it and found no water neither closed nor open flag for reconnaissance propagate as unknown act on it act on it go and look Collapsing nodata into zero deletes the third column — and with it the only signal that would have sent someone. Which is why the mask is not metadata. It is one third of the data.

Validation Checklist

  • Every read that feeds a write uses masked=True, or copies the source profile wholesale.
  • Every write states nodata= explicitly rather than inheriting it.
  • The nodata sentinel is outside the physically valid range for the quantity and is not NaN.
  • Resampling uses a mask-aware path, or resamples the mask with nearest and reapplies it.
  • Masked-cell counts are logged before and after every step, and asserted equal for non-resampling steps.
  • Any export to PNG, JPEG or an eight-bit render is tagged as display-only and cannot feed routing or assessment.
  • A downstream consumer treats nodata as a third state, not as zero — the routing layer flags those segments for reconnaissance.
  • The smoke test includes a fixture with a deliberate unconverged block, and asserts it is still nodata after the full pipeline.

Edge Cases and Gotchas

  • NaN as the sentinel. It cannot be compared with ==, several readers round-trip it inconsistently, and integer bands cannot hold it at all. Use an explicit out-of-range sentinel and reject NaN at the writer.
  • A valid value that collides with the sentinel. -9999 is safe for depth and unsafe for elevation in a region below sea level. Choose the sentinel against the quantity’s real range, not by convention.
  • Alpha bands standing in for a mask. Some pipelines carry validity as a fourth band. That works within the pipeline and is invisible to any consumer expecting nodata, so convert it to a real nodata declaration at the publication boundary.
  • Nodata surviving the raster and dying in the vectorisation. Polygonising a depth grid to produce flood extents will happily produce a polygon whose interior includes unconverged cells. Mask before polygonising, and carry the unconverged area as a separate “unassessed” polygon rather than folding it into either class.
  • Overviews built across the mask. The same averaging problem, one level up: an overview pyramid built without mask awareness contaminates every zoomed-out view, which is the view most people look at. Rebuild overviews after any mask repair.

Frequently Asked Questions

Why is losing the nodata mask more dangerous than losing other metadata? Because it changes an assertion rather than a description. A raster that loses its units is confusing and recoverable; a raster that loses its nodata mask converts every unmodelled cell into an ordinary value, and for a depth grid that value is usually zero, which means dry. The map then states that ground is passable in exactly the places where the model declined to produce a result — typically where the mesh failed to converge, which is often the hydraulically complex ground that most needed modelling. Nothing in the rendered output distinguishes it, because zero and dry share a colour.

What is the most common way the mask gets dropped in Python? Reading a band into a plain array with src.read(1) and writing a new file from it. The mask is a property of the dataset, not of the array, so the sentinel value arrives in the array as an ordinary number and is written back as one. The fix is src.read(1, masked=True), which returns a masked array carrying its own fill value, or copying the source profile wholesale and setting nodata explicitly on the write. Any function whose signature takes and returns a bare ndarray is a place where the mask cannot survive.

Can average resampling corrupt a hazard raster even when nodata is preserved? Yes, and this is the subtler failure. An average kernel over a window that straddles the mask boundary blends real depths with the sentinel, producing an output cell that is neither nodata nor a genuine value — a plausible number with no provenance that will pass any range check applied to it. Resample the values and the mask separately, using nearest for the mask, then reapply it. The same problem affects overview pyramids built without mask awareness, which contaminates the zoomed-out view most people actually look at.

Up: Raster Hazard Layers & Cloud-Optimized GeoTIFF

Other guides in Raster Hazard Layers & Cloud-Optimized GeoTIFF

Related
Serving Hazard Rasters Over Range Requests
→
Related
Sizing COG Overviews for Field Display Scales
→