🔧 Fix Natural hazards & risk Raster & rasterio Natural Hazards & Risk · step 20 Notebook By · Updated

A hazard raster's nodata is read as zero

On this page

Problem statement

The mean flood depth comes out at a quarter of a metre across a floodplain whose mapped water averages a metre deep. A DEM has a lake at −32,768 m. A warped hazard grid has lost a quarter of its flooded cells. All three are the same bug: a nodata value that something in the chain treated as a number.

It bites hardest on hazard data because hazard rasters are mostly nodata. The JRC European flood maps store only the wet cells — in a window south of Valencia, 76% of the file is nodata, and the sentinel is −3.4 × 10³⁸, the most negative float32 there is. One unguarded mean() on that array returns -inf.

Quick answer

Read masked, and decide explicitly what the mask means for this product:

import numpy as np, rasterio
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds

BBOX = [-0.46, 39.25, -0.33, 39.45]
with rasterio.open("flood/floodmap_EFAS_RP100_C.tif") as src:
    print(f"dtype {src.dtypes[0]}, nodata {src.nodata:.3g}")
    win = from_bounds(*transform_bounds(4326, src.crs, *BBOX), src.transform)
    raw = src.read(1, window=win)                      # nodata is just a number here
    masked = src.read(1, window=win, masked=True)      # nodata is masked

print(f"nodata cells: {masked.mask.sum():,} of {masked.size:,} ({masked.mask.mean():.1%})")
print(f"raw   min {raw.min():.3g}   mean {raw.mean():.3g}")
print(f"masked min {masked.min():.2f}  mean {masked.mean():.3f} m over the wet cells")
dtype float32, nodata -3.4e+38
nodata cells: 24,791 of 32,665 (75.9%)
raw   min -3.4e+38   mean -inf
masked min 0.10  mean 1.068 m over the wet cells

masked=True costs nothing and turns a silent wrong answer into a correct one. What it cannot tell you is whether those masked cells are dry land or unsurveyed ground — that is in the product documentation, not the file.

What nodata means in four hazard products and what may be substituted for it. the same mask, four meaningsproductnodata meansfill with zero?flood depth griddry, or outside the modelled reachyesDEMunsurveyed, or void-filled elsewhereno — invents a valleyshaking or severity surfaceoutside the grid, unknownno — invents calmland coveroutside the mapno — 0 is a classthe file cannot tell you which of these it is — the documentation canRead masked first, then decide what the mask means spatialworkflow.io
The same mask, four meanings: only the documentation distinguishes them.

Read the nodata, mask on read, decide what gaps mean

1. Look at src.nodata before reading anything

If it is None, the file declares nothing and you must find the convention yourself — check the minimum and the histogram for a suspicious spike.

2. Read masked, or convert to NaN immediately

read(1, masked=True) for float data, then .filled(np.nan) if you prefer arrays. For integer data there is no NaN, so keep the mask.

3. Decide what the gap means for this product

Depth grids usually mean "dry", so zero is right. Elevation, intensity and susceptibility usually mean "unknown", so zero is a lie. Land cover means "outside the map".

4. Never let the sentinel through a resampling

Bilinear or cubic resampling averages neighbours. One cell of −3.4 × 10³⁸ destroys every cell it touches, which is how a warp loses flooded area around every edge.

5. Use NaN-aware statistics

np.nanmean, np.nanpercentile, and a masked array where integers are involved. A plain mean() over a zero-filled hazard grid answers a different question — usually "per unit of study area" rather than "where it is hazardous".

6. Set nodata on everything you write

A derived raster with no declared nodata pushes this whole problem onto the next person, and they will not have your notebook.

A row of depth cells beside nodata cells and what bilinear resampling does to them. one sentinel in the kernel is enoughsource1.20.80.4−3.4e38−3.4e38−3.4e38bilinear1.00.6−8.5e37−1.7e38−3.4e38−3.4e38the two cells next to the gap were real water and are now nonsensedeclare src_nodata, or resample with nearest spatialworkflow.io
One sentinel in the kernel is enough: bilinear resampling spreads it to every neighbour.

What a warp loses and a mean misstates

The data used on this page

The European Flood Awareness System's pan-European flood hazard maps from the JRC (data catalogue): water depth in metres on a 100 m grid for the 1-in-100-year flood, one zipped GeoTIFF per return period, saved under flood/.

import os
import zipfile

import requests

EFAS = "https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/FLOODS/EuropeanMaps/floodMap_RP{}.zip"
os.makedirs("flood", exist_ok=True)
for rp in ("100",):
    tif = f"floodmap_EFAS_RP{rp}_C.tif"
    if not os.path.exists(f"flood/{tif}"):                      # about 180 MB per return period
        with requests.get(EFAS.format(rp), stream=True, timeout=600) as r:
            r.raise_for_status()
            with open("efas.zip", "wb") as f:
                for chunk in r.iter_content(1 << 20):
                    f.write(chunk)
        zipfile.ZipFile("efas.zip").extract(tif, "flood")
        os.remove("efas.zip")
print(sorted(os.listdir("flood")))
['floodmap_EFAS_RP100_C.tif']

Example 1 — what the warp loses

from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT

with rasterio.open("flood/floodmap_EFAS_RP100_C.tif") as src:
    win = from_bounds(*transform_bounds(4326, src.crs, *BBOX), src.transform)
    grid = dict(crs=src.crs, transform=src.window_transform(win),
                width=int(win.width), height=int(win.height), resampling=Resampling.bilinear)
    with WarpedVRT(src, **grid) as vrt:                       # source nodata respected
        good = vrt.read(1, masked=True).filled(np.nan)
    with WarpedVRT(src, src_nodata=None, **grid) as vrt:      # sentinel treated as data
        bad = vrt.read(1)

