Automatic Inundation Mapping Using Sentinel-2 Data Applicable to Both Camargue and Doñana Biosphere Reserves

Abstract

Flooding periodicity is crucial for biomass production and ecosystem functions in wetland areas. Local monitoring networks may be enriched by spaceborne derived products with a temporal resolution of a few days. Unsupervised computer vision techniques are preferred, since human interference and the use of training data may be kept to a minimum. Recently, a novel automatic local thresholding unsupervised methodology for separating inundated areas from non-inundated ones led to successful results for the Doñana Biosphere Reserve. This study examines the applicability of this approach to Camarque Biosphere Reserve, and proposes alternatives to the original approach to enhance accuracy and applicability for both Camargue and Doñana wetlands in a scientific quest for methods that may serve accurately biomes at both protected areas. In particular, it examines alternative inputs for automatically estimating thresholds while applying various algorithms for estimating the splitting thresholds. Reference maps for Camargue are provided by local authorities, and generated using Sentinel-2 Band 8A (NIR) and Band 12 (SWIR-2). The alternative approaches examined led to high inundation mapping accuracy. In particular, for the Camargue study area and 39 different dates, the alternative approach with the highest overall Kappa coefficient is 0.84, while, for the Doñana Biosphere Reserve and Doñana marshland (a subset of Doñana Reserve) and 7 different dates, is 0.85 and 0.94, respectively. Moreover, there are alternative approaches with high overall Kappa for all areas, i.e., at 0.79 for Camargue, over 0.91 for Doñana marshland, and over 0.82 for Doñana Reserve. Additionally, this study identifies the alternative approaches that perform better when the study area is extensively covered by temporary flooded and emergent vegetation areas (i.e., Camargue Reserve and Doñana marshland) or when it contains a large percentage of dry areas (i.e., Doñana Reserve). The development of credible automatic thresholding techniques that can be applied to different wetlands could lead to a higher degree of automation for map production, while enhancing service utilization by non-trained personnel.

1. Introduction

Wetlands, which constitute unique habitats for many different plant and animal species, are important for their water-related ecosystem services, such as food provision, water filtration, and protection against soil erosion [1]. Additionally, they provide important recreational and leisure activities, such as bird watching, fishing, and hiking [2]. Wetlands are in danger of rapid decline in both quantity and quality due to impacts related to climate change and human pressures [3]. Therefore, monitoring the variability of the surface water extent across time is important for taking actions to increase resilience. Satellite data can provide a cost-effective solution for frequent and accurate monitoring of surface water extent.

Numerous approaches, utilizing optical or radar data for estimating water surface areas, have been proposed in the literature. The advantage of radar-based approaches is that they can operate under nearly all-weather and day-night conditions. However, emergent vegetation [4], waves [5,6], sand [7], and radar shadows produced by terrain features [8] impede the efficient delineation between water and land. On the other hand, the extraction of the water surface from optical imagery is generally more straightforward than radar imagery [9,10], since the rich spectral information of optical data allows for the reliable detection of the water presence by utilizing various indices and bands. The main limitation of using optical data is cloud presence, which prohibits the observation of the earth’s surface [10]. Several approaches utilize both optical and radar data to deal with the lack of optical data during extended periods of cloud cover and overcome the limitations of radar data [11,12]. The methodology presented in this work relies on optical data. Thus, literature information focuses on this category of approaches.

Thresholding approaches detect water-covered areas by applying thresholds to one or more spectral bands or indices [13,14,15,16,17,18,19]. Commonly used indices include the Normalized Difference Water Index (NDWI) [15,16,19], Modified NDWI [14,18,19,20], and Automated Water Extraction Index [17,18,19,21]. Several approaches use information from Shortwave infrared (SWIR) spectral ranges to identify shallow inundated wetland areas, since it is less sensitive to sediment-filled waters and, hence, more efficient for registering the boundaries between water and dry areas in shallow wetlands [13,22,23,24]. Automatic thresholding approaches can be applied to different areas and are computationally inexpensive, but they may wrongly classify dark objects (i.e., shadows and buildings) as water when their spectral characteristics are similar [25]. Automatic thresholding approaches are distinguished into: (a) global approaches [15,17,19,20,26], which estimate thresholds based on the histogram analysis of the complete image, and (b) local thresholding approaches [23], which estimate local thresholds for image subsets containing high percentages of pixels belonging to the water and non-water classes, and then may take into consideration subsets’ thresholds to estimate an overall threshold. Local thresholding approaches overcome the incapability of global approaches to estimate an optimal histogram threshold when the class proportions within the image are imbalanced [27]. Various algorithms (such as Otsu’s [28], and Kittler and Illingworth’s [29]) have been used for estimating thresholds separating inundated and non-inundated pixels inside an image.

