Generating TauDEM Stream Networks

This workflow outlines the process of generating riverscapes using only a high-resolution DEM. It is particularly useful for working with LiDAR datasets or other detailed elevation models to create riverscapes. Developed by Philip Bailey in March 2025, this approach builds upon Jordan Gilbert's work with TauDEM for Riverscapes tools.

The motivation is to be able to build a stream network that can be used within the Riverscapes tools waterfall for generating riverscapes polygons. With nothing more than a DEM, this workflow is the first step to running the Channel Area, tauDEM and VBET portion of the tools waterfall.

Preparation

1. Launch a GitHub CodeSpace in the Riverscapes Tools repository. Experiment with different CodeSpace sizes to find the best performance. This workflow was developed using the largest size available and all tools ran in under 1 minute on 350Mb rasters.
2. Create a data processing folder at /workspace/data and open a new Visual Studio Code window in that folder.
3. Upload the high resolution DEM to the /workspace/data folder. (Simply drag the file into the Visual Studio Code window at the appropriate location.)
4. Open a terminal in the /workspace/data folder.

This workflow was developed on a LiDAR DEM for Cole Creek in southwestern Wyoming. The raster was provided by the Bureau of Land Management. The cartographic projection is EPSG:26913, UTM zone 13N.

Step 1. Pitfill the DEM

The following command will pitfill the DEM, ensuring that any internal sinks are filled. This is a critical step for ensuring that the flow direction algorithm works correctly.

Step 2. Generate Flow Direction And Slope

This process uses the D8 flow direction. tauDEM also has a D-infinity flow direction algorithm, but I couldn't get this to work with the final network and watersheds generation steps.

Step 3. Generate Contributing Area

It is critical during this step, to understand if all cells in the DEM slow towards a single outlet, or if areas of the DEM flow outwards in multiple directions. It's common when running LiDAR images for partial watersheds that some areas of the DEM will flow outwards in multiple directions. tauDEM includes the -nc switch for handling these situations and ignoring the areas of the DEM that flow outwards in multiple directions.

Step 4. Stream Network Raster

The stream network raster is generated by thresholding the flow accumulation raster. The threshold value is determined by the user and should be based on the size of the watershed and the resolution of the DEM. It can take some trial and error to find the right threshold value.

Step 5. Stream Network And watersheds

This is the critical step that generates the stream network and watersheds vector datasets.

Note that this tauDEM command took several iterations to make work in the CodeSpace. Omitting some of the output files caused a segmentation fault, perhaps because tauDEM attempts to create these outputs anyway and has trouble within permissions in the path it chooses. Once all the outputs were included, the command ran successfully.

Step 6. Height Above Nearest Drainage (HAND)

The height above nearest drainage (HAND) raster is generated by computing the height of each cell above the nearest stream. This raster can be useful for identifying areas of the landscape that are far from the stream network.

Optional Steps

The script below also includes the steps to generate a Height Above Nearest Drainage (HAND) raster, and then combine this with the other rasters generated above into a single multiband raster. This multiband raster is intended for image classification experiments to find valley bottoms automatically.

Script

#!/bin/bash

# Number of processing cores (adjust as needed)
CORES=4

# Input DEM file
DEM="dem.tif"

# Output files
FILLED_DEM="dem_filled.tif"
D8_FLOW="d8_flow.tif"
D8_SLOPE="d8_slp.tif"
D8_CONT_AREA="d8_cont_area.tif"
STREAM_RASTER="stream_raster.tif"
STREAM_ORDER="stream_order.tif"
STREAM_NETWORK="stream_network.gpkg"
WATERSHED="watershed.tif"
SUBWATERSHEDS="subwatersheds.tif"
HAND="hand.tif"

# 1. Pit Remove (Fill Sinks)
# Removes depressions from the DEM to ensure continuous flow.
mpirun -np $CORES pitremove -z $DEM -fel $FILLED_DEM

# 2a. D∞ Flow Directions & Slope
# mpirun -np $CORES dinfflowdir -fel $FILLED_DEM -ang $DINF_ANGLE -slp $DINF_SLOPE

# 2b. D8 Flow Directions & Slope
mpirun -np $CORES d8flowdir -fel $FILLED_DEM -p $D8_FLOW -sd8 $D8_SLOPE

# 3a. D∞ Contributing Area
# mpirun -np $CORES areadinf -ang $DINF_ANGLE -sca $DINF_CONTRIB -nc

# 3b. D8 Contributing Area
# Computes the contributing area based on the D8 flow direction.
mpirun -np $CORES aread8 -p $D8_FLOW -ad8 $D8_CONT_AREA -nc

# 4. Stream Raster (Thresholding)
# Manually adjust the threshold value based on inspection.
THRESHOLD=500000
# Generates a binary raster of streams using the selected threshold.
mpirun -np $CORES threshold -ssa $DINF_CONTRIB -src $STREAM_RASTER -thresh $THRESHOLD

# 5. Stream Network Extraction
mpirun -np $CORES streamnet -p $D8_FLOW -fel $FILLED_DEM -ad8 $D8_CONT_AREA -src $STREAM_RASTER \
  -net $STREAM_NETWORK -netlyr network -ord $STREAM_ORDER -tree stream_tree.dat -coord stream_coord.dat -w $SUBWATERSHEDS

# 6. Optional Height Above Nearest Drainage (HAND)
# mpirun -np $CORES dinfdistdown -fel $FILLED_DEM -ang $DINF_ANGLE -src $STREAM_RASTER -dd $HAND -m ave v

# gdalbuildvrt -separate $MERGEVRT $D8_FLOW $D8_SLOPE $D8_CONT_AREA $STREAM_RASTER $STREAM_ORDER $HAND

# gdal_translate -ot Float32 -co COMPRESS=LZW $MERGEVRT $MERGED_RASTER

Known issues

The output stream network can have several issues. The most common is that the stream network is not continuous, but almost as common are small tributaries in flat areas (see below). The network features are also extremely fragmented with lots of small geometries. These issues can be addressed by manually editing the stream network in a GIS. Additional research is needed to determine if these issues can be addressed in tauDEM itself.