Showcase - how much worse is a rarer flood, and how your window changes the answer¶

The JRC European Flood Hazard Map gives river-flood water depth for a chosen return period: how deep the water gets in a 1-in-10, 1-in-100 or 1-in-500 year event. The obvious expectation is that a rarer flood covers more ground.

Over the Dutch Rhine that holds - the flooded area grows by about a third from RP10 to RP500, and the water standing on it gets deeper as well.

The second half is a caution. Crop to a small box over the Merwede / Biesbosch and the flooded area barely grows at all - about seven per cent. That is not a fact about dikes. It is a fact about the box: it is already 89% under water at RP10, so its area cannot grow much whatever the water does. Widen the box and the growth returns to the whole-river figure.

This notebook reads both scales from the same three source rasters, so the comparison is like for like, and measures the saturation that explains the gap.

Setup¶

Each return period is a single whole-Europe GeoTIFF of roughly 23 GB. Nothing is downloaded whole: the backend opens each file lazily over /vsicurl and reads only the pixel window the area of interest needs.

In [ ]:

# NBVAL_SKIP
from IPython.display import HTML, display

# Jupyter scales an output image down to the cell width, so a larger `figsize`
# adds detail but not displayed size. Widening the container is what makes the
# maps physically bigger on screen.
display(
    HTML(
        "<style>"
        ".jp-Cell, .jp-Notebook, .cell, #notebook-container { max-width: 100% !important; }"
        ".jp-OutputArea-output img, .output_area img { max-width: 100% !important;"
        " width: 100% !important; height: auto !important; }"
        "</style>"
    )
)

# NBVAL_SKIP from IPython.display import HTML, display # Jupyter scales an output image down to the cell width, so a larger `figsize` # adds detail but not displayed size. Widening the container is what makes the # maps physically bigger on screen. display( HTML( "<style>" ".jp-Cell, .jp-Notebook, .cell, #notebook-container { max-width: 100% !important; }" ".jp-OutputArea-output img, .output_area img { max-width: 100% !important;" " width: 100% !important; height: auto !important; }" "</style>" ) )

In [ ]:

# NBVAL_SKIP
import math
import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
from cleopatra.styling.colors import DATA_STYLES

# This figure is a composite -- several glyphs' panels assembled in one
# plt.figure, so no single glyph owns it. cleopatra 0.40 added the public
# stamp_mark_on() as the supported route for exactly that case
# (serapeum-org/cleopatra#371); the local alias keeps the call sites unchanged.
from cleopatra.styling.watermark import stamp_mark_on as stamp_mark
from matplotlib.patches import Rectangle
from pyramids.dataset import Dataset
from pyramids.plot import DataStyle

from earthlens.core import EarthLens

# `oslo` (Crameri) reversed: a single-hue blue ramp running light -> dark, so
# deeper water reads as more ink. Chosen by measurement -- cleopatra scores its
# perceptual uniformity at 0.09 against matplotlib's `Blues` at 0.27.
DEPTH_CMAP = DATA_STYLES["oslo"]["oslo"]["cmap"].reversed()

# The stacked navy lockup from the repo brand kit. It is a solid card rather
# than a transparent mark, so it reads as a small brand panel in the corner
# instead of a wash over the maps -- which is why it goes in the figure margin
# and not on a panel. Paths in a notebook are relative to the notebook itself.
LOGO = Path("../../_images/branding/earthlens-brand-kit/logo")
LOGO = LOGO / "earthlens-lockup-stacked-navy.png"
_mark = plt.imread(LOGO)
LOGO_ASPECT = _mark.shape[0] / _mark.shape[1]

# The two extents are looked up by name in later cells, so the names live
# here rather than being retyped as dict keys each time.
RHINE = "Dutch Rhine"
DELTA = "Merwede / Biesbosch"

PERIODS = [10, 100, 500]

# One area of interest: Lobith at the German border across to the North Sea
# mouths. The zoom below is a sub-box of it, cropped from the same rasters
# rather than fetched again -- so the two scales cannot disagree.
LAT_LIM, LON_LIM = [51.6, 52.1], [4.0, 6.2]
# The Waal continuing as the Merwede, through the Biesbosch.
ZOOM = [4.6, 51.7, 5.1, 52.0]  # west, south, east, north

out = Path(tempfile.mkdtemp(prefix="efhm-rhine-"))

# NBVAL_SKIP import math import tempfile from pathlib import Path import matplotlib.pyplot as plt from cleopatra.styling.colors import DATA_STYLES # This figure is a composite -- several glyphs' panels assembled in one # plt.figure, so no single glyph owns it. cleopatra 0.40 added the public # stamp_mark_on() as the supported route for exactly that case # (serapeum-org/cleopatra#371); the local alias keeps the call sites unchanged. from cleopatra.styling.watermark import stamp_mark_on as stamp_mark from matplotlib.patches import Rectangle from pyramids.dataset import Dataset from pyramids.plot import DataStyle from earthlens.core import EarthLens # `oslo` (Crameri) reversed: a single-hue blue ramp running light -> dark, so # deeper water reads as more ink. Chosen by measurement -- cleopatra scores its # perceptual uniformity at 0.09 against matplotlib's `Blues` at 0.27. DEPTH_CMAP = DATA_STYLES["oslo"]["oslo"]["cmap"].reversed() # The stacked navy lockup from the repo brand kit. It is a solid card rather # than a transparent mark, so it reads as a small brand panel in the corner # instead of a wash over the maps -- which is why it goes in the figure margin # and not on a panel. Paths in a notebook are relative to the notebook itself. LOGO = Path("../../_images/branding/earthlens-brand-kit/logo") LOGO = LOGO / "earthlens-lockup-stacked-navy.png" _mark = plt.imread(LOGO) LOGO_ASPECT = _mark.shape[0] / _mark.shape[1] # The two extents are looked up by name in later cells, so the names live # here rather than being retyped as dict keys each time. RHINE = "Dutch Rhine" DELTA = "Merwede / Biesbosch" PERIODS = [10, 100, 500] # One area of interest: Lobith at the German border across to the North Sea # mouths. The zoom below is a sub-box of it, cropped from the same rasters # rather than fetched again -- so the two scales cannot disagree. LAT_LIM, LON_LIM = [51.6, 52.1], [4.0, 6.2] # The Waal continuing as the Merwede, through the Biesbosch. ZOOM = [4.6, 51.7, 5.1, 52.0] # west, south, east, north out = Path(tempfile.mkdtemp(prefix="efhm-rhine-"))