Machine-learning algorithms, including supervised and unsupervised ones, have also been used for detecting water bodies from multispectral imagery. Supervised approaches applied for water detection include random forests [30,31], neural networks [32], support vector machines [33,34], deep neural networks [35,36], and decision trees [13,22,37], while frequently used unsupervised approaches for water mapping include K-means [38] and Iterative Self-Organizing Data Analysis Technique (ISODATA) [39]. Even though machine-learning approaches usually exhibit higher accuracy than thresholding ones [30], they have several limitations. (a) In case reference data is not available, supervised approaches require collection of training samples, which is a time-consuming and tedious task that requires expert knowledge and/or validation in the field [25,40]. (b) Supervised methods may meet problems when mapping water bodies over large scale areas [41]. (c) Unsupervised methods need expert knowledge to select the initial class and iteration parameters [42], and may need post processing of the results to combine adjacent regions into larger regions corresponding to water bodies [30].

Most of the previously mentioned approaches aim to detect water in open-water bodies such as rivers, lakes, reservoirs, and watersheds [17,19,20,30,33,37,43], while a part of them focuses on wetland areas [13,22,23,24].

This study focuses on wetland areas and examines alternative approaches of the original automated local thresholding approach presented in Reference [23] with the objective to suggest alternative approaches, which may produce credible results for both Camargue and Doñana wetlands, i.e., be identified as possibly applicable to further wetland areas and biomes. Each examined alternative approach relies on a specific band or band combination, acknowledged as effective by the underlying physics, and a specific approach for estimating splitting thresholds. The different Sentinel-2 (S2) based inputs examined for estimating thresholds include: (i) Band 11 (SWIR-1), (ii) product (result of the multiplication) of Band 12 (SWIR-2) and Band 8A (NIR), and (iii) product of SWIR-1 and NIR. The different methods for estimating splitting thresholds include: (i) minimum entropy thresholding, and (ii) Otsu’s algorithm. The results of the alternative approaches are compared against reference maps, provided for Doñana and Camargue by local research institutes, based on locally developed water detection models [22,24].

2. Materials and Methods

2.1. Study Areas

2.1.1. Camargue

The Camargue (Figure 1) is a polderized delta created by the Rhône River that covers an area of 145,300 ha rarely exceeding 5 m in altitude. It comprises a high diversity of wetland types according to a water-salinity gradient such as lagoons, salt meadows, dense halophilous scrubs and steppes, brackish/freshwater marshes with tall emergent or permanent aquatic vegetation, and temporary pools. The climate is Mediterranean, being characterized by mild and wet winters and hot and dry summers. With mean annual precipitations of 600 mm, mainly concentrated from autumn to early spring, and a mean evapotranspiration of 1400 mm, the Camargue is characterized by high water deficits, especially in the summer [44,45]. As a result, 730 millions of cubic meters of water are pumped from the Rhône on average each year to compensate for river embankment, avoid soil salinization, and enhance primary production. This water, primarily pumped from March to September, is distributed through a complex network of channels for irrigating crops and pastures, as well as for flooding marshes, which support nature conservation, wildfowl hunting, and reed harvest. Accordingly, some Camargue wetlands are flooded year round or during most of the year as a result of artificial irrigation (e.g., hunting reed marsh). Other wetlands become naturally dry during the summer season (e.g., harvested reed marsh, bulrush marsh) and, lastly, a few are flooded only during a period of high rainfalls (e.g., halophilous scrubs, temporary pools) [46]. Within the context of climate change and increasing human pressures, the development of a remotely-sensed tool to monitor water seasonal variation in these wetlands is essential [47].

2.1.2. Doñana

