Tip

For an interactive online version click here:

Flow direction upscaling

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 geopandas as gpd
import numpy as np
import rasterio
from matplotlib import cm, colors

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

import pyflwdir

# read and parse flow direction 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]]
    prof = src.profile
    flw = pyflwdir.from_array(
        flwdir,
        ftype="d8",
        transform=src.transform,
        latlon=crs.is_geographic,
        cache=True,
    )

# vectorize streams for visualization
uparea = flw.upstream_area()
feats0 = flw.streams(uparea > 100, uparea=uparea)
# base color and labels on log10 of upstream area
gdf_stream0 = gpd.GeoDataFrame.from_features(feats0, crs=crs)
gdf_stream0["logupa"] = np.floor(np.log10(gdf_stream0["uparea"])).astype(int)
labels = {2: "1e2-1e3", 3: "1e3-1e4", 4: "1e4-1e5", 5: "1e5-1e6"}
gdf_stream0["loglabs"] = [labels[k] for k in gdf_stream0["logupa"]]
# keyword arguments for GeoDataFrame.plot()
gdf_plt_kwds = {
    "column": "loglabs",
    "cmap": colors.ListedColormap(cm.Blues(np.linspace(0.5, 1, 7))),
    "categorical": True,
    "legend": True,
    "legend_kwds": {"title": "Upstream area [km2]"},
}
title = "Original flow directions (upstream area > 100 km2)"
ax = quickplot(gdfs=[(gdf_stream0, gdf_plt_kwds)], title=title, filename="flw_original")

Flow direction upscaling

Models often have coarser resolution than the elevation data used to build their flow directions. Rather than deriving flow directions from upscaled elevation data, upscale the flow-direction data directly. The upscale() method implements the Iterative Hydrography Upscaling (IHU) algorithm (Eilander et al. 2020). It uses high-resolution flow directions and an upstream-area grid to iteratively select the best stream segment for each upscaled cell. It then traces each selected segment to the next downstream upscaled cell to determine the upscaled flow directions. See the referenced paper for details.

# upscale using scale_factor s
s = 10
flw1, idxs_out = flw.upscale(scale_factor=s, uparea=uparea, method="ihu")

The methods select an outlet pixel, the most downstream high-resolution pixel, from each upscaled cell. These outlet locations can be used with the ucat_area() method to calculate the contributing area of each upscaled cell.

# get the contributing unit area to each upscaled cell & accumulate
subareas = flw.ucat_area(idxs_out=idxs_out, unit="km2")[1]
uparea1 = flw1.accuflux(subareas)

# assess the quality of the upscaling
flwerr = flw.upscale_error(flw1, idxs_out)
percentage_error = np.sum(flwerr == 0) / np.sum(flwerr != 255) * 100
print(f"upscaling error in {percentage_error:.2f}% of cells")

upscaling error in 1.32% of cells

# vectorize streams for visualization
feats1 = flw1.streams(uparea1 > 100, uparea=uparea1)
# base color and labels on log10 of upstream area
gdf_stream1 = gpd.GeoDataFrame.from_features(feats1, crs=crs)
gdf_stream1["logupa"] = np.floor(np.log10(gdf_stream1["uparea"])).astype(int)
gdf_stream1["loglabs"] = [labels[k] for k in gdf_stream1["logupa"]]
# plot
title = f"IHU Upscaled flow directions ({s}x)"
ax = quickplot(
    gdfs=[(gdf_stream1, gdf_plt_kwds)], title=title, filename=f"flw_upscale{s:2d}"
)

Save upscaled flow directions to file

# update the profile and write the upscaled flow direction to a new raster file
prof.update(
    width=flw1.shape[1],
    height=flw1.shape[0],
    transform=flw1.transform,
    nodata=247,
)
with rasterio.open(f"rhine_d8_upscale{s}.tif", "w", **prof) as src:
    src.write(flw1.to_array("d8"), 1)

# save upscaled network to a vector file
gdf_stream1.to_file(f"rhine_d8_upscale{s}.gpkg", driver="GPKG")