Automated Extraction and Mapping for Desert Wadis from Landsat Imagery in Arid West Asia

Wadis, ephemeral dry rivers in arid desert regions that contain water in the rainy season, are often manifested as braided linear channels and are of vital importance for local hydrological environments and regional hydrological management. Conventional methods for effectively delineating wadis from heterogeneous backgrounds are limited for the following reasons: (1) the occurrence of numerous morphological irregularities which disqualify methods based on physical shape; (2) inconspicuous spectral contrast with backgrounds, resulting in frequent false alarms; and (3) the extreme complexity of wadi systems, with numerous tiny tributaries characterized by spectral anisotropy, resulting in a conflict between global and local accuracy. To overcome these difficulties, an automated method for extracting wadis (AMEW) from Landsat-8 Operational Land Imagery (OLI) was developed in order to take advantage of the complementarity between Water Indices (WIs), which is a technique of mathematically combining different bands to enhance water bodies and suppress backgrounds, and image processing technologies in the morphological field involving multi-scale Gaussian matched filtering and a local adaptive threshold segmentation. Evaluation of the AMEW was carried out in representative areas deliberately selected from Jordan, SW Arabian Peninsula in order to ensure a rigorous assessment. Experimental results indicate that the AMEW achieved considerably higher accuracy than other effective extraction methods in terms of visual inspection and statistical comparison, with an overall accuracy of up to 95.05% for the entire area. In addition, the AMEW (based on the New Water Index (NWI)) achieved higher accuracy than other methods (the maximum likelihood classifier and the support vector machine classifier) used for bulk wadi extraction.


Introduction
Arid regions are crucial areas for the study of global environmental change and sustainable development, and contain environment-specific habitats for both humans and natural species [1][2][3].However, desertification and water shortages caused by continuously decreasing precipitation in arid regions are a serious problem leading to a cascade of ecological and environmental problems [4][5][6][7].The term "wadi," in Arabic, refers to a valley-like landscape, commonly found in and around desert regions.They are ephemerally filled with water from heavy rainfall in the rainy season, resulting in a river that occurs intermittently in arid and semi-arid regions such as the Sahara Desert and the Arabian Peninsula.As one of the main river types in arid regions, the distribution and variation of wadis in these dry regions not only play a vital role in maintaining the health of the ecological environment and biological diversity, but they also have a significant impact on the sustainable development of human society and economy.In addition, considering the intermittently available sub-surface water contained within them, wadis are to some extent associated with the centers of local population [8].In addition, since the seasonal vegetation existing in wadis is regarded as a major food source in the desert regions, nomads and pastoralists rely on the distribution of wadis to design and locate their specific transhumance routes [9].Owing to the valuable water resources, they provide to sustain local human life, wadis therefore give rise to the distinctive scientific field of wadi hydrology, which attempts to explore the maximum potential for the utilization and planning of wadi systems [10,11].As a consequence, accurate wadi extraction is of great importance for the investigation of local water resources, construction of irrigation facilities, evaluation of geological engineering locations, and environmental protection.
Conventional in situ measurements are capable of providing an accurate mapping of wadis by means of field data; however, the method is limited to local coverage, and it does not provide timely, consistent, and systematic knowledge of wadis due to their short and erratic occurrence.The advent of remote sensing techniques has proved to be a useful solution for providing dynamic and synoptic monitoring of various landscapes by means of satellite imagery [12,13].Numerous research studies have attempted to extract various channel landscape features, e.g., river networks, which are morphologically similar to wadis, from satellite images [14,15].Common techniques for extraction of channel features can be primarily categorized into three main types: (i) a direct visual interpretation from original images by humans, a precise but time-and labor-consuming strategy which is restricted in terms of repetitive monitoring, large-scale wadi detection, as well as ongoing observation; (ii) a single threshold applied to spectral bands through spectral density-slicing on a single band or a synthesized band, i.e., the mathematical combination of several bands (despite its ease of implementation and extensive application [16,17], it is frequently constrained from a stable/optimal threshold perspective with respect to regionally diverse environments, which significantly limits a high level of automation in large regions [18]); and (iii) the same dilemma also occurs in the case of a classical remotely sensed method, i.e., the supervised classification method, a strategy in which relevant expert knowledge of satellite image characteristics is a prerequisite for further improving classification accuracy.In order to increase the level of automation and the extraction accuracy of channel-like features, further studies [19] have been carried out, and satisfactory results have been obtained.For example, an automatic method which used supervised classification based on vector data for creating training areas has been presented [20].Despite an acceptable degree of success, the poor availability of ancillary vector data limited its practical application.Other improvements [21][22][23] have focused on the combination of the classical classification and processing methods, represented by object-oriented analysis (OOA).The abundant spectral characteristics of landscapes provide a large amount of functional spectral indicators in the application of this type of method.Based on the guidance of image segmentation ideology, a classical supervised classification method, i.e., Support Vector Machine (SVM), was applied [21] in order to extract river features after preliminary enhancement of linear features by morphological processing.In addition, the OOA technique has recently been used extensively in the application of the gray-level co-occurrence matrix (GLCM) method [24].However, it is difficult to obtain satisfactory results in wadi regions using the aforementioned methods for the following reasons.First, because of the intermittent nature of seasonal flow activity, the wadis often exhibit morphological irregularities [11], leading to the failure of an OOA which is based on consistent morphological characteristics of a landscape feature.Subsequently, due to the absence of frequent inundation and the abundance of sediments derived from a large area, wadis are usually characterized by various braided stream patterns with a dense, crisscrossing drainage network, resulting in a more problematic extraction in comparison with regularly distributed channel-like features, such as rivers and tidal creeks [25,26].Furthermore, wadis, as seasonal rivers, appear erratically and their spectral characteristics are weakly developed and inconspicuous, especially in the dry season.In addition, since most rivers in wadi systems are very narrow (widths of less than 5 pixels) and shallow, the use of optical images to extract these tiny objects is likely to be affected by the problem of mixed pixels and adjacent pixels, resulting in inconsistent spectral discrimination between wadis and heterogeneous backgrounds.As a result, the successful extraction and mapping of wadis and sequential analyses associated with them necessitates a novel strategy which not only takes full advantage of the wadis' morphological features, but which also improves the spectral separability of the water in wadis and the adjacent pixels.
To overcome the aforementioned difficulties, we developed an automated method for extracting wadis (AMEW) (Table A1) from Landsat imagery, a series of freely available medium spatial resolution satellite images, which is capable of performing observation and monitoring of long time series.By systematically synthesizing several common and effective Water Indices (WIs), we developed a novel water index to effectively suppress the spectral mixing issue and to maximize the spectral separability between the water in wadis and the adjacent pixels.Subsequently, a multi-scale Gaussian matched filtering and local adaptive threshold segmentation strategy was constructed to automatically extract and map the wadi river distribution with satisfactory accuracy, with which the generalizability and transferability of the proposed method was also considered.

