A Method of Watershed Delineation for Flat Terrain Using Sentinel-2A Imagery and DEM: A Case Study of the Taihu Basin

Abstract

Accurate watershed delineation is a precondition for runoff and water quality simulation. Traditional digital elevation model (DEM) may not generate realistic drainage networks due to large depressions and subtle elevation differences in local-scale plains. In this study, we propose a new method for solving the problem of watershed delineation, using the Taihu Basin as a case study. Rivers, lakes, and reservoirs were obtained from Sentinel-2A images with the Canny algorithm on Google Earth Engine (GEE), rather than from DEM, to compose the drainage network. Catchments were delineated by modifying the flow direction of rivers, lakes, reservoirs, and overland flow, instead of using DEM values. A watershed was divided into the following three types: Lake, reservoir, and overland catchment. A total of 2291 river segments, seven lakes, eight reservoirs, and 2306 subwatersheds were retained in this study. Compared with results from HydroSHEDS and Arc Hydro, the proposed method retains crisscross structures in the topology and prevented erroneous streamlines in large lakes. High-resolution Sentinel-2A images available on the GEE have relatively greater merits than DEMs for precisely representing drainage networks and catchments, especially in the plains area. Because of the higher accuracy, this method can be used as a new solution for watershed division in the plains area.

1. Introduction

Drainage networks are essential to geospatial assessment, basin analysis, and applications for catchment delineation, flow statistics, flood risk assessment, and climate modeling [1,2]. They also play an important role in estimating the transmission of pollution and nutrients [3]. Drainage networks are a fundamental condition for watershed delineation and an indispensable component in hydrological modeling. Therefore, reasonable and accurate watershed delineation is a precondition for runoff, sediment transport, water quality simulation, and basin resource management [4].

With the arrival of DEM, there has been a growth spurt in studying DEM-based extraction algorithms for drainage networks [5,6,7]. DEM, which represents the continuous variation of relative elevation values in each pixel, can be used to identify a river direction or flow path on the ground surface. Many algorithms, including the drainage networks algorithm [8] and the delineating watersheds algorithm [9], are developed to automatically extract essential hydrological features with DEM. Those algorithms require filling DEM to avoid depressions on the slope and generate a reasonable stream network in the preliminary stage. The stream network process usually starts by allocating flow directions for each DEM cell, then analyzes the flow accumulation, and finally selects those cells that have a total higher threshold of flow accumulation than a defined value [10]. Many factors, including the DEM spatial resolution, the calculated algorithm, and the physical characteristics of the basin, can directly affect the accuracy of drainage networks derived from DEM data [11,12]. Using DEM with an algorithm has obvious advantages, because the watershed processing takes less time and is not influenced by human subjectivity.

Although DEM is widely used, many problems are inherent in its practical application. DEMs usually include shallow depressions, i.e., areas surrounded by higher elevation that could occur naturally or artificially in a river plain landscape. For the determination of stream networks, an effective method of eliminating pits or depressions is the stream burning algorithm [13]. The stream burning algorithm often identifies river channels or lakes that are invisible in the DEM, preventing serious errors in the streaming. However, many features, such as hardened roads, low dams, artificial rivers, and large lakes could result in significant depressions and limit the accuracy of the results, especially in low-relief areas [14]. In addition, the most freely available open-source DEM datasets have a medium resolution of 30 m, which limits the spatial accuracy. Therefore, traditional DEM processing can generate unrealistic drainage network models due to large depressions and subtle elevation differences in local-scale areas, especially on river plains [15]. Automated watershed delineation with medium resolution DEM is not as accurate as high-resolution raster images of drainage networks for a plains area.