print(f"nodata respected: {int(np.nansum(good > 0)):,} wet cells, "
      f"mean depth {np.nanmean(good[good > 0]):.3f} m, max {np.nanmax(good):.2f} m")
print(f"nodata as data:   {int((bad > 0).sum()):,} wet cells, "
      f"min {bad.min():.3g}, max {bad.max():.2f} m")
print(f"cells dragged to nonsense by the kernel: {int((bad < -1000).sum()):,}")
nodata respected: 7,736 wet cells, mean depth 1.059 m, max 12.60 m
nodata as data:   5,535 wet cells, min -3.4e+38, max 12.60 m
cells dragged to nonsense by the kernel: 26,757

Nearly 30% of the flooded cells vanish, all of them around the edges of the flood extent — which is exactly where exposure sits.

Example 2 — the mean that answers a different question

depth = good
zero_filled = np.where(np.isfinite(depth), depth, 0.0)
wet = depth > 0

print(f"mean over the whole window, nodata as 0 : {zero_filled.mean():.4f} m")
print(f"mean over wet cells only                : {np.nanmean(depth[wet]):.3f} m")
print(f"wet area {wet.sum() * 0.01:,.1f} km² of a {depth.size * 0.01:,.1f} km² window")
print(f"95th percentile of depth where wet      : {np.nanpercentile(depth[wet], 95):.2f} m")
mean over the whole window, nodata as 0 : 0.2537 m
mean over wet cells only                : 1.059 m
wet area 77.4 km² of a 322.9 km² window
95th percentile of depth where wet      : 2.32 m

Both means are arithmetically correct and they answer different questions. The first is a depth averaged over dry land, which is not a quantity anybody wants; quote it as "mean flood depth" and every comparison with another study is broken.

Example 3 — write it back so the next reader is safe

profile = dict(driver="GTiff", dtype="float32", count=1, crs=grid["crs"],
               transform=grid["transform"], width=grid["width"], height=grid["height"],
               nodata=np.nan, compress="deflate", tiled=True)
with rasterio.open("depth_rp100_window.tif", "w", **profile) as dst:
    dst.write(np.where(wet, depth, np.nan).astype("float32"), 1)
    dst.update_tags(nodata_means="dry or outside the modelled floodplain",
                    source="JRC EFAS flood hazard map, RP100")

with rasterio.open("depth_rp100_window.tif") as src:
    back = src.read(1, masked=True)
    print(f"declared nodata {src.nodata}, masked cells {back.mask.sum():,}")
    print(f"round-trip mean over wet cells {back.mean():.3f} m")
    print(f"tag: {src.tags()['nodata_means']}")
declared nodata nan, masked cells 24,556
round-trip mean over wet cells 1.059 m
tag: dry or outside the modelled floodplain

Why hazard rasters are mostly nodata

Why hazard rasters are mostly nodata

A flood map covers a continent and floods a river corridor. A burn severity map covers a scene and burns one perimeter. Storing only the affected cells keeps these files usable, so the fraction of nodata is high by design — and the fraction of your array that a mistake corrupts is correspondingly high.

Why the sentinel is so destructive in float32

-3.4e38 is the most negative representable float32. Averaged with anything it stays enormous; summed into an accumulator it overflows to -inf; compared with > 0 it is simply False, so a wet/dry test still works and the bug survives until someone takes a mean. That is why the failure often appears far downstream of the read.

Why zero is sometimes the right substitution

For a depth grid, "not mapped" and "dry" are the same thing for most purposes, and filling with zero is correct and convenient. For a DEM, an intensity surface or a susceptibility index, the missing cells are unknown and filling with zero invents a valley, a quiet spot or a safe area. The rule is not "never fill with zero" but "know which of the two products you have, from the documentation".

Wet cell counts and mean depths under correct and incorrect nodata handling. one file, two readingswet cells, masked7,736 cellswet cells, unmasked5,535 — 28% lost at the edgesmean depth, wet1.059 mmean depth, window0.254 m — averaged over dry land040008000depth bars ×2000 to share the scaleJRC EFAS RP100 depth map, window south of Valencia spatialworkflow.io
Two reads of one file: the flooded area and the mean depth both move.

Integer rasters, missing tags and alpha bands

  • Integer rasters have no NaN. Keep the masked array or carry the mask separately; astype(float) first if you need NaN.
  • nodata=None in the profile means "not declared", not "no nodata cells". Look at the data.
  • Multiple sentinels exist in some products (−9999 and 0 both meaning missing); the histogram shows them.
  • Alpha and mask bands: a GeoTIFF can carry a per-pixel mask instead of a sentinel — src.read_masks(1) reads it.
  • Compression does not care, but a nodata value that is not exactly representable in the dtype can round and stop matching.
  • Zonal statistics libraries usually respect nodata; check, because the default varies.

FAQ

Why is my mean flood depth so low?

Because nodata cells were counted as zero depth. Average over the wet cells, or say the average is over the whole area.

What is -3.4e38 in my raster?

The float32 minimum, used as a nodata sentinel. It is not a measurement.

Can I just fill nodata with zero?

For depth or damage grids, usually yes. For elevation, intensity or index surfaces, no — zero is a value there.

Why did reprojection lose flooded cells?

Bilinear resampling averaged the sentinel into real cells. Declare the source nodata, or resample with nearest.

How do I read a raster safely?

src.read(1, masked=True), then decide what the mask means before filling it.

Should I set nodata when writing?

Always. An undeclared nodata value is the same bug, passed on.