The Doñana wetlands (Figure 2), covering 108,429 ha, lie within and around the delta of the Guadalquivir River in Southwest Spain and contain two main habitat types: seasonal marshes and adjacent eolian sands holding aquifer-fed dune ponds. These are surrounded by scrublands, pine forests, and cultivated areas. The marshland area is comprised of seasonal open marsh with emergent plants, temporary pools with annual plant species, and scattered halophilous scrubs [48]. The Doñana climate is Mediterranean sub-humid with hot and dry summers and mild and wet winters. The annual precipitation, with an average of 550 mm, occurs mainly between October and April and is almost absent between May and September. The highest monthly rainfall usually occurs in November and maximum water levels are reached during February. Marshes dry up slowly in late spring and most of their surface gets completely dry by the end of July [49]. Marshland’s depth, turbidity, and vegetation cover varies depending on the amount and seasonal pattern of precipitation [24]. Marshes are breeding ground and host many species of migratory birds during the winter [49]. In parallel, rice-paddies, aquaculture ponds, and salt pans constitute an extensive system of artificial wetlands in Doñana that also provide a habitat for water birds and other species, especially when natural wetlands dry up seasonally [50].

2.2. Dataset

2.2.1. Camargue

Thirty-nine cloud-free and atmospherically corrected S2 Level-2A (L2A) products of the Camargue from 12 June, 2017 to 19 June, 2018 were downloaded from the Copernicus European Space Agency (ESA) hub (see

Table 1. Dates of cloud free S2 acquisitions over Camargue.

Cycle	Sept.	Oct.	Nov.	Dec.	Jan.	Feb.	Mar.	Apr.	May	June	July	Aug.
2016–2017										12,19	4,7 12,14 17	3,13 18,21
2017–2018	5,7 20,27	5,7 10,12 27,30	14,16 19,21	6,9 16,24	23	4,22 27	14	18,20	20,25	19

). The shapefile containing the boundaries of Camargue comprises two different tiles (namely 31TEJ, 31TFJ). The tile products are mosaicked and clipped to the extent of the shapefile.

The reference maps for Camargue were obtained by dichotomous partitioning of reflectance values encoded as 1 for water presence and 0 for water absence based on ground-truth (n = 1229) and optical-space derived (n = 2603) reference points covering the whole Biosphere Reserve area and all the main habitat types [22]. Ground-truth data refer to water level measures in different wetland types, focusing on those with a dense vegetation cover. For optical-space data, a water formula developed with SPOT-5 (Satellite Pour l’Observation de la Terre - 5) [13] was used to generate a water mask for dates when Landsat-8 scenes were also available. Fifty random points from each of the main land cover types in the Camargue (n = 17) were selected on the SPOT-5 water maps and transferred to Landsat-8 scenes with similar dates. These points were then used as training data along with ground-truth measures for creating a Landsat-8 water formula. The same procedure was repeated to develop a water formula with S2 from Landsat-8. Data mining was performed with the Rpart package in the R software using eight bands of S2, as well as various spectral indices found in the literature for explaining the variables (for more details, see References [13,22]). A random selection of 30% of all points was excluded from the sample and used for validation. The best model selected used the near (NIR) and short-wave (SWIR-2) infrared wavelengths. NIR was useful for discriminating areas that are completely dry, while SWIR-2 was efficient for water detection. This model provided an overall accuracy of 94% for predicting water presence/absence with Kappa coefficients of 0.82 on both the training and validation samples [22]. Wetlands characterized by a dense cover of vegetation were correctly classified at 89%.

2.2.2. Doñana

Seven cloud-free and atmospherically corrected S2 L2A products of the Doñana from 1 June, 2017 to 17 April, 2018 were downloaded from the ESA hub, so that they timely overlap with Landsat cloud-free parallel data (see

Table 2. Dates of cloud free S2 acquisitions over Doñana, coinciding with Landsat ones, with the exception of one date. Landsat acquisition was one day earlier than the S2 acquisition on 21 February, 2018.

Date

01/06/2017

11/07/2017

20/08/2017

08/11/2017

27/01/2018

21/02/2018

17/04/2018

). The shapefile containing the boundaries of Doñana Biosphere Reserve comprises three different tiles (namely 29SQA, 29SQB, and 29SPB). The tile products were mosaicked and clipped to the extent of the shapefile.

