An inundation map floods inland areas that cannot flood

On this page

Problem statement

The sea-level map shows water in a quarry twenty kilometres inland, in a drained polder behind a dike, and in every road cutting and gravel pit in the county. All of them are below the level. None of them can flood from the sea.

The cause is that a threshold has no concept of connectivity: dem <= level is true for every cell below the level regardless of whether water can reach it. On a coarse test over the Netherlands, thresholding at 0 m marked 20,944 cells as flooded, of which 580 cells in 85 separate pockets had no path to the edge of the grid.

A related and less fixable problem sits behind it: even the connected cells include everything behind a defence the terrain model does not resolve.

Quick answer

Label the below-threshold mask and keep only the components connected to the sea:

import numpy as np
from scipy import ndimage

below = (dem <= level).filled(False)
labels, n = ndimage.label(below)                      # 4-connectivity by default

edge = set(labels[0, :]) | set(labels[-1, :]) | set(labels[:, 0]) | set(labels[:, -1])
edge.discard(0)
connected = np.isin(labels, list(edge))

print(f"bathtub {below.sum():,}, connected {connected.sum():,}, "
      f"isolated {below.sum() - connected.sum():,} in {n - len(edge):,} pockets")

bathtub 20,944, connected 20,364, isolated 580 in 85 pockets

Seeding from the array edge assumes the sea reaches it. Where it does not — an inland tile, a study area cut short of the coast — seed from a polygon of known open water instead, or every component is isolated and the map is empty.

Five causes; only the first is fixed by a connectivity test.

Confirm it is connectivity, then filter from a sea seed

1. Confirm it is a connectivity problem

Label the mask and look at where the isolated components are. Quarries, pits, polders and cuttings with no path to the sea are the signature. If the flooded area is a continuous sheet reaching inland along a valley, the problem is something else.

2. Apply the connected-component filter

One ndimage.label call and a set of seed labels. It is the cheapest correction in coastal mapping and it should be in every pipeline.

3. Seed from the sea, not from the array edge

The edge shortcut works only when the sea touches the raster boundary. A polygon of open water — a coastline layer, an ocean mask, a manually drawn seed — is correct everywhere.

4. Use 4-connectivity for the water

Under 8-connectivity, water passes diagonally between cells that touch only at a corner, so a one-cell barrier with a diagonal step in it does not hold. Water does not flow through a corner.

5. Burn in the defences

This is the part a connectivity test cannot fix. A dike or seawall two to three metres wide is under-represented in a DEM produced at one to two metres with smoothing, so the terrain shows a gap where a wall stands. Rasterise the defence lines at their crest height, buffered by at least one cell, and take the maximum with the DEM.

6. Check the tidal limit up rivers

The sea reaches inland along estuaries, and a threshold carries it all the way to the head of the catchment. Cut the water surface at the tidal limit, or use a sloping water surface rather than a flat one.

7. Report the isolated pockets separately

They are genuinely below the level. They may flood from rainfall or groundwater, and dropping them without saying so hides a real result.

The edge seed is a convenience that fails silently on any tile the sea does not reach.

A connectivity filter with a proper seed

The data used on this page

The terrain is ETOPO 2022 at 60 arc-seconds (about 1.85 km), cut to the Netherlands and its coast, 3.0–7.5° E, 50.7–53.7° N, as a masked array dem with its transform. The sea seed sea_poly is Natural Earth's 1:10m ocean polygon for the same box. defences are the OpenStreetMap lines tagged man_made=dyke in the box; OpenStreetMap rarely records their height, so each is given a nominal 5 m crest — a stated assumption to show the mechanism, not the Dutch defence standard. level is 0 m.

import os
import subprocess

import geopandas as gpd
import numpy as np
import requests
import xarray as xr
from rasterio.transform import from_origin
from shapely.geometry import LineString, box

URL = ("https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/"
       "60s_surface_elev_netcdf/ETOPO_2022_v1_60s_N90W180_surface.nc")
if not os.path.exists("etopo_2022_60s.nc"):
    subprocess.run(["curl", "-s", "-C", "-", "--retry", "3", "-o", "etopo_2022_60s.nc", URL], check=True)
W, S, E, N = 3.0, 50.7, 7.5, 53.7
z = xr.open_dataset("etopo_2022_60s.nc").z.sel(lon=slice(W, E), lat=slice(S, N))
z = z.isel(lat=slice(None, None, -1))                    # north-up for rasterio
dem = np.ma.masked_invalid(z.values.astype("float32"))
res = float(abs(z.lon[1] - z.lon[0]))
transform = from_origin(float(z.lon[0]) - res / 2, float(z.lat[0]) + res / 2, res, res)