Researchers have made great efforts to improve the accuracy. The manual division of watershed boundaries from topographical maps is a very intense, time-consuming process, especially for large-scale studies, but is considered to be a very effective and accurate method [16]. Similar extraction and comparison of drainage networks could be obtained automatically from LiDAR-derived DEM or manually by drawing rivers from aerial photographs [17,18]. Automated delineation from LiDAR-derived DEM in plains areas have many problems. Stream networks extracted from satellite images or aerial photography are held to be truthful channels [19], primarily because the more detailed hydrological features in aerial photographs guarantee better drainage networks as compared with results obtained from DEM. Although high-resolution aerial photographs are ideal for extracting precise drainage networks and catchments, they are very expensive and not available on a large scale.

Due to their high spatial resolution and multispectral bands, optical remote sensing images, such as TM, SPOT, and Sentinel-2 images, can be used for water extraction. There are two common methods of extracting waterbodies. The first method which includes single threshold segmentation and an Otsu threshold segmentation algorithm, is to process the waterbody index threshold on the image such as normalized difference water index (NDWI), the modified NDWI (MNDWI), or the enhanced water index (EWI) index. Another method uses supervised classification, which employs training samples to select different trains, identify the characteristic parameters, and establish a discriminant function to distinguish different classifications. The second method usually involves extensive calculations, manual intervention, and a long time to obtain training samples.

For this situation, high-resolution satellite images, freely available on GEE, are a better choice for extracting drainage networks and catchments, especially for plains areas [20,21]. We develop a method to correct drainage networks by extracting rivers, lakes, and reservoirs. In addition to DEM, high-resolution Sentinel-2A images are used as input data. First, we extract drainage networks from Sentinel-2A images of the Taihu Basin on the GEE Platform. Second, we obtain results that show improved quality of watershed delineation as compared with catchments derived from the DEM. Compared to medium-resolution DEM of drainage networks and catchments, the objective of this study is to demonstrate the relative advantages of the proposed method for making precise and accurate watershed delineation in plains areas using Sentinel-2A images available on the GEE, taking the Taihu Basin as a case study.

2. Materials and Methods

2.1. Study Area

As shown in Figure 1, the Taihu Basin, which occupies an area of about 36,500 km², is situated in the downstream area of the Yangtze River basin, surrounded by water on three sides comprised of the Yangtze River in the north, the Hangzhou Bay in the south, and the East China Sea in the east. It extends across four provinces of China which are Jiangsu, Zhejiang, Anhui, and Shanghai. The Taihu Basin is divided into four geomorphologic types, high mountains, gentle hills, alluvial fans, and aggraded plains which are mainly dominated by gently sloping plains covering more than 70 percent of the total area.

Its overall terrain is distinctly variable, with elevations varying from 0 to 1600 m. The Taihu Basin has dense artificial rivers and numerous lakes, the largest of which has an area of 2338 km². The river network in the plains area is dense, crisscrossed, and surprisingly similar to an urban road network in its overall pattern. The average annual rainfall is 1177 mm for the total basin, the spatial distribution of which gradually decreases from 1408 mm in the southwest Tianmu mountains to less than 1100 mm in the eastern coastal area and northern plains [22]. Since it has one of the highest gross domestic products (GDP) and the most urbanized basin area in China, many projects have been allocated for watershed-based research since the 2000s and earmarked until 2025. Due to serious historical water pollution and flooding in the plains area, precise watershed delineation is of great significance for water conservation in the region [23,24].

2.2. Data Source

2.2.1. Sentinel-2 Images

GEE provides an online coding environment that can guarantee relatively rapid, distributed-cloud geospatial processing and analysis of large satellite image datasets that are convenient for everyone. It provides access to publicly available Sentinel-2A images and large-scale remote sensing analysis algorithms [25]. Therefore, GEE is used to code, analyze, and generalize the principal rivers, lakes, and reservoirs that could be derived from Sentinel-2A imagery of the Taihu Basin.