Inundation maps of the Doñana area with 30-m pixel resolution, which are generated from Landsat satellite data and provided by the Doñana Biological Station (EBD), are used as ground truth data. These reference maps were obtained by dichotomous partitioning of reflectance values on Landsat Thematic Mapper (TM) and Enhanced TM (ETM) band 5 (SWIR) based on 6005 ground-truth field sampling points, mostly collected within the marshland area of Doñana [24]. This model provided an overall accuracy of 93% for predicting water presence/absence with a Kappa coefficient of 0.65.

2.3. Methodology

The work presented in Reference [23] introduces an unsupervised approach, detecting automatically thresholds on the SWIR-1 band and on a Modified-Normalized Difference Vegetation Index (MNDVI) (estimated as the normalized difference of Band 7 and Band 5 of S2), to estimate open-water and water-vegetation subclasses. This approach demonstrated high classification accuracy for Doñana, with an overall Kappa coefficient reaching 0.94 and 0.88 for the marshland and the complete area (i.e., the Biosphere Reserve), respectively. This paper examines alternatives of the original thresholding approach driven by the findings for Camargue. Instead of the SWIR-1 band, other algebraic band combinations are also examined. Additionally, both the minimum cross entropy thresholding algorithm (MCET) [51] and Otsu’s algorithm [28] are used for estimating thresholds partitioning inundated and non-inundated pixels inside an image subset based on their class distribution. The methodology’s steps are presented in Figure 3.

2.3.1. Segmentation of the Satellite Image

As the first step, the satellite image is segmented into non-overlapping regions by utilizing the mean-shift segmentation algorithm [52]. The resulting segmentation map is utilized for selecting segments with a high percentage of inundated pixels.

The input to the segmentation algorithm is a false color image composed of Band 2 (BLUE), Band 3 (GREEN), and Band 4 (RED), which have been normalized to the range [0,255], by relying on minimum and maximum values representing the intensity percentile range from 1% to 99% per band. These bands are selected due to their 10-m resolution, which allows for a more accurate segmentation of the satellite image, compared to selection of other bands with a lower resolution [23].

2.3.2. Mapping of the Open-Water Subclass

Three input data alternatives are examined for estimating an initial threshold T_init separating inundated from non-inundated pixels: (i) SWIR-1 band, which is denoted as Alt1, (ii) product (per pixel multiplication) of SWIR-2 and NIR, which is denoted as Alt2, and (iii) product of SWIR-1 band and NIR, which is denoted as Alt3. The first alternative Alt1 is the one proposed in the original approach [23]. The second alternative Alt2 is inspired by the approach suggested in Reference [22] for Camargue, which applies strict thresholds to SWIR-2 and NIR (see paragraph 2.2.1), and the third alternative Alt3 is a variant of the Atl2 where SWIR-1 is used instead of SWIR-2. Each of Alt1, Alt2, and Alt3 is normalized in the range of [0,255], by relying on minimum and maximum values representing the intensity percentile range from 1% to 99%. The term “Inp”, used in the following, is a generic term for the input data, corresponding to either Alt1, Alt2, or Alt3.

The histogram of Inp is used to estimate an initial threshold T_init separating inundated from non-inundated pixels. T_init is the Inp value for which the first deep valley of this histogram is detected (Figure 4 shows histograms using Alt2 as input). If a pixel p has Inp(p) < T_init, it is denoted as inundated. Otherwise, it is denoted as non-inundated. Therefore, an initial inundation map is generated. Based on this inundation map, segments G_m, where m = 1,2,…,M, having a large percentage of inundated pixels (i.e., over 70% of a segment’s pixels are inundated) are selected and their corresponding centroids C_m are estimated.

Square patches ${\mathrm{P}}_{\mathrm{m}}^{\mathrm{k}}$ of expanding window size are centered around each C_m. The window size in pixels is given by [20∙k × 20∙k], where k = 1,2,…,20. The splitting thresholds $f_{\mathrm{m}}^k$ of patches with bimodal histogram (i.e., two distinctive classes of intensity values corresponding to inundated and non-inundated classes appear in the histogram) are estimated based on the histogram thresholding algorithms (MCET or Otsu’s algorithm). Their median is assumed to be the optimal threshold $f_{opt}^m$ corresponding to the segment G_m. The usefulness of using expanding patches is to calculate the optimal threshold per segment more robustly, since its estimation is based on multiple splitting thresholds. Then, the median of selected segments’ optimal thresholds is estimated as: $M_{opt}=\underset{m=1,2,\dots ,M}{\mathrm{median}}(f_{opt}^m)$.