if not os.path.exists("ne_10m_ocean.zip"):
    open("ne_10m_ocean.zip", "wb").write(
        requests.get("https://naciscdn.org/naturalearth/10m/physical/ne_10m_ocean.zip", timeout=300).content)
sea_poly = gpd.read_file("ne_10m_ocean.zip").clip(box(W, S, E, N)).union_all()

if not os.path.exists("dykes.gpkg"):                    # one Overpass query, geometry inline
    q = f'[out:json][timeout:800];way["man_made"="dyke"]({S},{W},{N},{E});out geom;'
    ways = requests.post("https://maps.mail.ru/osm/tools/overpass/api/interpreter",
                         data={"data": q}, timeout=900).json()["elements"]
    gpd.GeoDataFrame(geometry=[LineString([(p["lon"], p["lat"]) for p in w["geometry"]])
                               for w in ways if len(w.get("geometry", [])) > 1],
                     crs=4326).to_file("dykes.gpkg")
defences = gpd.read_file("dykes.gpkg").assign(crest_m=5.0)
level = 0
print(f"dem {dem.shape[1]} x {dem.shape[0]} cells of {res * 60:.0f}'; "
      f"{len(defences):,} dyke lines")

dem 270 x 180 cells of 1'; 1,250 dyke lines

Example 1 — a connectivity filter with a proper seed

import numpy as np, rasterio, geopandas as gpd
from rasterio.features import rasterize
from scipy import ndimage

def connected_inundation(dem, transform, shape, level, sea=None, connectivity=4):
    below = np.asarray(dem <= level)
    structure = None if connectivity == 4 else np.ones((3, 3), dtype=int)
    labels, n = ndimage.label(below, structure=structure)

    if sea is not None:
        seed_mask = rasterize([sea], out_shape=shape, transform=transform,
                              fill=0, default_value=1, dtype="uint8").astype(bool)
        seeds = set(np.unique(labels[seed_mask & below]))
    else:
        seeds = (set(labels[0, :]) | set(labels[-1, :])
                 | set(labels[:, 0]) | set(labels[:, -1]))
    seeds.discard(0)

    if not seeds:
        raise ValueError("no flooded component touches the sea — check the seed")

    connected = np.isin(labels, list(seeds))
    return connected, {
        "components": int(n), "seeded_components": len(seeds),
        "bathtub_cells": int(below.sum()), "connected_cells": int(connected.sum()),
        "isolated_cells": int(below.sum() - connected.sum()),
        "isolated_pockets": int(n - len(seeds)),
    }


connected, stats = connected_inundation(dem, transform, dem.shape, level, sea=sea_poly)
print(stats)

{'components': 89, 'seeded_components': 18, 'bathtub_cells': 20944, 'connected_cells': 20402, 'isolated_cells': 542, 'isolated_pockets': 71}

Raising when no component touches the seed is the important line: an empty map is a failure that otherwise looks like "nothing floods".

Example 2 — how much connectivity matters at different levels

import numpy as np
from scipy import ndimage

for level in (0, 1, 2, 5):
    below = dem <= level
    labels, n = ndimage.label(below)
    edge = (set(labels[0, :]) | set(labels[-1, :])
            | set(labels[:, 0]) | set(labels[:, -1]))
    edge.discard(0)
    connected = np.isin(labels, list(edge))
    print(f"level {level:>2} m: bathtub {below.sum():6,} ({below.mean():5.1%}), "
          f"connected {connected.sum():6,} ({connected.mean():5.1%}), "
          f"isolated {below.sum()-connected.sum():5,} in {n-len(edge):,} pockets")

level  0 m: bathtub 20,944 (43.1%), connected 20,364 (41.9%), isolated   580 in 85 pockets
level  1 m: bathtub 23,688 (48.7%), connected 23,521 (48.4%), isolated   167 in 48 pockets
level  2 m: bathtub 25,296 (52.0%), connected 25,186 (51.8%), isolated   110 in 30 pockets
level  5 m: bathtub 27,559 (56.7%), connected 27,495 (56.6%), isolated    64 in 23 pockets

The isolated count falls as the level rises, because higher water connects pockets to each other and to the sea. The grid here is one arc-minute — about 1.85 km — which is far too coarse for a real inundation study; it shows the mechanism, not the answer.

Example 3 — burn in the defences and see what changes

