Automated extraction of tidal creeks from airborne laser altimetry data

Journal of Hydrology 527 (2015) 1006–1020

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Automated extraction of tidal creeks from airborne laser altimetry data Yongxue Liu a,b,d,⇑, Minxi Zhou a, Saishuai Zhao a, Wenfeng Zhan c, Kang Yang a, Manchun Li a,b,⇑ a

Department of Geographic Information Science, Nanjing University, Nanjing, Jiangsu Province 210023, PR China Jiangsu Provincial Key Laboratory of Geographic Information Science and Technology, Nanjing University, Nanjing, Jiangsu Province 210023, PR China c International Institute for Earth System Science, Nanjing University, Nanjing, Jiangsu Province 210023, PR China d Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application, Nanjing, Jiangsu Province 210023, PR China b

a r t i c l e

i n f o

Article history: Received 10 December 2014 Received in revised form 4 May 2015 Accepted 30 May 2015 Available online 5 June 2015 This manuscript was handled by Konstantine P. Georgakakos, Editor-in-Chief, with the assistance of Marco Toffolon, Associate Editor Keywords: Tidal creeks Airborne LiDAR Channel extraction Geomorphology

s u m m a r y Tidal creeks (TCs) are transitional waterways between terrestrial and marine environments. Extracting geometric information for tidal creek networks (TCN) geometry from remote sensing is essential to understanding their characteristics, formation, and evolution. Currently, the major obstacles to automated recognition using digital elevation models (DEMs) derived from airborne light detection and ranging (LiDAR) data are low relief, varying widths, high density, strong anisotropy, and complicated patterns. Conventional methods, such as the optimal-elevation threshold method, the optimal-curvature threshold method, and the D8 method, cannot achieve satisfactory performance under these conditions. We propose an automated method for extracting tidal creeks (AMETC) using topographic features detected from LiDAR. Specifically, a multi-window median neighborhood analysis was designed to enhance depressions both in mudflat and marsh environments; a multi-scale and multi-directional Gaussian-matched filtering method was incorporated to enhance width-variant TCs; and a two-stage adaptive thresholding algorithm was implemented to segment low-contrast TCs. The AMETC was tested on two large LiDAR datasets of the Jiangsu coast with different resolutions. The quantitative assessments show that AMETC successfully extracted both small and large TCs from our study areas. The true positive extraction rate reached 95%, outperforming conventional methods. The AMETC is robust and weakly dependent on scale, and rarely requires manual intervention. Further applications suggests that the AMETC has potential for extraction of other types of channel features (e.g., badland networks and ravines). Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Tidal creeks (TCs) are linear/striped depressions embedded in the tidally dominated coastal landscapes. They consist of intricate systems of bifurcated channels, ultimately resulting in patterns that are among the most striking observed in natural environments. Functioning as channels for exchange sediments, nutrients, and pollutants (Pestrong, 1965; French and Stoddart, 1992), TCs propagate tidal waves, support fragile habitats (Kirwan and Mudd, 2012), and maintain delicate balances between sedimentary processes and hydrodynamics (Rinaldo et al., 1999; Allen, 2000; Lanzoni and Seminara, 2002; D’Alpaos et al., 2005; Marani et al., 2007; Mudd et al., 2010; Coco et al., 2013). Centuries of overexploitation and habitat transformation along coastal zones ⇑ Corresponding authors at: Department of Geographic Information Science, Nanjing University, Nanjing, Jiangsu Province 210023, PR China. Tel./fax: +86 25 89881181 (Y. Liu). E-mail addresses: [email protected] (Y. Liu), [email protected] (M. Zhou), [email protected] (S. Zhao), [email protected] (W. Zhan), [email protected] (K. Yang), [email protected] (M. Li). http://dx.doi.org/10.1016/j.jhydrol.2015.05.058 0022-1694/Ó 2015 Elsevier B.V. All rights reserved.

(Lotze et al., 2006), have jeopardized this equilibrium, resulting in a dramatically morphological adaptations of TCs. It is estimated that a total area of 13,380 km2 of coastal zones in China have been reclaimed between 1950 and 2008, and a further 2469 km2 of coastal reclamation is planned for 2012–2020 (Wang et al., 2014). As a consequence, salt marshes and upper intertidal zone are absence in many tidally dominated coastal zones. The morphological responses of TCs to intensifying anthropogenic disturbances, including reactivation, lateral migration, headward erosion, and bank erosion, may rework deposits, obliterating records of past environments, and posing more vulnerability to the coastal defense from waves and storms (Gabet, 1998; Hughes et al., 2009). To record changes in morphology, quantify network expansion and reduction rates, and understand future evolution trends, it is of critical importance to delineate tidal creeks precisely and objectively. Field measurements and visual interpretation of aerial images are the two conventional methods of mapping tidal creek networks (TCNs). However, the former is greatly restricted by poor accessibility and limited tidal exposure and is limited to local coverage;

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

