EFHM return periods — how flood depth grows with rarity¶

A return period expresses how rare a flood is: a 1-in-10-year event is common and shallow, a 1-in-500-year event is rare and deep. This notebook fetches three return periods for the same area and compares their depths — one windowed read per period.

In [ ]:

import math
import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
from cleopatra.styling.colors import DATA_STYLES
from pyramids.dataset import Dataset

from earthlens.core import EarthLens

# `oslo` from cleopatra's Crameri palettes, reversed so the ramp runs light ->
# dark and deeper water reads as more ink. Picked by measurement, not taste:
# `cleopatra.styling.perceptual.perceptual_uniformity` scores it 0.09 against
# matplotlib's `Blues` at 0.27 (0 == perfectly even steps), and it is the only
# candidate that stays a single blue hue end to end.
DEPTH_CMAP = DATA_STYLES["oslo"]["oslo"]["cmap"].reversed()

out = Path(tempfile.mkdtemp(prefix="efhm-rp-"))
lat_lim, lon_lim = [51.7, 52.0], [4.6, 5.1]
periods = [10, 100, 500]
paths = EarthLens(
    data_source="jrc",
    lat_lim=lat_lim,
    lon_lim=lon_lim,
    return_periods=periods,
    path=out,
).download()
[p.name for p in paths]

import math import tempfile from pathlib import Path import matplotlib.pyplot as plt from cleopatra.styling.colors import DATA_STYLES from pyramids.dataset import Dataset from earthlens.core import EarthLens # `oslo` from cleopatra's Crameri palettes, reversed so the ramp runs light -> # dark and deeper water reads as more ink. Picked by measurement, not taste: # `cleopatra.styling.perceptual.perceptual_uniformity` scores it 0.09 against # matplotlib's `Blues` at 0.27 (0 == perfectly even steps), and it is the only # candidate that stays a single blue hue end to end. DEPTH_CMAP = DATA_STYLES["oslo"]["oslo"]["cmap"].reversed() out = Path(tempfile.mkdtemp(prefix="efhm-rp-")) lat_lim, lon_lim = [51.7, 52.0], [4.6, 5.1] periods = [10, 100, 500] paths = EarthLens( data_source="jrc", lat_lim=lat_lim, lon_lim=lon_lim, return_periods=periods, path=out, ).download() [p.name for p in paths]

Side-by-side depth maps¶

In [ ]:

grids = [Dataset.read_file(p) for p in paths]

# One shared upper bound so the three return periods are visually comparable.
vmax = max(grid.stats(approx_ok=False)["max"].iloc[0] for grid in grids)

# Two rows: two maps on top, the third centred beneath them. A 4-column grid
# makes that exact: each top panel spans two columns, and the bottom one spans
# the middle two, so it sits centred rather than left-aligned under the first.
# `constrained_layout` sizes the margins from the real text extents instead of
# matplotlib's fixed fractions; the panels are aspect-locked and bound by
# width, so that reclaimed margin is what actually enlarges them.
fig = plt.figure(figsize=(13, 9), constrained_layout=True)
gs = fig.add_gridspec(2, 4)
axes = [
    fig.add_subplot(gs[0, 0:2]),
    fig.add_subplot(gs[0, 2:4]),
    fig.add_subplot(gs[1, 1:3]),
]

# The bar below is shared, so hold its mappable deliberately rather than
# reading whatever `glyph` the loop leaves bound. Every panel shares the cmap
# and the clim set below, so the first one speaks for all of them.
shared_im = None
for ax, rp, grid in zip(axes, periods, grids):
    glyph = grid.plot(
        fig=fig,
        ax=ax,
        cmap=DEPTH_CMAP,
        title=f"RP{rp} (1-in-{rp}-year)",
        title_size=15,
        colorbar=False,
    )
    glyph.im.set_clim(0, vmax)
    shared_im = shared_im or glyph.im
    ax.xaxis.set_ticks_position("bottom")
    ax.set_xlabel("longitude", fontsize=12)
    ax.set_ylabel("latitude", fontsize=12)
    ax.tick_params(labelsize=11)

# One horizontal colour bar along the bottom of the whole figure rather than
# one hung off a single panel, which would shrink only that panel and leave the
# three unequal. All three share `vmax`, so a single bar describes every map;
# whole-metre ticks replace the default labels, which land on the data's own
# min/max (0.100 / 1.601).
bar = fig.colorbar(
    shared_im,
    ax=axes,
    orientation="horizontal",
    location="bottom",
    shrink=0.5,
    pad=0.02,
    aspect=40,
)
bar.set_label("depth (m)", size=12)
bar.set_ticks(range(0, math.floor(vmax) + 1))
bar.ax.tick_params(labelsize=11)