Final threshold T_final (see Figure 4), discriminating the open-water subclass, is estimated as max(M_opt,T_init). Pixels p with Inp(p) < T_final are assumed to belong to the open-water subclass, comprising of water or water with sparse vegetation pixels.

When MCET or Otsu is used for estimating the splitting thresholds, T_final is denoted as $T_{final}^{\mathrm{MCET}}$ or $T_{final}^{\mathrm{Otsu}}$, respectively. Average (abbreviated as “avg”): $T_{final}^{\mathrm{avg}}=\left(T_{final}^{\mathrm{MCET}}+T_{final}^{\mathrm{Otsu}}\right)/\mathbf{2}$ is also examined as the final threshold.

2.3.3. Mapping of the Water-Vegetation Subclass

Inundated areas may be covered by emergent vegetation. The Inp value of the pixels in these areas is higher compared to the Inp values of the pixels with open water or water with sparse vegetation belonging to the open-water class. Therefore, a threshold T_upper, which is the Inp value for which the deep valley right after the first deep valley is detected in the Inp histogram, is assumed to represent the upper threshold for pixels comprising areas of water covered by dense vegetation (see Figure 4b). At the same time, areas with dense emergent vegetation are expected to have high MNDVI value. Thus, a threshold T_MNDVI, which should be more than 0.4, is selected. In particular, T_MNDVI is set equal to the MNDVI value, for which the first valley in the part of the histogram with MNDVI values over 0.4 is detected (see Figure 5b) following the findings in Reference [23].

A pixel p is assumed to belong to the water-vegetation subclass if the following two conditions are met: T_final < Inp(p) < T_upper & MNDVI(p) > T_MNDVI, presuming that T_upper and T_MNDVI can be detected in the histograms of Inp and MNDVI, respectively. These latter conditions are valid for Doñana on 20 August, 2017 (see Figure 4b and Figure 5b), but not for Camargue on 18 August, 2017 (see Figure 4a and Figure 5a).

3. Results

3.1. Comparison of Automatic Thresholding Results Against Reference Map of Camargue

The three input alternatives (Alt1 or Alt2, or Alt3) and the three different ways of estimating the final threshold ($T_{final}^{\mathrm{MCET}}$ or $T_{final}^{\mathrm{Otsu}}$ or $T_{final}^{\mathrm{avg}}$) form nine different alternatives, with each one combining an input alternative and a way of estimating the final threshold. The alternative ways of estimating final threshold are referred to as “MCET” (standing for $T_{final}^{\mathrm{MCET}}$), “OTSU” (standing for $T_{final}^{\mathrm{Otsu}}$), and “avg” (standing for $T_{final}^{\mathrm{avg}}$). For example, if “Alt2” is the input and “avg” is the way of estimating the final threshold, then the alternative will be named as “Alt2 (avg)”.

Regarding Camargue, each of the alternatives is applied to the 39 S2 products (see Table 1), and the results are compared to the reference maps. Pixels corresponding to the coastal waters are not taken into consideration for the comparison to obtain unbiased classification results. The results are presented using Kappa coefficient, which is estimated in terms of the observed agreement (p_o) and the agreement expected by chance (p_e) as: Kappa= (p_o − p_e) / (1-p_e) [53].

The curves presented in Figure 6 show how the Kappa coefficient varies for the time period between 12 June, 2017 and 19 June, 2018 for five out of nine alternatives with the best agreement with the reference map. The best agreement is given by “Alt3 (OTSU)”. “Alt3 (avg)” and “Alt2 (avg)” rank second and third, respectively, and have similar results. Lastly, “Alt3 (MCET)” and “Alt2 (MCET)” rank fourth and fifth, respectively, and give very close results. The previously mentioned ranking is based on the overall Kappa accuracy provided in Section 3.3. For all curves, it is seen that agreement decreases in the winter months from XII to II. For this period, Kappa is below 0.7 for most of the dates, while, for the rest of the year, the Kappa value is generally more than 0.7.