import numpy as np, geopandas as gpd
from rasterio.features import rasterize

def burn_defences(dem, transform, shape, defences, crest_col="crest_m", buffer_cells=1):
    cell = abs(transform.a)
    burn = rasterize(
        ((g.buffer(buffer_cells * cell), h)
         for g, h in zip(defences.geometry, defences[crest_col])),
        out_shape=shape, transform=transform, fill=np.nan, dtype="float32")
    return np.fmax(np.asarray(dem), burn)

level = 2                                    # below the 5 m nominal crest
before, _ = connected_inundation(dem, transform, dem.shape, level, sea=sea_poly)
after, _ = connected_inundation(
    burn_defences(dem, transform, dem.shape, defences), transform, dem.shape,
    level, sea=sea_poly)

print(f"flooded cells before defences {before.sum():,}, after {after.sum():,} "
      f"({after.sum()/before.sum()-1:+.1%})")

flooded cells before defences 25,185, after 23,847 (-5.3%)

Here the difference was 5.3%: 1,338 cells that the terrain alone floods at 2 m are held back once 1,250 OpenStreetMap dyke lines stand at a nominal 5 m. At 1.85 km a cell rarely shows the gap a dike fills, so the correction is modest; on a one- or two-metre DEM, where the dike is narrower than a cell but the gap is not, it is the single biggest correction available after connectivity.

A threshold has no memory of the sea

Why a threshold has no memory of the sea

dem <= level is evaluated per cell with no reference to any other cell. Every quarry, pit, polder and cutting below the level satisfies it. Connectivity is the missing constraint, and expressing it as a connected-component label is both the standard solution and a single function call.

Why the seed choice is the part that goes wrong

Seeding from the array edge is convenient and wrong for any tile the sea does not reach — an inland study area, a clipped extract, a tile in the middle of a coastline mosaic. The symptom is an empty map, which is easy to misread as "nothing floods at this level". Seeding from an explicit water polygon makes the assumption visible.

Why 4-connectivity is the conservative choice

With 8-connectivity, two cells that share only a corner are connected, so water passes through a diagonal gap in a one-cell wall. Physically it does not. The standard asymmetry in flood modelling is 4-connectivity for the water and 8-connectivity for the barriers, so a diagonal step in a dike still holds.

Why defences are the harder half

A connectivity test uses only the DEM, and the DEM does not contain the wall. Adding the defences needs a separate layer with crest heights, which many places do not publish, and even with it a still-water map cannot represent overtopping, breaching or drainage. That is the honest limit: a connected fill answers "what is below the level and reachable from the sea", not "what floods".

Counter-intuitively, the isolated share shrinks as the water rises.

Culverts, nodata barriers and tidal limits

Culverts and sluices connect where the terrain does not. Burn them as low cells.
Nodata acts as a barrier. Mask it, or the flood stops at the tile edge.
A DSM floods around buildings. Use a DTM.
The tidal limit up a river is not the same as the level everywhere.
Isolated pockets may still flood from rainfall or groundwater.
A diagonal gap in a dike holds under 4-connectivity and not under 8.
Report the level, the datum and the method with the extent.
Present a band, not a line. DEM error is a large fraction of a half-metre scenario.

Internal links

How to map coastal inundation from a DEM in Python — the full workflow
Sea level rise data explained — assembling the level
Tidal datums explained: which shoreline is the shoreline — the datum the level is in
How to fill sinks and compute flow direction in Python — depressions from the hydrological side
Height above nearest drainage explained — a better proxy for fluvial flooding
How to extract a coastline from a raster in Python — the same connectivity rule
How to rasterize a vector layer in Python — burning in defences
Digital elevation models explained — DTM versus DSM

FAQ

Why does my sea-level map flood a quarry inland?

Because a threshold has no connectivity. Label the below-threshold mask and keep only components that reach the sea.

How do I seed the connectivity test?

From a polygon of known open water. Seeding from the array edge works only if the sea touches the raster boundary.

Why is my connected map empty?

Because no below-threshold component touched the seed — usually an inland tile with an edge-based seed. Raise rather than returning an empty map.

Should I use 4- or 8-connectivity?

Four for the water. Under 8-connectivity, water passes diagonally through a one-cell gap in a wall.

Why does water still appear behind the sea wall?

Because the wall is narrower than a DEM cell and was smoothed away. Rasterise the defences at their crest height, buffered by at least one cell.

What about the isolated pockets?

Report them. They are genuinely below the level and may flood from rainfall or groundwater, just not from the sea.