Corresponding to the latest temporal resolution and highest spatial resolution data available freely in GEE, Sentinel-2A OLI images that completely covered the whole Taihu Basin were collected from 1 May 2018 to 31 December 2018. Sentinel-2A images, obtained from the GEE platform throughout 2018, were considered to be Level 2 products and archived as top-of-atmosphere reflectance collections. The preprocessing work, including stitching images and masking clouds, was coded and executed through the developed application programming interface (API) in the GEE platform [26]. To ensure complete coverage with little cloud cover, less than 10% of scenarios were selected from the image sequences, using each specific path and row shown in Figure 2. To reduce the influence of seasonal variability on the identification of waterbodies, satellite images in summer were preferred, because plentiful water contributed to the extraction of rivers and lakes during this period. A total of three rows and four paths in the Sentinel-2A images, and approximately 15 scenarios in 2018, were obtained and used to satisfy the basic coverage needs of the Taihu Basin.

2.2.2. DEM Data

With approximately 30 m resolution and high vertical accuracy (less than 8 m), the SRTM 1 arc-second DEM data, acquired from the United States Geological Survey, is produced for the year 2000 and provided the GCS_WGS_84 geographic coordinates for the Taihu Basin. The SRTM data is easily downloadable from https://earthexplorer.usgs.gov/. The entire DEM is composed of eight scenario small DEM images of the Taihu Basin, whose respective names are n30_e119_1arc_v3, n30_e120_1arc_v3, n31-e119_1arc_v3, n31_e120_1arc_v3, n31_e121_1arc_v3, n32_e119_1arc_v3, and n32_e120_1arc_v3.

2.3. Processing Methods

The overall approach used in this study for the watershed delineation of the Taihu Basin is shown in Figure 3. Significant characteristics of the proposed method include the following: First, GEE, a distributed-cloud satellite image processing platform, was applied to process a number of Sentinel-2A images covering the Taihu Basin within a short timeframe. It is a relatively automatic, precise means of extracting waterbodies and determining the drainage network mainly from Sentinel-2A images, rather than from the DEM data. Second, it delineated catchments for each hydrological element by considering the modifying flow direction of each river, lake, and reservoir, in the plains area. These steps are discussed in detail in the following sections.

First, we obtained the actual rivers, lakes, and reservoirs from the Sentinel-2A images to make up a drainage network, as opposed to results that were extracted from the DEM. Secondly, for the key point of watershed delineation we acquired a drainage network and modified the flow direction of the hydrological features to obtain accurate catchments in the Taihu Basin. Therefore, the major process, shown in Figure 3, was as follows: (a) Calculate the MNDIW index with GEE; (b) use the Canny edge detector to extract water bodies; (c) determine the flow direction for each river, lake, and reservoir; (d) determine the flow direction of the overland flow; and (e) conduct subwatershed delineation for each river, lake, and reservoir.

2.3.1. MNDWI Index Introduction

Most methods that detect waterbodies from multispectral satellite images rely on the fact that the absorbing radiation effect of waterbodies is evident in the near-infrared and even longer bands. This principle enables them to highlight waterbodies clearly with raster band calculation, using the NDWI index to represent differences in the spectral curve and highlight the waterbody information in the image. Due to its use of the near-infrared band to replace the shortwave infrared band, the MNDWI index appears to be more sensitive than the NDWI index, resulting in a greater ability to highlight waterbody information [27,28]. We processed Sentinel-2A images to determine the surface reflectance and followed formula (1) to compute the MNDWI for surface water detection on GEE. In fact, the spectral reflection spectrum curve of water usually varies, owing to different subtle objects in the water. This makes it hard to adopt a spectral index with a fixed threshold to detect waterbodies on the ground. Numerous regional errors could be generated if a zero threshold is applied to detect waterbodies from water index images. There are regional gaps in the spectral properties of different waterbodies.

$$MNDWI=\frac{{\rho }_{Green}-{\rho }_{SWIR}}{{\rho }_{Green}+{\rho }_{SWIR}}$$

(1)

Green indicates a green band (band 3, 10 m) and SWIR indicates a short infrared band (band 11, 20 m) in

Table 1. An introduction to the bands used in the Sentinel-2A images.