Study Area
Jordan is located in the northwestern part of the Arabian Peninsula.Much of the country has a sub-tropical desert climate with an average temperature of 30 ˝C in summer and 13 ˝C in winter.The annual rainfall is approximately 100 mm, most of which occurs in winter [27].Due to the extreme shortage of water resources and serious desertification, the desert area of Jordan exceeds 80% of the national territorial area.The abrupt but infrequent heavy rainfall events in the desert environment result in a well-developed system of wadi landforms.The drainage patterns in the study area are extremely complex, and the river networks are fairly dense, making it highly suited for extraction and mapping; moreover, the products of the research are potentially of significant value for local decision makers in performing hydrological constructions and planning.In order to validate the effectiveness and robustness of the proposed AMEW, we carefully selected three regions from the southeastern part of Jordan as the test sites, with each study site covering an area of 3600 km 2 , in order to ensure the stability of the AMEW.Study Area 1 is located in the middle of Ma'am Governorate; Study Area 2 is in the south of Ma'am Governorate; and Study Area 3 is located in center of Jordan, crossed by Azraq, Jizeh, and Ma'am Governorates (Figure 1).Given that the geomorphological features of these complex landscapes are spatially diverse, the three test sites were selected in accordance with the variation of latitude and location-specific characteristics of wadis in order to ensure both a comprehensive and systematic analysis.

Satellite Data
Satellite images of varying spatial and temporal resolution have proven their outstanding effectiveness in the detection and mapping of changes in land cover [28].High temporal resolution datasets (e.g., moderate-resolution imaging spectroradiometer (MODIS) data) have been extensively applied to land cover mapping and change detection in a short repetition cycle [29].However, the coarse spatial resolution (generally larger than 500 m) has limited their practical application at a regional scale for tiny objects such as wadis.Conversely, a widespread application of high-resolution images for monitoring is restricted due to their limited spatial coverage and high expense.Landsat satellite series data, with its free accessibility, appropriate spatial resolution (30 m for optical bands), global data coverage, and long temporal continuity (over 40 years), is an ideal data resource for delineating the distribution and analyzing the spatial variation of wadi rivers [30,31].The Landsat-8 Operational Land Imager (OLI) launched recently has already provided data free to users, with increased spatial coverage and improvements with respect to signal-to-noise ratio and radiometric performance.Thus, the Landsat-8 OLI imagery was chosen in the present study to conduct the relevant experiments and to extend the application of Landsat series images in the features extraction field.Related Landsat 8 OLI images processed by Level 1T Product Generation System (LPGS) were downloaded from USGS EarthExplorer [32].
shortage of water resources and serious desertification, the desert area of Jordan exceeds 80% of the national territorial area.The abrupt but infrequent heavy rainfall events in the desert environment result in a well-developed system of wadi landforms.The drainage patterns in the study area are extremely complex, and the river networks are fairly dense, making it highly suited for extraction and mapping; moreover, the products of the research are potentially of significant value for local decision makers in performing hydrological constructions and planning.In order to validate the effectiveness and robustness of the proposed AMEW, we carefully selected three regions from the southeastern part of Jordan as the test sites, with each study site covering an area of 3600 km 2 , in order to ensure the stability of the AMEW.Study Area 1 is located in the middle of Ma'am Governorate; Study Area 2 is in the south of Ma'am Governorate; and Study Area 3 is located in center of Jordan, crossed by Azraq, Jizeh, and Ma'am Governorates (Figure 1).Given that the geomorphological features of these complex landscapes are spatially diverse, the three test sites were selected in accordance with the variation of latitude and location-specific characteristics of wadis in order to ensure both a comprehensive and systematic analysis.In order to avoid the influence of prevalent clouds, we selected two images, obtained in September 2014, as the experimental data (path/row: 173/38 and 173/39).The two selected images, stored as GeoTIFF format and Universal Transverse Mercator (UTM) map projections, were all geo-referenced within a precision of 0.4 pixels.Necessary prerequisite processing was carried out via ENVI 5.1 [33] to determine the real reflectance of ground objects and to maintain the transferability of the methodology to other images [18].Radiometric calibration consisted of two major steps: (i) Top of the atmosphere (TOA) radiances R λ (W/(m 2 ¨sr¨µm)) were derived from a linear transformation process (Equation ( 1)).Corresponding calibration coefficients (i.e., the gain factor G λ and the offset O λ ) were provided in the LPGS metadata file; (ii) Subsequently, TOA reflectance ρ TOA was calculated from available provided calibration coefficients (Equation (2)).
Remote Sens. 2016, 8, 246 5 of 23 Subsequently, the atmospheric correction was conducted on the obtained ρ TOA via the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) module in ENVI v.5.1 [34], the coefficients of which were obtained from a look-up table on the official website.

Reference Data
The reference maps, i.e., the "true value" in accuracy evaluation, were manually digitized onscreen from the original Landsat images based on visual interpretation by three experts.Final reference maps were synthesized from the three digitization results, and the high-resolution Google Earth images were used as an additional reference to help distinguish inconspicuous wadis from the background.Given the inevitable time span between Landsat OLI data with Google Earth data, human expertise was necessarily involved to prevent false judgments and unqualified reference maps.However, given the fact that there were too many branches to be extracted manually, errors may be inevitable.Therefore, different experienced interpreters were employed as a cross validation to delineate the outlines of wadis to minimize the omission of small channels.Note that images from Google Earth were only used as a reference for ambiguous regions.

Methods
As mentioned in the Introduction, the extraction of wadis in the current study required a systematic improvement in both spectral separability and morphological utility.Therefore, the AMEW proposed in this study involves three main steps (Figure 2): (i) calculation of various general WIs and construction of a new WI (i.e., various WIs originally proposed to enhance water features from heterogeneous backgrounds were synthesized to improve the separation of the spectral features of wadis); (ii) Gaussian matched filtering to enhance the channel-like wadis; (iii) local adaptive threshold segmentation to minimize the manual interventions and to maximize the automation of the proposed method.
Remote Sens. 2016, 8, 246 5 of 23 maps.However, given the fact that there were too many branches to be extracted manually, errors may be inevitable.Therefore, different experienced interpreters were employed as a cross validation to delineate the outlines of wadis to minimize the omission of small channels.Note that images from Google Earth were only used as a reference for ambiguous regions.

Methods
As mentioned in the Introduction, the extraction of wadis in the current study required a systematic improvement in both spectral separability and morphological utility.Therefore, the AMEW proposed in this study involves three main steps (Figure 2): (i) calculation of various general WIs and construction of a new WI (i.e., various WIs originally proposed to enhance water features from heterogeneous backgrounds were synthesized to improve the separation of the spectral features of wadis); (ii) Gaussian matched filtering to enhance the channel-like wadis; (iii) local adaptive threshold segmentation to minimize the manual interventions and to maximize the automation of the proposed method.

Water Index Technique to Enhance the Spectral Separability of Wadi Water Bodies
The WI techniques were originally designed to enhance the discrimination between water bodies and land.The WI was derived from a mathematical calculation between several optical bands of optical satellite images, which was aimed to enhance the distinction of the spectral signals for targeted