1 - Fetch once¶

Three rasters, one per return period, over the whole river. The zoom used later is a sub-box of this extent, so it is cropped from these same rasters rather than fetched a second time -- half the requests, and the two scales are guaranteed to be the same pixels.

In [ ]:

# NBVAL_SKIP
paths = EarthLens(
    data_source="jrc",
    return_periods=PERIODS,
    lat_lim=LAT_LIM,
    lon_lim=LON_LIM,
    path=out,
).download()

whole = [Dataset.read_file(p) for p in paths]
size = sum(p.stat().st_size for p in paths) / 1e6
print(
    f"{whole[0].columns} x {whole[0].rows} px   {size:.1f} MB for {len(paths)} rasters"
)

# The zoom is cut from the rasters already in hand. `touch=False` keeps only
# pixels fully inside the box; the default `touch=True` would also take the ones
# the boundary grazes, which is ~0.9% more cells and would shift the numbers.
zoomed = [ds.crop(bbox=ZOOM, touch=False) for ds in whole]
print(
    f"zoom: {zoomed[0].columns} x {zoomed[0].rows} px, cropped -- not downloaded again"
)

# NBVAL_SKIP paths = EarthLens( data_source="jrc", return_periods=PERIODS, lat_lim=LAT_LIM, lon_lim=LON_LIM, path=out, ).download() whole = [Dataset.read_file(p) for p in paths] size = sum(p.stat().st_size for p in paths) / 1e6 print( f"{whole[0].columns} x {whole[0].rows} px {size:.1f} MB for {len(paths)} rasters" ) # The zoom is cut from the rasters already in hand. `touch=False` keeps only # pixels fully inside the box; the default `touch=True` would also take the ones # the boundary grazes, which is ~0.9% more cells and would shift the numbers. zoomed = [ds.crop(bbox=ZOOM, touch=False) for ds in whole] print( f"zoom: {zoomed[0].columns} x {zoomed[0].rows} px, cropped -- not downloaded again" )

2 - What changes between return periods¶

Two numbers per raster: how many cells are flooded at all, and how deep the water is on average. stats() excludes the nodata, and read_array(masked=True) gives a MaskedArray built from the band's own nodata, so counting valid cells needs no sentinel comparison.

In [ ]:

# NBVAL_SKIP
scales = {RHINE: whole, DELTA: zoomed}

measured = {}
for name, grids in scales.items():
    rows = []
    for rp, grid in zip(PERIODS, grids):
        stats = grid.stats(approx_ok=False).iloc[0]
        rows.append(
            {
                "rp": rp,
                "cells": int(grid.read_array(masked=True).count()),
                "mean": float(stats["mean"]),
                "max": float(stats["max"]),
            }
        )
    measured[name] = rows

for name, rows in measured.items():
    first = rows[0]
    print(name)
    for row in rows:
        d_area = 100 * (row["cells"] / first["cells"] - 1)
        d_deep = 100 * (row["mean"] / first["mean"] - 1)
        print(
            f"  RP{row['rp']:<4} {row['cells']:>10,} cells  mean {row['mean']:.2f} m"
            f"  max {row['max']:.2f} m   area {d_area:+5.1f}%  depth {d_deep:+5.1f}%"
        )
    print()

# NBVAL_SKIP scales = {RHINE: whole, DELTA: zoomed} measured = {} for name, grids in scales.items(): rows = [] for rp, grid in zip(PERIODS, grids): stats = grid.stats(approx_ok=False).iloc[0] rows.append( { "rp": rp, "cells": int(grid.read_array(masked=True).count()), "mean": float(stats["mean"]), "max": float(stats["max"]), } ) measured[name] = rows for name, rows in measured.items(): first = rows[0] print(name) for row in rows: d_area = 100 * (row["cells"] / first["cells"] - 1) d_deep = 100 * (row["mean"] / first["mean"] - 1) print( f" RP{row['rp']:<4} {row['cells']:>10,} cells mean {row['mean']:.2f} m" f" max {row['max']:.2f} m area {d_area:+5.1f}% depth {d_deep:+5.1f}%" ) print()

The two scales disagree, and the disagreement is the lesson. Across the whole Dutch Rhine the flooded area grows by about a third from RP10 to RP500. Inside the small Merwede / Biesbosch box it grows by about seven per cent, while the mean depth still rises by roughly half.

It is tempting to read that second number as embankments holding the water in. The measurement does not support it. The box is already about 89% under water at RP10, so there is almost no dry ground left for a rarer flood to take - the ceiling is the crop, not the dikes.

3 - Upstream, Lobith to Gorinchem¶