Sentinel-2 Bands	Wavelength (µm)	Resolution (m)
Band 2—Blue	0.490	10
Band 3—Green	0.560	10
Band 4—Red	0.665	10
Band 8—NIR	0.842	10
Band 11—SWIR	1.610	20

. Due to automatically and dynamically projecting the 20 m (band 11) to the 10 m (band 3) resolution in GEE, we were able to directly calculate the MNDWI between the green and SWIR bands without considering the resolution difference of the Sentinel-2A images. According to previous studies, varying thresholds in the MNDWI index have generally required manual threshold adjustment and limited their applicability. The challenge was to establish a varying optimal threshold method, allowing for more fully automated waterbody detection. On the basis of a statistical histogram of all the MNDWI values in a satellite image, histogram-based threshold detection methods help to overcome this problem [29].

2.3.2. Canny Edge Detector for Extracting Water Polygons

The Canny edge detector is a widely used method applied in visualization computer systems, which allows for accurate edge detection in images [30]. The scope of this study employed histogram-based threshold methods, using a Canny edge filter to detect the boundaries between water and nonwater pixels. The Canny edge algorithm, which is applied to a MNDWI image with a threshold, revealed water edges that are in rapidly changing areas of an index value [31,32]. In terms of the MNDWI index, this situation generally involves obvious value changes in image pixels when shifting from one neighborhood to another. Waterbody and overland pixels are distinguished by applying buffer analysis to count the MNDWI values, generate a histogram, and determine a segmentation threshold. In a real situation, the MNDWI histogram presents bimodal statistical distribution, and therefore an obvious distinction could be made between waterbody and overland pixels.

In this study, we used binary images where the boundaries between water and nonwater pixels were clearly defined. The algorithm operated as follows: (1) For the image smoothing, generally, edge detectors are prone to noise, therefore, the image was smoothed with a square-sized Gaussian structural element, usually of 5 × 5 size; (2) for the gradient intensity calculation, the gradient magnitudes and directions were computed, and the gradient magnitude defined the edge, with a high gradient magnitude generating a rapid change in color (an edge), in contrast to a low gradient magnitude. The gradient direction, as the name implies, defined the orientation of the edge; and (3) for the hysteresis thresholding, this stage removed small pixel noises based on two threshold values of the intensity gradient, and discarded pixels below a certain threshold.

GEE supplies a built-in interpretation of the Canny edge detection algorithm, which is utilized to compute the edges [33]. To estimate the boundaries of a river under cloudy circumstances, multiple historical Sentinel-2A satellite images were used to extract composite waterbodies, representing the boundaries of rivers and lakes. Multitemporal composite results were used to identify segments of a river that were initially unknown due to the presence of clouds. As shown in Figure 4, the process of extracting waterbodies comprised the following five steps: (a) calculate the MNDWI index from the Sentinel-2A images, (b) detect sharp edges with the Canny edge detector, (c) count the threshold values, (d) extract the waterbodies using the statistical threshold in the MNDWI image, and (e) compose waterbodies from multitemporal Sentinel-2A images. The detailed code can be found in Appendix A.

After extracting the waterbodies from the GEE platform, the rivers, lakes, and reservoirs were generalized in the ArcMap software. First, each river segment was simplified into a single centerline and the width calculated using the “Collapse Dual Lines to Centerline” tool in ArcMap. Second, the lakes were manually distinguished from reservoirs based on the fact that lakes in the plains area were usually less than 10 m deep or even below zero, while reservoirs in the mountain area were higher than 20 m. Finally, waterbody polygons were preliminarily processed so that they could be analyzed and organized as linear rivers, polygonal lakes, and polygonal reservoirs based on their hydrological features in ArcMap.

2.3.3. Flow Direction Determination for Each River, Lake, and Reservoir

