Tip

For an interactive online version click here:

Flow direction data

The FlwdirRaster object is at the core of pyflwdir. It stores gridded flow-direction data in a common, actionable format, where each cell is represented by the linear index of its next downstream cell.

We currently support two local flow-direction (D8) formats: ArcGIS D8 and PCRaster LDD (see figure). We also support one global flow-direction format: CaMa-Flood NEXTXY. Local formats describe the downstream cell by its relative direction to a neighboring cell, while the global format describes it by row and column indices.

We read the flow direction raster data, including meta-data, using rasterio

import rasterio

with rasterio.open("rhine_d8.tif", "r") as src:
    flwdir = src.read(1)
    transform = src.transform
    crs = src.crs
    latlon = crs.to_epsg() == 4326

Next, we parse this data to a FlwdirRaster object, the core object to work with flow direction data. In this step the D8 data is parsed to an actionable format.

NOTE: The first call to most methods may be slow because Numba compiles the code just in time. Subsequent calls to the same method, even with different arguments, are usually much faster!

import pyflwdir

flw = pyflwdir.from_array(
    flwdir, ftype="d8", transform=transform, latlon=latlon, cache=True
)

# When printing the FlwdirRaster instance we see its attributes.
print(flw)

{'ftype': 'd8',
 'idxs_ds': array([-1, -1, -1, ..., -1, -1, -1], shape=(679954,), dtype=int32),
 'idxs_pit': array([20994], dtype=int32),
 'idxs_seq': None,
 'latlon': True,
 'nnodes': 349847,
 'shape': (682, 997),
 'transform': Affine(0.008333333333325754, 0.0, 3.5666666664997138,
       0.0, -0.008333333333339965, 52.00833333330708)}

We can then use the many methods available on the FlwdirRaster object. See the FlwdirRaster API.

To visualize the flow directions, we derive a vector stream network using the streams() method. Each line represents a stream segment with a Strahler order of at least min_sto, as computed by stream_order(). The line features are converted to a GeoDataFrame for visualization.

import geopandas as gpd

feats = flw.streams(min_sto=4)
gdf = gpd.GeoDataFrame.from_features(feats, crs=crs)
gdf.head()

	geometry	idx	idx_ds	pit	strord
0	LINESTRING (8.7125 46.65417, 8.72083 46.6625, ...	640691	635723	False	4
1	LINESTRING (8.8375 46.62917, 8.84583 46.6375, ...	643697	635723	False	4
2	LINESTRING (8.17917 46.57083, 8.1875 46.57083,...	650597	634651	False	4
3	LINESTRING (8.85417 46.69583, 8.8625 46.69583,...	635723	632737	False	5
4	LINESTRING (8.87917 46.7625, 8.8875 46.7625, 8...	627750	632737	False	4

import numpy as np
from matplotlib import cm, colors

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

# keyword arguments passed to GeoDataFrame.plot()
gdf_plot_kwds = {
    "column": "strord",
    "cmap": colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))),
}
# plot streams with hillshade from elevation data (see utils.py)
ax = quickplot(gdfs=[(gdf, gdf_plot_kwds)], title="Streams")