while the latter is cumbersome, subjective, and rather impractical for repetitive observations or inspection of large areas (Mason et al., 2006). In contrast, airborne laser altimetry (Light Detection And Ranging, LiDAR) can represent TC geometry more accurately by providing high-resolution topography data (Klemas, 2011). Recent advances in LiDAR mapping have improved the efficiency of channel feature extraction. This technique has been applied to interdisciplinary research, including detection of typical channel-like TCs (Lohani and Mason, 2001; Lohani et al., 2006; Mason et al., 2006). Other applications include the recognition of channel heads (Clubb et al., 2014), moraines (Rutzinger et al., 2006), gullies (James et al., 2007; Baruch and Filin, 2011), channel-bed morphology (Cavalli et al., 2008), streams (Cho et al., 2011), coastal structural lines (Brzank et al., 2008), and surface fissures (Stumpf et al., 2013). The methods developed in interdisciplinary studies are all pertinent to the extraction of TCs, so we considered methods designed for both TCs and other channel-like features. In general, the relevant methods can be divided into the following four categories: (1) Flow-routing models. River networks can be extracted from digital elevation models (DEMs) according to their flow accumulation grids, calculated using by the classic D8 method (Ocallaghan and Mark, 1984), the Rho8 method (Fairfield and Leymarie, 1991), the FD8 method (Freeman, 1991), or the D-infinity method (Tarboton, 1997). (2) Thresholding methods. Fagherazzi et al. (1999) extracted tidal creeks from low-resolution DEMs based on combined elevation and curvature thresholds. Lohani and Mason (2001) introduced an adaptive elevation threshold method to locate TC fragments. Rutzinger et al. (2006) used predefined windows and thresholds on curve segments for moraine extraction. (3) Methods based on mathematical morphology operations. Cho et al. (2011) detected forested streams using a bottom-hat operation. (4) High-level image processing (HLIP) methods. Mason et al. (2006) presented a semi-automatic multi-level knowledge-based approach. High-gradient pixels were found using edge detectors, and this was followed by edge association, centerline generation, network repair, and channel expansion. Brzank et al. (2008) extracted structural lines in the Wadden Sea by fitting a pre-defined surface model—a two-dimensional (2D) hyperbolic tangent curve). Stumpf et al. (2013) detected surface fissures with a combination of Gaussian filters, mathematical morphology and object-oriented image analysis. Other methods beyond those designed for TCs may also be suitable for use. TCs are characterized by low relief, variable widths, strong anisotropy of waterway directions, and widespread distribution in inter-tidal zones. TCs are often of high density near tidal estuaries. All of the aforementioned characteristics present challenges to accurate extraction. Flow-routing models only achieve moderate performance in low-relief intertidal zones (Lohani and Mason, 2001). Thresholding and mathematical morphology methods use pre-defined threshold/mask/structure-element that reduce their adaptability and may limit them to use with low-resolution data (Baruch and Filin, 2011). HLIP methods use the geometry and geomorphology of linear features and yield reasonably good results for uncomplicated areas (Cho et al., 2011). Although HLIP is relatively promising, its automation and application require further improvements (Lohani et al., 2006). An automated method for extracting tidal creeks (AMETC) from LiDAR DEMs based on HLIP techniques is proposed in this paper. The method will help coastal geomorphologists detect and map TCs and better understand their formation and evolution. The AMETC was developed to be weakly dependent on scale, robust, and automatic. Weak dependence on scale helps in the detection of TCs of different widths across DEMs of different resolutions. The efficient detection of high-density TCs at any orientation across vast tidal flats makes the method robust. Automation is

1007

