A Review Of Upscaling Algorithms For Flow Direction Rasters

INTRODUCTION

The shape of the Earth's surface can be shown by a regular grid, which can be rectangular, triangular, hexagonal, or any other shape. Flow direction is a morphometric parameter that describes this shape. For a grid cell, flow direction is a value representing the direction (or multiple directions) water will flow out of the cell. There are different schemes (algorithm) of flow direction assignment; they can route flow from the cell to the one (D8 (O’Callaghan and Mark 1984)), or many (Dinf (Tarboton 1997), MFD (Freeman 1991; Quinn et al. 1991), or one of its derivatives (Qin et al. 2017)). Depending on the applied algorithm, flow direction could be presented as actual azimuth (e.g. for Dinf) or a coded value substituting one or many of the allowed directions. For example, the output of the D8 flow direction computation is presented as integer values in the domain {0; 2ⁿ, n Є [0, 7]}, where each code stands for a direction multiple of 45°.

Global Earth system models, macroscale hydrologic models, and many other applications require river networks at coarse resolutions. These models typically use gridded data as an input, so river networks are usually formalized as flow directions. Typical resolution for macroscale land surface models ranges from 0.5° to 1° (Decharme et al. 2010; Oki and Sud 1998; Pappenberger et al. 2010). The same resolution is applied in global hydrological models (Sood and Smakhtin 2015), although there are more recent examples of direct implementation of a finer resolution river network data (Lin et al. 2019).

The concept of flow directions was initially introduced in morphometric analysis of digital elevation models (DEMs) (Jenson and Domingue 1988; O’Callaghan and Mark 1984). According to the D8 algorithm, the flow direction of the cell is the direction of the steepest slope line between the central cell and the one of its 8 neighbors. This modeling principle only considers the movement of water down the slope due to gravity. If this principle is applied to a global scale and coarse grid resolutions, it can lead to a substantial error because the typical gradient between cells is relatively lower than in fine resolutions, so the ‘formally’ calculated slope between neighboring cells will not always reflect real flow behavior (Döll and Lehner 2002; Fekete et al. 2001; O’Donnell et al. 1999; Vörösmarty et al., 2000). However, relatively simple methods of aggregating more detailed (fine-resolution) data, such as calculating mean or modal value, are not suitable for flow directions because they are strictly related to the grid geometry (Decharme et al. 2010). Therefore, several upscaling methods were developed to produce flow direction coverage on coarse-resolution grids.

It should be noted that different research fields have opposite interpretations of the terms ‘upscaling’ and ‘downscaling’. In Earth sciences, downscaling is usually understood as increasing spatial resolution, i.e. decreasing cell size, while upscaling means decreasing spatial resolution and increasing cell size (Eilander et al. 2021; Peng et al. 2017; Wu et al. 2011). At the same time, in the practice of digital image processing, upscaling is a synonym for increasing image resolution (Panda and Meher 2024). In this paper, we follow the tradition established in the Earth sciences and understand upscaling as the production of coarse-resolution grids based on more detailed data.

The literature presents many approaches to creating low spatial resolution flow direction grids. These approaches could be divided into the following groups:

1) manual flow direction assignment based on map image analysis (e.g. Miller et al. 1994);

2) resampling high-resolution flow direction grids and/or their derivatives to a coarser resolution, often with additional coverages (e.g. Fekete et al. 2001; O’Donnell et al. 1999; Reed 2003);

3) tracing vector stream network to derive raster flow direction coverage (Mayorga et al. 2005; Olivera and Raina 2003).

Manual data processing is very labor-intensive, and the use of methods based on the vector representation of the hydrographic network is limited by the availability of initial global coverage data. Therefore, in this paper, we will focus on methods for the second and third groups.

FLOW DIRECTION UPSCALING ALGORITHMS

Grid-based tracing