Figure 7 shows an example of inundation maps obtained from an S2 scene on 27 February, 2018 when using the alternative approaches of Alt3 input as compared to the reference map. Green and pink colors indicate classification differences between the reference map and the Alt3 maps. The map of “Alt3 (OTSU)” (Figure 7b) agrees better with the reference map (Figure 7a), compared to “Alt3 (MCET)” (Figure 7d), which largely underestimates water presence, and “Alt3 (avg)” (Figure 7c), which shows intermediate results. However, while “Atl3 (OTSU)” is able to detect water in vegetated areas more efficiently, such as the ones enclosed in the yellow squares, it tends to overestimate water presence in agricultural areas.

3.2. Comparison of Automatic and Thresholding Results Against Landsat Reference Maps of Doñana

In order to account for the uncertainty caused by the lower spatial resolution of the Landsat-derived reference maps, pixels in the transition zones between inundated and non-inundated areas (i.e., pixels including in their eight closest-neighbor pixels at least one pixel of a different class) in the S2 maps were excluded from the accuracy estimation (see Reference [23] for more information). Pixels in the area corresponding to sea coastal waters were also excluded. The curves, presented in Figure 8, show how Kappa coefficient varies for seven different dates regarding the five top ranking alternatives at the Doñana marshland (Figure 8a) and Doñana complete area (Figure 8b).

In Doñana marshland, “Alt1 (OTSU)” provides the highest accuracy, followed by “Alt2 (MCET)” and “Alt1 (Avg)” (see Section 3.3 for more details on the ranking that is based on overall Kappa).

In the Doñana complete area, “Alt3 (MCET)” provides the highest accuracy, followed by “Alt1 (MCET)” and “Alt2 (MCET)” (see Section 3.3 for more details on the ranking that is based on overall Kappa). “Alt1 (avg)” provides the most consistent Kappa values across the study period.

Figure 9 shows an example of the Landsat derived inundation reference map on 27 January, 2018 and the inundation maps estimated based on the alternative approaches using Alt3 input. Green and pink colors indicate the classification differences between the reference map and the alternative maps. “Alt3 (MCET)” (Figure 9d) agrees well with the reference map (Figure 9a), while “Alt3 (OTSU)” (Figure 9b) largely overestimates the water presence and “Alt3 (avg)” (Figure 9c) shows intermediate results.

3.3. Overall Kappa Per Approach and Examined Areas

Figure 10 provides overall Kappa coefficients, which are estimated when the number of true positive pixels of the inundated class, false positive pixels of the inundated class, true positive pixels of the non-inundated class and false positive pixels of the non-inundated class are added for dates per approach and examined study. The combined Overall Accuracy of each approach per study case is also given as a number close the point corresponding to its overall Kappa. Kappa values can be classified according to Reference [53] as follows: “Moderate” when 0.40 < Kappa ≤ 0.60 (shortened as “Mod”), “substantial” when 0.60 < Kappa ≤ 0.80 (shortened as “Sub”) or “almost perfect” when 0.80 < Kappa ≤ 1 shortened as (“Alm Perf”).

“Alt2 (avg)” and “Alt3 (avg)” (see yellow arrows in Figure 10) have both “Alm Perf” overall Kappa when considering Camargue and Doñana marshland study areas. While, for the Doñana complete area, their Kappa is “Sub”, since overall Kappa for “Alt2 (avg)” and “Alt3 (avg)” is 0.63 and 0.73, respectively. On the other hand, “Alt2 (MCET)” and “Alt3 (MCET)” (see blue arrows in Figure 10) exhibit more consistent performance across sites, since their Kappa is “Alm Perf” for both Doñana marshland and Doñana complete area, while, for Camargue, their overall Kappa is very close to “Alm Perf.”

4. Discussion

The aim of this paper is to examine the performance of various alternative automatic thresholding approaches for inundation mapping in two different wetland areas, which include Camargue and Doñana. These are considered as testbeds to evaluate the transferability of the different approaches and select the ones that favor the acquisition of consistently credible results for both wetlands. In this study, the local thresholding alternative approaches utilize multispectral satellite data, while the vast majority of local thresholding approaches have been designed for utilizing radar data (e.g., [27,54,55,56,57,58]). Moreover, the estimation of local thresholds relies on patches of expanding size, contrary to most of the approaches using image subsets of fixed size [55,56,57,58]. This helps to increase robustness, when input data have different spatial resolution such as when inputs from different satellites are used.

