Topographic Position based Stream definition (TPS)

Python code for stream definition in a Digital Elevation Model (DEM) conditioned by Topographic Position Index (TPI).

Reference:

Rafael Barbedo, Vinícius Siqueira & Walter Collischonn (2022). Topographic Position-based Stream definition (TPS): A simple method to address spatial variability of drainage density in stream networks. Hydrological Sciences Journal.

Requirements

Python > 3.5
Numpy
GDAL > 3
Shapely
Fiona

Inputs required

Digital Elevation Model (DEM)
Flow direction (FDR) - must be in HydroSHEDS format (1=E, 2=SE, 4=S, 8=SW, 16=W, 32=NW, 64=N, 128=NE)
Flow accumulation (FAC) - in Area OR cells
Topographic Position Index (TPI) - from Weiss (2001); a radius of 900m works well for this method with a 90m resolution DEM (~10 DEM cells)

All these maps are easily obtained using a GIS processing tool in a DEM (e.g. QGIS). You can also download them via Google Earth Engine. Just click here to check the GEE code (using MERIT Hydro database).

Parameters

TPI threshold (V)
Min Area threshold (Amin)
Max Area threshold (Amax)

Stream is defined where FAC > Amin AND TPI < V OR FAC > Amax. Then, gaps are filled using FDR.

Example usage

Run the tps_functions.py file in your python IDLE
Load the input rasters

# Source file to copy saving format
region = 'Corrente'
srcfn = 'GIS/'+region+'_dem_merit_90.tif'

# Load rasters
dem = gdal.Open(srcfn).ReadAsArray()
fdr = gdal.Open('GIS/'+region+'_fdr_merit_90.tif').ReadAsArray()
fac = gdal.Open('GIS/'+region+'_fac_merit_90.tif').ReadAsArray()
tpi = gdal.Open('GIS/'+region+'_tpi_merit_90.tif').ReadAsArray()

Define parameters

# V (m)
v = -8 

# Amin (km²)
amin = 1.6

# Amax (km²)
amax = 320

Compute stream network with TPS method

str_mask = get_str_mask(fdr, fac, tpi, amax, amin, v)

# Save tif file
array2tif('GIS/'+region+'_str_mask.tif', srcfn, str_mask)

(Optional) Convert stream mask to line vector

str_line = array2gdf(rasref, mask_tpi, fdr, fac)

# Save shp file
str_line.to_file('GIS/'+region+'_str_line.shp')