All algorithms for generalizing flow directions by tracing flow on fine-resolution grids are based on similar principles. Firstly, it is implicitly assumed that the geometry of the target grid is aligned to the geometry of the source grid. The target grid’s cell size is an integer multiple of the source grid’s cell size, and source and target grids’ reference points, which are usually in the lower and left corners, are either the same or are moved by an integer number of the source grid cell size. In most cases, both grids are constructed in a geographic coordinate system (latitude-longitude), although there are exceptions to this rule. Then, based on the flow direction raster, we obtain a catchment area (contributing area, upstream area) raster before performing the actual upscaling procedure. Finally, the procedure is usually based on overlaying a source (fine-resolution) grid and a target (coarse-resolution) grid; for each target grid cell, an ‘outlet’ source grid cell is determined. This ‘outlet’ cell has the highest catchment area value among other fine cells within the coarse cell; it is also often (but not always) adjacent to the edge of the coarse cell. The outlet is later used as a starting point to trace flow downstream along the initial (fine-resolution) drainage graph. Tracing is performed until a certain condition is met. The resulting flow direction for the current target cell is determined from the position of the cell where tracing was stopped.

It should also be noted that a certain terminology has developed in the works devoted to this problem. In particular, the term “cell” usually refers to a coarse-resolution grid cell, while fine-resolution grid cells are denoted as “pixels”. Hereinafter we, while possible, do not follow this terminology and always use ‘coarse cell’ and ‘fine cell’, respectively.

One of the first implementations of the principle described above was presented in (Döll and Lehner 2002). First, it finds outlet cells for each target cell. Then, it simply assigns flow direction for the target cell according to the direction of the outlet cell. This procedure itself leads to an uneven distribution of resulting flow directions: cardinal directions (north, east, south, and west) prevail over diagonal directions (northeast, southeast, southwest, and northwest). In the work of Döll and Lehner (2002), this procedure was used as a first step in an iterative procedure, which included manual review and correction. It is also stated that the ‘automatic’ part of this method is used in the mRM routing model (Thober et al. 2019).

A more complicated technique was presented earlier by O’Donnell et al. (1999). They expand the concept of corner cells, including not only a fine-resolution grid cell located exactly at the corner of a coarse-resolution target cell but also its fine-resolution neighbors. The exact number of neighbors is a parameter of the algorithm. For non-corner cells, the procedure is the same; for corner cells, the flow is traced within a 3×3 vicinity of the central target cell, and the water is routed to the coarse cell when the tracing is stopped.

The COTAT method (Cell Outlet Tracing with an Area Threshold) applies tracing of the source flow directions from all outlet cells (Reed 2003). It also utilizes different stopping conditions. For each coarse grid cell, the procedure starts with its fine-resolution outlet; the catchment area value of the outlet is stored for further processing. Next, flow is traced from each outlet cell along the fine-resolution grid. During tracing, the following conditions are checked: 1) catchment area difference between the current source grid cell and the outlet cell should be lower than the user-defined threshold; 2) the current source grid cell should be within the 3×3 vicinity of the central target cell. If one of the conditions is violated, tracing stops, and the coarse cell where the tracing ended is identified. The identified coarse cell receives the central coarse grid cell direction. If the outlet cell has the highest catchment area value among all its fine-resolution neighbors, it practically means that this cell represents a local sink, and flow direction for the target cell is set to zero, which stands for no outward flow. The basic principle of the COTAT procedure is shown in Fig. 1.

Upon completion of the routing, COTAT implements a post-processing procedure to prevent instances of crossing paths. When the flow directions of adjacent coarse cells intersect, a special reclassification table reassigns one of them to the non-crossing pattern.

An upgraded version of COTAT was presented in (Paz et al. 2006); this method is referred to as COTAT+ (Davies and Bell 2009) or the effective area method (EAM) (Eilander et al. 2021). This modification has three differences from the original COTAT procedure. First, it introduces one more limiting parameter: minimum distance along the flow path. If the threshold catchment area difference is exceeded, but the distance of the constructed path is less than the assigned minimum, the tracing continues until the minimum is reached. The second modification affects cases when no valid outlet fine-resolution cell could be defined; in this case, the last visited cell is assigned as the draining cell. The third modification is introduced for crossing flow directions: instead of a reclassification table, a more sophisticated decision rule is used. According to the authors (Paz et al. 2006), these modifications improve flow direction assignment over flat areas with large parallel flows, e.g., in the Amazon River basin.