The overall Kappa of the alternative approaches varies from 0.73 to 0.84 for Camargue, from 0.74 to 0.94 for Doñana marshland, and from 0.45 to 0.85 for the complete area of Doñana. The experimental results demonstrate that there is not a common approach across all study cases achieving absolute top accuracy. “Alt3 (OTSU)” is best for Camargue, “Alt1 (OTSU)” is best for Doñana marshland, and “Alt3 (MCET)” is best for Doñana as a whole. The use of “Alt2” and “Alt3” input alternatives leads to the estimation of additional inundated areas compared to “Alt1”, but sometimes leads to overestimating the water presence in dry areas. “OTSU” algorithm identifies the final threshold at a relatively high value leading to a possible overestimation of the inundated areas, while the “MCET” algorithm identifies a final threshold at a lower value, which leads to a possible underestimation of the inundated areas, especially in emergent vegetation areas. As a sequence, the use of “avg” helps to balance between the overestimation and underestimation of inundated areas. “Alt2 (avg)” and “Alt3 (avg)” are the most accurate approaches when considering, at the same time, Camargue and Doñana marshland, with both approaches exhibiting overall Kappa over 0.81 and 0.86 for the two study areas, respectively. These areas are dominated by emergent vegetation and temporary flooded areas, and, thus, “Alt2 (avg)” and “Alt3 (avg)” could be utilized for other wetlands with similar characteristics. However, the overall Kappa is below 0.74 for both approaches when considering the Doñana complete area. The complete Doñana area is comprised of larger parts that are permanently dry. In this case, the utilization of MCET achieves a higher accuracy. Overall, “Alt2 (MCET)” and “Alt3 (MCET)” showcase a consistent overall Kappa exceeding 0.79 for all three study areas, since they avert the erroneous detection of water in dry areas. Thus, these two alternatives should be preferred for wetlands comprising large parts of dry areas.

The fact that different approaches seem to operate best in one or the other study area may be related to the models used for generating the reference maps, as well as the original field data used to build the local water presence models. In Doñana as in Camargue, the reference maps are from formulas developed to optimize water detection in the wetland types that represent very effective local conditions. Doñana marshland mostly includes seasonal marshes with relative sparse and low emergent vegetation due to the short flooding period. Tall emergent plants are rare and they are generally grazed by cattle when present [48]. Hence, the formula originally developed in Doñana does not have to perform well under dense vegetation cover. On the other hand, tall emergent plants cover large areas of semi-permanent marshes in Camargue, and the formula originally developed was specifically meant to detect water under dense and tall vegetation cover [22]. The sensors and spectral bands used to develop the original formula at each site could also influence the performance of each alternative approach tested. In Doñana, the estimation of the reference maps relied on band 5 of Landsat TM and ETM, which is similar to the SWIR-1 band of S2. The best performing alternative at this site is associated with the SWIR-1 band (e.g., Alt1 data input). In Camargue, the original formula used a combination of SWIR-2 and NIR S2 bands [22], and the three best performing alternatives use NIR in combination with SWIR-1 (e.g., Alt2 data input) or SWIR-2 (e.g., Alt3 data input) bands.

In Camargue, the accuracy of all alternative approaches decreases during winter (from December to February). This systematic error is largely due to the misclassification of the water presence under dense cover of scrubby vegetation (Salicornia marshes). This habitat, which is flooded by rainfall during this specific time of year [22], contributes to 25% of water pixels that were misclassified as dry in Figure 7. Another example is given in Figure 11, where the comparison between (a) and (b) shows that “Alt3 (OTSU)” mainly detects open water areas, and neglects the water presence under Salicornia salt marshes classified as flooded in the reference map.

During the study period, S2 satellites were passing over Camargue during late morning (between 10.10 and 10.40 CET). As a consequence, shadows, which are mainly observed in the town of Arles in the north of Camargue area, appear in S2 data. Longer shadows, which are observed from the end of autumn to early spring and mainly during winter, are misclassified as inundated areas to a larger degree with the reference map compared to the automatic thresholding alternatives. Therefore, shadow misclassification related to time of the satellite passage factor can also contribute to the reduction of agreement between the reference map and automatic thresholding approaches in the winter. Another example of disagreement is related to deep waters in a large lagoon that were misclassified as non-inundated areas by the reference map in November (14 November 2017) but not by “Alt2 (avg)” (Figure 12). This misclassification (see within the red square of Figure 12b) was attributed to the presence of waves caused by strong winds during the S2 image acquisition. Presumably, other sources of “false misclassifications” could arise from the reference maps of Doñana and Camargue because they were built from models that were not 100% accurate.