The eastern half of the Dutch Rhine, from the German border down to Gorinchem. Here the river still has floodplain to fill, so a rarer return period reads as pale blue reaching further from the channel.

In [ ]:

# NBVAL_SKIP
# All three figures below share this colour range, so a shade means the same
# depth in every one. That is the whole point of the comparison, so it is
# computed once here rather than per figure.
vmax = max(row["max"] for rows in measured.values() for row in rows)

BASEMAP = "Esri.WorldImagery"
OVERLAY = DataStyle(alpha=0.85)

# Titles live here because each figure is rendered twice -- once at print scale
# and once, further down, with the fonts enlarged for posting.
UP_TITLE = (
    "Upstream Dutch Rhine \u00b7 Lobith to Gorinchem\n"
    "JRC European Flood Hazard Map \u00b7 Lobith to Gorinchem"
)
DOWN_TITLE = (
    "Downstream Dutch Rhine \u00b7 Gorinchem to the North Sea\n"
    "JRC European Flood Hazard Map \u00b7 Gorinchem to the North Sea"
)
ZOOM_TITLE = (
    "A closer window \u00b7 its area barely grows because the box is already full\n"
    "the red box above, at the same colour scale"
)


def growth(rows, index):
    """Return the area / depth change of one return period against RP10."""
    first, row = rows[0], rows[index]
    return (
        100 * (row["cells"] / first["cells"] - 1),
        100 * (row["mean"] / first["mean"] - 1),
    )


def add_depth_bar(fig, mappable, font_scale=1.0):
    """Put the shared depth scale along the bottom of a figure.

    Built with `fig.colorbar` rather than through the glyph's own `.cbar`: each
    panel's glyph owns one bar for one axes, and what these figures need is a
    single bar spanning every panel so one shade means one depth across all of
    them. The mappable handed over is still the glyph's `.im`, so the bar tracks
    the range set by `set_clim`.
    """
    bar = fig.colorbar(
        mappable,
        ax=fig.axes,
        orientation="horizontal",
        location="bottom",
        shrink=0.45,
        pad=0.01,
        aspect=50,
    )
    bar.set_label("river-flood water depth (m)", size=26 * font_scale)
    bar.set_ticks(range(0, math.floor(vmax) + 1))
    bar.ax.tick_params(labelsize=22 * font_scale)
    return bar


def stamp_logo(fig, ax, height_frac=0.22, pad_in=0.35):
    """Stamp the earthlens mark inside a map panel's top-right corner."""
    # A mark in the figure margin is an orphan; a watermark belongs on the map.
    # `stamp_mark` anchors to a FIGURE corner, so the offset that lands it inside
    # `ax` is measured from the figure edge across to the panel edge. Panel
    # positions are only final after the layout pass, hence `execute` first.
    fig.get_layout_engine().execute(fig)
    fig_w, fig_h = fig.get_size_inches()
    box = ax.get_position()
    # `frac` is read as the mark's longer on-figure side, judged on FRACTIONS
    # rather than inches -- so derive both candidates and hand over the larger.
    height = box.height * height_frac
    width = height * fig_h / (fig_w * LOGO_ASPECT)
    return stamp_mark(
        fig,
        LOGO,
        frac=max(width, height),
        corner="upper right",
        margin=(1 - box.x1 + pad_in / fig_w, 1 - box.y1 + pad_in / fig_h),
        shadow=False,
    )


# Map panels are aspect-locked, so at full width their height is fixed by their
# own shape -- the whole river is 4.4:1, and no `figsize` can make that strip
# taller without stretching the geography. Cutting the river in half at Gorinchem
# halves the aspect to 2.2:1, which is what doubles the height each map gets.
# The two halves then get a figure each, so neither is a mile-long scroll.
MID = 5.1
UPSTREAM = [MID, LAT_LIM[0], LON_LIM[1], LAT_LIM[1]]
DOWNSTREAM = [LON_LIM[0], LAT_LIM[0], MID, LAT_LIM[1]]