The next method is FLOW, proposed by Yamazaki et al. (2009). In this method, flow is traced down from the outlet cell until the next outlet cell (of the neighbor target cell) is reached, and the length of the flow path is measured. If the length is lower than the user-defined threshold, the second outlet cell is rejected as an outlet cell. If so, the outlet cell for the neighbor target cell is reassigned: a new outlet is selected among the cells allocated on the border of a (neighbor) target cell, excluding the one previously rejected. The procedure repeats until the flow path length becomes longer than the threshold value for each target grid cell. After that, flow directions of a coarse-resolution grid are assigned following the constructed paths. A distinctive feature of the method is that the resulting flow directions do not follow the D8 pattern. Instead, flow from the "central" target cell may be directed to a cell that does not belong to its immediate neighborhood. According to Yamazaki et al. (2009), this allows for more realistic flow patterns and better catchment area matching between source and target datasets.

The method developed by Lucas-Picher et al. (2003) for the Canadian Regional Climate Model (CRCM) is significantly different from the above methods because it does not require grid alignment and is specifically designed for the situation where the coordinate systems of the source and target datasets are not coincident. The method consists of the following steps: 1) derive spatial extent of the target grid cell in the domain of the source grid; 2) select all the source grid cells which intersect the computed extent; 3) from the cells selected in step 2, select cells whose flow direction is oriented outside the extent; 4) from the cells selected in step 3, find a cell with the highest catchment area value; 5) identify the coarse grid cell that covers the fine cell selected in step 4; 6) assign flow direction from the central coarse cell to the coarse cell identified in step 5. The process is shown in Fig. 2.

The Dominant River Tracing (DRT) method was proposed by Wu et al. (2011). It also relies on tracing the flow over a regular grid, but this is done in a specially defined hierarchical order. DRT applies the concept of the longest effective dominant river (LEDR) segment–a river segment that dominates the local drainage of the cell. For a coarse cell, the LEDR is identified as the river segment that has a relatively large (but not necessarily the largest) catchment area value compared to other river segments in this cell and is longer than a minimum length threshold set by the user within the cell. The method calculates flow directions consequently, analyzing one basin at a time. The basins are sorted and ordered according to their respective catchment areas.

DRT starts by finding the dominant basin and river of the study area and assigns flow directions for cells along dominant rivers beginning from headwater cells to basin, also called sink cells. Subdominant rivers and tributary flow paths are then identified and ordered according to their respective catchment areas. The priority for assigning flow directions is assigned to successively higher-order rivers until all cells in the most dominant basin have assigned flow directions. The DRT then selects progressively smaller, less dominant basins and rivers and assigns flow directions in a similar manner until all cells in the given study area have been assigned flow directions.

The assignment of the flow direction itself for a cell in the DRT is performed as follows. First, the cell is divided into eight π/4 sectors radiating from the cell center. Then, an LEDR from the cell is traced downstream for a specified length (0.6 of cell size for cardinal sectors and 0.8 of cell size for diagonal sectors). The cell where the tracing stopped receives the flow direction from the target cell. DRT also applies two supplementary procedures to preserve subdominant rivers and sinuous flow paths.

DRT was applied to create a set of global coverages of flow directions and related parameters at 2°, 1°, 30' (1/2°), 15'(1/4°), 7.5' (1/8°), 3.75' (1/16°) resolutions (Wu et al. 2012)¹. HYDRO1K² and HydroSHEDS (Lehner et al. 2008) were used as input data for these coverages..

