A fast and simple algorithm for calculating flow accumulation matrices from raster digital elevation models

s of the International Cartographic Association, 1, 2019. © Authors 2019. CC BY 4.0 License. 29th International Cartographic Conference (ICC 2019), 15–20 July 2019, Tokyo, Japan | https://doi.org/10.5194/ica-abs-1-434-2019


Introduction
The automatic extraction of drainage networks from raster digital elevation models (DEMs) is required in many scenarios such as soil erosion modeling, hydrological process simulation, and geomorphological analysis (Fu et al., 2011;Nobre et al., 2011;Yamazaki et al., 2012;Buchanan et al., 2014;Bai et al., 2015). A widely used method for extracting drainage networks from DEMs is based on the simulation of surface flow (Jenson and Domingue, 1988;Wang and Liu, 2006;Zhou et al., 2016). The method is composed of multiple steps which include removing depressions, assigning flow directions, and calculating the flow accumulation matrix. Among these procedures, calculating the flow accumulation matrix is an important step. The flow accumulation of a cell is equal to the number of cells that drain to it (O'Callaghan and Mark, 1984). Flow accumulation is an essential input for many hydrological and topographic analyses such as stream channel extraction, stream channel ordering, and sub-watershed delineation (Bai et al., 2015;Su et al., 2015;Barnes, 2017).
Some algorithms derive the flow accumulation matrix directly from a DEM (Arge et al., 2003;Bai et al., 2015). These algorithms generally require cells to be sorted based on their elevation values and have O(NlogN) time complexity. They start from the highest cells and gradually move to lower cells. The algorithms encounter problems in flat areas, where flow directions cannot be determined based solely on the values of the neighboring cells. These algorithms access cells in an order based on their elevations and can result in random scattered data swapping between memory and the hard drive when they are applied to massive DEMs that do not fit in the main memory (Su et al., 2015).
It is more common to derive the flow accumulation matrix from a flow direction matrix rather than directly from a DEM. There are two methods for flow direction determination from a DEM: the single-flow direction method and the multiple-flow direction method. In the single-flow direction method, each cell only drains to one neighboring cell. The D8 method, which uses the direction of steepest descent as the flow direction of a cell, is the most widely adopted single-flow direction method (O'Callaghan and Mark, 1984;Garbrecht and Martz, 1997;Nardi et al., 2008;Barnes et al., 2014). In the multiple-flow direction method, each cell can flow to more than one neighboring cell (Freeman, 1991;Quinn et al., 1991;Qin and Zhan, 2012). Because the D8 method is the most widely used method for determining flow direction, this study focuses on calculating the flow accumulation matrix from the flow direction matrix that is derived using the single-flow D8 method. The assignment of flow directions in depressions and flat areas in a DEM must be treated with special algorithms when the D8 method is used (Garbrecht and Martz, 1997;Wang and Liu, 2006;Barnes et al., 2014).
With the advent of airborne LiDAR (Light Detection and Ranging) technology, DEMs have become increasingly large. It is common for a DEM to contain billions of cells. The time required to calculate flow accumulation matrices of massive DEMs using conventional methods is becoming prohibitively long. In recent years, new algorithms have been proposed for calculating flow accumulation matrices. In this study, we propose a fast and simple algorithm to calculate the flow accumulation matrix. The remainder of the paper is organized as follows. Section 2 gives an 4 overview of the algorithms for flow accumulation calculation from single-flow flow direction matrices. Our proposed algorithm is presented in Section 3. Section 4 presents the experimental results and compares the proposed algorithm with existing algorithms. Section 5 concludes the paper.
2 Overview of algorithms for flow accumulation calculation In this section, we give an overview of the algorithms for calculating flow accumulation from single-flow flow direction matrices. All the algorithms produce the same flow accumulation matrix for the same input flow direction matrix. For convenience, Table 1 lists a group of symbols that are used in the pseudocode of the algorithms. Among the symbols, the symbol of NextCell(c) represents a function that returns a Boolean value and is used for tracing the immediate downstream cell of input cell c. If the input cell c drains towards the outside of the DEM or to a NODATA cell, the function returns a false value. Otherwise, the function returns a true value and cell c is updated to point to the downstream cell to which it drains.  The NIDP of a cell c is the number of neighboring cells that drain to c. Cells with an NIDP of zero are usually located on ridges, and cells with an NIDP greater than one are the intersection cells of more than one drainage path. The pseudocode for calculating the NIDP matrix from a flow direction matrix is shown in Algorithm 1 (Fig. 1). 2.3 BTI-based algorithm Su et al. (2015) propose an algorithm to use basin tree indices (BTI) to guide the calculation of the flow accumulation matrix. Their algorithm starts from the outlet cells that drain to the outside of the DEM and builds basin trees by tracing the drainage paths from the outlet cells to the source cells of each basin. The flow accumulation matrix is then calculated by tracing the trees from the leaves to the roots. The time complexity of the algorithm is O(N). The pseudocode of this algorithm is shown in Algorithm 3 (Fig. 3). To avoid the repeated reallocation of the arrays for 7 storing the trees, an array of size N is pre-allocated for storing all of the basin trees in our implementation of the algorithm in Section 4.

Recursive algorithm
This algorithm computes the accumulation value of a cell c by recursively computing the accumulation values of all of its neighboring cells that drain to it (Freeman, 1991;Choi, 2012).
This is the first algorithm for flow accumulation computation that has a time complexity of O(N).
The pseudocode of this algorithm is shown in Algorithm 4 (Fig. 4). Usually, the program that implements the recursive algorithm needs to determine the maximum size for the call stack during compilation and this size must be sufficiently large to process input data, which is a challenging issue to deal with when the size of the input data is unknown. In our algorithm, we define three types of cells within the flow direction matrix: source cells, interior cells, and intersection cells. A source cell does not have any neighboring cells that drain to it and its NIDP value is zero. An interior cell has only one neighboring cell that drains to it and its NIDP value is one. An intersection cell has more than one neighboring cell that drains to it and its NIDP value is greater than one.
The  A worked example of the proposed algorithm is given in Fig. 6. A synthetic DEM with a size of 3 rows by 4 columns, including its flow direction matrix, is shown in Fig. 6(a). The algorithm first computes the NIDP matrix ( Fig. 6(b)) and initializes the flow accumulation matrix with 1 ( Fig.   6(c)). The algorithm then traverses the NIDP matrix from top to bottom and from left to right. It 10 encounters the first source cell H. The algorithm starts the first round of tracing from H and all downstream cells of H including D, C, and F, are traced ( Fig. 6(d)). The accumulation value of each traced cell is increased by the accumulation of its immediate upstream cell. Because F is an intersection cell, the tracing stops and the NIDP value of F is decreased by 1. F is treated as an interior cell hereafter. The algorithm keeps traversing the remaining cells in the NIDP matrix and encounters the second source cell J. The algorithm starts the second tracing from J and all downstream cells including I, E, and A are traced (Fig. 6(e)). After the second tracing is done, the NIDP value of A is decreased by 1 and A is treated as an interior cell hereafter. In Fig. 6(f), the third tracing starts from cell L and all downstream cells including K, G, F, B, and A are traced.
Note that cell F is being treated as an interior cell and that the third tracing does not stop when F is encountered. After the third tracing is done, no source cells remain and the traversing process is complete. The final flow accumulation matrix is obtained and shown in Fig. 6(f). Unlike the BTI-based algorithm, however, the recursive algorithm does not build the basin trees explicitly. The running time of a recursive algorithm is highly dependent on the optimization applied to it by the compiler. Different compilers may generate machine codes with considerable differences in running times.

Memory requirement analysis
All algorithms require an input flow direction matrix and an output flow accumulation matrix.
Suppose that a DEM has N cells in it. Our algorithm requires an NIDP matrix. Because the NIDP value varies between 0 and 8, the NIDP matrix requires N bytes. Wang's algorithm also requires an NIDP matrix. In addition, it requires a queue Q to hold source cells. Each cell that is pushed into Q has its row index and column index information, which requires 8 bytes for storage in the most general cases, where the number of rows and columns may exceed 70,000 and a 4-byte unsigned integer is needed to store the row index and column index separately. For our test DEMs in Section 4, about 33.71% of the total number of the cells averaged across all the test DEMs are source cells. On average, the initial space required by the queue is more than 2.5N bytes. Jiang's algorithm does not use the NIDP matrix. Instead it uses a stack and requires at least 2.5N bytes as well. The BTI-based algorithm requires additional memory space to store the basin trees. Similar to Wang's algorithm, each cell needs 8 bytes for its row index and column index. The additional total memory space required by the BTI-based algorithm is at least 8N bytes. The recursive algorithm requires the call stack to track the calls. The memory requirement of a recursive algorithm depends on many factors, including the compiler and the platform on which the algorithm runs, and it is difficult to make a quantitative estimate of its memory requirement. All five algorithms produce the same flow accumulation matrix for each tested DEM. Each flow direction matrix is completely loaded into the main memory for processing and the loading time is excluded from the total running time. Each algorithm is run multiple times and the average running time is used for analysis. Figures 7 and 8 Su et al. (2015), who report that the BTI-based algorithm is much faster than Wang's algorithm, which is called the improved flow number matrix algorithm in their study. The difference in the running times may be caused partly by different details of the implementations.
For example, Wang et al. (2011) originally use two collection structures for storing the source cells in two consecutive iteration steps and copy elements from one collection structure to the other structure. For the BTI-based algorithm, we pre-allocate a matrix of the same size as the input 13 DEM to avoid the repeated reallocation of the dynamic arrays for storing the basin trees. Third, in terms of the absolute amount of running times, with the exception of the Wang's algorithm, other four algorithms generally run slower on the Linux system than on the Windows system. This may partially be attributed to the fact that the main frequency of the processor on which the Linux system runs is lower than that of the processor on which the Windows system runs.
It is clear that the running times of the algorithms for calculating flow accumulations are subject to such settings as the hardware configurations, operating systems, compilation optimization options, and implementation details. Because our proposed algorithm has a smaller constant coefficient before the time complexity of O(N) than other four algorithms, our proposed algorithm is expected to run the fastest as long as similar settings are adopted for all algorithms. time complexity, our algorithm runs substantially faster, requires less space than all non-recursive algorithms, does not require a collection structure, and is easy to understand and implement.
Our proposed algorithm is only applicable to single-flow flow direction matrices. In future work, we will adapt the algorithm to make it applicable to multiple-flow flow direction matrices.