Water Index Technique to Enhance the Spectral Separability of Wadi Water Bodies
The WI techniques were originally designed to enhance the discrimination between water bodies and land.The WI was derived from a mathematical calculation between several optical bands of optical satellite images, which was aimed to enhance the distinction of the spectral signals for targeted water bodies.However, the different mathematical formulations make them inherently sensitive to certain types of noise [35].This study investigated four commonly utilized WIs, including the Normalized Difference Water Index (NDWI) [36], the Modified Normalized Difference Water Index (MNDWI) [37], and the Automated Water Extraction Index (AWEI sh and AWEI nsh ) [38].Calculated by the normalized ratio of green and near infrared (NIR) bands, the NDWI can effectively remove vegetation and soil information, but it may be unsuitable for use in certain environments such as built-up areas, deserts, and ice sheet surfaces [39].The MNDWI, refined from the NDWI, can be better adjusted to various landscapes by substituting the mid-infrared band (SWIR1 of OLI images) for the near-infrared band.The AWEI sh (with shadow) and the AWEI nsh (no shadow), derived from a band combination of blue, green, NIR, and short-wave infrared (SWIR) bands, can not only enhance water information, but it can also reduce errors caused by shadows and dark regions.These four WIs have been extensively applied in recent years to address the water extraction problem (Table 1).Because the WIs have dissimilar enhancement effects in different regions due to the aforementioned intrinsic water characteristics and land cover complexity, we calculated the four WIs (Figure 3a-e) in order to comprehensively assess their reflectance separability and practical performance for the test sites.For convenient inspection and comparison, the pixel value of WI images was normalized to a range from ´1 and 1. Wadis represented low values while the backgrounds were characterized by higher values.In addition, 800 samples, representing wadis, bare land, vegetation, and built-up land, were carefully selected from the four WI images of three study areas to systematically compare the discrimination between wadis and their heterogeneous environmental backgrounds (Figure 3a-d).The more the value in the box plots deviates from the center (middle line), the less the reflectance is concentrated.In reality, the first and third quartile of reflectance can be used to demonstrate the most concentrated reflectance of wadis with the other three landscapes.The box plots enable us to visually assess how the degree of overlapping reflectance varies between the four landscape types.To quantify the spectral separability between the wadis and the other three types, we have calculated the Jeffries-Matusita distance (JMD) which quantifies the average distance between the density functions of two compared land cover types i and j in the multidimensional spectral feature space (Equation (3)) [40].
where B is the Bhattacharyya distance (Equation ( 4)) between landscape i and j, m i and m j are their average reflectance, and ř i and ř j are corresponding covariance matrices, respectively.
In reality, the JMD ranges from 0 to 2, the value of which corresponds to spectrally identical classes and perfect spectrally separable classes.Separability was considered to be excellent for JMD values were larger than 1.9, favorable when JMD values were between 1.4 and 1.9, and poor when they were lower than 1.4.From both visual inspection and the quantified results (Figure 3), we could determine that the New Water Index (NWI) improves the reflectance separability between wadis and the others, followed by the MNDWI and the AWEI combined with two efficient indexes with respect to different situations.Due to the severe desertification, sparseness of vegetation, and general absence of buildings in the study areas, the fundamental issue for wadi extraction is to distinguish water from bare land.Figure 4a illustrates the fact that wadis in NDWI images are often indistinct from bare land.Wadi mainstreams can be significantly enhanced by the MNDWI, whereas the AWEIsh and the AWEInsh images are better able to distinguish several small branches missed by the MNDWI.Figure 4a,d demonstrates that wadis in the AWEIsh and the AWEInsh images are more clearly distinguished from bare land.The MNDWI was less effective and the NDWI performed the worst.To find an optimum WI which can effectively enhance both main streams and tributaries, we developed an NWI using a compound procedure.Firstly, we produced a new composite color image with the four WI images, where the MNDWI was applied as the red band, the AWEIsh as the green band, and the AWEInsh as the blue band.Subsequently, we applied the HIS (Hue, Intensity, and Saturation combination) transformation [41] to the new color image and then used the I component as the NWI (Figure 4f).After a careful sampling process (Figure 4g), we found that the NWI, taking advantage of the merits of all the three WIs, not only enhanced the main streams and branches at the same time, but it also enhanced the discrimination between wadis and adjacent objects.To further determine the optimal WI for wadi extraction, we also used the MNDWI, the AWEIsh, the AWEInsh, and the NWI for extraction.From preliminary experiments and previous relevant research [18,38], the performance of the NDWI, which is capable of discriminating water from vegetation, was found to be inconsistent when applied to complex and heterogeneous aquatic environments.In order to enhance the transferability and generalizability of the proposed method under other circumstances, the NDWI was previously excluded.Due to the severe desertification, sparseness of vegetation, and general absence of buildings in the study areas, the fundamental issue for wadi extraction is to distinguish water from bare land.Figure 4a illustrates the fact that wadis in NDWI images are often indistinct from bare land.Wadi mainstreams can be significantly enhanced by the MNDWI, whereas the AWEI sh and the AWEI nsh images are better able to distinguish several small branches missed by the MNDWI.Figure 4a,d demonstrates that wadis in the AWEI sh and the AWEI nsh images are more clearly distinguished from bare land.The MNDWI was less effective and the NDWI performed the worst.To find an optimum WI which can effectively enhance both main streams and tributaries, we developed an NWI using a compound procedure.Firstly, we produced a new composite color image with the four WI images, where the MNDWI was applied as the red band, the AWEI sh as the green band, and the AWEI nsh as the blue band.Subsequently, we applied the HIS (Hue, Intensity, and Saturation combination) transformation [41] to the new color image and then used the I component as the NWI (Figure 4f).After a careful sampling process (Figure 4g), we found that the NWI, taking advantage of the merits of all the three WIs, not only enhanced the main streams and branches at the same time, but it also enhanced the discrimination between wadis and adjacent objects.To further determine the optimal WI for wadi extraction, we also used the MNDWI, the AWEI sh , the AWEI nsh , and the NWI for extraction.From preliminary experiments and previous relevant research [18,38], the performance of the NDWI, which is capable of discriminating water from vegetation, was found to be inconsistent when applied to complex and heterogeneous aquatic environments.In order to enhance the transferability and generalizability of the proposed method under other circumstances, the NDWI was previously excluded.

Wadi Enhancement Based on Gaussian Matched Filtering
Rivers of sufficient width, specifically the main streams, can be successfully delineated from satellite imagery using a simple threshold method with satisfactory accuracy, after improvement of their spectral contrast with backgrounds by means of the WI technique.However, most of the rivers in wadi systems are narrow and in the form of tiny linear features with strong anisotropy.Moreover, the morphological features of river tributaries, such as orientation and concavities, vary continuously along the mainstreams, while their widths also vary at the same time.These characteristics mean that further enhancement and subsequent extraction of these tiny irregular features should be based on their specific physical morphologies rather than on inconspicuous spectral features alone.Matched filtering technology in image processing fields is commonly used to enhance linear features such as retinal blood vessels, road networks, and gullies and fissures, given that their cross-sections generally appear as Gaussian-shaped curves [42,43].In many cases, the successful application of this image processing approach to extracting terrestrial targets is hindered by massive noise and false alarms in the background.Spectral representation of much noise can be improved with enhancements of the targeted objects, and noise can be compensated for by the pre-processing of the WI technique (Figure 5a).Considering the aforementioned similarity in morphological characteristics between wadis and linear features, and the complementarity of the method with the WI technique, we adopted the use of Gaussian matched filtering for identifying wadis.In addition, after the neighborhood analysis processed in the original NWI, different wadis were found to be more similar in terms of morphological characteristics, which facilitated the subsequent parameter tuning process for the Gaussian matched filtering.