Another upscaling procedure based on flow tracing was presented in (Eilander et al. 2021). The method, called Iterative Hydrography Upscaling (IHU), consists of four consecutive stages (iterations). The first stage involves computing initial flow directions for coarse grid cells. The following stages are needed to correct initial directions to fix errors (stage 2), optimize in-between outlet distances (stage 3), and minimize the error caused by erroneous flow directions (stage 4). At the first stage, like all methods listed above, IHU is based on identifying outlet fine-resolution cells for each coarse-resolution cell and tracing flow downstream; but unlike the above approaches, here the output cell is not located at the edge of the target cell, but rather in the inner region. The inner region is defined as a rounded rectangle excluding edge cells and some corner cells adjacent to them. The flow is traced from the outlet cell down to the next outlet cell. The second step is repairing erroneous flow directions. The authors recognize direction as erroneous if the first outlet cell downstream (from fine-resolution grid) is not coincident with the coarse-resolution downstream cell. An iterative procedure identifies and assigns potential alternatives at this step. Not all erroneous situations could be repaired; if there are any, they are further addressed in stage 4. The third stage aims to optimize the distance between outlet pixels, measured along the fine-resolution flow directions. If this distance is short, one of the outlet pixels can potentially be removed in favor of another river segment within the same cell. The fourth stage addresses erroneous directions that were not corrected at stage 2; for these cells, tracing is performed again but now flow path length is considered. The IHU procedure was implemented in the open-source PyFlowDir package³. The authors also prepared a generalized datasets based on MERIT Hydro data (Yamazaki et al. 2019) with the resolution of 30'', 5' and 15'; these datasets are called MERIT Hydro IHU⁴.

Deriving flow directions from catchment area grids

Two methods allow producing coarse-resolution flow direction grids directly from fine-resolution catchment area grids without flow tracing. These methods are NSA (Network Scaling Algorithm) and DMM (Double Maximum Method). Both methods have the same requirements for source and target grids as flow tracing methods: the geometry of grids should be aligned, and the cell size of the target grid should be an integer multiple of the cell size of the source grid.

The NSA proposed by Fekete et al. (2001) works as follows. First, it aggregates (summarizes) catchment area values from a source grid to the target grid. Then, the flow direction for each target grid cell is assigned to the neighbor cell with the largest aggregated catchment area value. The procedure is shown in Fig. 3.

The DMM procedure proposed by Olivera et al. (2002) is more complicated. To implement it, two grids with coarse spatial resolution are constructed (primary and auxiliary). The auxiliary grid’s origin is shifted relative to the origin of the primary grid by a half of the coarse grid cell size. First, an outlet (fine) cell is identified for a cell of the primary grid. Second, an auxiliary coarse cell is identified for the outlet cell found in step 1. Third, an outlet cell is found for the auxiliary grid cell identified in step 2. Fourth, a primary grid cell is identified for the outlet found in step 3. This primary grid cell is a target cell for the ‘starting’ primary grid cell, which was processed in step 1. The final flow direction is assigned from the ‘starting’ cell to the ‘target’ cell. An overview of the process is shown in Fig. 4.

Tracing flow along vector hydrographic network

Another two methods were proposed to derive flow directions from a vector representation of the hydrographic network. Network Tracing Method (NTM) performs an overlay between vector streamlines and the target grid and traces the flow along vector lines (Olivera and Raina 2003). The procedure begins with a preprocessing step: a geometric intersection of the input river network lines graph with the boundaries of the target grid cells is performed. As a result, the vector network is divided into separate segments; each segment lies entirely within one cell. Then, endpoints of the lines lying on the cell boundaries are identified. After that, for each cell, its output point is determined as one of the previously identified endpoints of the lines through which rivers leave the cell, which has the greatest flow path length. The flow path length of a point is determined as the longest possible path upwards along the graph of the hydrographic network. The algorithm then traces network down along the graph, starting from the output point. The conditions for stopping tracing are reaching a threshold of length increase, crossing more than two coarse cells, and falling beyond the 3×3 neighborhood of the cell. Comparing the increment of the upstream flow length with the threshold value when choosing the flow direction allows skipping cells during tracing, through which the flow passes too short a path. This parameter is also used to balance orthogonal and diagonal flow directions. The higher this threshold, the more often the algorithm will choose diagonal directions over orthogonal ones.