plt.show()
print(f"shared colour range: 0 to {vmax:.2f} m")

grids = [Dataset.read_file(p) for p in paths] # One shared upper bound so the three return periods are visually comparable. vmax = max(grid.stats(approx_ok=False)["max"].iloc[0] for grid in grids) # Two rows: two maps on top, the third centred beneath them. A 4-column grid # makes that exact: each top panel spans two columns, and the bottom one spans # the middle two, so it sits centred rather than left-aligned under the first. # `constrained_layout` sizes the margins from the real text extents instead of # matplotlib's fixed fractions; the panels are aspect-locked and bound by # width, so that reclaimed margin is what actually enlarges them. fig = plt.figure(figsize=(13, 9), constrained_layout=True) gs = fig.add_gridspec(2, 4) axes = [ fig.add_subplot(gs[0, 0:2]), fig.add_subplot(gs[0, 2:4]), fig.add_subplot(gs[1, 1:3]), ] # The bar below is shared, so hold its mappable deliberately rather than # reading whatever `glyph` the loop leaves bound. Every panel shares the cmap # and the clim set below, so the first one speaks for all of them. shared_im = None for ax, rp, grid in zip(axes, periods, grids): glyph = grid.plot( fig=fig, ax=ax, cmap=DEPTH_CMAP, title=f"RP{rp} (1-in-{rp}-year)", title_size=15, colorbar=False, ) glyph.im.set_clim(0, vmax) shared_im = shared_im or glyph.im ax.xaxis.set_ticks_position("bottom") ax.set_xlabel("longitude", fontsize=12) ax.set_ylabel("latitude", fontsize=12) ax.tick_params(labelsize=11) # One horizontal colour bar along the bottom of the whole figure rather than # one hung off a single panel, which would shrink only that panel and leave the # three unequal. All three share `vmax`, so a single bar describes every map; # whole-metre ticks replace the default labels, which land on the data's own # min/max (0.100 / 1.601). bar = fig.colorbar( shared_im, ax=axes, orientation="horizontal", location="bottom", shrink=0.5, pad=0.02, aspect=40, ) bar.set_label("depth (m)", size=12) bar.set_ticks(range(0, math.floor(vmax) + 1)) bar.ax.tick_params(labelsize=11) plt.show() print(f"shared colour range: 0 to {vmax:.2f} m")

Flooded area and depth grow with the return period¶

In [ ]:

# `read_array(masked=True)` returns a MaskedArray built from the band's own
# nodata, so the flooded-cell count needs no sentinel comparison.
for rp, grid in zip(periods, grids):
    stats = grid.stats(approx_ok=False).iloc[0]
    flooded = grid.read_array(masked=True).count()
    print(
        f"RP{rp:>3}: {flooded:>7,} flooded cells | "
        f"mean {stats['mean']:.2f} m | max {stats['max']:.2f} m"
    )

# `read_array(masked=True)` returns a MaskedArray built from the band's own # nodata, so the flooded-cell count needs no sentinel comparison. for rp, grid in zip(periods, grids): stats = grid.stats(approx_ok=False).iloc[0] flooded = grid.read_array(masked=True).count() print( f"RP{rp:>3}: {flooded:>7,} flooded cells | " f"mean {stats['mean']:.2f} m | max {stats['max']:.2f} m" )

Takeaway¶

Rarer floods inundate more cells and reach greater depths. Selecting the right return period lets you match a hazard layer to a design standard (e.g. RP100 for many flood defences, RP500 for critical infrastructure).

Three ways to name a return period¶

return_periods accepts an int, a plain string, or the RP-prefixed form that matches the source filenames. They are equivalent, so a value read from a config file or a filename does not need converting first.

In [ ]:

for spec in ([100], ["100"], ["RP100"]):
    picked = EarthLens(
        data_source="jrc",
        lat_lim=lat_lim,
        lon_lim=lon_lim,
        return_periods=spec,
        path=out,
    )
    print(f"{str(spec):10s} -> {picked.count()} product(s)")

for spec in ([100], ["100"], ["RP100"]): picked = EarthLens( data_source="jrc", lat_lim=lat_lim, lon_lim=lon_lim, return_periods=spec, path=out, ) print(f"{str(spec):10s} -> {picked.count()} product(s)")