Wadi Enhancement Based on Gaussian Matched Filtering
Rivers of sufficient width, specifically the main streams, can be successfully delineated from satellite imagery using a simple threshold method with satisfactory accuracy, after improvement of their spectral contrast with backgrounds by means of the WI technique.However, most of the rivers in wadi systems are narrow and in the form of tiny linear features with strong anisotropy.Moreover, the morphological features of river tributaries, such as orientation and concavities, vary continuously along the mainstreams, while their widths also vary at the same time.These characteristics mean that further enhancement and subsequent extraction of these tiny irregular features should be based on their specific physical morphologies rather than on inconspicuous spectral features alone.Matched filtering technology in image processing fields is commonly used to enhance linear features such as retinal blood vessels, road networks, and gullies and fissures, given that their cross-sections generally appear as Gaussian-shaped curves [42,43].In many cases, the successful application of this image processing approach to extracting terrestrial targets is hindered by massive noise and false alarms in the background.Spectral representation of much noise can be improved with enhancements of the targeted objects, and noise can be compensated for by the pre-processing of the WI technique (Figure 5a).Considering the aforementioned similarity in morphological characteristics between wadis and linear features, and the complementarity of the method with the WI technique, we adopted the use of Gaussian matched filtering for identifying wadis.In addition, after the neighborhood analysis processed in the original NWI, different wadis were found to be more similar in terms of morphological characteristics, which facilitated the subsequent parameter tuning process for the Gaussian matched filtering.It is very difficult to generate satisfactory results based on original WI images due to the inconspicuous and weak spectral signals of wadis.Therefore, we first utilized a neighborhood analysis and median filtering prior to the Gaussian matched filtering in order to transform the absolute spectral difference to a relative one, thereby enhancing the contrast between wadis and their neighboring backgrounds.The response of neighborhood analysis is strongly affected by the size of the utilized windows.Narrow wadis can be effectively enhanced by small windows, while wide wadis may be omitted, and vice versa.Hence, a neighborhood analysis with progressively increasing window sizes was designed and utilized to achieve efficient enhancement of wadis with varying widths.Changing the size of windows is an important strategy for optimum matching of objects of different scale and has been used extensively in calculating linear object characteristics such as width and curvature [44,45].Equation ( 4) is the expression for neighborhood analysis.F(x, y) is the WI image, Median (x, y) is the result of median filtering, and N(x, y) is the result of a neighborhood analysis.Figure 6 illustrates the results of multi-window neighborhood analysis, from which we were able to determine the effectiveness of the progressively increasing window sizes in relation to different wadis.

N(x, y) = Median(x, y) -F(x, y)
(5) After the improvement provided by the neighborhood analysis, the wadi cross sections appear as more distinct Gaussian-shaped curves (Figure 5b); therefore, the use of Gaussian matched filtering It is very difficult to generate satisfactory results based on original WI images due to the inconspicuous and weak spectral signals of wadis.Therefore, we first utilized a neighborhood analysis and median filtering prior to the Gaussian matched filtering in order to transform the absolute spectral difference to a relative one, thereby enhancing the contrast between wadis and their neighboring backgrounds.The response of neighborhood analysis is strongly affected by the size of the utilized windows.Narrow wadis can be effectively enhanced by small windows, while wide wadis may be omitted, and vice versa.Hence, a neighborhood analysis with progressively increasing window sizes was designed and utilized to achieve efficient enhancement of wadis with varying widths.Changing the size of windows is an important strategy for optimum matching of objects of different scale and has been used extensively in calculating linear object characteristics such as width and curvature [44,45].Equation ( 4) is the expression for neighborhood analysis.F(x, y) is the WI image, Median (x, y) is the result of median filtering, and N(x, y) is the result of a neighborhood analysis.Figure 6 illustrates the results of multi-window neighborhood analysis, from which we were able to determine the effectiveness of the progressively increasing window sizes in relation to different wadis.
Npx, yq " Medianpx, yq ´Fpx, yq It is very difficult to generate satisfactory results based on original WI images due to the inconspicuous and weak spectral signals of wadis.Therefore, we first utilized a neighborhood analysis and median filtering prior to the Gaussian matched filtering in order to transform the absolute spectral difference to a relative one, thereby enhancing the contrast between wadis and their neighboring backgrounds.The response of neighborhood analysis is strongly affected by the size of the utilized windows.Narrow wadis can be effectively enhanced by small windows, while wide wadis may be omitted, and vice versa.Hence, a neighborhood analysis with progressively increasing window sizes was designed and utilized to achieve efficient enhancement of wadis with varying widths.Changing the size of windows is an important strategy for optimum matching of objects of different scale and has been used extensively in calculating linear object characteristics such as width and curvature [44,45].Equation ( 4) is the expression for neighborhood analysis.F(x, y) is the WI image, Median (x, y) is the result of median filtering, and N(x, y) is the result of a neighborhood analysis.Figure 6 illustrates the results of multi-window neighborhood analysis, from which we were able to determine the effectiveness of the progressively increasing window sizes in relation to different wadis.

N(x, y) = Median(x, y) -F(x, y)
(5) After the improvement provided by the neighborhood analysis, the wadi cross sections appear as more distinct Gaussian-shaped curves (Figure 5b); therefore, the use of Gaussian matched filtering After the improvement provided by the neighborhood analysis, the wadi cross sections appear as more distinct Gaussian-shaped curves (Figure 5b); therefore, the use of Gaussian matched filtering to further enhance the wadis is effective and reliable.The 2-D Gaussian matched filtering function (Equation ( 6)) used in this paper is generated by repeating the 1-D Gaussian function in a neighborhood along the y-axis (Equation ( 5)), and the 1-D function is defined along the x-axis [46].
In Equations ( 5) and ( 6), σ represents the scale of the filtering corresponding to the width of wadis; L is the length of the neighborhood along the y-axis to smooth noise, and M(x, y) is the filtering result.
Since the orientations of wadis are unknown in advance, the kernel of Gaussian matched filtering must cover every orientation in order to capture the wadis.Subsequently, the responses from every orientation were comprehensively compared, and the orientation with the strongest response was considered as the correct orientation of the wadis (Equation ( 7)) [25].
where b represents the filtering process and R(x, y) is the result of matched filtering.The kernel of M (x, y) was rotated according to θ with an angle interval of 10 ˝to construct 18 different kernels that spanned all possible orientations.In addition, only if σ matches the width of wadis can they provide the strongest response.However, the entire wadi system is highly complex, and a single filtering scale can only enhance wadis of the same width.In order to ensure that wadis of all widths were effectively enhanced, we utilized multi-scale Gaussian matched filtering.As σ increases, all of the wadis are enhanced.

Local Adaptive Threshold Segmentation
Threshold segmentation is an essential final step in extracting wadis from the background following the application of Gaussian matched filtering.According to the literature, there are a large number of image threshold segmentation methods, typically including histogram shape-based methods, clustering-based methods, and entropy-based methods [47][48][49].All of these methods are based on a global threshold in regard to the whole segmentation object.Since Gaussian matched filtering clearly results in a unimodal distribution, a global optimum threshold cannot accommodate both large rivers and very small rivers at a same time, since it would lead to many omissions of regional tiny rivers, especially.Therefore, in the present study the local adaptive threshold segmentation method was used to delineate wadis of widely varying size.

Narrow Wadi Extraction
Firstly, we built a window of wˆw pixels for the Gaussian matched filtering results derived from the preceding step, i.e., R(x, y), and calculated the local threshold T(x, y) according to Equation (8).In this equation, µ(x, y) is the local mean value for a specific analysis window while ε(x, y) is the local standard deviation.During this process, pixels with values greater than a local threshold (defined below) would be classified as wadis and others would be classified as backgrounds, and vice versa (Equation ( 9)).Subsequently, the process was completed by sliding the window over the entire image to ensure that each of the local narrow wadis could be extracted with an appropriate threshold.
T px, yq " µ px, yq `k ˆε px, yq W 1 px, yq " t In principle, the aforementioned local threshold was determined by the local mean and standard deviation values in a neighboring window; however, this threshold varied considerably and was likely to be influenced by the local maximum and minimum values.Therefore, wadis with pixel values close to the heterogeneous background were likely to be misclassified, producing incomplete and discontinuous results.To suppress the local effect caused by an extremum value, we replaced the wadi pixels detected in the first round of threshold segmentation with the local mean values, and then repeated the threshold segmentation.Finally, we combined the first round results (w 1 ) with the second round results (w 2 ) using a mathematical union (Equation ( 10)) process to yield the complete assembly of narrow wadis.