An important drawback of NTM is high requirements to the input data. To obtain a correct result, the initial hydrographic network graph must be planar, have a tree structure, and have all its edges oriented strictly downstream (Olivera and Raina 2003). To satisfy the planarity condition, all rivers must be divided into separate segments at confluence points. The tree structure implies that bifurcation is forbidden. The orientation of the vector lines downstream is important since it determines the scanning direction along the graph. In addition, the density of the initial network is also important for the NTM algorithm: it must be dense enough for modeling at a desired resolution. To satisfy this condition, at least one network segment should intersect each cell boundary.

Another method was proposed in (Mayorga et al. 2005). This method is limited to single basins and requires an outlet point for the processing area. It first converts the vector stream network into a flow direction grid and then applies an iterative procedure to correct the resulting flow directions. The method imposes very strict requirements on the input vector data: the entire network must be interconnected; network segments should form a tree-like structure (channels that form loops or are split into parallel braided channels must be simplified to single channels following single, well-defined directions); no polygons are allowed; and all segments should be presented as polylines.

INPUT DATA FOR UPSCALING

The methods discussed above use different types of data as input datasets. Grid-oriented upscaling procedures require catchment area coverage, or a catchment area with a flow direction grid. There are three global datasets that provide these coverages: a relatively old HYDRO1K, and more recent HydroSHEDS (Lehner et al. 2008; Lehner and Grill 2013) and MERIT Hydro (Yamazaki et al. 2019). HydroSHEDS grids were utilized in the work of (Wu et al. 2011, Wu et al. 2012). MERIT Hydro coverage is implied in (Eilander et al. 2021).

The issue of input data for vector-based approaches appears to be more complex. The work of Mayorga et al. (2005) relied on Digital Chart of the World (DCW) data (Danko 1991), but the authors emphasize that this data requires careful and time-consuming preprocessing. Another vector-based approach (Olivera and Raina 2003) used a stream network derived from the HYDRO1K dataset. HydroSHEDS and MERIT Hydro also include vector river networks derived from flow direction rasters. However, if raster coverages are the main sources for these datasets, it begs the question of whether vector-based methods are necessary or whether grid-based methods are better.

EVALUATION OF UPSCALING RESULTS

All the above methods produce reliable flow direction distributions at low spatial resolution and global coverage, but the computational results still differ in detail. Fig. 5 shows upscaling results according to three different approaches (COTAT, DRT, NTM) for a limited area in the lower Mississippi River basin. The overall distribution of flow directions is similar between all three results, but the details are noticeably different: for example, the main course of the Mississippi in the result obtained by the DRЕ algorithm is oriented strictly from north to south, while other algorithms show some sinuosity for the same river.

Since the upscaling results vary, quantitative criteria for their evaluation are needed. Several criteria were proposed to evaluate flow direction distribution, including:

• visual comparison of the generated network with reference data (Döll and Lehner 2002; O’Donnell et al. 1999);
• estimation of basin and sub-basin areas for the selected points (Döll and Lehner 2002; Fekete et al. 2001; O’Donnell et al. 1999; Olivera et al. 2002; Vörösmarty et al. 2000; Wu et al. 2011);
• comparison of the overall river segment distribution and statistics, including length and number of segments (Vörösmarty et al. 2000) ;
• calculation and evaluation of flow direction distribution statistics itself (Olivera et al. 2002; Reed 2003)

The comparison is usually made relative to some reference distribution, which is based on higher resolution flow directions. Basin and sub-basin areas are estimated directly through overlay with reference polygons or indirectly via comparison of catchment areas. In the latter case, a linear regression between the upscaled and reference results is calculated, and typical statistics (R², RMSE, etc.) are computed. The process of comparing river segments follows a similar way. For flow direction distribution, Olivera et al. (2002) proposed a following criterion, called ‘side-to corner ratio’: 59% of the cells should have orthogonal flow direction (north, east, south, west), and 41% of the cells should have diagonal flow direction. This criterion is based on a theoretical assumption and has not been tested for fine-resolution datasets.

