WhereTheWaterFlows

This package implements the D8 flow routing algorithm [1] as well as a basin-filling algorithm, also by [1]. This allows to calculate water pathways on a digital elevation model (DEM).

So far little efforts have been made to make this fast or memory efficient. The algorithm seems to be of order Q(n) where n is the number of grid points (provided the number of pits is constant). For a 1000x1000 map with 8 pits, it runs in 5s on my laptop from 2012.

Example of upslope area calculated in below example.

Manual

using WhereTheWaterFlows, PyPlot
const WWF = WhereTheWaterFlows

"Synthtic DEM with a few maxs and mins"
function peaks2(n=100, randfac=0.05)
    coords = range(-pi, pi, length=n)
    return coords, coords, sin.(coords) .* cos.(coords') .-
        0.7*(sin.(coords.+1) .* cos.(coords')).^8 .+
        randfac*randn(n,n)
end
x,y,dem = peaks2(200)
area, slen, dir, nout, nin, pits, c, bnds = waterflows(dem)

# log-upslope area as well as pits (sinks)
plotarea(x, y, area, pits)

# log-upslope area over contours of the dem
plotarea_dem(x, y, dem, area, pits)

# catchments
figure()
WWF.heatmap(x,y,c)

# stream length
figure()
WWF.heatmap(x,y,slen)

demf = fill_dem(dem, pits, dir)
# "lake-depth"
figure()
WWF.heatmap(x,y,demf.-dem)

References

[1] O’Callaghan, J. and Mark, D.: The extraction of drainage networks from digital elevation data, Comput. Vision Graph., 28, 323–344, 1984. (behind a pay-wall)