Wide Wadi Extraction
In the case of wide wadis, the 2-D Gaussian kernel could only provide peak responses close to the edges, leaving the central parts unenhanced (Figure 7b) and resulting in further misclassifications.However, applying local adaptive threshold segmentation directly to the neighborhood analysis outcomes can achieve a satisfactory performance.Therefore, for wide wadis, the outcomes from a large-window neighborhood analysis were used for local threshold segmentation instead of a second round by mean value.Note that, in order to guarantee the completeness of the final results, we not only fused the results of both wide and narrow wadis, but also merged the outcomes from different window sizes.likely to be influenced by the local maximum and minimum values.Therefore, wadis with pixel values close to the heterogeneous background were likely to be misclassified, producing incomplete and discontinuous results.To suppress the local effect caused by an extremum value, we replaced the wadi pixels detected in the first round of threshold segmentation with the local mean values, and then repeated the threshold segmentation.Finally, we combined the first round results (w1) with the second round results (w2) using a mathematical union (Equation ( 10)) process to yield the complete assembly of narrow wadis.

Wide Wadi Extraction
In the case of wide wadis, the 2-D Gaussian kernel could only provide peak responses close to the edges, leaving the central parts unenhanced (Figure 7b) and resulting in further misclassifications.However, applying local adaptive threshold segmentation directly to the neighborhood analysis outcomes can achieve a satisfactory performance.Therefore, for wide wadis, the outcomes from a large-window neighborhood analysis were used for local threshold segmentation instead of a second round by mean value.Note that, in order to guarantee the completeness of the final results, we not only fused the results of both wide and narrow wadis, but also merged the outcomes from different window sizes.