The well-known D8 algorithm was used to determine the flow direction from each cell to its downslope neighbor or neighbors on a topographic surface [34]. This method can track and represent the flow from each pixel to neighboring pixels by using eight discrete flow angle values. Due to the large numbers of ring-shaped and well-shaped river structures in the plains area, a traditional D8 algorithm, used to express dendritic nodes, could not properly digitize the flow direction of river nodes connecting more than three rivers [35]. In order to digitize multiple directions uniformly into a raster value, river nodes connected to more than three river segments had to be recorded. Therefore, the multi-D8 algorithm, which is an extension of the traditional D8 algorithm, was used to express multiple flows from a crossed node. The multi-D8 algorithm could express multiple flow directions from a node cell in a binary exponential manner, so there was no confusion when multistreams were encoded and decoded in a digital representation [36]. As shown in Figure 5, there was a detailed explanation for 30 (2¹ + 2² + 2³ + 2⁴) equaling four flow directions in a node cell.

Because the spatial resolution of Sentinel-2A images was 10 m, the principal rivers represented by image pixels in the GEE were usually wider than 100 m in the results. In a sense, this was a generalization of rivers from satellite imagery, neglecting narrow rivers of the same scale. In this study, a variety of means were used to comprehensively assess and analyze the flow direction of each river, such as hydrological observation data, DEM height differences, and node directions for each river segment. First, hydrological observation data from observation stations directly determined the flow direction of each generalized river, but it could not cover all the rivers in the Taihu Basin. Second, node junctions consisting of rivers, lakes, or reservoirs in a hydrological network were useful for determining the flow direction of undefined river segments. If a node junction was composed of three main rivers, an undefined river flow could be estimated from the flow directions of two other rivers [37]. Third, the flow direction of each main river was estimated according to the average elevation difference between the start point and end point values, which could be extracted from the DEM data of a river segment [38]. Following the above three steps, we determined the flow direction of each main river in the Taihu Basin. The hydrological observation data was provided by the Taihu Lake Basin Hydrological Information Service System on the internet, which was maintained by the local water conservancy bureau and can be visited at http://218.1.102.107:8100/indexWater.html. The Taihu Basin Hydrological Information Service System covered the 98 hydrological observation sites that collected daily hydrological data for the study area, including flow rate and water flow direction.

Since the catchment of a lake or reservoir usually had a specific shape, it was not necessary to determine flows for a cell overlapped by itself. Instead, a unique identifier was assigned to each cell that was overlapped by an integer value of 500 for a lake and 600 for a reservoir. As shown in Figure 6, the flow direction of a lake or reservoir was represented by its edge line with an arrow.

2.3.4. Flow Direction Determination for the Overland Flow

All the flow direction grids comprised river segments, lakes, reservoirs, and overland flows in the Taihu Basin. The flow direction of the overland flow derived from the overland DEM is shown in the Figure 7. First, the overland DEM need to be pretreated by masking waterbody grid cells that include the main rivers, large reservoirs, and lakes, then filling the pits or depressions of the overland DEM in ArcMap. Second, the flow directions of the overland flow can be digitized through the D8 algorithm and calculated from the overland DEM in ArcMap. Therefore, we cut the original DEM into pieces of overland polygons, which excluded the geometry of rivers and lakes; calculated the D8 flow direction for each overland flow; and merged the flow direction grids using the hydrological tools in ArcMap. Compared to calculating the flow direction directly from the whole DEM, this method prevented a large number of straight flow directions in depression areas and made streams more precise in the overland flows.

2.3.5. Subwatershed Division

The Taihu Basin can be divided into three kinds of catchments: Reservoirs, lakes, and onshore catchments. The sub-basin partitioning method was distinct from traditional methods that made use of flow direction, flow accumulation from DEM, and outlets to the complete catchment division in the Taihu Basin. This method extracted waterbodies based on the high-resolution Sentinel-2A satellite images and composed a drainage network for the plains area. The catchment area was delineated for each actual main river in this section.

