I was just asking ChatGPT about Python libraries that can calculate D8 flow directions from a digital elevation model, and was pleased to see it suggest some convincing looking code for Landlab! I haven’t checked it works, but at a glance it looks convincing (though why didn’t it use the FlowDirectorD8 component directly?). Any Landlab gurus want to comment?
import numpy as np
from landlab import RasterModelGrid
from landlab.components import FlowAccumulator
import rasterio
# Load DEM from a GeoTIFF
with rasterio.open("input_dem.tif") as src:
dem = src.read(1, masked=True) # Read as masked array
transform = src.transform
dx = transform[0] # cell size in x direction
dy = -transform[4] # cell size in y direction
# Create a Landlab grid (assumes square cells)
nrows, ncols = dem.shape
grid = RasterModelGrid((nrows, ncols), xy_spacing=(dy, dx)) # (y, x)
# Add the topographic data to the grid
grid.add_field("topographic__elevation", dem.filled(dem.mean()), at="node")
grid.at_node["topographic__elevation"][grid.core_nodes] = dem.data[1:-1, 1:-1].flatten()
# Run D8 flow direction + accumulation
flow_accumulator = FlowAccumulator(grid, flow_director="D8")
flow_accumulator.run_one_step()
# D8 directions are stored in grid.at_node['flow__receiver_node']
receiver_nodes = grid.at_node["flow__receiver_node"]
drainage_area = grid.at_node["drainage_area"]
# If you want to export flow directions, you can map them to D8 codes manually
# D8 codes for 8 directions (like richdem): E=1, SE=2, S=4, SW=8, etc.
# Optional: convert node IDs to D8 directions (advanced, not built-in)
I really should learn Landlab…