Results
The automated method for extracting wadis (AMEW) used in this study was implemented under Matlab (R2012b 8.0.0783).The window size of neighborhood analysis varied from 8 to 80 pixels with an interval of 10 pixels.This window size interval was decided by a balance between accuracy and computing time and by the specific extent of the test sites, which could be changed manually according to the practical situation.The scales of Gaussian matched filtering varied from 1.0 to 1.4 in 0.2 intervals.There were three parameters in the local adaptive threshold segmentation procedure: w was the window size of local calculation, which varied according to the neighborhood analysis; k1 was the parameter for narrow wadi threshold segmentation; and k2 was the parameter for wide wadi threshold segmentation.For all three regions with four different WIs, k1 was set to a value range between 0.05 to 0.1 in terms of diverse aquatic environments, with a mean value of 0.07 and a subtle standard deviation of 0.02; and k2 was set to a value range between 0.10 to 0.35, with a mean value of 0.18 and a standard deviation of 0.08.To evaluate the results of each wadi extraction method, we superimposed the extraction results upon the reference map, pixel by pixel, and made a comparison with the reference map, which was achieved by a laborious manual digitization procedure discussed in Section 2.2 in order to determine the omissions and commission errors.There were four possibilities in the comparison results: true positive (TP), true negative (TN), false negative (FN), and false positive (FN) [50].Five metrics were used to quantitatively assess the extraction results: TPR (the ratio of correct wadi results to all true wadis: TPR = TP/(TP + FN)); FPR (the ratio of erroneous wadi results to the background of the reference map: FPR = FP/(FP + TN)); EC (the ratio of erroneous

Results
The automated method for extracting wadis (AMEW) used in this study was implemented under Matlab (R2012b 8.0.0783).The window size of neighborhood analysis varied from 8 to 80 pixels with an interval of 10 pixels.This window size interval was decided by a balance between accuracy and computing time and by the specific extent of the test sites, which could be changed manually according to the practical situation.The scales of Gaussian matched filtering varied from 1.0 to 1.4 in 0.2 intervals.There were three parameters in the local adaptive threshold segmentation procedure: w was the window size of local calculation, which varied according to the neighborhood analysis; k 1 was the parameter for narrow wadi threshold segmentation; and k 2 was the parameter for wide wadi threshold segmentation.For all three regions with four different WIs, k 1 was set to a value range between 0.05 to 0.1 in terms of diverse aquatic environments, with a mean value of 0.07 and a subtle standard deviation of 0.02; and k 2 was set to a value range between 0.10 to 0.35, with a mean value of 0.18 and a standard deviation of 0.08.To evaluate the results of each wadi extraction method, we superimposed the extraction results upon the reference map, pixel by pixel, and made a comparison with the reference map, which was achieved by a laborious manual digitization procedure discussed in Section 2.2 in order to determine the omissions and commission errors.There were four possibilities in the comparison results: true positive (TP), true negative (TN), false negative (FN), and false positive (FN) [50].Five metrics were used to quantitatively assess the extraction results: TPR (the ratio of correct wadi results to all true wadis: TPR = TP/(TP + FN)); FPR (the ratio of erroneous wadi results to the background of the reference map: FPR = FP/(FP + TN)); EC (the ratio of erroneous wadi results to all true wadis: EC = FP/(TP + FN)); EO (the ratio of omissions to all true wadis: EO = FN/(TP + FN)); and OA (the ratio of correctly classified pixels to the total number of pixels: OA = (TP +TN)/N).The TPR represents the pixels correctly classified as wadi pixels, while the FPR is the fraction of pixels incorrectly identified as wadi.EO shows how many wadi pixels were omitted by the methodology, and the EC indicates the misclassified backgrounds.In fact, the overall accuracy represents the entire extraction results and should be regarded as a generalized standard to evaluate the methodology; the TPR represents the reliability of the extracted wadis; and the FPR demonstrates the performance of a method to suppress the influence from backgrounds.The OA reflects the global performance of the method, with high OA indicating high accuracy.The practical significance of the TPR and FPR should be assessed according to the specific water environments.

Comparison of Various WIs and the Newly Synthesized WI
As discussed in Section 3.1, a comprehensive comparison of different original WIs was first carried out to explore the practical effect accomplished with respect to the complex wadi system.Preliminary results obtained by four water indexes, i.e., the MNDWI, the AWEI sh , the AWEI nsh , and the NWI, are illustrated in Figure 8 for three representative areas.Each study area covered a vast area (over than 2000 ˆ2000 pixels) to ensure the comprehensiveness of the method assessment.Thresholds in this step were determined by an iterative experiment to find the optimal threshold for each WI.Due to the high degree of complexity and high density of the entire wadi system, six blocks with typical morphological features were selected to illustrate the performance as well as the robustness of different WIs with clear perspectives.Note that for clarity only the results for two blocks are illustrated in Figure 8; further illustrations of the remaining four blocks with the same expression are given in Figures S1 and S2 in the supplementary online materials.
Figure 8 reveals that the extraction effects of four WIs are visually dissimilar.From visual inspection, we found that both the false negative rate and the false positive rate of the proposed NWI were significantly less than the other three, demonstrating the considerable improvement provided by the NWI.The AWEI nsh achieved poorer results compared to the NWI, followed by the AWEI sh and the MNDWI, the errors for which were located in tiny wadi tributaries as expected.For wadi extraction in these two blocks, the false positive rate of the AWEI nsh was relatively low, benefiting from the intrinsic robust performance of this AWEI with regard to aquatic environments without a significant amount of shadow.The enhancement of main streams by the AWEI nsh and the AWEI sh were relatively poor, which resulted in many omissions (refer to Block 4, Block 5 in the supplementary online materials).With regard to the MNDWI, the main streams were well enhanced, whereas many narrow wadis were ignored, resulting in many omissions occurring around narrow wadis.In addition, the drainage density of narrow wadis was generally high, and there were numerous mixed-pixels in their vicinity, significantly increasing the difficulty of the extraction by numerous commissions (see Block 6 in the supplementary online materials).Table 2 lists the statistical accuracy assessments of the four WIs in the three study areas.The results demonstrate that the NWI performed best among the four WIs, with the highest TPR (88.40%-94.29%)and lowest FPR (2.53%-4.50%);OA, EC, and EO varied from 93.26%-95.62%,3.63%-10.60%,and 5.71%-11.60%,respectively.The performance of the three other WIs varied over different study areas.In Study Area 1, the performance of the AWEI nsh was better according to its relatively high TPR (79.87%-87.03%)and low FPR (4.7%%-6.16%).the AWEI sh performed worst on the basis of its low TPR (about 70%) and high FPR (11.08%-11.26%),which means that most wadis were omitted, and certain backgrounds were unexpectedly identified as wadis.In Study Area 2, the MNDWI performed worst.Its TPR only reached about 67%, and its EC (24.69%-27.58%)and EO (29.84%-37.63%)were significantly higher than the others.The AWEI nsh and the AWEI sh performed slightly better with their OA, varying from 85.36%-86.5% and 86.5%-83.91%.In Study Area 3, the MNDWI performed better, with TPR and FPR of 77.23% and 9.69%, respectively, compared with 74.32% and 10.7%, and 63.6% and 13.46% for other WIs.According to the assessments, the NWI outperformed other WIs based on both visual inspections as well as on statistical quantification.Therefore, it was selected as the optimal WI and was used for pre-processing for the subsequent wadi extraction work.Hereafter, the AMEW refers to the method combined with the NWI.

Comparison of the AMEW and Other Conventional Extraction Methods
To further validate the robustness and adaptability of the AMEW, the maximum likelihood (ML) classifier and the SVM classifier were compared with the AMEW.Supervised classifiers were used as comparison groups mainly because of their high classification accuracy compared to unsupervised classifiers [51].Several previous studies have used this type of method for a reference and comparison with proposed methods [18,43].Among various supervised methods, the ML classifier is one of the most representative parameter classifiers, while the SVM classifier is one of the most representative non-parameter classifiers.We selected a substantial number of training samples to describe the spectral characteristics of each class in order to ensure the effectiveness of the supervised method, in accordance with the assumption that a small number of intelligent training samples should be sufficient to ensure an accurate image classification [52].The radial basis function (RBF) was used as the kernel function of the SVM.The highest accuracy derived from various kernel parameters was used to represent the performance of the supervised classifier.Landscapes in the were categorized into four classes (i.e., wadi, vegetation, bare land, and built-up land).The extraction results for the three methods, i.e., AMEW, ML, and SVM, are illustrated in in Figure 9.To better display the relevant results, we again used specific blocks selected from the three study areas to reveal the local performance.Four possibilities in each set of extraction results (i.e., TP, TN, FN, and FP) and five metrics, including TPR, FPR, EC, EO, and OA, were calculated to assess the performance of each method and were derived from error matrices through ENVI [34].The principle used was to calculate the specific pixel number correctly classified as wadi or incorrectly classified as background by comparing the extracted result with the reference map.The statistical accuracy assessment results are listed in Table 3.
the kernel function of the SVM.The highest accuracy derived from various kernel parameters was used to represent the performance of the supervised classifier.Landscapes in the were categorized into four classes (i.e., wadi, vegetation, bare land, and built-up land).The extraction results for the three methods, i.e., AMEW, ML, and SVM, are illustrated in in Figure 9.To better display the relevant results, we again used specific blocks selected from the three study areas to reveal the local performance.Four possibilities in each set of extraction results (i.e., TP, TN, FN, and FP) and five metrics, including TPR, FPR, EC, EO, and OA, were calculated to assess the performance of each method and were derived from error matrices through ENVI [34].The principle used was to calculate the specific pixel number correctly classified as wadi or incorrectly classified as background by comparing the extracted result with the reference map.The statistical accuracy assessment results are listed in Table 3.   Visual inspection revealed that the ML classifier and the SVM classifier performed poorly for wadi extraction in terms of both main streams and tributaries.In principle, they could only extract wide wadis but they ignored the majority of narrow wadis that displayed extremely weak contrast with their heterogeneous background (refer to Block 2 in Figure 9).Worse still, they usually misclassified background with similar spectral characteristics to wadis as rivers (refer to Block 6 in Figure S4 of the supplementary online materials).The AMEW performed significantly better than the other two with significantly fewer omissions and commissions.The rare omissions mainly occurred in: (i) specific segments of the wide wadis (discussed in Section 5.1); and (ii) areas that could not be significantly enhanced by WIs.Commissions mainly occurred in: (i) narrow wadis with very low spectral contrast with backgrounds; (ii) non-wadis features that were enhanced by WIs by mistake, such as roads and mountain shadows; and (iii) linear false alarms in backgrounds.Because the practical effects of the AMEW were highly reliant on the enhancement effect of WIs, there is no doubt that areas erroneously enhanced by WIs will be identified as wadis by subsequent image processing steps.Moreover, since WIs could only enhance parts of wadis given their high complexity, commissions are likely to appear in areas where the discrimination between wadis and background was limited.It was quite possible to produce extraction results with large commission errors if the network of tributaries was too dense.Note that the reference map was manually digitized based on human interpretation, and random errors were inevitably generated even when we tried to minimize them using a synthesis of results from three experienced operators.According to Table 3, there are no obvious differences between the ML classifier and the SVM classifier.Their TPR and OA reached only around 75% while the TPR and OA of the AMEW were generally above 88% and 93%, respectively.In addition, the ML and SVM classifiers demonstrated very high FPR, i.e., averages greater than 20%.Their average EC and EO were above 40% and 25%, respectively, which exceeded that of the AMEW (with an EC of 3.63%-10.60%and an EO of 5.71%-11.60%).These statistics were basically consistent with the visual illustrations in Figure 9.

Limitations of the Proposed AMEW
Although we have demonstrated an effective performance and robustness of the proposed AMEW, certain limitations are inevitable in the use of the method given the high level of complexity of the wadis networks.In this study, 2-D Gaussian matched filtering was used to enhance wadis based on their Gaussian-curve-shaped spectral characteristics.However, because a 2-D Gaussian kernel could not detect the central part of wide wadis, wide wadis and narrow wadis had to be extracted separately.Thus, we applied the local adaptive threshold segmentation on the outcome of large window-sized neighborhood analysis to obtain the wide wadis, and the complete results were generated by a fusion of narrow and wide wadis.The results indicate that the AMEW is capable of identifying most of the narrow wadis but performed comparatively poorly in the case of wide wadis.In general, it produced many gaps in the extraction results, and sometimes only the edges of wadis detected (Figure 10).large window-sized neighborhood analysis to obtain the wide wadis, and the complete results were generated by a fusion of narrow and wide wadis.The results indicate that the AMEW is capable of identifying most of the narrow wadis but performed comparatively poorly in the case of wide wadis.In general, it produced many gaps in the extraction results, and sometimes only the edges of wadis detected (Figure 10).Tests demonstrated that the larger window size (used in the neighborhood analysis of wide wadi extraction) yielded a better extraction performance.As the window size increased, the wide wadis became more and more complete (Figure 10a-d), and the TPR rates also increased (Table 4).The influence of window size is common in research using a sliding-window technique given that the extraction objects are characterized by different width and curvatures.A small window size always seems to be sufficient for recognizing the objects of interest, while a large window size often omits the linear objects in satellite images.Therefore, an iterative process for different window sizes is regarded as an effective solution [53,54].The parameter k2 of the local adaptive threshold segmentation in wide wadi extraction also affected the wide wadis performance.We found that the lower the k2, the greater the TPR (Figure 11a-d, Table 4).It should be noted that the wadi systems were extremely dense and small wadis of very low contrast made up the majority of the wadi systems.Too large a window size and too low a k2 value invariably resulted in an unacceptably crumbly extraction (Figure 10e-h and Figure 11e-h), greatly increasing the FPR as a result (Table 4, Table 5).For the majority of wadis, a medium-sized window (150 pixels) and intermediate k2 value were used to optimize the extraction effects.Tests demonstrated that the larger window size (used in the neighborhood analysis of wide wadi extraction) yielded a better extraction performance.As the window size increased, the wide wadis became more and more complete (Figure 10a-d), and the TPR rates also increased (Table 4).The influence of window size is common in research using a sliding-window technique given that the extraction objects are characterized by different width and curvatures.A small window size always seems to be sufficient for recognizing the objects of interest, while a large window size often omits the linear objects in satellite images.Therefore, an iterative process for different window sizes is regarded as an effective solution [53,54].The parameter k 2 of the local adaptive threshold segmentation in wide wadi extraction also affected the wide wadis performance.We found that the lower the k 2 , the greater the TPR (Figure 11a-d, Table 4).It should be noted that the wadi systems were extremely dense and small wadis of very low contrast made up the majority of the wadi systems.Too large a window size and too low a k 2 value invariably resulted in an unacceptably crumbly extraction (Figure 10e-h and Figure 11e-h), greatly increasing the FPR as a result (Table 4, Table 5).For the majority of wadis, a medium-sized window (150 pixels) and intermediate k 2 value were used to optimize the extraction effects.

Preliminary Effects of WIs on the AMEW
Since the wadi extraction process of the AMEW is based on the WI calculation, the extraction performance of the AMEW was highly dependent on the preliminary enhancement effect of WIs.Extraction results indicated large differences with the application of different WIs.Although the NWI outperformed the NDWI, the MNDWI, the AWEInsh and the AWEIsh, it still could not achieve a completely satisfactory enhancement.Omissions and commission errors frequently occurred in places where the discrimination between wadis and backgrounds was weak and inconspicuous with an ambiguous separability (Figure 12).Moreover, linear features that were erroneously enhanced by WIs would undoubtedly be misclassified as wadis, thereby increasing the FPR (Figure 12).In addition, the pixel values of roads in some areas were rather high in WI images, which formed distinctive contrasts to the surrounding wadis and backgrounds.Since the main task of neighborhood analysis and Gaussian-matched filtering was to detect dark linear features, the dark background around the roads was frequently misclassified as wadis, thereby generating two parallel false wadis on both sides of the road and increasing the FPR (Figure 13).Note that this error was inevitable since the method is based on morphological features.Linear features such as roads were therefore extracted by mistake.However, this would not be an irretrievable error since these obvious non-wadi linear features were not significant in number and could easily be excluded by manual inspection during refinement of the extraction results.

Preliminary Effects of WIs on the AMEW
Since the wadi extraction process of the AMEW is based on the WI calculation, the extraction performance of the AMEW was highly dependent on the preliminary enhancement effect of WIs.Extraction results indicated large differences with the application of different WIs.Although the NWI outperformed the NDWI, the MNDWI, the AWEI nsh and the AWEI sh , it still could not achieve a completely satisfactory enhancement.Omissions and commission errors frequently occurred in places where the discrimination between wadis and backgrounds was weak and inconspicuous with an ambiguous separability (Figure 12).Moreover, linear features that were erroneously enhanced by WIs would undoubtedly be misclassified as wadis, thereby increasing the FPR (Figure 12).In addition, the pixel values of roads in some areas were rather high in WI images, which formed distinctive contrasts to the surrounding wadis and backgrounds.Since the main task of neighborhood analysis and Gaussian-matched filtering was to detect dark linear features, the dark background around the roads was frequently misclassified as wadis, thereby generating two parallel false wadis on both sides of the road and increasing the FPR (Figure 13).Note that this error was inevitable since the method is based on morphological features.Linear features such as roads were therefore extracted by mistake.However, this would not be an irretrievable error since these obvious non-wadi linear features were not significant in number and could easily be excluded by manual inspection during refinement of the extraction results.

Extended Applications of the AMEW by Means of Other Optical Satellite Data
As mentioned in Section 5.2, the extraction results were dependent on the preliminary enhancement effect of the WIs.As a consequence, there was no strict requirement for satellite images with similar spectral features.The AMEW performed well as long as the preliminary enhancement of the WIs was carried out effectively.The aforementioned tests indicated that, using Landsat 8 OLI images, the AMEW is capable of precisely mapping wadis with a satisfactory degree of accuracy.Thus, it is assumed that optical data with similar wavelength features to Landsat OLI images could also be used.This point is of great importance since the extended application of Landsat OLI from spatial and temporal resolution perspectives can be used to accomplish more meaningful research, such as time series change detection, combining OLI imagery with previous Landsat series imagery such as TM and ETM+ of high quality.Moreover, based on the same principle, it can also be assumed that the AMEW can be extended to other series of optical imagery with different spatiotemporal resolutions, such as the visible infrared imaging radiometer suite (VIIRS) and the MODIS [18], which is important for research on a continental scale as well as for near-real time monitoring of valuable natural resources.

Conclusions
Wadis are a valuable water resource in desert regions, and their accurate mapping and monitoring is of great importance for regional hydrological management and utilization.While currently available river extraction methods are reported to be successful in specific aquatic environments, they have limited applicability to wadi extraction given their irregular drainage patterns, complex network systems, and inconspicuous spectral characteristics.In this study, an automated method for extracting wadis (AMEW) was developed in order to delineate wadis networks from Landsat-8 OLI imagery.The approach combines WIs and subsequent image processing technologies from their morphological features, including multi-scale Gaussian matched filtering and local adaptive threshold segmentation.Experimental results demonstrate that the AMEW with high accuracy, automation, and promising applicability can be used for wadi extraction on a large scale.
We produced a composite of three classical WIs (MNDWI, the AWEI nsh , and the AWEI sh ) by HIS transformation in order to construct a new WI (NWI) to preliminarily suppress the influence from massive noise and false alarms, in an attempt to take advantage of various WIs.Experimental results demonstrated that the NWI, combining the merits of other WIs, achieved a significantly increased performance and was the optimal WI for suppressing heterogeneous noise.In addition, extended application to other optical data with similar spectral wavelengths demonstrates that the proposed AMEW could be easily transferred to different data with varying spatiotemporal resolution and should be useful for extraction and mapping of other land cover types.

Figure 1 .
Figure 1.(a) Locations of the three study areas in Jordan; (b) true-color image (432 band for OLI) of Study Area 1; (c) true-color image of Study Area 2; and (d) true-color image of Study Area 3.

Figure 3 .
Figure 3. (a-e) WI value distributions for wadi, bare land, vegetation, and built-up land.In the box plots, the top and bottom of the vertical lines represent the maximum and minimum; the top and bottom of the boxes represent the first quartile and the third quartile; the lines in the boxes represent the median.The red line shows the reflectance separability of the third quartile of wadis reflectance with other landscapes, indicating possible reflectance confusion.The JMD is the average value for wadis and other three types.

Figure 3 .
Figure 3. (a-e) WI value distributions for wadi, bare land, vegetation, and built-up land.In the box plots, the top and bottom of the vertical lines represent the maximum and minimum; the top and bottom of the boxes represent the first quartile and the third quartile; the lines in the boxes represent the median.The red line shows the reflectance separability of the third quartile of wadis reflectance with other landscapes, indicating possible reflectance confusion.The JMD is the average value for wadis and other three types.

23 Figure 5 .
Figure 5. (a) WI value profiles of the cross sections of several wadis in a WI image; (b) pixel value profiles of the cross-sections of several wadis in neighborhood analysis results.

Figure 6 .
Figure 6.Results of neighborhood analysis: (a) NWI image; (b) window size = 16 pixels; (c) window size = 48 pixels; (d) window size = 80 pixels.The window sizes shown here are only used to provide a conceptual illustration of the application of the technique.

Figure 5 .
Figure 5. (a) WI value profiles of the cross sections of several wadis in a WI image; (b) pixel value profiles of the cross-sections of several wadis in neighborhood analysis results.

) 23 Figure 5 .
Figure 5. (a) WI value profiles of the cross sections of several wadis in a WI image; (b) pixel value profiles of the cross-sections of several wadis in neighborhood analysis results.