On the basis of the above modified flow directions of each river, node, lake, reservoir, and overland flow, an entire flow direction raster grid covering the Taihu Basin was generated by the overlapping process of ArcMap, as shown in Figure 8. First, the sub-basin division searched the inflow paths of each grid in the upward direction, from the end to the start point of a river segment. Second, an inflow path for each grid in a river segment linked all the same inflow adjacent grids of the overland flow into a chain. Third, all inflow paths flowing into a unique river were marked with the same integer identification (ID). The sub-basin of each river segment consisted of grids with the same ID. Finally, a few sub-basins with an area of less than 10 km² threshold were combined into neighboring sub-basins to meet the requirements. In this way, we could accurately divide the two pavement catchments of each river, the catchment of each lake, and the catchment of each reservoir.

It was clear that each river segment had only two catchments, on the left and right. Each river node was connected to at least three catchments, which had the same directions as their river segments. The catchments’ topology, which was important for the calculation and analysis of the distributed hydrological modeling, was based on the flow directions in the river network and organized according to hydrographic nodes. Following the determined flow directions of the stream network, the relationships between upstream and downstream nodes became clear. Therefore, it was easily possible to derive the up–down relationship of the catchments from the downstream and upstream nodes in the river network. In this way, we established a watershed delineation with clear upstream and downstream relationships and provided accurate catchment data for watershed application in the Taihu Basin.

3. Results

3.1. Accuracy Analysis for Waterbody Extraction

We compare the effects of waterbody extraction among rivers, lakes, and reservoirs in different terrains. Although some very narrow rivers have broken, the overall structure of rivers can be shown in the river plain. The effects make it better in reservoirs and lakes.The effect is shown in Figure 9.

The Global Land Cover 30 (GLC30) is produced using multi-seasonal Landsat data based on the random forest-based mapping framework of Tsinghua University, for which the overall accuracy of the global land cover mapping is 73% at 30 m spatial resolution. The classification types include cropland, forest, grassland, shrub land, wetland, water, impervious surface, bare land, and snow. To validate the effect of waterbody extraction, this study selected validation points from GLC30, which can be downloaded from http://0-data-ess-tsinghua-edu-cn.brum.beds.ac.uk/ [39]. A random sampling method was used to select 500 water surface points and 500 nonwater surface points from GLC30, in which the water and nonwater verification points were evenly distributed across the Taihu Basin. The sample points that were extracted from the GLC30 data were checked and modified manually using the Google Earth software.

$$OA=\frac{1}{N}{\displaystyle \sum _{i=1}^M{n}_{ii}}$$

(2)

$$P_c=\frac{1}{N^2}{\displaystyle \sum _{i=1}^M{n}_{i+}*n_{+i}}$$

(3)

$$Kappa=\frac{P_o-P_c}{1-P_c}$$

(4)

The overall accuracy (OA) was calculated from the confusion matrix using Equation (2), which essentially represented the proportion that was classified correctly from all the reference sites. Cohen’s kappa coefficient was generated from a statistical test to evaluate the accuracy of a classification, which was calculated from Equations (3) and (4). In Equation (4), Po represented the ratio of the total sample, i.e., the ratio of the correct classification of each category to the total sample, which meant that Po was the overall accuracy. The user’s accuracy was calculated by taking the total number of correct classifications for a particular class and dividing it by the row total. In this study, the user’s accuracy was calculated from the above error matrix. As shown in

Table 2. The confusion matrix and accuracy evaluation of waterbodies in the Taihu Basin.

Waterbody (2018)	Water	No-Water	Sum Pixels	Accuracy
Water	441	47	488
No Water	59	453	512	OA = 89.3%
Sum reference pixels	500	500		Kappa = 79.3%

, the user’s accuracy of water was 90.3%, the overall accuracy was 89.3%, and the Kappa coefficient was 0.793.

3.2. Drainage Networks and Catchment of the Taihu Basin