improved because the parameters used in the AMETC should be stable enough to minimize the need for manual intervention. 2. Study area and datasets 2.1. Study area The Jiangsu coast is located in the west side of the Southern Yellow Sea (SYS) (Fig. 1a). Its offshore area is characterized by the Southern Yellow Sea Radial Sand Ridges (SYSRSR)—the largest offshore tidal ridges on the Chinese continental shelf. Over the last few centuries, the Jiangsu coast has received large sediment loads from both the Yellow River and the Yangtze River (Wang et al., 2012), resulting in broad tidal flats (>5000 km2 above the lowest normal low water). These tidal flats are located in a macro-tidal environment, with an average tidal range of 3.9–5.5 m (Xing et al., 2012). Two tidal wave systems, i.e., a progressive Poincaré wave from the East China Sea and an amphidromic system in the Yellow Sea, converge near the center part of Jiangsu coast, i.e., Jianggang (Ni et al., 2014). Under these hydrodynamic conditions, the mean tidal range reaches a maximum in the Jianggang area (up to 9.28 m), decreasing southward and northward (Gong et al., 2012). In the broad Jiangsu coastal zone, intricate TCNs (Fig. 1b–c) have been sculpted by the tidal currents (especially ebb currents from tidal flats). The mean tidal creek density of the Jiangsu coast is approximately 3.61 km/km2, and can reach a density of 13.27 km/km2 in some tidal estuaries (e.g., Fig. 1b). A variety of drainage patterns occur in the region, including dendritic, parallel, radial, deranged, and pinnate pattern. During the past 40 years, over 2000 km2 of tidelands along the Jiangsu coast (including 1444 km2 of native salt marshes and 229.2 km2 of Spartina alterniflora salt marshes) have been reclaimed (Zhang et al., 2004; Zhao et al., 2015). The remaining salt marshes are no more than 213.21 km2 in total area, according to our interpretation from synchronous Landsat images (Fig. 1b–c). The mean intertidal flat width has decreased from 7.3 km to 2.9 km, and the mean elevation of the tidal flats has decreased from 1.54 to 0.31 m. As a consequence, these intricate networks are undergoing dramatic morphological changes in response to the changing physical conditions. 2.2. Datasets Two Airborne LiDAR DEM datasets of different spatial resolutions were used for our study. Both sets were acquired at low tide. 2.2.1. Dataset A (5400  6700 pixels, 2.5-m resolution The dataset covers the southern section of Jiangsu (Region A in Fig. 1a, enlarged in Fig. 1b), and the study area is approximately 226 km2. On July 17, 2009, the LiDAR data were acquired in first-return mode by an Optech ALTM 3100 at an altitude of 1500 m and a speed of 240 km h1. The laser wavelength was 1064 nm, the laser pulse rate was 85 kHz, the scan rate was 70 Hz, and the laser footprint at nadir was approximately 30 cm. The spacing of the LiDAR points was 2 m, resulting in a gridded DEM with a cell size of 2.5 m. The root mean square error (RMSE) of the DEMs was approximately 13.5 cm (GPS-RTK synchronous validation). 2.2.2. Dataset B (11,000  19,300 pixels, 5-m resolution) The dataset covers the central region of the Jiangsu coast (Region B in Fig. 1a, enlarged in Fig. 1c), and the study area is approximately 5300 km2. The LiDAR datasets were produced by the Jiangsu Provincial Bureau of Surveying, Mapping and

1008

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

Fig. 1. Location maps of the study area: (a) HJ-1 CCD image of the Jiangsu coast (R4G3B2, GMT 02:41, December 22, 2013); (b) tidal creeks in Region A; (c) tidal creeks in Region B.

Geo-information and were collected in April and May of 2006. The spacing of the original LiDAR points was 1 m, and the LiDAR DEMs used in this study were resampled to a cell size of 5  5 m using the bilinear interpolation method.

3. Methodology Extracting TCs from DEMs with image processing technology presents a number of challenges, which were addressed through a variety of techniques. Field investigations show that the tidal creeks of the Jiangsu coast exhibit the following common characteristics (Fig. 2a–f): (1) Continuous concavity. The TCs are embedded within the coastal zone, and creek depths vary from centimeters to meters with a V-shape/trapezoidal cross section. (2) Linear structure. The TCNs typically exhibit linear or striped features, although their general geometry can be complicated. (3) Varying width. The creek widths may vary from a dozen centimeters in the upper intertidal zone, to several kilometers in the lower

intertidal zone. (4) Strong anisotropy. Different orders of TC segmentation can exhibit preferred orientations that deviate from a random network. Based on the aforementioned characteristics, we extracted TCs from DEMs using the following methods (Fig. 3): (1) Median neighborhood analysis (MNA) was used to generate an effective indicator (i.e., residual topography) to enhance depressions (Section 3.1). (2) Gaussian-matched filtering (GMF) was adopted to enhance linear structures with Gaussian-like shape (Li et al., 2012) (Section 3.2). (3) The kernel of the GMF was rotated to better match TCs with all possible orientations (Section 3.2). (4) Width variation was addressed with multi-window and multi-scale strategies in the MNA and GMF. (5) Automation was accomplished using two-stage adaptive thresholding (TAT) that automatically segmented TCs (Section 3.3). Notably that the width variation of TCs is too large to not be simultaneously extracted under a unified standard. Our coastal field investigation showed that TCs less than 50 m wide are broadly and densely distributed but difficult to extract. Therefore, we defined tidal creeks less than 20 pixels wide

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

1009

Fig. 2. Tidal creeks on the Jiangsu coast: (a) a TC in a Spartina alterniflora salt marsh, (b) a TC in the transition zone between a Spartina alterniflora salt marsh and a bare tidal flat, (c) a pinnate TC drainage in the upper tidal flat zone, (d) a meander in the middle tidal flat zone, (e) an extra-large TC in the lower tidal flat zone, (f) a deep gulley in a large offshore sandbank.

as ‘‘small’’ TCs (i.e., 50 m wide in Dataset A) and those greater than 20 pixels as ‘‘large’’ creeks. The small and large TCs were processed separately. 3.1. Multi-window Median Neighborhood Analysis (MNA) A TC is a type of depression embedded in a coastal zone (Fig. 2a–f). Hence, elevation thresholds may not be reliable indicators for use in delineating TCs. While elevations of TCs are lower than the average elevation of their neighboring areas. Therefore, a median neighborhood analysis (MNA) was utilized to enhance depressions (Eq. (1)). The rationale for this analysis is illustrated in Fig. 4a.

Rw ¼ Median ðf ðx  u; y  v Þ þ wðu; v ÞÞ  f ðx; yÞ u;v

Fig. 3. Flow diagram for tidal creeks extraction from LiDAR DEM.

ð1Þ

where f is the DEM; w is the window used to calculate the median value; (x,y) and (u,v) are the pixel coordinates of f and w, respectively; Medianu;v is the median operation for a given window w; and Rw denotes the residual topography between the median and the DEM. The MNA converts depressions (tidal creeks, pans, etc.) to positive values and vice versa. To simplify further extraction, all negative values (i.e., positive relief) were automatically set to zero, and only values greater than zero were retained (Fig. 4a). The MNA response to TCs of different widths varies with the window size used. Small TCs can be effectively enhanced with a small

1010

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

Fig. 4. Multi-window MNA: (a) Conceptual graph of for MNA. (b) Block 1 from Dataset A (1000  1000 pixels, lighter = higher); (c) window size (W) = 9 pixels; (d) W = 49 pixels; (e) W = 99 pixels; and (f–i) the four profiles (I–IV) in (b) approximated with Gaussian curves (r is the standard deviation of the Gaussian function).

window size, while large TCs may go undetected (Fig. 4a), so a progressively increasing multi-window strategy was adopted to enhance width-variant TCs. Fig. 4c–e illustrate how TCs of different widths are gradually enhanced using increasing window sizes from 9 to 49 to 99 pixels. In general, large TCs can be easily enhanced using a large window size and can be segmented using a simple thresholding method. However, small tidal creeks are more common along the coast and are more difficult to delineate, so further adaptation (see Section 3.2) of the methods was necessary. 3.2. Small TC enhancement based on Gaussian-matched filtering (GMF) The concave cross sections of TCs usually exhibit shapes from trapezoidal to V-shaped, which can be approximated well using Gaussian-shaped curves (Fig. 4e–f). Therefore, two-dimensional (2-D) Gaussian-matched filtering (GMF) was used to enhance small TCs (Zhang et al., 2010; Fraz et al., 2012; Li et al., 2012). In this study, the 2-D matched filter was defined as a Gaussian function along the x-axis with the function repeated in a neighborhood along the y-axis (Li et al., 2012). For ease of analysis, the inverted second derivative of the Gaussian g 00r was used for the matched filter. Eq. (2) defines the one-dimensional matched filter, and Eq. (3) defines the 2-D matched filter.

x 2  r 2  x2 g 00r ¼ pffiffiffiffiffiffiffi e 2r2 2pr5

ð2Þ

 2  x  r2  x2 MF r ðx; yÞ ¼  pffiffiffiffiffiffiffi e 2r2 ; 2pr5

where jxj  3r;

jy j 

L 2

ð3Þ

where r is the standard deviation of the second derivative of the Gaussian function, which is related to the width of the targeted feature, and L is the length of the neighborhood along the y-axis to smooth noise. The matched filter given by Eq. (3) only has a peak response at certain alignment angles: h  p=2 (Zhang et al., 2010). Because the TC orientation was unknown, the matched filter was rotated from 0 to p at certain intervals to better capture tidal creeks with different orientations. Eq. (4) is the clockwise rotation of the matched filter for a given angle h.

8 h 0 0 > < MF r ðx ; y Þ ¼ MF r ðx; yÞ; 0 x ¼ xcosðhÞ þ ysinðhÞ; > : 0 y ¼ ycosðhÞ  xsinðhÞ:

where h 2 ½0; p

ð4Þ

The maximum response of the input image to the matched filter at a given scale r among all possible orientations h can be calculated as follows:

  g wr ðx; yÞ ¼ max Rw ðx; yÞ  MF hr 0hp

ð5Þ

1011

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

Fig. 5. Small-TC GMF results at different scales: (a–c) GMF results for Fig. 4b–d at r = 1.5 (lighter = higher); (d–f) GMF results for Fig. 4 b–d at r = 2.5; (g–i) GMF results for Fig. 4b–d at r = 3.5; (j–l) GMF results for Fig. 4b–d at r = 4.5. Notably the large TC (e.g., green circles) are gradually enhanced with the increasing scale (r). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

where Rw ðx; yÞ is the MNA result for a given window w, MF hr is the matched filter at a given scale r with an orientation h,  denotes the convolution operation, and gwr denotes the maximum response for all orientations. In this study, we used an angle interval of 10° to construct 18 different kernels. We compared the corresponding responses of all 18 orientations and only retained the maximum as the final response. In a previous study, the scale, r, of the matched filter was usually a fixed value selected on the basis of experience (Zhang et al., 2010). However, given the great variation in small TC widths (several pixels to dozens of pixels in a 2.5-m resolution DEM), a single scale may not be effective (Fig. 4e–h). A progressive increasing multi-scale strategy was adopted to enhance small TCs of different widths. Fig. 5a–l show the response of tidal creeks to the different window sizes (9, 49, and 99 pixels) in the MNA and different r values (1.5, 2.5, 3.5, and 4.5) used in the GMF (Fig. 4b–d). 3.3. TC Segmentation using two-stage adaptive thresholding (TAT) Global optima or user-defined thresholds are often used to extract water bodies from background data (Frazier and Page, 2000; Ji et al., 2009; Verpoorter et al., 2012). However, many indistinct tidal creeks would be ignored if a single optimal value was

used (Fig. 6b). We propose a two-stage adaptive thresholding (TAT) algorithm that consists of the following steps: 3.3.1. First-round thresholding Based on the GMF results gwr, the threshold surface T is calculated according to the local mean m(x, y) and the local standard deviation d(x, y) within a sliding w  w window (Eq. (6)). The 1

first-round segmented results bwr are illustrated in Fig. 6c. 1

bwr ðx; yÞ ¼



1; IF g wr > T ðx;yÞ 0; Otherwise

; where T ðx; yÞ ¼ mðx;yÞ þ k  dðx; yÞ ð6Þ

3.3.2. Neighboring effect removal The local threshold surface T is closely related to the local mean value, which is greatly influenced by the local maximum in the window. Strong local maximums can greatly increase the local mean, resulting in a high local threshold. Faint TCs neighboring high pixels may be misclassified as background, resulting in a spatial discontinuity (Fig. 6c). To suppress the neighboring effect induced by local maximums, pixels detected in the first round are replaced with the local mean m(x, y) (Eq. (7)).

1012

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

Fig. 6. Two-stage adaptive thresholding: (a) GMF results (r = 1.5) of Fig. 4c; (b) OTSU thresholding results; (c) first-round thresholding (w = 49; k = 0.2); (d) second-round thresholding (w = 49; k = 0.2). The small isolated segments (labeled in red) can be sieved by 30 pixels. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

( g 0wr ðx; yÞ

¼

mðx; yÞ;

1

IF bwr ðx; yÞ ¼ 1

g wr ðx; yÞ; Otherwise

ð7Þ

3.3.3. Second-round thresholding As in the first round, the local adaptive threshold algorithm is applied to the modified g 0wr , and the consequent binary results are merged with first-round results by a logical OR operation (Fig. 6d).

Through TAT, small TCs with low contrast can be more easily segmented from the background (Fig. 6d). According to visual interpretation, the TAT result is significantly better than results of the OTSU method (Fig. 6b). Moreover, the original DEM was not particularly de-noised, so many small isolated segments can be seen in Fig. 6d (labeled in red). These segments are shorter or smaller than extracted tidal creeks, allowing filtering of areas smaller than the minimum grain size (t). Fig. 7a–l shows the segmentation results for Fig. 5a–l.

Fig. 7. Small-TC segmentation results for Fig. 5a–l. The window size used in the TAT is equal to the window size used in the MNA, and k = 0.2.

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

1013

Fig. 8. The fusion strategy in the AMETC: (a) small-TC fusion results for Fig. 7a, d, g, j; (b) small-TC fusion results for Fig. 7b, e, h, k; (c) small-TC fusion results for Fig. 7c, f, i, l; (d) small-TC fusion results; (e) LiDAR DEM in Block 1; (f) MNA results at window size = 99; (g) large-TC extraction results for f (w = 99, k = 0.2); (h) the final TC extraction results.

3.4. Fusion strategy for the TC segmentation results Because of scale variance, tidal creeks of different widths are processed differently by the multi-window MNA and multi-scale GMF (Figs. 4b–d and 5a–l). Therefore, the multi-window and multi-scale binary results (Fig. 7a–l) need to be combined as follows:

Eðx; yÞ ¼

[

[

a1
bwr ðx; yÞ

ð8Þ

where bwr is the final TAT result at the specific window size (w) used in the MNA and scale (r) used in the GMF. For large TCs, the second-derivative Gaussian kernel using a small scale can only enhance their edges, leaving the centers unenhanced and undetected (green1 circle in Fig. 8d). Large TCs can be efficiently enhanced in the MNA using a large window (Fig. 8f) and effectively detected by the proposed TAT (Fig. 8g). To guarantee completeness, the results of the two methods were combined with the logical OR operation (Fig. 8h). 4. Results 4.1. Validation strategy 4.1.1. Validation areas To validate the robustness and adaptability of the AMETC, two datasets were applied: Dataset A (2.5-m resolution, 5400  6700 pixels, 226 km2, Fig. 1b) and Dataset B (5-m resolution, 11,000  19,300 pixels, 5300 km2, Fig. 1c). The global extracted results are illustrated in Fig. 1b–c. Due to space limitations, we selected five blocks (each 1000  1000 pixels) from Dataset A (Fig. 1b) and six blocks (each 2000  2000 pixels) from 1 For interpretation of color in Figs. 8–10 and 12, the reader is referred to the web version of this article.

Dataset B (Fig. 1c) to illustrate details of the extracted results. The eleven blocks were selected for the following reasons: (1) Blocks 1–8 are all near tidal estuaries where TCs are highly developed. In these eight blocks, S. alterniflora salt marshes are generally distributed along the seawall, whereas the remaining regions are dominated by bared tidal flats. (2) Blocks 9–11 (Fig. 1b) are located in the two sandbanks (i.e., the Dongsha sandbank and the Gaoni sandbank). These two offshore sandbanks are unvegetated, and the TCs widths here are generally larger than the widths of those distributed along the coastal zone. 4.1.2. Data truthing To generate a ‘‘true’’ TC reference map, the TCs in Blocks 1–5 (2.5-m resolution DEM, 1000  1000 pixels) were manually traced at a scale of 1: 5000. As noted by Mason et al. (2006), TCs are not always easy for an operator to define because of irregular banks and nested channels. In the interpretation, two discriminants were mainly used: (1) the elevation discrepancy of tidal creeks and their surrounds; and (2) the continuity of a TC. Three experienced interpreters were employed to delineate the outlines of TCs. The three independent manual digitizations were the merged by a logical OR operation to avoid omission of small branches. The interpretation was then converted to a binary image with a cell size of 2.5 m for use in the subsequent accuracy assessment. Blocks 6–11 (5-m resolution, 2000  2000 pixels) have four times more area than Blocks 1–5. Production of the reference maps for all six blocks was exceptionally laborious, therefore, only Block 6 was used to demonstrate applicability. The extraction accuracy of the other five blocks (Blocks 7–11) was evaluated by visual interpretation. 4.1.3. Validation metrics The performance of the AMETC was assessed quantitatively by pixel-to-pixel comparison of the reference maps and extraction results. Comparisons yielded four cases: true positive (TP), true negative (TN), false negative (FN), and false positive (FP)

1014

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

(Fawcett, 2006). In this study, five metrics were used: (i) the true positive rate (TPR), defined as the rate of pixels correctly identified as TCs: TPR = TP/(TP + FN); (ii) the false positive rate (FPR), defined as the rate of pixels erroneously identified as TCs: FPR = FP/(FP + TN); (iii) the error of commission (EC), defined as the ratio of erroneous TCs to all true TCs: EC = FP/(TP + FN); (iv) the error of omission (EO), defined as the ratio of omitted TCs to all true TCs: EO = FN/(TP + FN); and (v) Overall accuracy (OA), defined as the ratio of correctly classified pixels to the total number of pixels in the image OA = (TP + TN)/(N). 4.1.4. Comparison between methods Three widely used methods were used for comparison: the optimal-elevation threshold method (OETM), the optimal-curvature threshold method (OCTM), and the D8 method. The optimal value in the OETM and OCTM was selected by the OTSU method (Otsu, 1979). The D8 method is available in the ArcGIS software. It should be noted that the extracted results obtained with the D8 method were single-pixel width in most cases, while the reference map was non-single-pixel; therefore, the commissions (those segments located outside of the reference area) can be identified by a computer, while the omissions require manual delineation. 4.2. Extraction results and validation The AMETC was implemented in MATLAB R2013a (8.1.0.604). The key experimental parameters were as follows: (1) 10 windows, varying from 9 to 99 pixels in 10-pixel intervals, were used in the multi-window MNA (refer to Eq. (1)); (2) four scales (1.5 6 r 6 4.5, with 1.0-r intervals) were used in the GMF (refer to Eq. (5)); (3) the window size for calculating the local mean in the TAT was set to be equal to the window size in the MNA, and k was fixed at 0.2 (refer to Eq. (6)). To filter out small isolated segments, the minimum grain size (t) was set to 30 pixels. These parameters remained constant for Datasets A and B. Visual inspection showed that the AMETC successfully extracted most of the TCs and in particular extracted all of the small TCs in all of the blocks (Figs. 9i–l and 10g–l). A few omissions occurred in: (1) several segments of extra-large TCs (purple circles in Figs. 9i, l, 10g–l), and (2) a few upstream segments of small TCs. Small segments are clearly discontinuous in the LiDAR DEMs and were effectively detected in all enhancement procedures but were filtered out in the fusion procedure. Over-extraction mainly occurred in: (1) some man-made linear reclamation depressions along outer seawalls (green circles in Figs. 9i, k–l and 10g), and (2) some bed ripples formed on offshore sandbanks (green circles in Figs. 10k–l). These two linear depression types form differently than TCs, but are similar in geometry and can only be distinguished by experienced interpreters. Moreover, a few un-extracted or over-extracted pixels were distributed along extracted TCs because the manual reference maps were subject to errors of interpretation and screen resolution. The OETM, OCTM, and D8 methods were employed for comparisons. Strong threshold dependency caused poor performance on the two full datasets, so results (Fig. 9m–x) were extracted from the individual subset of eleven blocks (Fig. 9a–d). The OETM generally detected only large TCs, while almost all small TCs were un-extracted (Fig. 9m–p). Moderately sized TCs were not detected well in low-contrast regions (Fig. 9p). The OCTM performed better than the OETM (Fig. 9q–t). However, many TCs were not differentiated from their neighboring regions and the shapes of the TCs were not continuous. Moreover, only the edges of large TCs were extracted. The D8 method detected almost all of the small TCs (Fig. 9u–x), but many ‘‘false creeks’’ were found in very plain tidal

flats (yellow in Fig. 9u–x) or ‘‘straight lines’’ within large TCs (straight lines in Fig. 9x). The pixel-to-pixel accuracy assessments using the five metrics for Blocks 1–6 are presented in Table 1. The performance of the four methods differed over the five blocks: (1) Block 1–4, and 6. The TPR, EO, and the OA of the AMETC were the best among the four methods, varying between 95–97.63%, 5–2.37%, and 96.33– 98.31%, respectively. The FPR and EC of the AMETC were minimal, varying between 1.45–3.07%, and 4.12–6.59%, respectively. The TPR of the D8 method was competitive, varying between 89.97% and 94.79%; however, the EC was rather high, varying between 18.60% and 49.43%. In other words, it would be very laborious to distinguish true TCs from a high error-rate set in post-processing for the D8 method. The extracted TCs were only a small subset of the referenced TCs using OETM and the OCTM, so their TPR was significantly lower (generally no more than 40%). Although the FPR of the OETM performed the best (varying from 0.31% to 2.78%), it was meaningless because it would be very cumbersome to delineate numerous omitted TCs (Figs. 9m–p) in post-processing. (2) Block 5. The AMETC had a TPR (40.93%) significantly worse than that of D8 (98.15%) because (i) the extra-large TC in Block 5 (width > 800 m, i.e., >320 pixels) was not extracted well by the AMETC (purple circles in Fig. 9l), and (ii) too many ‘‘straight lines’’ were produced by the D8 within the extra-large TC, resulting in a high but not indicative TPR. Moreover, the error caused by non-extraction of the extra-large TC in the AMETC results can be easily corrected during post-processing.

5. Discussion 5.1. Parameter tuning In application of the AMETC, multi-window, multi-scale and multi-direction image analyses were performed to handle complicated TCs of various widths, directions, and patterns. Nine simple parameters were utilized that can be categorized into three groups: (1) the size and number of windows used in the MNA; (2) the length L, width (r and X domain), and orientation h of the matched filter used in the GMF; and (3) the window size, k value, and minimum grain size (t) used in the TAT. The size and number of windows used in the MNA affected the separability of TCs of different widths in low-contrast DEMs. We evaluated how the window properties, including the window size and number, influenced the extraction accuracy (Fig. 11a) and computational complexity (Fig. 11b). The evaluations for a single window used 11 windows from 9 to 109 pixels in size, at intervals of 10-pixels. The evaluations for multiple windows used 11 window combinations (i.e., window sizes = [9], [9, 19], [9, 19, 29], . . ., [9, 19, 29. . ., 109]). We made the following observations: (1) with increasing window/window combinations, the increasing TPR based on a single- and multi-window strategy was high at first but gradually declined (Fig. 11a). (2) TCs of varying widths were enhanced more efficiently with multi-window rather than single-window strategies (Fig. 11a). For example, the TPR was 62.09% with only a single window (window size = 9), while the TPR for the two-window combination (window size = [9, 19]) increased to 79.32%. (3) Although the accuracy increased with larger window/window-combinations, the computing time increased exponentially. The median operation in the MNA is very sensitive to the size and number of windows. For instance, the computing time for multiple windows using the combination [9, 19, 29. . ., 99] was 644 s, and it increased to 829 s for the larger combination of [9, 19, 29. . ., 109] (Fig. 11b). The computing time cloud become prohibitively large if the multi-window strategy were used for a large dataset (e.g., 10,000  10,000 pixels). To balance accuracy

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

1015

Fig. 9. TC extraction results from Dataset A (2.5-m resolution): (a–d) are LiDAR DEMs of Block 1–4 (each 1000  1000 pixels, lighter = higher); (e–h) are reference maps of Block 1–4; (i–l) are extraction results from the AMETC; (m–p) from the OETM; (q–t) from the OCTM, and (u–x) from the D8.

and computing time, no more than ten windows were used in this study, and the maximum window size was limited to 100 pixels. This limitation makes it difficult to effectively enhance, detect,

and extract extra-large TCs (e.g., the TC > 800 m in Fig. 9d). However, extra-large TCs are not common and are easy to correct manually.

1016

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

Fig. 10. TC extraction results from Dataset B (5-m resolution): (a–f) are LiDAR DEMs of Block 6–11 (each 2000  2000 pixels, lighter = higher); (g–l) are extraction results by the AMETC.

The properties of the matched filter in the GMF, namely (1) the length L of the kernels, (2) the width of the kernels (determined by the scale r and the X domain), and (3) the number of orientations of the kernels, directly influence the extraction of small TCs. In previous research on retinal blood vessel extraction, a single scale was often adopted, and the L was usually fixed to obtain square filters. The X domain was usually set to 3r because more than 99% of the area under the Gaussian curve lies within the range [3r, 3r], and r was fixed at 2, based on experience (Zhang et al., 2010; Fraz et al., 2012; Li et al., 2012), allowing enhancement of a vessel less than

12 pixels wide. In contrast, TCs in LiDAR DEMs are typically more complicated. Widths vary greatly both upstream and downstream, and between main channels and branches. As mentioned earlier, we replicated the second derivative of the Gaussian function in the y-direction to generate a 2-D kernel, which means that the width of the targeted feature needs to be uniform to obtain a peak response. The width of a TC is reasonably consistent along very short segments, so L was fixed at 2 to obtain a 3-pixel length (Eq. (3)). The X domain is the half width of the kernel. In our study area, most TCs are less than 50 meters wide (i.e., 20 pixels in a DEM with

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

1017

Fig. 10 (continued)

a 2.5-m resolution), so r was set to no more than 4.5, and the X domain was set to 3r to construct a 2-D kernel (27  3) to enhance TCs less than 67.5 m wide (i.e., 27 pixels). Our experiments suggest that setting the number of orientations to 18 is sufficient. To handle the varying widths between main channels and small branches, a multi-scale approach was adopted. However, the larger the X domain is, the more orientations there are, and the more scales that are used in the GMF, the heavier calculation burden is. In this study, four scales, 1.5 6 r 6 4.5, with 1.0-r intervals, were used to enhance TCs of different widths.

Finally, the parameters in the TAT, i.e., the window size for local adaptive thresholding, the coefficient k, and the minimum grain size (t), significantly impacted the final extraction accuracy. (1) Window sizes from 6 to 106 pixels at a two-pixel interval were used in the TAT to determine the optimal parameters for local adaptive thresholding. The results show that TPR reaches its maximum only when the window size for the TAT matches the window size for the MNA (Fig. 12a). (2) Similarly, values of the coefficient k from 0.01 to 1.00 at intervals of 0.01, were tested. The results indicate that a higher k may identify fewer

1018

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

Table 1 Performance of the AMETC, D8, OETM, and OCTM methods, compared with manual digitization. Area

Method

Accuracy (%) EC

EO

OA

Block 1

AMETC D8 method OETM OCTM

97.63 89.97 31.40 36.59

1.45 – 0.31 4.67

4.12 18.60 0.87 13.02

2.37 10.03 68.60 63.41

98.31 – 81.68 79.84

Block 2

AMETC D8 method OETM OCTM

95.00 90.09 32.76 23.07

3.07 – 2.78 2.97

6.86 18.99 6.21 6.63

5.00 9.91 67.24 76.93

96.33 – 77.26 74.13

Block 3

AMETC D8 method OETM OCTM

97.03 94.79 0.89 22.55

1.55 – 0.51 1.65

6.59 49.43 2.18 7.03

2.97 5.21 99.11 77.45

98.18 – 80.76 83.95

Block 4

AMETC D8 method OETM OCTM

96.32 92.24 3.59 24.88

2.13 – 0.83 2.51

6.57 36.69 2.57 7.73

3.68 7.76 96.41 75.12

97.49 – 75.76 79.71

Block 5

AMETC D8 method OETM OCTM

40.93 98.15 86.97 9.68

2.50 – 27.67 2.41

4.34 23.70 48.08 4.19

59.07 1.85 13.03 90.32

76.83 – 77.68 65.48

Block 6

AMETC

96.41

2.06

6.43

3.59

97.55

TPR

FPR

Fig. 11. The relation between window properties and extraction accuracy/computation complexity (data are calculated from Block1): (a) The relation between TPR and the window sizes and numbers used in the MNA; (b) the relation between computing time and the window sizes and numbers.

Fig. 12. (a) The relation between k value and errors of commission and omission; (b) the relation between the TPR and the window size in the TAT.

TCs (i.e., higher omission error, indicated by red points in Fig. 12b). A lower k may identify many actual TCs but could cause false positives (i.e., higher commission error, indicated by yellow

points in Fig. 12b). The two curves intersect near k = 0.21. For convenience, k was set to 0.2. (3) The minimum grain size (t) was set to 30 pixels.

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

1019

Fig. 13. (a) A ravine landform DEM; (b) ravine extraction results by the D8 method and AMETC.

5.2. Transferability of the AMETC for terrestrial ravine network extraction Terrestrial ravines typically have linear depressions with V-shaped or U-shaped cross sections, which suggests that the AMETC might be useful for terrestrial ravine network extraction. We applied the AMETC to a ravine landform region in the northern Loess Plateau, Hengshan County, Shaanxi Province, China. In this area, the ravine density is as great as 16 km/km2. The tested DEM (Fig. 13a) was 459  485 pixels with a low resolution of 25  25 m. Four windows from 9 to 39 pixels at 10-pixels intervals were used in the multi-window MNA. Three scales from 1.5 to 2.5 at intervals of 0.5 were used in the GMF. The minimum grain size (t) in the TAT procedure was set to 10 pixels, and the other parameters were the same as those in the previous TC extraction. The preliminary experiments indicated that the AMETC could effectively extract ravines in this high-relief terrain from a low-resolution DEM (Fig. 13b). Visual interpretations show that the AMETC performed better than the D8 method, particularly on the upstream ravines (Fig. 13b). However, more applications are required to assess the feasibility and robustness of the AMETC on DEMs of different resolutions over different types of terrain.

thresholding (TAT). The AMETC was successfully applied to a LiDAR DEM covering more than 5000 km2, and the results showed that the proposed method is feasible for tidal creeks of different widths at different resolutions. Our quantitative assessments demonstrated that the AMETC generally outperformed conventional methods, especially for extraction of small tidal creeks. The parameters of the AMETC were sufficiently stable to ensure good automation and the method rarely required manual intervention when applied to the vast study area. Future work should address the following needs: (1) Validation of the AMETC using very-high-resolution LiDAR datasets or aerial images with resolutions up to 0.03 m. (2) Development of an effective algorithm for enhancing extra-large tidal creeks. Extra-large tidal creeks (>100 pixels) were not detected effectively in this study because of the heavy computational burden of the MNA when large windows. (3) Development of an algorithm for determining the thalwegs of tidal creeks (with single-pixel width). Producing single-pixel tidal creek networks is important to effectively representing key features (e.g., the width/depth ratio, network order, and watershed) that can give geomorphologists a better understanding of the network evolution. Acknowledgements

6. Conclusions To help geomorphologists better understand the characteristics, formation, and evolution of tidal creek networks, an automated method for extracting tidal creeks based on LiDAR DEMs is needed. The distinctive features of tidal creek networks such as low relief, variable widths, high density, strong anisotropy, and complicated drainage patterns present a challenge to accurate extraction. Using common features of tidal creeks, particularly their continuous linear traces and concave profiles, we developed an automated extraction method (AMTEC) that employs several image processing technologies, including multi-window median neighborhood analysis (MNA), multi-scale and multi-direction Gaussian-matched filtering (GMF), and two-stage adaptive

The MATLAB codes mentioned in this paper are available in the GitHub (https://github.com/yongxuenju/tidal-creeks-extraction). The LiDAR DEM dataset is owned by the Jiangsu Provincial Bureau of Surveying, Mapping and Geo-information (JPBSMG, http://www.jschj.gov.cn). The LiDAR DEM cannot be accessed online. However, it can be accessed, in accordance with the national spatial data sharing policy of China, by contacting JPBSMG ([email protected]). This research was supported by the National Natural Science Foundation of China (Nos. 41471068, 41171325, 41230751, and J1103408), the Program for New Century Excellent Talents in University (NCET-12-0264), the Fundamental Research Funds for the Central Universities, and the National Key Project of Scientific and Technical Supporting

1020

Y. Liu et al. / Journal of Hydrology 527 (2015) 1006–1020

Programs funded by the Ministry of Science & Technology of China (No. 2012BAH28B04). Any errors or shortcomings in the paper are the responsibility of the authors. References Allen, J.R.L., 2000. Morphodynamics of Holocene Salt marshes: a review sketch from the Atlantic and Southern North Sea coasts of Europe. Quatern. Sci. Rev. 19, 1155–1231. Baruch, A., Filin, S., 2011. Detection of gullies in roughly textured terrain using airborne laser scanning data. ISPRS J. Photogramm. Remote Sens. 66 (5), 564– 578. Brzank, A., Heipke, C., Goepfert, J., Soergel, U., 2008. Aspects of generating precise digital terrain models in the Wadden Sea from lidar-water classification and structure line extraction. ISPRS J. Photogramm. Remote Sens. 63 (5), 510–528. Cavalli, M., Tarolli, P., Marchi, L., Fontana, G.D., 2008. The effectiveness of airborne LiDAR data in the recognition of channel-bed morphology. CATENA 73, 249– 260. Cho, H.C., Slatton, K.C., Krekeler, C.R., Cheung, S., 2011. Morphology-based approaches for detecting stream channels from ALSM data. Int. J. Remote Sens. 32 (24), 9571–9597. Clubb, F.J., Mudd, S.M., Milodowski, D.T., Hurst, M.D., Slater, L.J., 2014. Objective extraction of channel heads from high-resolution topographic data. Water Resour. Res. 50 (5), 4283–4304. Coco, G., Zhou, Z., van Maanen, B., Olabarrieta, M., Tinoco, R., Townend, I., 2013. Morphodynamics of tidal networks: advances and challenges. Mar. Geol. 346, 1–16. D’Alpaos, A., Lanzoni, S., Marani, M., Fagherazzi, S., Rinaldo, A., 2005. Tidal network ontogeny: channel initiation and early development. J. Geophys. Res.-Earth Surface 110 (F2), Unsp f02001. Fagherazzi, S., Bortoluzzi, A., Dietrich, W.E., Adami, A., Lanzoni, S., Marani, M., Rinaldo, A., 1999. Tidal networks 1. Automatic network extraction and preliminary scaling features from digital terrain maps. Water Resour. Res. 35 (12), 3891–3904. Fairfield, J., Leymarie, P., 1991. Drainage networks from grid digital elevation models. Water Resour. Res. 27 (5), 709–717. Fawcett, T., 2006. An introduction to ROC analysis. Pattern Recogn. Lett. 27 (8), 861– 874. Fraz, M.M., Remagnino, P., Hoppe, A., Uyyanonvara, B., Rudnicka, A.R., Owen, C.G., Barman, S.A., 2012. Blood vessel segmentation methodologies in retinal images: a survey. Comput. Methods Programs Biomed. 108 (1), 407–433. Frazier, P.S., Page, K.J., 2000. Water body detection and delineation with Landsat TM data. Photogramm. Eng. Remote Sens. 66 (12), 1461–1467. Freeman, T.G., 1991. Calculating catchment-area with divergent flow based on a regular grid. Comput. Geosci. 17 (3), 413–422. French, J.R., Stoddart, D.R., 1992. Hydrodynamics of salt-marsh creek systems: implications for marsh morphological development and material exchange. Earth Surf. Proc. Land. 17 (3), 235–252. Gabet, E.J., 1998. Lateral migration and bank erosion in a saltmarsh tidal channel in San Francisco Bay, California. Estuaries 21 (4B), 745–753. Gong, Z., Wang, Z.B., Stive, M.J.F., Zhang, C.K., Chu, A., 2012. Process-based morphodynamic modeling of a schematized mudflat dominated by a longshore tidal current at the central Jiangsu coast, China. J. Coastal Res. 28 (6), 1381–1392. Hughes, Z.J., FitzGerald, D.M., Wilson, C.A., Pennings, S.C., Wieski, K., Mahadevan, A., 2009. Rapid headward erosion of marsh creeks in response to relative sea level rise. Geophys. Res. Lett. 36 (3), L03602. James, L.A., Watson, D.G., Hansen, W.F., 2007. Using LiDAR data to map gullies and headwater streams under forest canopy: South Carolina, USA. CATENA 71, 132– 144. Ji, L., Zhang, L., Wylie, B., 2009. Analysis of dynamic thresholds for the normalized difference water index. Photogramm. Eng. Remote Sens. 75 (11), 1307–1317. Kirwan, M.L., Mudd, S.M., 2012. Response of salt-marsh carbon accumulation to climate change. Nature 489, 550. Klemas, V., 2011. Beach profiling and LIDAR bathymetry: an overview with case studies. J. Coastal Res. 27 (6), 1019–1028.

Lanzoni, S., Seminara, G., 2002. Long-term evolution and morphodynamic equilibrium of tidal channels. J. Geophys. Res.-Oceans 107 (C1), 1–13. Li, Q., You, J., Zhang, D., 2012. Vessel segmentation and width estimation in retinal images using multiscale production of matched filter responses. Expert Syst. Appl. 39 (9), 7600–7610. Lohani, B., Mason, D.C., 2001. Application of airborne scanning laser altimetry to the study of tidal channel geomorphology. ISPRS J. Photogramm. Remote Sens. 56 (2), 100–120. Lohani, B., Mason, D.C., Scott, T.R., Sreenivas, B., 2006. Extraction of tidal channel networks from aerial photographs alone and combined with laser altimetry. Int. J. Remote Sens. 27 (1), 5–25. Lotze, H.K., Lenihan, H.S., Bourque, B.J., Bradbury, R.H., Cooke, R.G., Kay, M.C., Kidwell, S.M., Kirby, M.X., Peterson, C.H., Jackson, J.B.C., 2006. Depletion, degradation, and recovery potential of estuaries and coastal seas. Science 312 (5781), 1806–1809. Marani, M., D’Alpaos, A., Lanzoni, S., Carniello, L., Rinaldo, A., 2007. Biologicallycontrolled multiple equilibria of tidal landforms and the fate of the Venice lagoon. Geophys. Res. Lett. 34 (11), L11402. Mason, D.C., Scott, T.R., Wang, H.J., 2006. Extraction of tidal channel networks from airborne scanning laser altimetry. ISPRS J. Photogramm. Remote Sens. 61 (2), 67–83. Mudd, S.M., D’Alpaos, A., Morris, J.T., 2010. How does vegetation affect sedimentation on tidal marshes? Investigating particle capture and hydrodynamic controls on biologically mediated sedimentation. J. Geophys. Res.-Earth Surface 115, F03029. Ni, W.F., Wang, Y.P., Zou, X.Q., Zhang, J.C., Gao, J.H., 2014. Sediment dynamics in an offshore tidal channel in the southern Yellow Sea. Int. J. Sedim. Res. 29 (2), 246– 259. Ocallaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks from digital elevation data. Comput. Vis. Graph. Image Process. 28 (3), 323–344. Otsu, N., 1979. Threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern. 9 (1), 62–66. Pestrong, R., 1965. The Development of Drainage Patterns in Tidal Marshes, vol. 10. Stanford University Publications in Earth Science, pp. 1–87. Rinaldo, A., Fagherazzi, S., Lanzoni, S., Marani, M., Dietrich, W.E., 1999. Tidal networks 2. Watershed delineation and comparative network morphology. Water Resour. Res. 35 (12), 3905–3917. Rutzinger, M., Höfle, B., Pfeifer, T., Geist, T., Stötter, J., 2006. Object based analysis of airborne laser scanning data for natural hazard purposes using open source components. J. Photogramm. Remote Sens. 36 (Pt. 4/C42). Stumpf, A., Malet, J.P., Kerle, N., Niethammer, U., Rothmund, S., 2013. Image-based mapping of surface fissures for the investigation of landslide dynamics. Geomorphology 186, 12–27. Tarboton, D.G., 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resour. Res. 33 (2), 309– 319. Verpoorter, C., Kutser, T., Tranvik, L., 2012. Automated mapping of water bodies using Landsat multispectral data. Limnol. Oceanogr.-Methods 10, 1037–1050. Wang, W., Liu, H., Li, Y., Su, J., 2014. Development and management of land reclamation in China. Ocean Coast. Manage. 102, 415–425. Wang, Y.P., Gao, S., Jia, J.J., Thompson, C.E.L., Gao, J.H., Yang, Y., 2012. Sediment transport over an accretional intertidal flat with influences of reclamation, Jiangsu coast, China. Mar. Geol. 291, 147–161. Xing, F., Wang, Y.P., Wang, H.V., 2012. Tidal hydrodynamics and fine-grained sediment transport on the radial sand ridge system in the southern Yellow Sea. Mar. Geol. 291, 192–210. Zhang, B., Zhang, L., Zhang, L., Karray, F., 2010. Retinal vessel extraction by matched filter with first-order derivative of Gaussian. Comput. Biol. Med. 40 (4), 438– 445. Zhang, R.S., Shen, Y.M., Lu, L.Y., Yan, S.G., Wang, Y.H., Li, J.L., Zhang, Z.L., 2004. Formation of Spartina alterniflora salt marshes on the coast of Jiangsu Province, China. Ecol. Eng. 23 (2), 95–105. Zhao, S.S., Liu, Y.X., Li, M.C., Sun, C., Zhou, M.X., Zhang, H.X., 2015. Analysis of Jiangsu tidal flats reclamation from 1974 to 2012 using remote sensing. China Ocean Eng. 29 (1), 143–154.