def draw_rhine_half(
    bbox,
    half_label,
    title,
    png_name,
    zoom_box=False,
    font_scale=1.0,
    ticks=True,
    height_in=60,
    logo_frac=0.22,
):
    """Draw the three return periods over one half of the river, and save it.

    `font_scale` enlarges every label without changing the pixel count, which is
    what makes a copy legible once a feed scales the image down to its column
    width -- see the posting-size export at the end of the notebook.
    """
    # 60in holds three 2.2:1 strips at full width; past that the panels stop
    # growing (they are width-bound) and the extra only becomes gap. Enlarged
    # text needs more room for itself, hence the taller `height_in` there.
    fig = plt.figure(figsize=(40, height_in), dpi=70, constrained_layout=True)
    gs = fig.add_gridspec(3, 1)
    panels = []
    # The bar is shared, so hold its mappable deliberately rather than reading
    # whatever `glyph` the loop happens to leave bound. Every panel gets the
    # same cmap and clim below, so the first one speaks for all of them.
    shared_im = None
    for rp_i, (rp, grid) in enumerate(zip(PERIODS, whole)):
        piece = grid.crop(bbox=bbox, touch=False)
        ax = fig.add_subplot(gs[rp_i, 0])
        panels.append(ax)
        glyph = piece.plot(
            fig=fig,
            ax=ax,
            cmap=DEPTH_CMAP,
            colorbar=False,
            basemap=BASEMAP,
            data_style=OVERLAY,
        )
        glyph.im.set_clim(0, vmax)
        shared_im = shared_im or glyph.im
        ax.set_title(
            f"{RHINE}  \u00b7  RP{rp}  \u00b7  {half_label}",
            fontsize=30 * font_scale,
        )
        if ticks:
            ax.xaxis.set_ticks_position("bottom")
            ax.tick_params(labelsize=20 * font_scale)
            ax.set_ylabel("latitude", fontsize=22 * font_scale)
        else:
            # Tick LABELS are unreadable once a feed scales the image down, so
            # at posting size they are clutter rather than information. Only the
            # labels go: dropping the locators with `set_xticks([])` would throw
            # away the georeferencing that makes these axes lon/lat.
            ax.tick_params(labelbottom=False, labelleft=False)
        # The numbers are for the whole river, not this half, so say so rather
        # than letting them read as a measurement of the panel above.
        stat = measured[RHINE][rp_i]
        d_area, d_deep = growth(measured[RHINE], rp_i)
        change = (
            "RP10 baseline"
            if rp_i == 0
            else f"area {d_area:+.0f}%   depth {d_deep:+.0f}%  against RP10"
        )
        ax.set_xlabel(
            f"whole river at RP{rp}: {stat['cells']:,} cells  |  "
            f"mean {stat['mean']:.2f} m  |  {change}",
            fontsize=22 * font_scale,
        )
        if zoom_box:
            ax.add_patch(
                Rectangle(
                    (ZOOM[0], ZOOM[1]),
                    ZOOM[2] - ZOOM[0],
                    ZOOM[3] - ZOOM[1],
                    fill=False,
                    edgecolor="#c1121f",
                    # A hairline survives print but vanishes when the image is
                    # scaled to a feed column, so it scales with the text.
                    linewidth=3.2 * font_scale,
                    zorder=5,
                )
            )
            if rp_i == 0:
                # `text`, not `annotate`: an annotation anchors on `xy`, and at
                # lon 5.1 that point falls just outside the cropped half, which
                # makes matplotlib clip the whole label away. Satellite tiles are
                # dark and busy, so the label needs its own ground.
                ax.text(
                    ZOOM[0] + 0.01,
                    ZOOM[3] + 0.012,
                    "the last figure covers this box",
                    color="#c1121f",
                    fontsize=24 * font_scale,
                    fontweight="bold",
                    va="bottom",
                    ha="left",
                    bbox={
                        "boxstyle": "round,pad=0.3",
                        "facecolor": "white",
                        "alpha": 0.85,
                        "edgecolor": "#c1121f",
                    },
                )
    add_depth_bar(fig, shared_im, font_scale)
    fig.suptitle(title, fontsize=38 * font_scale)
    # Stamped after the title so the layout it measures is the final one. Saved
    # without `bbox_inches="tight"`, which would rescale the mark's margin.
    stamp_logo(fig, panels[0], height_frac=logo_frac)
    path = out / png_name
    fig.savefig(path, dpi=110)
    print(f"written to {path}")
    return fig


fig_upstream = draw_rhine_half(
    UPSTREAM,
    "upstream, Lobith to Gorinchem",
    UP_TITLE,
    "dutch_rhine_upstream.png",
)
plt.show()
# These figures are large; release each once shown so a run of this
# notebook does not carry all three at once.
plt.close(fig_upstream)