This section introduces the results of watershed division analysis through the proposed method in the study area, which occupies a large plain in southeastern China. The northwestern and southern parts of the Taihu Basin are dominated by undulating hills and low mountains, accounting for approximately 20% of the total area, while the central and eastern parts are flatlands. The wide plain area is filled with dense river channels, large lakes, and a considerable number of reservoirs, which makes the hydrological elements relatively complex in the Taihu Basin. The hydrological elements were obtained and identified from Sentinel-2 images using the Canny algorithm on the GEE platform. Considering the resolution limitation of the DEM and Sentinel-2 images, some small lakes and reservoirs with an area of less than 10 km² were replaced by the centerlines of each waterbody. It was important to maintain the actual topological connections between each linear river, polygonal lake, and polygonal reservoir, which were vital for the pretreatment processing of the drainage network.

As shown in Figure 10, almost 2291 river segments with a width greater than 100 m were extracted, filtered, and generalized from the Sentinel-2 images. Seven lakes and eight reservoirs, which had clearly topological entrances and outlets to river segments, were extracted. The final catchments in the Taihu Basin were divided into the following three types using the proposed approach: Lake catchments, reservoir catchments, and overland catchments. The flow directions in a stream were expressed by arrows along the centerline of a river, or the edge line of a lake or reservoir. Because the streams used in the catchment division were the real drainage networks from the Sentinel-2A images, they retained the original crossing river structures of the Taihu Basin. Each river segment had two overland catchments on the left and right. The divided subwatersheds were detailed, comprising a total of 2306 catchments across the entire basin. These actual catchments with topology connectivity provided more accurate fundamental spatial data for distributed hydrological analysis.

3.3. Comparison of the Stream Burning Methods in Arc Hydro and HydroSHEDS

Simultaneously, this paper compared the sub-basin division results of the stream burning methods in the Arc Hydro and HydroSHEDS tools. In order to contrast this way with the stream burning method frequently used in the plains area, we used rivers and DEM to divide the Taihu Basin through Arc Hydro tools with an integrated stream burning method. The stream burning method in Arc Hydro mainly includes the following steps: (1) generating filled and flawless DEM, (2) using the DEM reconditioning tool to reduce a certain elevation artificially for vector rivers and burn the rivers into 30 m DEM, (3) generating flow direction grids and extracting streamlines based on a flow accumulation threshold, and (4) using streamlines and outlets to extract catchments using the Watershed tool. According to several processing experiments, a good result was obtained using the burning algorithm when the flow accumulation was set to 4000 ha in Arc Hydro. This result brought the source of the drainage network closer to the actual situation. The results of the stream burning method are shown in Figure 11a.

HydroSHEDS, a hydrological product that offer local- to global-scale catchment data, is mainly generated from SRTM DEM data at three arc-second spatial resolution [40]. Several procedures, such as eliminating small pits from elevation data (the sink filling method), combining digital river maps (the stream burning method), and calculating flow directions, are used in the data generation process and have been well documented in the literature. HydroSHEDS generates high-precision drainage network and subwatershed models, because it uses the above methods to produce hydrologically conditioned elevation and reduce errors [41]. As shown in Figure 11b, the streamlines derived by the proposed method were contrasted with HydroSHEDS datasets to evaluate the accuracy of the streamlines, which were downloaded from the official website https://www.hydrosheds.org.

The comparative results in Figure 11 showed that the streamlines obtained by stream burning in Arc Hydro and HydroSHEDS were generally similar in their geometric shape and spatial distribution. However, compared to actual river streamlines drawn from the Sentinel-2A images, HydroSHEDS did not solve the problem of ring-shaped and well-shaped river networks in the plains area. Meanwhile, it can clearly be seen from Figure 11 that streamlines obtained by stream burning in Arc Hydro and HydroSHEDS could not eliminate the large numbers of incorrect streamlines in some large lakes in the Taihu Basin. These fallacious streamlines, which were inconsistent with the actual situation, could introduce serious errors into hydrological applications used in the Taihu Basin.

4. Discussion

4.1. Differences from the Stream Burning Method

