Algorithms – such as D8 (O’Callaghan & Mark 1984) or MD8 (Quinn P. et al. 2006) – indicate flow directions using an elevation raster. Flow directions are calculated by evaluating elevation of neighbouring cells. However, they don’t seem to represent flat surfaces with an accumulation of cells with similar elevation – such as lakes – appropriately.
This post proposes a simple algorithm to represent the flow within a lake by implementing a gravity flow towards the lake outflow. In contrary to the algorithm by Garbrecht & Martz (1997) (improved by Barnes et al. 2014), this solution does not increment elevations within the flat surfaces or lakes. Rather it calculates a “downnode” (i.e. the neighbouring raster cell where the water flows) for each cell with gravity towards the outflow.
- It is assumed that all the nodes (raster cells) and the outflow of a lake are known
- It is assumed that the outflow does not flow back into the lake
- It is assumed that lakes / flat surfaces don’t have multiple outflows
Description of the Algorithm
The algorithm starts at the outflow of the lake (Figure 1, a). It then searches for neighbours that are within the lake and sets their downnode to the outflow. From there the nearest neighbour to the outflow is chosen and again its neighbours are identified (Figure 1, b). Then again, the nearest node to the outflow is chosen and its neighbours are identified. The algorithm continues until every cells downnode is reset.
Figure 1: Explanation of lake flow algorithm. The rastercell with the star in (a) represents the lake outflow. Four different states are marked; Not touched (bright blue), Downnode reset (darker blue), Checknode (yellow star, current node to be checked), Checked (red cross, cells that have already been checknodes).
checknode = outflow tocheck = [outflow] while not every downnode is set: neighbours=getNeighbours(checknode) for neighbour in neighbours if neighbour in lake and downnode not set: neighbour.setDownnode(lake) tocheck.add(neighbour) tocheck.remove(checknode) #get nearest node from outflow for node in tocheck: distance = outflow.distance(node) if distance<distanceNearest: distanceNearest=distance nodeNearest=node checknode = nodeNearest
outflow – a raster cell representing the outflow of a lake
getNeigbours – gets the 8 neighbours of a cell/node
lake – represents all the nodes/cells within a lake
outflow.distance(node) – returns the Euclidian distance between the outflow and the node, an improved version of the algorithm would calculate the shortest path distance to the outflow instead ot the Euclidian distance
Barnes, R., Lehman, C. & Mulla, D. 2014. An efficient assignment of drainage direction over flat surfaces in raster digital elevation models. Computers & Geosciences, 62, 128–135, 10.1016/j.cageo.2013.01.009.
Garbrecht, J. & Martz, L.W. 1997. The assignment of drainage direction over flat surfaces in raster digital elevation models. Journal of Hydrology, 193, 204–213, 10.1016/S0022-1694(96)03138-1.
O’Callaghan, J.F. & Mark, D.M. 1984. The extraction of drainage networks from digital elevation data. Computer vision, graphics, and image processing, 28, 323–344.
Quinn P., Beven K., Chevallier P. & Planchon O. 2006. The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models. Hydrological Processes, 5, 59–79, 10.1002/hyp.3360050106.
Thanks to Manuel Bär for the inspiration!