Figure 6 .
Figure 6.Results of neighborhood analysis: (a) NWI image; (b) window size = 16 pixels; (c) window size = 48 pixels; (d) window size = 80 pixels.The window sizes shown here are only used to provide a conceptual illustration of the application of the technique.

Figure 6 .
Figure 6.Results of neighborhood analysis: (a) NWI image; (b) window size = 16 pixels; (c) window size = 48 pixels; (d) window size = 80 pixels.The window sizes shown here are only used to provide a conceptual illustration of the application of the technique.
Remote Sens. 2016, 8, 246 13 of 23was used for pre-processing for the subsequent wadi extraction work.Hereafter, the AMEW refers to the method combined with the NWI.

Figure 8 .
Figure 8. (a-d) WI images of Study Area 1; (e-h) wadi extraction results for Study Area 1; (i-l) WI images of Block 1; (m-p) wadi extraction results for Block 1; (q-t) WI images for Block 2; (u-x) wadi extraction results for Block 2.

Figure 8 .
Figure 8. (a-d) WI images of Study Area 1; (e-h) wadi extraction results for Study Area 1; (i-l) WI images of Block 1; (m-p) wadi extraction results for Block 1; (q-t) WI images for Block 2; (u-x) wadi extraction results for Block 2.

Figure 9 .
Figure 9. (a) NWI image of Study Area 1; (b-d) extraction results for Study Area 1; (e) NWI image of Block 1; (f-h) extraction results for Block 1; (i) NWI image for Block 2; (j-l) extraction results for Block 2. For further illustrations of results, see Figures S3 and S4 in the supplementary online materials.