From the perspective of data precision, we find that there are some obvious differences between the catchments extracted using a traditional method and those obtained by using the proposed method in the Taihu Basin. The traditional method depends on 30 m DEM to define the flow accumulation and extract the stream network using a critical threshold value in Arc Hydro [42,43]. However, a digital hydrological network, which is composed of rivers, lakes, and reservoirs, is precisely extracted from 10 m high-resolution Sentinel-2A images covering the entire Taihu Basin. The spatial resolution of the satellite images and actual local hydrological elements have a significant impact [44]. The accurate and realistic rivers and lakes, which are important components of hydrological networks, are extracted from the Sentinel-2A images on the GEE and prevent many cumbersome manual vectorizations.

Compared with previous studies, the stream network extracted by the DEM in Arc Hydro is a simple dendritic structure, with flow direction represented by the standard D8 algorithm for each cell [45,46]. In fact, the usually circular and well-shaped river structures of the southern plains in the Taihu Basin are ignored. Meanwhile, multi-flow directions in the crossed nodes of rivers are redefined and digitized by the multi-D8 algorithm in the flow direction determination process of this study. The proposed method directly adjusts the flow direction of each river, lake, and reservoir in the real stream network, hence, it could avoid many unrealistic river channels in the lake area and retain the complete structure of the plains area. With the help of high-resolution Sentinel-2A satellite images, this automated method proves conducive to the timely and accurate division of the catchment area.

4.2. Catchments for Hydrological Simulation

The detailed catchment results are derived from rivers and lakes in the Sentinel-2A images and retain their crossed and annular structure. From the perspective of hydrological simulation, a simple dendritic structure is usually needed. In fact, the rivers could simply be changed to a dendritic structure by interruption and used for the hydrological slope and sub-basin processes in the SWAT model for a smaller area [47,48,49]. Simultaneously, the detailed extracted catchments that corresponded to these rivers could be merged and used in the hydrological runoff simulation. This finer division of catchments in the local plains may provide improved accuracy of rainfall runoff simulations.

The drainage network results could also be used for the main river channel calculation in the HEC-RAS model, which has the ability to calculate crossed and annular rivers in the plains area [50,51]. This offers the possibility of performing two joint models to improve flood simulation accuracy based on the extracted drainage network and catchments [52].

4.3. Further Improvement

Some aspects of the proposed method require further improvement in order to more readily obtain accurate catchments in the future. The finer the degree of the stream network, the greater the workload required to preprocess the flow direction cells. The flow direction of rivers that are assessed using hydrological data from observatory stations may be the most accurate method, and therefore more extensive data should be collected and utilized to correct the flow direction of rivers in the future. In addition, the consideration of artificial highways and railways, which form relatively high-terrain pseudo-topographic valleys in the plains area, could make the division of overland flow catchments more reliable and reasonable [53,54]. Water bodies derived from satellite images at different resolutions is corresponded to catchment areas of different scales. In addition, the multiresolution satellite images used in the GEE could be further considered and tested for their usefulness in extracting multiscale drainage network models and watershed delineation in the future.

5. Conclusions

Due to large depressions and subtle differences of DEM in the plains area, catchments division are difficult. The traditional method, using stream burning to delineate the watershed, is not sufficiently accurate for actual local hydrological purposes. This study proposes a method to delineate the watershed using high-resolution Sentinel-2A images as the data source for the plains area. On the basis of the Canny edge detection algorithm, rivers, lakes, and reservoirs are automatically extracted from Sentinel-2A satellite images, then combined to form the drainage network. Catchments are extracted for each hydrological element by directly modifying the flow direction of each river, lake, reservoir, and overland flow. A watershed is delineated into three types of catchment, i.e., lake, reservoir, and overland flow catchment in the Taihu Basin. The proposed method is found to be suitable for the Taihu Basin plains study area. Compared to the traditional method, the new method is more consistent with the actual precise watershed delineation. Therefore, this study provides a new method for extracting the subwatershed from Sentinel-2A images and providing more accurate basic data for subsequent hydrological applications, especially in the plains area.