A small number of papers compare different upscaling methods with each other. Davies and Bell (2009) compared COTAT, COTAT+ (EAM) and NTM methods, suggesting that COTAT+ shows the best results among tested approaches. Wu et al. (2011) evaluated their DRT algorithm against the approach of Döll and Lehner (2002). The results of the IHU application were compared with COTAT+ (EAM) and DMM (Eilander et al. 2021). Cohesive testing of all modern upscaling methods in global coverage has not been conducted to date.

DISCUSSION

There are many algorithms that have been developed for upscaling flow directions and/or creating coarse-resolution flow direction coverage. Most of them use procedures based on flow tracing over a grid of lower spatial resolution. The two most novel algorithms of this kind (Dominant River Tracing, DRT, and Iterative Hydrography Upscaling, IHU) also belong to this group. Alternative approaches include the derivation of flow directions from fine-resolution catchment area coverages (e.g. Double Maximum Method, DMM) and the calculation of flow directions from vector hydrography network data (e.g. Network Tracing Method, NTM). The authors of DRT and IHU also made sets of global coarse-resolution flow direction coverages; typical cell size within these datasets varies from 30 arc seconds to 1 degree. These datasets have been shown to have acceptable accuracy for application in modeling hydrological and climatic processes in global coverage.

Higher-resolution flow direction and catchment area coverage, such as HydroSHEDS and MERIT Hydro, are most often used as input data for generating low-resolution flow directions. In some cases, vector representations of the hydrographic network based on raster coverage are also used. Alternative sources of vector representations of the hydrographic network, such as Digital Chart of the World, are used less frequently because existing algorithms have high input data requirements that these sources do not possess. In particular, river segments should form a continuous network with a tree-like structure; branching of the channel is not allowed. Data sets created by digitizing geographic maps do not have these properties.

There are several ways to evaluate the results of upscaling. They include visual assessment of the result, calculation of catchment area statistics, estimation of length distributions and number of river network segments, and calculation of direct flow direction statistics. While all the above methods allow characterization of the distribution as a whole (and are used accordingly), their ability to detect local differences on coarse-resolution grids still needs to be assessed. This is especially the case for the criterion proposed in (Olivera et al. 2002) to evaluate the distribution of flow directions (59/41 ratio). It is necessary to find out whether such a ratio is observed in the raw data used for upscaling. It also seems important to automate the estimation of river network shape generated from coarse-resolution flow directions; to date, no effective solution to this problem has been proposed.

CONCLUSIONS

The paper considers contemporary methods for generating flow direction grids (rasters) at coarse spatial resolution (about 1°). Most of these methods can be seen as specific algorithms of generalization, or in this case, the upscaling of flow direction grids with relatively fine resolution. These methods are based on flow tracing in fine resolution and further derivation of directions at coarse resolutions. The most advanced methods of this type are Dominant River Tracing (DRT) and Iterative Hydrography Upscaling (IHU). DRT output data (resolutions 2°, 1°, 1/2°, 1/4°, 1/8°, 1/16°). These coverages are generated from fine-resolution datasets: HYDRO1K, HydroSHEDS, and MERIT Hydro, respectively. Alternative approaches include the generation of the desired datasets directly from catchment area grids (without flow tracing) and the aggregation of vector stream networks. There are no publicly available ready-made datasets based on these methods, making them less frequently used.

Assessing the validity of the resulting distributions is still an underdeveloped issue. Although several estimation methods have been proposed, all of them mostly reveal global characteristics of the obtained distributions, while local features of flow direction distribution remain unclear. It is necessary to develop estimation methods that would allow efficient (and automatic) analysis of river network configurations derived from coarse-resolution flow directions and to apply these approaches to generation results based on all modern upscaling methods.