The isotropic organization of DEM structure and extraction of valley lines using hexagonal grid

The automatic extraction of valley lines (VLs) from digital elevation models (DEMs) has had a long history in the GIS and hydrology fields. The quality of the extracted results relies on the geometrical shape, spatial tessellation, and placement of the grids in the DEM structure. The traditional DEM structure consists of square grids with an eight‐neighborhood relationship, where there is an inconsistent distance measurement between orthogonal neighborhoods and diagonal neighborhoods. The directional difference results in the extracted VLs by the D8 algorithm not guaranteeing isotropy characteristics. Alternatively, hexagonal grids have been proved to be advantageous over square grids due to their consistent connectivity, isotropy of local neighborhoods, higher symmetry, increased compactness, and more. Considering the merits above, this study develops an approach to VL extraction from DEMs based on hexagonal grids. First, the pre‐process phase contains the depression filling, flow direction calculation, and flow accumulation calculation based on the six‐neighborhood relationship. Then, the flow arcs are connected, followed by estimating the flow direction. Finally, the connected paths are organized into a tree structure. To explore the effectiveness of hexagonal grids, comparative experiments are implemented against traditional DEMs with square grids using three sample regions. By analyzing the results between these two grid structures via visual and quantitative comparison, we conclude that the hexagonal grid structure has an outstanding ability in maintaining the location accuracy and bending characteristics of extracted valley networks. That is to say, the DEM‐derived VLs based on hexagonal grids have better spatial agreement with mapped river systems and lower shape diversion under the same resolution representation. Therefore, the DEMs with hexagonal grids can extract finer valley networks with the same data volume relative to traditional DEM.