# NBVAL_SKIP # All three figures below share this colour range, so a shade means the same # depth in every one. That is the whole point of the comparison, so it is # computed once here rather than per figure. vmax = max(row["max"] for rows in measured.values() for row in rows) BASEMAP = "Esri.WorldImagery" OVERLAY = DataStyle(alpha=0.85) # Titles live here because each figure is rendered twice -- once at print scale # and once, further down, with the fonts enlarged for posting. UP_TITLE = ( "Upstream Dutch Rhine \u00b7 Lobith to Gorinchem\n" "JRC European Flood Hazard Map \u00b7 Lobith to Gorinchem" ) DOWN_TITLE = ( "Downstream Dutch Rhine \u00b7 Gorinchem to the North Sea\n" "JRC European Flood Hazard Map \u00b7 Gorinchem to the North Sea" ) ZOOM_TITLE = ( "A closer window \u00b7 its area barely grows because the box is already full\n" "the red box above, at the same colour scale" ) def growth(rows, index): """Return the area / depth change of one return period against RP10.""" first, row = rows[0], rows[index] return ( 100 * (row["cells"] / first["cells"] - 1), 100 * (row["mean"] / first["mean"] - 1), ) def add_depth_bar(fig, mappable, font_scale=1.0): """Put the shared depth scale along the bottom of a figure. Built with `fig.colorbar` rather than through the glyph's own `.cbar`: each panel's glyph owns one bar for one axes, and what these figures need is a single bar spanning every panel so one shade means one depth across all of them. The mappable handed over is still the glyph's `.im`, so the bar tracks the range set by `set_clim`. """ bar = fig.colorbar( mappable, ax=fig.axes, orientation="horizontal", location="bottom", shrink=0.45, pad=0.01, aspect=50, ) bar.set_label("river-flood water depth (m)", size=26 * font_scale) bar.set_ticks(range(0, math.floor(vmax) + 1)) bar.ax.tick_params(labelsize=22 * font_scale) return bar def stamp_logo(fig, ax, height_frac=0.22, pad_in=0.35): """Stamp the earthlens mark inside a map panel's top-right corner.""" # A mark in the figure margin is an orphan; a watermark belongs on the map. # `stamp_mark` anchors to a FIGURE corner, so the offset that lands it inside # `ax` is measured from the figure edge across to the panel edge. Panel # positions are only final after the layout pass, hence `execute` first. fig.get_layout_engine().execute(fig) fig_w, fig_h = fig.get_size_inches() box = ax.get_position() # `frac` is read as the mark's longer on-figure side, judged on FRACTIONS # rather than inches -- so derive both candidates and hand over the larger. height = box.height * height_frac width = height * fig_h / (fig_w * LOGO_ASPECT) return stamp_mark( fig, LOGO, frac=max(width, height), corner="upper right", margin=(1 - box.x1 + pad_in / fig_w, 1 - box.y1 + pad_in / fig_h), shadow=False, ) # Map panels are aspect-locked, so at full width their height is fixed by their # own shape -- the whole river is 4.4:1, and no `figsize` can make that strip # taller without stretching the geography. Cutting the river in half at Gorinchem # halves the aspect to 2.2:1, which is what doubles the height each map gets. # The two halves then get a figure each, so neither is a mile-long scroll. MID = 5.1 UPSTREAM = [MID, LAT_LIM[0], LON_LIM[1], LAT_LIM[1]] DOWNSTREAM = [LON_LIM[0], LAT_LIM[0], MID, LAT_LIM[1]] def draw_rhine_half( bbox, half_label, title, png_name, zoom_box=False, font_scale=1.0, ticks=True, height_in=60, logo_frac=0.22, ): """Draw the three return periods over one half of the river, and save it. `font_scale` enlarges every label without changing the pixel count, which is what makes a copy legible once a feed scales the image down to its column width -- see the posting-size export at the end of the notebook. """ # 60in holds three 2.2:1 strips at full width; past that the panels stop # growing (they are width-bound) and the extra only becomes gap. Enlarged # text needs more room for itself, hence the taller `height_in` there. fig = plt.figure(figsize=(40, height_in), dpi=70, constrained_layout=True) gs = fig.add_gridspec(3, 1) panels = [] # The bar is shared, so hold its mappable deliberately rather than reading # whatever `glyph` the loop happens to leave bound. Every panel gets the # same cmap and clim below, so the first one speaks for all of them. shared_im = None for rp_i, (rp, grid) in enumerate(zip(PERIODS, whole)): piece = grid.crop(bbox=bbox, touch=False) ax = fig.add_subplot(gs[rp_i, 0]) panels.append(ax) glyph = piece.plot( fig=fig, ax=ax, cmap=DEPTH_CMAP, colorbar=False, basemap=BASEMAP, data_style=OVERLAY, ) glyph.im.set_clim(0, vmax) shared_im = shared_im or glyph.im ax.set_title( f"{RHINE} \u00b7 RP{rp} \u00b7 {half_label}", fontsize=30 * font_scale, ) if ticks: ax.xaxis.set_ticks_position("bottom") ax.tick_params(labelsize=20 * font_scale) ax.set_ylabel("latitude", fontsize=22 * font_scale) else: # Tick LABELS are unreadable once a feed scales the image down, so # at posting size they are clutter rather than information. Only the # labels go: dropping the locators with `set_xticks([])` would throw # away the georeferencing that makes these axes lon/lat. ax.tick_params(labelbottom=False, labelleft=False) # The numbers are for the whole river, not this half, so say so rather # than letting them read as a measurement of the panel above. stat = measured[RHINE][rp_i] d_area, d_deep = growth(measured[RHINE], rp_i) change = ( "RP10 baseline" if rp_i == 0 else f"area {d_area:+.0f}% depth {d_deep:+.0f}% against RP10" ) ax.set_xlabel( f"whole river at RP{rp}: {stat['cells']:,} cells | " f"mean {stat['mean']:.2f} m | {change}", fontsize=22 * font_scale, ) if zoom_box: ax.add_patch( Rectangle( (ZOOM[0], ZOOM[1]), ZOOM[2] - ZOOM[0], ZOOM[3] - ZOOM[1], fill=False, edgecolor="#c1121f", # A hairline survives print but vanishes when the image is # scaled to a feed column, so it scales with the text. linewidth=3.2 * font_scale, zorder=5, ) ) if rp_i == 0: # `text`, not `annotate`: an annotation anchors on `xy`, and at # lon 5.1 that point falls just outside the cropped half, which # makes matplotlib clip the whole label away. Satellite tiles are # dark and busy, so the label needs its own ground. ax.text( ZOOM[0] + 0.01, ZOOM[3] + 0.012, "the last figure covers this box", color="#c1121f", fontsize=24 * font_scale, fontweight="bold", va="bottom", ha="left", bbox={ "boxstyle": "round,pad=0.3", "facecolor": "white", "alpha": 0.85, "edgecolor": "#c1121f", }, ) add_depth_bar(fig, shared_im, font_scale) fig.suptitle(title, fontsize=38 * font_scale) # Stamped after the title so the layout it measures is the final one. Saved # without `bbox_inches="tight"`, which would rescale the mark's margin. stamp_logo(fig, panels[0], height_frac=logo_frac) path = out / png_name fig.savefig(path, dpi=110) print(f"written to {path}") return fig fig_upstream = draw_rhine_half( UPSTREAM, "upstream, Lobith to Gorinchem", UP_TITLE, "dutch_rhine_upstream.png", ) plt.show() # These figures are large; release each once shown so a run of this # notebook does not carry all three at once. plt.close(fig_upstream)

4 - Downstream, Gorinchem to the North Sea¶

The western half, where the Rhine becomes the Merwede and fans out through the delta. The red box marks the stretch the last figure zooms into.

In [ ]:

# NBVAL_SKIP
fig_downstream = draw_rhine_half(
    DOWNSTREAM,
    "downstream, Gorinchem to the North Sea",
    DOWN_TITLE,
    "dutch_rhine_downstream.png",
    zoom_box=True,
)
plt.show()
# These figures are large; release each once shown so a run of this
# notebook does not carry all three at once.
plt.close(fig_downstream)

# NBVAL_SKIP fig_downstream = draw_rhine_half( DOWNSTREAM, "downstream, Gorinchem to the North Sea", DOWN_TITLE, "dutch_rhine_downstream.png", zoom_box=True, ) plt.show() # These figures are large; release each once shown so a run of this # notebook does not carry all three at once. plt.close(fig_downstream)

5 - A closer window, and why its number misleads¶

The same three return periods over the boxed stretch, on the same colour scale. The flooded outline barely moves - about seven per cent - while the whole scene darkens.

Before reading that as confinement, check how much of the box was dry to begin with. The cell below measures it, then widens the box until there is dry ground left to take. The depth increase is real; watch what happens to the area curve.

In [ ]:

# NBVAL_SKIP
# A flooded area can only grow into ground that is still dry, so a box that is
# nearly full caps its own growth no matter what the water does. Measured on the
# same rasters, widening the box restores the whole-river figure -- which is what
# tells you the flat curve belongs to the window, not to the dikes.
#
# Per-cell ground area on a geographic grid is a pyramids primitive
# (serapeum-org/pyramids#1085, in the 0.65 floor): `domain_area` sums the valid
# (non-nodata) cells -- the flooded ones here -- and `cell_area` gives every
# cell's area, so the whole box is their sum. Both integrate over the CRS's own
# ellipsoid, so the notebook no longer hand-rolls the spherical maths.


def box_coverage(grid):
    """Return (flooded km2, whole-box km2), each cell weighted by its latitude."""
    return grid.domain_area(unit="km2"), float(grid.cell_area(unit="km2").sum())


WIDENING = {
    "the zoom above": ZOOM,
    "wider": [4.3, 51.6, 5.4, 52.1],
    "widest": [4.0, 51.6, 5.5, 52.1],
    "whole area of interest": None,
}

print(f"{'box':24}{'% full at RP10':>16}{'area growth to RP500':>22}")
for label, bbox in WIDENING.items():
    first = whole[0].crop(bbox=bbox, touch=False) if bbox else whole[0]
    last = whole[-1].crop(bbox=bbox, touch=False) if bbox else whole[-1]
    wet_first, box_first = box_coverage(first)
    wet_last, _ = box_coverage(last)
    print(
        f"{label:24}{100 * wet_first / box_first:15.1f}%"
        f"{100 * (wet_last / wet_first - 1):+21.1f}%"
    )

# NBVAL_SKIP # A flooded area can only grow into ground that is still dry, so a box that is # nearly full caps its own growth no matter what the water does. Measured on the # same rasters, widening the box restores the whole-river figure -- which is what # tells you the flat curve belongs to the window, not to the dikes. # # Per-cell ground area on a geographic grid is a pyramids primitive # (serapeum-org/pyramids#1085, in the 0.65 floor): `domain_area` sums the valid # (non-nodata) cells -- the flooded ones here -- and `cell_area` gives every # cell's area, so the whole box is their sum. Both integrate over the CRS's own # ellipsoid, so the notebook no longer hand-rolls the spherical maths. def box_coverage(grid): """Return (flooded km2, whole-box km2), each cell weighted by its latitude.""" return grid.domain_area(unit="km2"), float(grid.cell_area(unit="km2").sum()) WIDENING = { "the zoom above": ZOOM, "wider": [4.3, 51.6, 5.4, 52.1], "widest": [4.0, 51.6, 5.5, 52.1], "whole area of interest": None, } print(f"{'box':24}{'% full at RP10':>16}{'area growth to RP500':>22}") for label, bbox in WIDENING.items(): first = whole[0].crop(bbox=bbox, touch=False) if bbox else whole[0] last = whole[-1].crop(bbox=bbox, touch=False) if bbox else whole[-1] wet_first, box_first = box_coverage(first) wet_last, _ = box_coverage(last) print( f"{label:24}{100 * wet_first / box_first:15.1f}%" f"{100 * (wet_last / wet_first - 1):+21.1f}%" )

In [ ]:

# NBVAL_SKIP
# Two panels share a row and the third sits centred beneath. Four columns is what
# makes that centring exact -- with three, the last panel could only sit under the
# first or the second. Each panel spans half the width, so it is ~11.7in tall.
ZOOM_SLOTS = [(0, slice(0, 2)), (0, slice(2, 4)), (1, slice(1, 3))]