Another finding is related with the rice fields within the study areas. The criteria for mapping the water-vegetation subclass (see paragraph 2.4.3) as defined for Doñana are not satisfied for Camargue for any of the dates or any of the alternative approaches. On the contrary, for Doñana, the criteria are satisfied for the dates falling within summer (e.g., 11 July, 2017 and 20 August, 2017) and the water-vegetation subclass is mainly detected in the upper east area where rice fields are located. This result complies with the results in Reference [23] and is in agreement with the growing cycle of rice, since the rice grows during the summer period in Doñana and, at the same time, the paddies are flooded [59]. The difference between Camargue and Doñana is that the thresholds examined in the criteria for detecting water-vegetation subclass can be detected in the Inp for Alt1 and Alt2 and MNDVI histograms of Doñana, but not of Camargue (see Figure 4 and Figure 5). This likely relates to the land cover synthesis of each area. In particular, for the rice paddies of Camargue, it is found that there are other land cover types (e.g., reed marsh) with similar spectral behavior (i.e., similar Inp and MNDVI values) to the rice paddies. This is evident in a much smaller degree for Doñana. Moreover, the average MNDVI value of the Camargue rice paddies is lower than the average MNDVI value of the Doñana rice paddies and, thus, the discrimination of Camargue rice paddies from other vegetated areas in the MNDVI histogram is impeded. Due to the previously mentioned reasons, emergent rice cannot be detected in Camargue from the middle of July to early September.

Furthermore, by comparing Alt1, Alt2, and Alt3, it is evident that Alt2 and Alt1 histograms allow the detection of T_upper (see Figure 4b and Figure 13a), while, in the histogram of Alt3, T_upper cannot be detected (see Figure 13b). As a consequence, when comparing the inundation maps presented in Figure 14, it is evident that the water-vegetation subclass, corresponding mainly to rice fields located in the upper right part of the inundation map, cannot be detected when using the “Alt3 (MCET)” approach.

This study proves that automatic thresholding can be applied to more than one study areas and achieve high inundation mapping accuracy, without the need for simultaneous ground truth data or user’s intervention. Machine learning approaches that are developed based on ground truth data derived from a specific site [13,30,33,35,37] may perform very accurately as well. However, a credible performance for other sites cannot be safeguarded. Several studies have identified that surface reflectance accuracy of S2 L2A products may vary for different dates [60,61]. Hence, the performance of machine learning approaches could be negatively affected for dates that there is a notable variance in the surface reflectance accuracy compared to the dates, for which the training samples were derived in order to train the classification models.

5. Conclusions

This study examines automatic thresholding alternative approaches for separating inundated class pixels from non-inundated class pixels by utilizing atmospherically corrected S2 data. The experimental results show that alternative approaches are able to achieve high classification accuracy for Camargue and Doñana study areas. Out of the nine alternative approaches tested, three, seven, and four approaches were “Alm Perf” for Camargue, Doñana marshland, and Doñana complete area, respectively. “Alt2 (avg)” and “Alt3 (avg)” provided “Alm Perf” results for both Camargue and Doñana marshland, while “Alt2 (MCET)” and Alt3 (MCET)” provided the most consistent results for all areas, including the Doñana complete area. Thus, “Alt2 (avg)” and “Alt3 (avg)” are suggested for wetlands extensively covered by temporary flooded and emergent vegetation areas, such as the Camargue and the Doñana marshland, while “Alt2 (MCET)” and “Alt3 (MCET)” are expected to give more consistent results for wetlands including a large portion of dry areas, such as the Doñana complete area.

Future steps could consider the exploitation of ancillary information, such as digital elevation models to improve water detection under emergent vegetation, by inferring the water presence based on detected adjacent water covered areas having similar elevation, and land cover information to correct areas erroneously classified as water covered, where water presence is not expected. Furthermore, S2 inundation maps of a site generated via the automatic thresholding alternative achieving top accuracy among other alternatives can be fused with S1 data in order to allow for inundation mapping during extended cloudy periods, based on the example of Reference [12].