| INTRODUC TI ON
The automatic extraction of valley lines (VLs) derived from digital elevation models (DEMs) with well-established digital terrain analysis (DTA) technology plays an important role in digital watersheds, hydrological models, soil erosion, and other environmental analyses (Dietrich, Wilson, Montgomery, & McKean, 1993;Gopinath, Swetha, & Ashitha, 2014). The DEMs used in this process can be categorized into three forms: triangular irregular network (TIN;de Azeredo Freitas, da Costa Freitas, Rosim, & de Freitas Oliveira, 2016), contour line model (Ai, 2007;Ai & Li, 2010), and regular grid model (Costa-Cabral & Burges, 1994;O'Callaghan & Mark, 1984). Among these forms, the regular grid model has been the most widely utilized DEM structure due to its simple data structure, ease, and efficiency of computation (Walker & Willgoose, 1999). This model is a mathematical representation and simulation of the ground surface via an ordered array of elevation values from the discrete sampling points. In addition to the precision of elevation, VLs deduced from the regular grid-based DEMs rely on the spatial relationship among neighborhood grids (Liu, 2002), which is determined by the DEM structure, including grid size, grid direction, as well as grid geometry.
The flow-routing algorithm plays a critical role in delineating VLs and is employed to simulate the overland flow over landscapes (Martz & Garbrecht, 1992). This algorithm assigns the transmission and fluid processes of surface materials (e.g., water, sediments, and nutrients) from higher to lower units. The results predicted by this method are represented as fully connected and unidirectional networks (Tribe, 1992), which depend not on the isolated grid but on the elevation difference between adjacent neighborhoods according to their spatial relationship (Zhou & Liu, 2002). Thus, the VLs derived from gridded DEMs can be influenced by the spatial tessellation and placement of grids in the gridded structure, resulting simultaneously from the grid size and geometry.
However, the majority of studies have focused only on the grid size based on a square DEM structure with an eight-neighborhood relationship (Ariza-Villaverde, Jiménez-Hornero, & De Ravé, 2015;Davies & Bell, 2009;Woodrow, Lindsay, & Berg, 2016;Yang et al., 2014). There are two kinds of neighboring grids in a square gridbased structure: four orthogonal neighborhoods (or "von Neumann neighborhoods") with edge connectivity and four diagonal neighborhoods (or "Moore neighborhoods") with corner connectivity. The unfixed connectivity results in a different distance measure; that is, the distance between adjacent grids in the diagonal direction is √ 2 times that in the orthogonal direction. Tang, Zhao, and Cao (2000) proved that diagonal neighborhoods have lower accuracy in terrain representation due to the larger distance measurement. To solve distance imbalance, for a given grid, the D8 algorithm (O'Callaghan & Mark, 1984) utilizes the inverse distance weighted (IDW) method to calculate the slope relative to each neighborhood, represented by the ratio of the elevation difference to distance. However, this causes insufficient isotropy of the terrain to change the degree between orthogonal and grid structures via visual and quantitative comparison, we conclude that the hexagonal grid structure has an outstanding ability in maintaining the location accuracy and bending characteristics of extracted valley networks. That is to say, the DEM-derived VLs based on hexagonal grids have better spatial agreement with mapped river systems and lower shape diversion under the same resolution representation.
Therefore, the DEMs with hexagonal grids can extract finer valley networks with the same data volume relative to traditional DEM.
diagonal directions simultaneously; namely, the calculated results express the terrain variation in a unit of length for orthogonal neighborhoods, while they express the terrain variation in √ 2 units of length for diagonal neighborhoods. This anisotropy is more apparent in a DEM with a larger grid size, resulting in the inaccurate representation of the ground surface in diagonal directions, further affecting the results of VL extraction. In view of the above question, it is necessary to find another regular grid instead of square grids to construct a DEM tessellation structure and to explore a compatible method for VL extraction.
With respect to the DEM structure, a regular grid should satisfy the characteristics of single shape, seamless splicing, standard coding, defined hierarchy, and ease of interconverting with geographic space. In addition to squares, a large number of scholars (Carr, Olsen, & White, 1992;Peuquet, 2002;Raposo, 2013) have suggested that the equilateral triangle and hexagon are also utilized as basic units to establish regular tessellation models for terrain surfaces without any gap or overlap, and they have examined the relative advantages or disadvantages for corresponding tessellation models. Equilateral triangles have two different orientations, 12 neighborhoods with two types of connectivity and three distance measurements (the ratio is 1 ), increasing the difficulty of the subsequent calculations and analysis.
In image processing, cartography, ecology, computer games, and other fields, equilateral hexagons have been proved to be advantageous over square grids due to their consistent connectivity, compact structure, isotropy of local neighborhoods, and visual advantages (Carr et al., 1992;Condat, Van De Ville, & Blu, 2005;Middleton & Sivaswamy, 2006;Raposo, 2013;Singh & Oberoi, 2015;White & Kiester, 2008). However, there are many mature and widespread applications based on square grids in the aforementioned fields (Dale, 1998;Gonzalez & Woods, 2008;Shen, Ai, & He, 2018;Shen, Ai, Wang, & Zhou, 2018). In the hexagonal tessellation model, one cell is connected to six neighbor cells via edges only, which simplifies the neighborhood relationship. The distance between the centroids of adjacent cells is the same, meaning that the hexagon also has a consistent distance measurement (Luczak & Rosenfeld, 1976). Therefore, the hexagonal grid is more suitable for neighborhood processes, especially for issues related to nearest neighborhood, movement, or connectivity (Birch, Oom, & Beecham, 2007;Wang & Kwan, 2018). In comparison with the square, the hexagon is closer to the circle in shape and has a shorter perimeter under equal area, which makes its structure more compact and reduces edge effects (Hales, 2001). Theoretically, this property could increase the representation accuracy for the spatial features of a land surface, and it can realize the same degree of reliability with fewer sampling mesh elements (Lin & Peng, 2001). In particular, the line elements in hexagonal tessellation are less distracting due to their increased isotropic quality and greater angular resolution of hexagons (Kamgar-Parsi & Sander, 1989;Poorthuis & Zook, 2015). It is apparent that the hexagon has six-fold symmetry, whereas the square only has four-fold symmetry, which is one reason explaining its isotropy characteristics. Another reason is that the ratio of grid distance to Euclidean distance has less sensitivity to direction in hexagonal tessellation because of its consistent distance measurement (Anneville, Cury, Le Page, & Treuil, 1998). In terms of visualization, the angular resolution of the hexagon is 60°, showing as the non-orthogonal arrangement between neighborhoods in the tessellation model. Hexagons only conform to the sensitivity to oblique angles by photoreceptors in the human eye (He & Jia, 2005). As a result, hexagons have attracted increasing attention in spatial modeling, and they have been efficiently utilized in various applications related to the fields of image processing, cartography, ecology, computer games, and so on.
Applying the above superiorities of hexagons to the flow direction analysis, relevant studies (de Sousa, Nery, Sousa, & Matos, 2006;Hyväluoma, 2017) believe that the grid-based DEMs of hexagons can reduce the grid orientation dependence of flow routing to a certain extent and that these models have the outstanding ability to maintain the flow direction from the original resolution at a coarser resolution. However, there is little application in DEM analysis because of the complicated composition and decomposition of hexagons. Therefore, whether hexagons have dominance in VL extraction has not been examined.
This article proposes an innovative approach to the automatic extraction of VLs from DEMs based on a hexagonal tessellation model. The remainder of the article is arranged as follows. In Section 2, we introduce the grid structure and the approach for generating a grid-based DEM with hexagons. Then, based on the hexagonal grid structure, the implementation of two crucial steps (including depression filling and flow direction calculation) in VL extraction is discussed in Section 3. Furthermore, the experiments, with three sample areas in the Yangtze River Basin, are presented in Section 4 to compare the DEM-derived VLs from hexagonal and square grids. We analyzed the advantages and disadvantages of the hexagonal grid in VL extraction from DEMs by means of visual and quantitative comparison. Finally, Section 5 offers some conclusions.

| Hexagonal grid structure
The grid structure is the digital grid space that the DEM uses to simulate the real geographic space. To satisfy the demand of DEM analysis, two notable points should be taken into account: 1. A suitable grid coordinate system for the computer to store, process, and transform geographic coordinates.
2. The establishment of a relationship between neighbors for DEM analysis.
Focusing on the spatial configuration of hexagons, there are five coordinate systems proposed for the hexagonal grid structure, including offset coordinates, axial coordinates, symmetric coordinates, doubled coordinates, and Gosper curve coordinates (Birch et al., 2007;Her, 1993Her, , 1995Lukić & Nagy, 2019;Nagy, 2003;Xin, Ai, & He, 2017). The offset coordinate system is the simplest way to label each hexagon with the number of rows and columns; therefore, it is the most similar to the orthogonal coordinate system of squares, as displayed in Figure 1a.
Moreover, it is compatible with the Cartesian coordinate system, making it easier to convert to and from geographic coordinates. Hence, it is the preferred choice for discretization of geographical entities with a hexagonal grid structure. However, the arrangement between two adjacent rows or columns has a "zigzag" misplacement.
As shown in Figure 1b, in the offset coordinate system, each hexagon is encoded with (i, j), where i and j are the number of columns and rows, respectively, and the centroid of the hexagon (0, 0) corresponds to the origin in the geographic coordinate system.

F I G U R E 1
The coordinate systems of squares and hexagons: (a) orthogonal coordinate system of squares; and (b) offset coordinate system of hexagons According to the "point raster" view (Tang, Liu, & Lü, 2005), the attribute of each cell centroid represents the replacement for the whole cell. Therefore, the calculation and analysis between hexagons is actually replaced by their centroids. With this in mind, the transformation between gridded coordinate systems and geographic coordinate systems is a basic step in representing the real planar space. In hexagonal structures, the ratio of grid spacing between a row and a column is 1: √ 3 2 , so this transformation can be defined as: where X (i,j) and Y (i,j) are geographic coordinates of the centroid of the hexagon (i, j); x 0 and y 0 are the geographic coordinates of the origin; and w is the horizontal width of the hexagon.
Six-connectivity is the basic neighboring relationship of the hexagonal grid structure applied for operation in DEM analysis. The neighbors of each hexagon include the hexagons to its left, upper left, upper right, right, lower right, and lower left. However, even-row hexagons have different patterns for neighboring grid coordinates from odd-row hexagons due to misplacement among columns, as shown by the hexagon (5, 2) with its neighborhoods identified in orange and the hexagon (3, 5) with its neighborhoods identified in green in Figure 1b. The distinctions are reflected in the neighborhoods to the upper left, upper right, lower right, and lower left, which are shown in red in Figure 2b. For even-row hexagons (i, j), the column numbers of these neighborhoods are i, i + 1, i, and i + 1, while they are i − 1, i, i -1, and i for odd-row hexagons (i, j). The relationship between neighborhoods can be established via this pattern.

| Generation of hexagon grid-based DEM
The existing data sources used to generate hexagonal grid-based DEMs include discrete sample points, contour line data (Wang & Ai, 2019), TIN data, square grid-based DEM products (e.g., Advanced Spaceborne Thermal

Emission and Reflection Radiometer Global DEM [ASTER GDEM] and Shuttle Radar Topography Mission [SRTM])
and more. This generation can be achieved with the interpolation methods proposed for traditional grid-based The pattern of neighboring grid coordinates: (a) the pattern for cell (i, j) in a square coordinate system; and (b) the pattern for even-row cell (i, j) and odd-column cell (i, j) in a hexagonal coordinate system, respectively DEMs, where the only difference is in calculating the geographic coordinates of the cell centroid. Given the availability of these sources, existing square grid-based DEMs are the most reasonable choice for this process.
The centroids of all pixels in the existing square grid-based DEMs are regarded as discrete sampling elevation points with a regular distribution. Then, the DEM establishment method based on these sampling points is utilized to build a hexagonal grid-based DEM.
The detailed procedure is as follows: 1. Divide the study region into a set of basic hexagons, where the size of the hexagon depends on the application goals of the DEM. Based on the basic hexagons and the topological relationship between these hexagons, the respective tessellation model is established.
2. For each hexagon, identify the four closest discrete sampling elevation points to its centroid as participating interpolation points. This is based on the correspondence between hexagons and the discrete sampling elevation points transformed from existing square grid-based DEMs, as shown in Figure 3.
3. Calculate the elevation value of each hexagon centroid as the elevation value of its whole cell based on the bilinear interpolation method (Schowengerdt, 2012). Taking the hexagon with central point P as an example, points A, B, C, and D are the participating interpolation points, as presented in Figure 3. The elevation value of point P can be expressed as: where z p is the elevation value of point P; (x p ,y p ) and (x a , y a ) are the geographical coordinates of points P and A, respectively; and z a , z b , z c , and z d are the elevation values of points A, B, C, and D, respectively. (3) The correspondence between hexagonal grids and discrete sampling elevation points

| E X TR AC TI ON OF VL S BA S ED ON HE X AG ON -G RIDDED S TRUC TURE
The method of VL extraction derives from the natural characteristics of surface runoff; that is, water always flows down the steepest slope (Fairfield & Leymarie, 1991;O'Callaghan & Mark, 1984). Generally, the automatic extraction of VLs can be summed up in five steps: (1) depression filling; (2)

| Depression filling
In gridded DEMs, depressions (or pits) refer to pixels with lower elevation values than their surrounding pixels, resulting in a lack of outlets to lower regions. These depressions are divided into two categories: natural depressions, showing real landscape features; and spurious depressions, which are introduced by man-made error.
Natural depressions are normally rare or absent in most terrain types, and are far less common than spurious depressions (Marks, 1984) caused by technical issues and interpolation errors during the collection and processing of a DEM (Martz & Garbrecht, 1998). The existence of these depressions fails to simulate water from one unit to its neighbor until it reaches the border of the DEM via the flow direction, making flow paths discontinuous. Thus, their identification and removal is a pre-processing step for automated extraction of VLs based on gridded DEMs.
Extensive algorithms have been proposed to fill depressions. Zhou, Sun, and Fu (2016) summarized three representative algorithms, the Jenson-Domingue algorithm (Fairfield & Leymarie, 1991), Planchon-Darboux algorithm (Planchon & Darboux, 2002), and priority-flood algorithm (Wang & Liu, 2006), respectively. Many studies have shown that the priority-flood algorithm represents a considerable improvement for its simplicity, speed, and versatility in comparison with the other two algorithms.
The priority-flood algorithm was first introduced by Ehlschlaeger (1989). In contrast to common operations with local 3 × 3 neighborhoods, this method identifies the location of general water flow and drainage basins by simulating the pattern recognition methodology and human behavior from a macro perspective. In this algorithm, the DEM cells are inundated from the outside to the inside, which is similar to the inundated process of natural terrain; that is, the water level rises gradually from the lowest point of the surface up to the whole terrain. Soille and Gratin (1994) first applied this algorithm to depression filling using DEM data with integer elevation value.
Then, Wang and Liu (2006) extended it to floating-point data, which is known as the W&L algorithm, and is a classic depression method based on the priority-flood algorithm. Therefore, the W&L algorithm is a good choice for implementing depression filling in hexagonal grid-based DEMs.
The W&L algorithm employs the minimum priority queue data structure and least-cost search technique to build the optimal spill path progressively based on the spill elevation of each cell.
1. Spill elevation: the minimum elevation value that the cell needs to be raised by to make water spill from this cell to the outlet.
2. Least-cost search technique: the method preferentially searches for the direction of the least cost in the search process. In the W&L algorithm, the least cost refers to spill elevation. Consequently, the optimal spill path is the lowest spill elevation path.
3. Minimum priority queue data structure: this structure is utilized to accomplish three basic operations for leastcost search: insert, find-minimum, and delete-minimum operations.
The lowest cell within the border of the DEM is regarded as an outlet. Then, in comparison with the elevation differences of neighborhoods, the elevation of each cell is raised by its spill elevation. In essence, the process of depression filling is carried out by reverse searching from potential outlets to the interior.
Analogously, the six-connected neighborhoods operation is introduced to the above reverse searching for the hexagonal grid structure. The specific steps are as follows: 1. Set the elevation of the cells contained in the boundary of the DEM as the spill elevation. Then, push these cells into a minimum priority queue and mark them as processed.
2. Obtain the head element and corresponding hexagon c in the minimum priority queue as the processing cell.
Search for unmarked cells among its six neighborhoods and calculate the spill elevation value for each unmarked cell with the formula below: where H n is the spill elevation value of the nth neighbor; E c and E n are the original elevation values of hexagon c and its nth neighbor, respectively.
3. Record the spill elevation of the unmarked grids calculated in step 2, push them into the minimum priority queue, and mark them as processed. Then, delete hexagon c.
4. Repeat steps 2 and 3 until the minimum priority queue is empty. The depression filling is then finished. In the second iteration, cell C is the head element. It is easy to obtain the unmarked cells F and G among the neighbors via searching. The calculated spill elevations of these two cells are both 40, which is higher than their original elevations. It means that cells F and G must be in a depression. Therefore, the elevation values need to be raised to 40. Then, push cells F and G into the priority queue and delete cell C. The depression-filling operations for cells F and G are finished. With the same operation, cells K and J are filled in the fourth and fifth iterations, respectively.
The entire process of depression filling is completed in 16 iterations. Consequently, the hexagon-based W&L algorithm is successful in depression filling with six-connected neighborhood operations.

| Deterministic six-node algorithm
Flow-routing algorithms are proposed to estimate the outflow from the current pixel to one or more adjacent pixels. According to the number of downslope pixels, they are divided into two kinds of generic methods: (1) the SFD method, which allocates the total amount of flow from the given cell to one of its eight neighborhoods, such The difference between the D8 and D6 algorithms results substantially from the neighborhood operation, which is manifested in two aspects. The first is shown in terms of distance measurement. Compared to the D8 algorithm, with a distance measurement relationship of 1: √ 2, the consistent distance between neighborhoods and the central hexagon simplifies the slope calculation step and results in the same flow arcs between adjacent cells.
Another aspect is presented in the angular resolution of the flow direction, where it is 60° for the D6 algorithm and 45° for the D8 algorithm. This variable determines not only the angle of the topological connection of the flow arcs between two algorithms in the local results, but also the angle at which branches drain into the trunk in a global result. In conclusion, the differences in neighborhood operations between these two algorithms influence the characteristics of the extracted VLs.

| Flow direction in flat region
With the operation of depression filling, there are three classes of estimated results by the D6 algorithm: (1) Only one neighbor with a maximal slope value. An example is cell L shown in the left panel of Figure 5b, which has a determinate flow direction value of 6; (2) More than one neighbor with a maximal slope value greater than 0. Under this condition, the water is forced to flow toward the first cell among neighbors with maximal slope value in coding order. As shown in the middle panel of Figure 5b, the flow direction of cell B is 1; and (3) More than one neighbor with a maximal slope value equal to 0. The hexagonal cells in this class are represented as flat areas whose flow direction cannot be determined directly. In a case shown in the right panel of Figure  b. If there is no such neighborhood n, as described in step (a), take these neighborhoods as the central cell c in sequence. Then, perform step (a).

c. Repeat steps (a) and (b) until the pending array is empty. The flow direction calculation of non-flat cells is
finished. Figure 6a shows a case of a flat area caused by depression filling. Hexagons C, D, F, G, J, and K are original flat cells in Figure 6b. The detailed process of estimation is presented in Figure 6c, in which hexagons with yellow background and red edges correspond to the processing cell c in step 2, hexagons with yellow background but black edges are in the pending array, hexagons with blue edges are the eligible neighborhood n in step 2(a), and the remainder of the hexagons with black edges have calculated flow direction. Taking flat cell C as an example, among its adjacent cells F and G with steepest downward slopes, no one has determined flow direction. Then, cells F and G are regarded successively as central cells with the same operation. The operation for neighbor F fails to find the path toward the outlet. In the operation for neighbor G, its neighbor L only adjoins outlet cell P. Therefore, the path from cell C to outlet cell P is C → G → L → P, which determines the flow direction of cells G and C as 6.
Applying a similar operation for cells D, F, J, and K, the calculated flow directions are 6, 1, 2, and 1, respectively.
The estimated result, as shown in Figure 6c, identifies that this rule is capable of effectively avoiding the flow vectors pointing toward each other in flat areas.
The following steps are presented in Figure 7.

| Study region and data
Three sample areas lying in the Yangtze River Basin in China are selected as study regions. They are subcatchments of the Jialing River system, Wujiang River system, and upper Yangtze River system, respectively, as shown in Figure 8. These three regions are chosen because of: (1) the distinct topographical characteristics among them, including elevation range, relief degree, and landform type, which are utilized to analyze the suitability of the VL extraction methods on different terrains; (2) the well-developed river systems, which provide reference data for evaluating the extracted VLs based on the D8 and D6 algorithms; and (3) the diverse structures of river systems in these areas, which are used to test the ability to simulate different shapes and bending characteristics of channels for the D8 and D6 algorithms. By implementing the algorithms in these three areas with square and hexagonal grid-based DEMs, we can objectively analyze the influence of grid geometry on VL extraction and explore the advantages and disadvantages of hexagonal grid-based DEMs in VL extraction.
The detailed descriptions of each region are as follows: 1. The sub-catchments of the Jialing River system (region 1). Region 1 lies between the Jialing River and the Peijiang River in Central Sichuan province and covers an area of 216 km 2 with an elevation from 269 to 507 m. In this area, the landform is the middle hill in Sichuan Basin, with undulating, low-elevation mountains, and the well-developed river system has a parallel distribution.

F I G U R E 7
The steps of flow accumulation calculation, determination for the threshold of flow accumulation, and connection among flow arcs based on calculated flow direction 2. The sub-catchments of the Wujiang River system (region 2). Region 2 falls southeast of Nanshan district,

Chongqing. It is in a transition area between the Yunnan-Guizhou Plateau and the southeastern edge of Sichuan
Basin, where the terrain is dominated by middle mountains and the rivers are distributed in a dendritic pattern.
Its area is 200 km 2 , and its elevation is between 594 and 1,993 m. However, the central and southeastern areas of region 2 both have a low-relief land surface (Zha, Cheng, & Zhao, 2014).
3. The sub-catchments of the upper Yangtze River system (region 3). Region 3 is located in eastern Chongqing.
The northeast part is near the southeast region of the Liangping District, where the landform is a deep hill with a gentle slope. The remaining region is located northwest of Zhong County, which comprises the center of the Three Gorges region. It is a typical hilly landform with low mountains and a dense network of rivers. The entire region has a total area of 369 km 2 with an elevation ranging from 267 to 1,061 m.
In the experiments, the original data source used to establish the grid-based DEMs with both squares and hexagons is the ASTER GDEM (version 2), with 30 m resolution, which is available from the Geospatial Data Cloud.
Mapped river systems are regarded as the reference data to estimate the quality of the estimated VLs, which are obtained from a 1:250,000 digital line graphic produced by the national surveying and mapping department.
F I G U R E 8 The location of three study regions

| Comparison of generated grid-based DEMs
To ensure the same sampling density in tessellation models of hexagons and squares for the DEM structure, the same resolution is defined as the same area. Theoretically, the ratio of the horizontal grid width of a hexagon to that of a square is � 2∕ √ 3 at the same resolution, which is approximated to 1:0.93 in experiments. Based on the generation method proposed in Section 2.2, the data sources are resampled to 10 resolutions for gridded DEMs with hexagons and squares. Table 1 shows the horizontal grid width of these two kinds of grids under each resolution, where the numbers 1-10 are the resolution ID. Figures 9-11 show the established results of the hexagonal grid-based DEMs and square grid-based DEMs with resolutions of 1, 6, and 10. From the perspective of visualization, for each region the results for the highest resolution between these two grid-based DEMs are similar, as observed in Figures 9a,d, 10a,d, and 11a,d. With the decrease in resolution, the differences between the two grid-based DEMs become more significant; that is, the regional outlines and the boundaries between different elevation values are clearer and smoother presented by hexagonal grids in the lowest-resolution results, especially in resolution 10 of regions 1 and 2, as shown in Figures   9c,f and 10c,f within the red outline. However, the above differences are insignificant in resolution 10 of region F I G U R E 9 The comparison of grid-based DEMs with hexagons and squares at resolutions 1, 6, and 10 for region 1 3, because this region has a larger area with more pixels relative to the other two regions. Thus, the smaller the number of pixels in the region, the greater the visual differences between the hexagonal grid-based DEM and the square grid-based DEM.

| Comparison of extracted VLs
Based on the generated grid-based DEMs with both hexagons and squares under varying resolutions in Section 4.2, the VLs are extracted by the D6 and D8 algorithms according to the aforementioned extraction steps. After repeating the experiments, the flow accumulation threshold is set to 50 to ensure that the VLs corresponding to the reference drainage networks can be extracted at coarser resolutions. We evaluated the effectiveness between the hexagon grid and the square grid in VL extraction from DEMs by means of visual and quantitative comparison.

| Visual comparison
This approach can directly judge the performance of either the hexagon-based method or the traditional squarebased method by overlaying estimated VLs on the real river system. Figures 12-14 show the results of delineated VLs from the square-based method and the hexagon-based method, respectively, for regions 1-3 under various resolutions. As the resolution decreases, the total length and number of valleys in the delineated results produced by the two methods both show a decreasing tendency. As shown in Table 2, when comparing the number of valleys at the same resolution, there is less difference between the two methods. However, this is evident from F I G U R E 1 0 The comparison of grid-based DEMs with hexagons and squares at resolutions 1, 6, and 10 for region 2 the total length of the VLs extracted by the hexagon-based method being mostly longer than that extracted by the square-based method, especially at higher resolutions. Only in the cases of region 1 with resolution 10 and region 3 with resolution 9 are the lengths of the VLs extracted by the hexagon-based method equal to 66.00 and 144.14 km, respectively, which are shorter than those extracted by the square-based method (69.27 and 147.32 km, respectively). Furthermore, the bending characteristic of valleys experiences a steady loss with a coarsening of resolution. These results demonstrate that grid size is a critical factor affecting the results of extraction, whether for hexagonal or square DEMs.
Notably, the hexagon grid exhibits superior capability in retaining location accuracy and shape characteristics at coarser resolutions. Comparing the comparative results of the three regions, the hexagon-based method has the most similar performance with the square-based method shown in region 1 (see Figure 12), where the extracted main valleys of the two methods have a slight offset, however, the remainder of the valleys coincide with each other to a great extent. Even so, compared with the reference channels with blue lines, the hexagon-based results show better agreement, and they have outstanding performance in maintaining bending characteristics at a lower resolution. To better display this advantage, we take the results within the dark blue box of Figure 12 as an example to give a detailed description. In Figure 15, we only reserve the flows corresponding to the mapped river system for a clear comparison, which is the same in Figures 16 and 17. At resolution 1, although the flows in the square-based results are similar to the "blue lines" channels in shape, they deviate more from the "blue lines" channels in terms of location than they do by the hexagon-based method. With the decrease in resolution, the shapes of the streams extracted by both algorithms are gradually simplified, but the differences in shape loss are not obvious in resolution 3. Up to resolution 8, the streams in the square-based results have a higher Additionally, they become discontinuous in the area within the green outline in Figure 12. In contrast, the corresponding streams identified by the hexagon-based method are continuous and smoother, and the meanders are preserved to some extent.
The significant distinctions between the results of the two methods mainly occur in areas with a flat terrain, as with the cases in the central part of region 2 (the dark blue outline in Figure 13) and the northeast part of region 3 (the dark blue outline in Figure 14)  In addition, the details of valleys within the dark blue outline of region 2 in Figure 16 show that the gridded DEMs of both squares and hexagons have difficulty in recognizing the detailed shape and accurate location of a natural river system for flat surfaces. However, it also confirms that the hexagon grid outputs a better simulation for

| Quantitative comparison
To examine the aforementioned performance of the two methods in the three regions, the mean separation distance (MD) is employed to calculate the degree of spatial agreement between the DEM-derived valley networks and reference data. Given that it is difficult to automatically detect the correct extracted valleys corresponding to reference channels due to the massive spurious streams, this article redefines the MD value as the average vertical distance between the vertices contained in the reference data nearest to the VLs, as displayed in Figure 18. It is calculated as follows: where MD is the mean separation distance between the reference data and extracted VLs; d i is the vertical distance from the ith vertex of the reference data to the nearest VL; and n is the total number of vertices. To calculate the MD value, the reference vector data are divided into 30 m segments. Then, these segments are converted into vertices in ArcGIS with the "Feature Vertices to Points" toolbox. Theoretically, the smaller the MD value, the higher the accuracy of the extracted VLs matching the mapped river system. Figure 19 compares the calculated MD values of the generated VLs in the three sample regions among different resolutions, where the horizontal axis represents the resolution ID, the vertical axis represents the estimated MD value, the black line represents the result using the square grid, and the red line represents the result using the hexagon grid. In all regions, as the resolution coarsens, the MD value using the square grid exhibits a trend that decreases at first and subsequently increases, whereas in the hexagon-based results, the MD value increases as a whole. This illustrates the effect of grid size on the quality of DEM-derived VLs. It is obvious that the MD values using the hexagon grid are smaller than those using the square grid, which demonstrates the superiority of the hexagon grid in retaining the location accuracy. The only exception occurs in the case of region 1 at resolution 7, where the MD value is 70.03 for the hexagon-based result and 68.67 for the square-based result.
Computation of MD between the reference data and extracted VLs Figure 20 shows the calculated reduction ratio of the MD value produced by the hexagon-based method relative to the square-based method to measure the differences in spatial agreement between them. In this figure, the horizontal axis represents the resolution ID; the vertical axis represents the estimated ratio; and the blue, orange, and gray bars represent the results of regions 1, 2, and 3, respectively. The ratios are almost more than 10%, with the exception of regions 1 and 2 at resolutions 7 and 8, respectively. This result strengthens the improved quality of VLs estimated by the hexagon-based method. The most significant distinctions are presented in region 2 at 7 out of 10 resolutions (excluding resolutions 7, 8, and 10), resulting from the flatness of the central area. For the other two areas, region 3 has a larger ratio than region 1 at 6 out of 10 resolutions (excluding resolutions 3, 5, 6, and 10), which is also caused by the flat terrain of the northeast part of region 3. Hence, we conclude that grid geometry is also a critical factor influencing the results of the DEM-derived VLs, and has an apparent impact in the areas with flat terrain.
Comparing the MD values of the hexagon-based results among the three regions at the same resolution, region 2 holds the smallest values at 7 out of 10 resolutions (excluding resolutions 6, 8, and 10). The MD values of region 1 are slightly lower than those of region 3, except for resolution 1, meaning that the estimated VLs of region 1 match the mapped data by a narrow margin. Therefore, in the middle-mountainous region 2, the DEM-derived VLs have the best agreement with mapped river systems as a whole, which is not influenced by the slight deviation in the central area with flatness. With respect to regions 1 and 3, the generated VLs have worse performance. It is inferred that the DEMs of the hexagon-gridded structure would produce more accurate VLs in regions with larger relief, which is in accordance with the performance of the traditional DEMs of square-gridded structure.
F I G U R E 1 9 Comparison of estimated MD values between the square-based and hexagon-based methods: (a) region 1; (b) region 2; and (c) region 3 F I G U R E 2 0 The reduction ratios of MD value produced by the hexagon-based method relative to the square-based method

| CON CLUS IONS
The VLs extracted from square grid-based DEMs may be affected by insufficient isotropy in the terrain representation among the eight directions due to the inherent distinctive distance measurement. By contrast, as a regular grid in the tessellation model that represents the real planar space, hexagons have more advantageous performance due to their consistent connectivity, compact structure, isotropy of local neighborhoods, and visual advantages. Hence, this article proposes a hexagon-based method for VL extraction from DEMs. This method is applied to three sample regions to measure its effectiveness. By comparing this method with the traditional method based on a square structure, we can draw several conclusions as follows: (1) The hexagon grid has superior capability in maintaining a detailed shape and location accuracy. As the grid size increases, the VLs obtained from the hexagon-gridded structure match the mapped river system with a smaller error, and they have a lower loss rate for meanders; (2) The results between the two methods are quite different in areas with flatness, resulting from the parallel flows in various directions due to the inconsistent angular resolutions of hexagons and squares.
However, they have better agreement with each other in areas with higher relief; (3) The hexagon-based method has less improvement for VL extraction in flat areas. The D6 algorithm introduced to calculate the flow direction in this method has an inherent defect in the SFD flow-routing algorithm, which is prone to simulate parallel channels in areas with low relief, leading to unrealistic results; and (4) This article demonstrates that grid geometry plays a major role in the quality of VLs obtained from grid-based DEMs, which provides novel thinking in DEM analysis.
Given the advantages and disadvantages of DEM-derived VLs based on hexagonal grids, future research will focus on two aspects: (1) improving the flow-routing algorithm with the MFD method or modifying DEMs with mapped topographic data to enhance the accuracy of the extracted valleys for flat areas; (2) applying the proposed method to hydrological models, such as watershed analysis and hydrological parameter extraction. In addition, the hexagon-gridded structure will be extended to the extraction of topographic characteristic lines, topographic factors, and contour lines extraction.