def draw_zoom(png_name, font_scale=1.0, height_in=28, logo_frac=0.22, ticks=True):
    """Draw the three return periods over the zoom box, and save it."""
    fig = plt.figure(figsize=(40, height_in), dpi=70, constrained_layout=True)
    gs = fig.add_gridspec(2, 4)
    panels = []
    # As above: the shared bar's mappable is chosen, not inherited from the last
    # iteration of the loop.
    shared_im = None
    for col, (rp, grid) in enumerate(zip(PERIODS, zoomed)):
        row, span = ZOOM_SLOTS[col]
        ax = fig.add_subplot(gs[row, span])
        panels.append(ax)
        glyph = grid.plot(
            fig=fig,
            ax=ax,
            cmap=DEPTH_CMAP,
            colorbar=False,
            basemap=BASEMAP,
            data_style=OVERLAY,
        )
        glyph.im.set_clim(0, vmax)
        shared_im = shared_im or glyph.im
        stat = measured[DELTA][col]
        d_area, d_deep = growth(measured[DELTA], col)
        change = (
            "RP10 baseline"
            if col == 0
            else f"area {d_area:+.0f}%   depth {d_deep:+.0f}%"
        )
        ax.set_title(f"{DELTA}  \u00b7  RP{rp}", fontsize=26 * font_scale)
        ax.set_xlabel(
            f"{stat['cells']:,} cells  |  mean {stat['mean']:.2f} m\n{change}",
            fontsize=20 * font_scale,
        )
        if ticks:
            ax.xaxis.set_ticks_position("bottom")
            ax.tick_params(labelsize=17 * font_scale)
            if span.start == 0 or row == 1:
                ax.set_ylabel("latitude", fontsize=19 * font_scale)
        else:
            # Tick LABELS are unreadable once a feed scales the image down, so
            # at posting size they are clutter rather than information. Only the
            # labels go: dropping the locators with `set_xticks([])` would throw
            # away the georeferencing that makes these axes lon/lat.
            ax.tick_params(labelbottom=False, labelleft=False)
    add_depth_bar(fig, shared_im, font_scale)
    fig.suptitle(ZOOM_TITLE, fontsize=31 * font_scale)
    stamp_logo(fig, panels[0], height_frac=logo_frac)
    path = out / png_name
    fig.savefig(path, dpi=110)
    print(f"written to {path}")
    return fig


fig_zoom = draw_zoom("dutch_rhine_zoom.png")
plt.show()
# These figures are large; release each once shown so a run of this
# notebook does not carry all three at once.
plt.close(fig_zoom)

# NBVAL_SKIP # Two panels share a row and the third sits centred beneath. Four columns is what # makes that centring exact -- with three, the last panel could only sit under the # first or the second. Each panel spans half the width, so it is ~11.7in tall. ZOOM_SLOTS = [(0, slice(0, 2)), (0, slice(2, 4)), (1, slice(1, 3))] def draw_zoom(png_name, font_scale=1.0, height_in=28, logo_frac=0.22, ticks=True): """Draw the three return periods over the zoom box, and save it.""" fig = plt.figure(figsize=(40, height_in), dpi=70, constrained_layout=True) gs = fig.add_gridspec(2, 4) panels = [] # As above: the shared bar's mappable is chosen, not inherited from the last # iteration of the loop. shared_im = None for col, (rp, grid) in enumerate(zip(PERIODS, zoomed)): row, span = ZOOM_SLOTS[col] ax = fig.add_subplot(gs[row, span]) panels.append(ax) glyph = grid.plot( fig=fig, ax=ax, cmap=DEPTH_CMAP, colorbar=False, basemap=BASEMAP, data_style=OVERLAY, ) glyph.im.set_clim(0, vmax) shared_im = shared_im or glyph.im stat = measured[DELTA][col] d_area, d_deep = growth(measured[DELTA], col) change = ( "RP10 baseline" if col == 0 else f"area {d_area:+.0f}% depth {d_deep:+.0f}%" ) ax.set_title(f"{DELTA} \u00b7 RP{rp}", fontsize=26 * font_scale) ax.set_xlabel( f"{stat['cells']:,} cells | mean {stat['mean']:.2f} m\n{change}", fontsize=20 * font_scale, ) if ticks: ax.xaxis.set_ticks_position("bottom") ax.tick_params(labelsize=17 * font_scale) if span.start == 0 or row == 1: ax.set_ylabel("latitude", fontsize=19 * font_scale) else: # Tick LABELS are unreadable once a feed scales the image down, so # at posting size they are clutter rather than information. Only the # labels go: dropping the locators with `set_xticks([])` would throw # away the georeferencing that makes these axes lon/lat. ax.tick_params(labelbottom=False, labelleft=False) add_depth_bar(fig, shared_im, font_scale) fig.suptitle(ZOOM_TITLE, fontsize=31 * font_scale) stamp_logo(fig, panels[0], height_frac=logo_frac) path = out / png_name fig.savefig(path, dpi=110) print(f"written to {path}") return fig fig_zoom = draw_zoom("dutch_rhine_zoom.png") plt.show() # These figures are large; release each once shown so a run of this # notebook does not carry all three at once. plt.close(fig_zoom)

6 - Posting-size export¶

The figures above are sized for print. This section re-renders all three for a social feed, writing a *_post.png beside each original: the type goes up 2.2x relative to the figure, and the lat/lon ticks come off - at feed size they are unreadable anyway, so dropping them buys back room for the maps. Nothing about the data changes, only the presentation.

In [ ]:

# NBVAL_SKIP
# Two changes make a copy read well once a feed scales it into its column, and
# neither is "more dpi" -- dpi changes the pixel count, not the proportions.
# First the type goes up 2.2x relative to the figure. Second the axis ticks come
# off: at that size they are illegible anyway, so they cost clutter and buy
# nothing, and dropping them lets the maps stay the largest thing in the frame.
# Chasing legible TICK text instead needs about 2.9x and looks shouty -- the
# titles and the numbers under each map are what a reader actually reads.
POST_SCALE = 2.2
FEED_PX = 552

