Tip

For an interactive online version click here:

Elevation indices

Here we assume that flow directions are known. We read the flow direction raster data, including meta-data, using rasterio and parse it to a pyflwdir FlwDirRaster object, see earlier examples for more background.

# import pyflwdir, some dependencies and convenience methods
import numpy as np
import rasterio
import pyflwdir

# local convenience methods (see utils.py script in notebooks folder)
from utils import quickplot, plt  # data specific quick plot method

# read and parse flow direciton data
with rasterio.open("rhine_d8.tif", "r") as src:
    flwdir = src.read(1)
    crs = src.crs
    extent = np.array(src.bounds)[[0, 2, 1, 3]]
    flw = pyflwdir.from_array(
        flwdir,
        ftype="d8",
        transform=src.transform,
        latlon=crs.is_geographic,
        cache=True,
    )
# read elevation data
with rasterio.open("rhine_elv0.tif", "r") as src:
    elevtn = src.read(1)

height above nearest drain (HAND)

The hand() method uses drainage-normalized topography and flowpaths to delineate the relative vertical distances (drop) to the nearest river (drain) as a proxy for the potential extent of flooding (Nobre et al. 2016). The pyflwdir implementation requires stream mask drain and elevation raster elevtn. The stream mask is typically determined based on a threshold on upstream_area() or stream_order(), but can also be set from rasterizing a vector stream file.

# first we derive the upstream area map
uparea = flw.upstream_area("km2")

# HAND based on streams defined by a minimal upstream area of 1000 km2
hand = flw.hand(drain=uparea > 1000, elevtn=elevtn)
# plot
ax = quickplot(title="Height above nearest drain (HAND)")
im = ax.imshow(
    np.ma.masked_equal(hand, -9999),
    extent=extent,
    cmap="gist_earth_r",
    alpha=0.5,
    vmin=0,
    vmax=150,
)
fig = plt.gcf()
cax = fig.add_axes([0.82, 0.37, 0.02, 0.12])
fig.colorbar(im, cax=cax, orientation="vertical")
cax.set_ylabel("HAND [m]")
plt.savefig("hand.png")

Floodplains

The floodplains() method delineates geomorphic floodplain boundaries based on a power-law relation between upstream area and a maximum HAND contour as developed by Nardi et al (2019). Here, streams are defined based on a minimum upstream area threshold upa_min and floodplains on the scaling parameter b of the power-law relationship.

floodplains = flw.floodplains(elevtn=elevtn, uparea=uparea, upa_min=1000)
# plot
floodmap = (floodplains, -1, dict(cmap="Blues", alpha=0.5, vmin=0))
ax = quickplot(
    raster=floodmap, title="Geomorphic floodplains", filename="flw_floodplain"
)