Figure 9 .
Figure 9. (a) NWI image of Study Area 1; (b-d) extraction results for Study Area 1; (e) NWI image of Block 1; (f-h) extraction results for Block 1; (i) NWI image for Block 2; (j-l) extraction results for Block 2. For further illustrations of results, see Figures S3 and S4 in the supplementary online materials.

Figure 10 .
Figure 10.Reference maps and wadi extraction results based on different window sizes.As the window size increases, wide wadis become more and more complete, but commission errors around small wadis become increasingly unacceptable; (a-d) different window sizes for block 1; (e-h) different window sizes for block 2.

Figure 11 .
Figure 11.Reference images and wadi extraction results based on different values of k2.As k2 decreases, wide wadis become more and more complete, but over-extractions around small wadis become more and more significant; (a) reference images for block 1; (b-d) the results for different k2 in block 1; (e) reference image for block 2; (f-h) the results of different k2 for block 2.

Figure 11 .
Figure 11.Reference images and wadi extraction results based on different values of k2.As k2 decreases, wide wadis become more and more complete, but over-extractions around small wadis become more and more significant; (a) reference images for block 1; (b-d) the results for different k2 in block 1; (e) reference image for block 2; (f-h) the results of different k2 for block 2.

Figure 12 .
Figure 12. (a,b) original image band and the commission and omission errors in the vicinity of linear false alarms in a local block; (c,d) the other original image band and the commission and omission errors in the vicinity of linear false alarms in another local block.

Figure 13 .
Figure 13.NWI images and over-extractions around roads (two parallel false wadis on both sides).(a,b) are the reference map and the extraction result for a local block; (c,d) are the reference map and the extraction map for another local block.

Figure 12 . 23 Figure 12 .
Figure 12. (a,b) original image band and the commission and omission errors in the vicinity of linear false alarms in a local block; (c,d) the other original image band and the commission and omission errors in the vicinity of linear false alarms in another local block.

Figure 13 .
Figure 13.NWI images and over-extractions around roads (two parallel false wadis on both sides).(a,b) are the reference map and the extraction result for a local block; (c,d) are the reference map and the extraction map for another local block.

Figure 13 .
Figure 13.NWI images and over-extractions around roads (two parallel false wadis on both sides).(a,b) are the reference map and the extraction result for a local block; (c,d) are the reference map and the extraction map for another local block.
) WI images of Study Area 2; (e-h) wadi extraction results of Study Area 2; (i-l) WI images of Block 3; (m-p) wadi extraction results of Block 1; (q-t) WI images of Block 4; (u-x) wadi extraction results of Block 4., Figure S2: (a-d) WI images of Study Area 2; (e-h) wadi extraction results of Study Area 2; (i-l) WI images of Block 5; (m-p) wadi extraction results of Block 5; (q-t) WI images of Block 6; (u-x) wadi extraction results of Block 6., Figure S3: (a) NWI image of Study Area 2; (b-d) the extraction results of Study Area 2; (e) NWI image of Block 3; (f-h) the extraction results of Block 3; (i) NWI image of Block 4; (j-l) the extraction results of Block 4, Figure S4: (a) NWI image of Study Area 3; (b-d) the extraction results of Study Area 3; (e) NWI image of Block 5; (f-h) the extraction results of Block 5; (i) NWI image of Block 6; (j-l) the extraction results of Block 6.

Table 1 .
Mathematical formulations of four common WIs.

Table 2 .
Performance of the AMEW based on the AWEInsh, the AWEIsh, the MNDWI and the NWI.

Table 3 .
Performance evaluation of the ML classifier, the SVM classifier and the AMEW.

Table 3 .
Performance evaluation of the ML classifier, the SVM classifier and the AMEW.

Table 4 .
Performance evaluation of wadi extraction results based on different window sizes.
Figure 10.Reference maps and wadi extraction results based on different window sizes.As the window size increases, wide wadis become more and more complete, but commission errors around small wadis become increasingly unacceptable; (a-d) different window sizes for block 1; (e-h) different window sizes for block 2.

Table 4 .
Performance evaluation of wadi extraction results based on different window sizes.

Table 5 .
Performance evaluation of wadi extraction results based on different values of k2.

Table 5 .
Performance evaluation of wadi extraction results based on different values of k2.