# Built one at a time. These are the largest renders in the notebook -- 62 in
# tall at 2.2x type -- so materialising all three before measuring any would
# undo the close-as-you-go the earlier cells go out of their way to do.
recipes = [
    (
        draw_rhine_half,
        (
            UPSTREAM,
            "upstream, Lobith to Gorinchem",
            UP_TITLE,
            "dutch_rhine_upstream_post.png",
        ),
        {"font_scale": POST_SCALE, "ticks": False, "height_in": 62, "logo_frac": 0.26},
        30,
    ),
    (
        draw_rhine_half,
        (
            DOWNSTREAM,
            "downstream, Gorinchem to the North Sea",
            DOWN_TITLE,
            "dutch_rhine_downstream_post.png",
        ),
        {
            "zoom_box": True,
            "font_scale": POST_SCALE,
            "ticks": False,
            "height_in": 62,
            "logo_frac": 0.26,
        },
        30,
    ),
    (
        draw_zoom,
        ("dutch_rhine_zoom_post.png",),
        {"font_scale": POST_SCALE, "ticks": False, "height_in": 29, "logo_frac": 0.26},
        26,
    ),
]

print()
# `base_pt` travels with its recipe: each builder sets its own panel-title size,
# and carrying it here beats inferring it from the figure's position in a list.
for build, args, kwargs, base_pt in recipes:
    fig = build(*args, **kwargs)
    w_in, h_in = fig.get_size_inches()
    # A feed fits a portrait image inside about 4:5, so anything taller than 1.25
    # is shown narrower than the full column -- which costs back some legibility.
    fitted = min(FEED_PX, FEED_PX * 1.25 / (h_in / w_in))
    title_px = base_pt * POST_SCALE / (w_in * 72) * fitted
    print(
        f"  aspect {h_in / w_in:.2f}  ->  shown {fitted:.0f}px wide  "
        f"->  panel title {title_px:.1f}px"
    )
    # Written to disk inside the builder, so close before the next one is drawn
    # and only one of these figures is ever resident.
    plt.close(fig)

# NBVAL_SKIP # Two changes make a copy read well once a feed scales it into its column, and # neither is "more dpi" -- dpi changes the pixel count, not the proportions. # First the type goes up 2.2x relative to the figure. Second the axis ticks come # off: at that size they are illegible anyway, so they cost clutter and buy # nothing, and dropping them lets the maps stay the largest thing in the frame. # Chasing legible TICK text instead needs about 2.9x and looks shouty -- the # titles and the numbers under each map are what a reader actually reads. POST_SCALE = 2.2 FEED_PX = 552 # Built one at a time. These are the largest renders in the notebook -- 62 in # tall at 2.2x type -- so materialising all three before measuring any would # undo the close-as-you-go the earlier cells go out of their way to do. recipes = [ ( draw_rhine_half, ( UPSTREAM, "upstream, Lobith to Gorinchem", UP_TITLE, "dutch_rhine_upstream_post.png", ), {"font_scale": POST_SCALE, "ticks": False, "height_in": 62, "logo_frac": 0.26}, 30, ), ( draw_rhine_half, ( DOWNSTREAM, "downstream, Gorinchem to the North Sea", DOWN_TITLE, "dutch_rhine_downstream_post.png", ), { "zoom_box": True, "font_scale": POST_SCALE, "ticks": False, "height_in": 62, "logo_frac": 0.26, }, 30, ), ( draw_zoom, ("dutch_rhine_zoom_post.png",), {"font_scale": POST_SCALE, "ticks": False, "height_in": 29, "logo_frac": 0.26}, 26, ), ] print() # `base_pt` travels with its recipe: each builder sets its own panel-title size, # and carrying it here beats inferring it from the figure's position in a list. for build, args, kwargs, base_pt in recipes: fig = build(*args, **kwargs) w_in, h_in = fig.get_size_inches() # A feed fits a portrait image inside about 4:5, so anything taller than 1.25 # is shown narrower than the full column -- which costs back some legibility. fitted = min(FEED_PX, FEED_PX * 1.25 / (h_in / w_in)) title_px = base_pt * POST_SCALE / (w_in * 72) * fitted print( f" aspect {h_in / w_in:.2f} -> shown {fitted:.0f}px wide " f"-> panel title {title_px:.1f}px" ) # Written to disk inside the builder, so close before the next one is drawn # and only one of these figures is ever resident. plt.close(fig)

Notes¶

• The basemap is context, not data. Tiles are Esri World Imagery, (C) Esri and the GIS User Community; the flood layer sits over them at 85% opacity. Fetching them needs network at render time - drop basemap / data_style from the plot calls to run the figures offline.
• Three figures, one colour range. They are separate so each can be read (and posted) on its own, but vmax is computed once across all of them, so a shade means the same depth in every one. The zoom never reaches the dark end - its deepest water is about 5.5 m against 9.2 m along the wider river - which leaves it looking pale. A per-figure scale would fix that and make them incomparable, which is the one thing this notebook exists to do.
• A small window can manufacture a finding. Measured on the zoom alone this would read as "a rarer flood does not spread", which is false for the river as a whole. The saturation cell in section 5 computes the numbers quoted here -- the zoom box is 88.9% full at RP10 and its flooded area grows +7.1% out to RP500, against +27%, +35% and +33% as the box widens. Being nearly full is what holds it down: a box already 88.9% covered cannot grow by more than about twelve per cent however rare the flood, so +7.1% is close to the most it could have shown. They are repeated in the prose because the docs site does not execute notebooks, so a reader who has not run this sees no output; run the cell and they should match. Check how much of an area of interest is already covered before reading a small change as a result.
• Licence. The EFHM is CC-BY-4.0. Cite Dottori, F., Alfieri, L., Bianchi, A., Skoulikaris, C., Salamon, P. (2020), River flood hazard maps for Europe and the Mediterranean Basin region, JRC / Copernicus Emergency Management Service.