remote sensing

Article
Automatic Inundation Mapping Using Sentinel-2
Data Applicable to Both Camargue and Doñana
Biosphere Reserves
Georgios A. Kordelas 1 , Ioannis Manakos 1, * , Gaëtan Lefebvre 2 and Brigitte Poulin 2
1 Information Technologies Institute, Centre for Research and Technology Hellas, Charilaou-Thermi Rd. 6th km,
57001 Thessaloniki, Greece; kordelas@iti.gr
2 Tour du Valat Research Institute for the conservation of Mediterranean wetlands, Le Sambuc, 13200 Arles,
France; lefebvre@tourduvalat.org (G.L.); poulin@tourduvalat.org (B.P.)
* Correspondence: imanakos@iti.gr; Tel.: +30-2311-257760




Received: 19 July 2019; Accepted: 25 September 2019; Published: 27 September 2019


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.

Keywords: inundation mapping; flood mapping; automatic thresholding; Sentinel-2; wetlands;

marshland; Camargue; Doñana

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

Remote Sens. 2019, 11, 2251; doi:10.3390/rs11192251 www.mdpi.com/journal/remotesensing

Remote Sens. 2019, 11, 2251 2 of 20

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–19]. Commonly used indices include the Normalized Difference Water
Index (NDWI) [15,16,19], Modified NDWI [14,18–20], and Automated Water Extraction Index [17–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–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–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.,
Remote Sens. 2019, 11, 2251 3 of 20

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].
Remote Sens. 2019, 11, 2251 4 of 20
Remote Sens. 2019, 11, x FOR PEER REVIEW 4 of 20

Figure 1. Map
Figure of the
1. Map Camargue
of the Camargue Biosphere Reserve
Biosphere Reserve located
located in south
in south France,
France, partially
partially fallingthe
falling within within
the Occitanie
Occitanie andandProvence-Alpes-Côte
Provence-Alpes-Côte d’Azur
d’Azur administrative
administrative regions,regions, with underlying
with underlying S2 RGB
S2 RGB (red-
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 20
green-blue) image
(red-green-blue) imageonon 1818
April 2018.
April RedRed
2018. line:line:
boundary of theof
boundary Biosphere Reserve.Reserve.
the Biosphere

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].

Figure 2. Map of the Doñana Biosphere Reserve located in Southwest Spain, falling within the Andalusia
Figure 2. Map of the Doñana Biosphere Reserve located in Southwest Spain, falling within the
administrative
Andalusia region, with underlying
administrative S2 underlying
region, with RGB imageS2 onRGB
21 February
image on2018. Red line:
21 February boundary
2018. of the
Red line:
Biosphere Reserve. Yellow shaded area: Doñana marshland wetland area.
boundary of the Biosphere Reserve. Yellow shaded area: Doñana marshland wetland area.

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
Remote Sens. 2019, 11, 2251 5 of 20

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). 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.

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

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

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). 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.

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

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.
Remote Sens. 2019, 11, 2251 6 of 20

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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 20
The methodology’s steps are presented in Figure 3.

Figure 3. Schematic flow diagram of the automatic thresholding methodology.

Figure 3. Schematic flow diagram of the automatic thresholding methodology.
2.3.1. Segmentation of the Satellite Image
2.3.1. Segmentation of the Satellite Image
As the first step, the satellite image is segmented into non-overlapping regions by utilizing the
As thesegmentation
mean-shift first step, the satellite
algorithm image
[52].is The
segmented
resultinginto non-overlapping
segmentation mapregions by utilizing
is utilized the
for selecting
mean-shift
segments segmentation
with algorithm
a high percentage [52]. The resulting
of inundated pixels. 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
The input
(GREEN), to the4 segmentation
and Band (RED), which algorithm
have beenisnormalized
a false color
toimage composed
the range [0,255],ofby
Band 2 (BLUE),
relying Band
on minimum
3 (GREEN),
and maximum and Bandrepresenting
values 4 (RED), which have been
the intensity normalized
percentile to the
range fromrange
1% to[0,255],
99% per byband.
relying on
These
minimum and maximum values representing the intensity percentile range from 1% to
bands are selected due to their 10-m resolution, which allows for a more accurate segmentation of the 99% per band.
These bands
satellite image,arecompared
selected due to their 10-m
to selection resolution,
of other which
bands with a allows for a more [23].
lower resolution 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 Tinit 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
Remote Sens. 2019, 11, 2251 7 of 20

2.3.2. Mapping of the Open-Water Subclass

Three input data alternatives are examined for estimating an initial threshold Tinit 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].
Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 20
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
Square patches P of expanding window size are centered around each Cm. The window size
Alt3 is a variant of the Atl2 where SWIR-1 is used instead of SWIR-2. Each of Alt1, Alt2, and Alt3 is
in pixels is given by [20∙k × 20∙k], where k = 1,2,…,20. The splitting thresholds 𝑓 𝐤 of patches with
normalized in the range of [0,255], by relying on minimum and maximum values representing the
bimodal histogram (i.e., two distinctive classes of intensity values corresponding to inundated and
intensity percentile range from 1% to 99%. The term “Inp”, used in the following, is a generic term for
non-inundated classes appear in the histogram) are estimated based on the histogram thresholding
the input data,
algorithms (MCET corresponding to either Alt1,
or Otsu’s algorithm). Alt2,
Their or Alt3.
median is assumed to be the optimal threshold 𝑓 𝐦
The histogram of Inp is used to estimate
corresponding to the segment Gm. The usefulness of using an initial threshold separating
Tinit patches
expanding inundated
is to calculate fromthe
non-inundated pixels. T init is the Inp value for which the first
optimal threshold per segment more robustly, since its estimation is based on multiple deep valley of this histogram
splitting is
detected (Figure
thresholds. 4 shows
Then, the histograms
median ofusing Alt2 as
selected input). If optimal
segments’ has Inp(p) <isTinit
a pixel p thresholds , it is denoted
estimated as:
as inundated. Otherwise, m it is denoted as non-inundated. Therefore, an initial inundation map is
= median(
MoptFinal threshold fopt ) . final (see Figure 4), discriminating the open-water subclass, is estimated as
generated. Based
m =1,2,..., M on T this inundation map, segments Gm , where m = 1,2, . . . , M, having a large
max(Mopt,Tof
percentage init). Pixels p with Inp(p) < Tfinal are assumed to belong to the open-water subclass,
inundated pixels (i.e., over 70% of a segment’s pixels are inundated) are selected and
comprising
their of watercentroids
corresponding or water C withare
sparse vegetation pixels.
estimated.
m

(a) (b)
Figure4.4. Histogram
Figure Histogram of ofthe
theAlt2
Alt2map
map forfor
(a)(a)
Camargue
Camargue on on18 August,
18 August 2017, and and
2017, (b) Doñana
(b) Doñana on 20,on
August 2017. Tinit is the value for which the first deep valley in the Alt2 histogram is detected (black
20 August 2017. Tinit is the value for which the first deep valley in the Alt2 histogram is detected (black
point).The
point). The red
red point
point on on the
the Alt2
Alt2histogram
histogramcorresponds
corresponds totothethe
final Alt2
final threshold
Alt2 threshold TfinalTestimated by
final estimated
relying on MCET. The deep valley right after the first valley in the Alt2 histogram,
by relying on MCET. The deep valley right after the first valley in the Alt2 histogram, detected only detected only on
the
on right
the histogram,
right histogram, corresponds
corresponds to threshold
to threshold Tupper (green
Tupper point).
(green Histogram
point). valuesvalues
Histogram reach reach
up to up255to
normalized intensity value. However, since there are a few pixels with an Alt2
255 normalized intensity value. However, since there are a few pixels with an Alt2 value over 200 and value over 200 and
theydo
they donot
notcause
causefluctuations
fluctuationsofofthe
theAlt2
Alt2histogram
histogramcurve,
curve,the
theupper
upperlimit
limitofofthe
thex-axis
x-axisisissetsettoto200
200for
for visualization
visualization reasons.reasons.

When MCET
Square patchesorPOtsu
k is used for estimating the splitting thresholds, Tfinal is denoted as 𝐓𝐟𝐢𝐧𝐚𝐥 size
or
m of expanding window size are centered around each Cm . The window
𝐓 , respectively. Average (abbreviated as “avg”): 𝐓
pixels is given by [20·k × 20·k], where k = 1,2, . . . , 20.
in 𝐟𝐢𝐧𝐚𝐥 𝐓𝐟𝐢𝐧𝐚𝐥
𝐟𝐢𝐧𝐚𝐥The splitting 𝐓𝐟𝐢𝐧𝐚𝐥 /𝟐
thresholds is also
k examined as
fm of patches with
the final threshold.
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
2.3.3. Mapping of the Water-Vegetation Subclass m
algorithms (MCET or Otsu’s algorithm). Their median is assumed to be the optimal threshold fopt
corresponding to the segment Gm . 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: Mopt = median ( fopt m ).
m=1,2,...,M
Remote Sens. 2019, 11, 2251 8 of 20

Final threshold Tfinal (see Figure 4), discriminating the open-water subclass, is estimated as
max(Mopt ,Tinit ). Pixels p with Inp(p) < Tfinal 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, Tfinal is denoted as TMCET final
or
Otsu (a) avg

MCET (b)
Otsu
Tfinal , respectively. Average (abbreviated as “avg”): Tfinal = Tfinal + Tfinal /2 is also examined as
the finalFigure 4. Histogram of the Alt2 map for (a) Camargue on 18 August, 2017, and (b) Doñana on 20,
threshold.
August 2017. Tinit is the value for which the first deep valley in the Alt2 histogram is detected (black
2.3.3. Mapping
point). Theof the
red Water-Vegetation
point Subclass
on the Alt2 histogram corresponds to the final Alt2 threshold Tfinal estimated by
relying on MCET. The deep valley right after the first valley in the Alt2 histogram, detected only on
Inundated areas may be covered by emergent vegetation. The Inp value of the pixels in these areas
the right histogram, corresponds to threshold Tupper (green point). Histogram values reach up to 255
is higher compared to the Inp values of the pixels with open water or water with sparse vegetation
normalized intensity value. However, since there are a few pixels with an Alt2 value over 200 and
belonging to the open-water class. Therefore, a threshold Tupper , which is the Inp value for which the
they do not cause fluctuations of the Alt2 histogram curve, the upper limit of the x-axis is set to 200
deep valley right afterreasons.
for visualization 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𝐓𝐟𝐢𝐧𝐚𝐥 ,same
respectively.
the When MCET orAverage
time, areas with
Otsu (abbreviated
dense
is used as “avg”):
foremergent
estimating 𝐓𝐟𝐢𝐧𝐚𝐥 arethresholds,
vegetation
the splitting 𝐓
expected
𝐟𝐢𝐧𝐚𝐥 𝐓𝐟𝐢𝐧𝐚𝐥
to have
Tfinal /𝟐 is alsoMNDVI
high
is denoted examined as
as 𝐓𝐟𝐢𝐧𝐚𝐥value.
or
the final
Thus, threshold.
a threshold TMNDVI , which should be more than 0.4, is selected. In particular, TMNDVI is set
equal to the MNDVI value, for which the first valley in the part of the histogram with MNDVI values
2.3.3.
over 0.4Mapping of the
is detected (seeWater-Vegetation Subclass
Figure 5b) following the findings in Reference [23].

(a) (b)

Figure 5. Histogram of the MNDVI map for (a) Camargue on 18 August, 2017, and (b) Doñana on
20 August, 2017. TMNDVI is the MNDVI value for which the first deep valley after 0.4 is detected in the
MNDVI histogram (red point). TMNDVI is detected only on the right histogram.

A pixel p is assumed to belong to the water-vegetation subclass if the following two conditions
are met: Tfinal < Inp(p) < Tupper & MNDVI(p) > TMNDVI , presuming that Tupper and TMNDVI can
be detected in the histograms of Inp and MNDVI, respectively. These latter conditions are valid
for Doñana on 20 August, 2017 (see Figures 4b and 5b), but not for Camargue on 18 August, 2017
(see Figures 4a and 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
avg
final threshold (TMCET
final
or TOtsu
final
or Tfinal ) 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 TMCET final
), “OTSU” (standing for TOtsu
final
), and “avg”
avg
(standing for Tfinal ). 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)”.
 are met: Tfinal < Inp(p) < Tupper & MNDVI(p) > TMNDVI, presuming that Tupper and TMNDVI 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).
Remote Sens. 2019, 11, 2251 9 of 20
3. Results

3.1.Regarding
ComparisonCamargue,
of Automaticeach of the alternatives
Thresholding Results Against is Reference
applied to Mapthe 39 S2 products (see Table 1),
of Camargue
and the results are compared to the reference maps. Pixels corresponding to the coastal waters are not
The three input alternatives (Alt1 or Alt2, or Alt3) and the three different ways of estimating the
taken into consideration for the comparison to obtain unbiased classification results. The results are
final threshold (𝐓𝐟𝐢𝐧𝐚𝐥 or 𝐓𝐟𝐢𝐧𝐚𝐥 or 𝐓𝐟𝐢𝐧𝐚𝐥 ) form nine different alternatives, with each one combining
presented using Kappa coefficient, which is estimated in terms of the observed agreement (po ) and the
an input alternative and a way of estimating the final threshold. The alternative ways of estimating
agreement expected
final threshold are by chance
referred to (p
ase )“MCET”
as: Kappa= (po − for
(standing / (1-pe),) [53].
pe ) 𝐓 “OTSU” (standing for 𝐓𝐟𝐢𝐧𝐚𝐥 ), and
𝐟𝐢𝐧𝐚𝐥
The(standing
“avg” curves presented
for 𝐓𝐟𝐢𝐧𝐚𝐥 ).in
ForFigure
example, 6 show howisthe
if “Alt2” the Kappa
input and coefficient varies
“avg” is the wayfor the time period
of estimating the
between 12 June, 2017 and 19 June, 2018 for five out of
final threshold, then the alternative will be named as “Alt2 (avg)”. nine alternatives with the best agreement with
the reference map.Camargue,
Regarding The best agreement
each of the is given by “Alt3
alternatives (OTSU)”.
is applied to the“Alt3 (avg)” and
39 S2 products “Alt2
(see (avg)”
Table rank
1), and
second and third,
the results respectively,
are compared to theand have similar
reference results.
maps. Pixels Lastly, “Alt3to(MCET)”
corresponding the coastalandwaters
“Alt2are
(MCET)”
not
rank fourth
taken into and fifth, respectively,
consideration and give very
for the comparison closeunbiased
to obtain results. The previously
classification mentioned
results. ranking
The results are is
based on theusing
presented overall Kappa
Kappa accuracywhich
coefficient, provided in Section
is estimated in 3.3.
terms Forofall
thecurves,
observedit isagreement
seen that agreement
(po) and
the agreement
decreases expected
in the winter by chance
months from(pXIIe) as:
to Kappa= (po period,
II. For this – pe) / (1-p e) [53].
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
Figure 6. 6.Kappa
Kappacoefficient
coefficientvariation
variation of
of five
five alternatives
alternativesininthe
thetime
timeperiod
periodbetween
between 1212
June, 2017
June, 2017andand
19 June, 2018. The months’ numbers are given as Latin numerals and are colored
19 June, 2018. The months’ numbers are given as Latin numerals and are colored according to according to thethe
season they belong, i.e., green for summer, blue for autumn, black for winter, and pink for spring.
season they belong, i.e., green for summer, blue for autumn, black for winter, and pink for spring.

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.
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
Remote efficiently,
Sens. 2019, 11, 2251 such as the ones enclosed in the yellow squares, it tends to overestimate 10
water
of 20
presence in agricultural areas.

Figure 7. Inundation
Inundationmap mapon on2727February,
February,2018 using
2018 (a)(a)
using anan
S2 S2
derived inundation
derived inundationreference map,
reference (b)
map,
“Alt3
(b) (OTSU)”,
“Alt3 (c) “Alt3
(OTSU)”, (avg)”,
(c) “Alt3 and (d)
(avg)”, and“Atl3
(d)(MCET)”. Blue andBlue
“Atl3 (MCET)”. gray and
colors correspond
gray to inundated
colors correspond to
and non-inundated
inundated areas, respectively.
and non-inundated The green color
areas, respectively. Thecorresponds
green color to areas classified
corresponds as inundated
to areas classified
in (a) but non-inundated
as inundated in (b), (c), (d),
in (a) but non-inundated in while
(b–d), pink
whilecolor
pink corresponds to areas
color corresponds classified
to areas as non-
classified as
inundated in (a)inbut
non-inundated (a) inundated
but inundatedin (b),
in (c), (d). The
(b–d). The yellow
yellow colored squares in (a) enclose emergent
vegetated areas. The The upper
upper yellow
yellow square
square encloses
encloses reed
reed beds
beds and
and the
the lower one encloses short
helophytes, submerged macrophytes,
macrophytes, halophilous
halophilous scrubs,
scrubs, and
and annual
annual and
and perennial
perennial herbs.
herbs.

3.2. Comparison
3.2. Comparison of
of Automatic
Automatic and
and Thresholding
Thresholding Results
Results Against
Against Landsat
Landsat Reference
Reference Maps
Maps of
of Doñana
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).
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
Remote Sens.how Kappa
2019, 11, 2251 coefficient varies for seven different dates regarding the five top ranking
11 of 20
alternatives at the Doñana marshland (Figure 8a) and Doñana complete area (Figure 8b).

Figure 8. Kappa
Kappacoefficient
coefficientvariation
variationofoffive
fivealternatives forfor
alternatives seven different
seven dates
different with
dates respect
with to (a)
respect to
Doñana
(a) marshland
Doñana andand
marshland (b)(b)
Doñana complete
Doñana complete area.
area.

In Doñana marshland, “Alt1 (OTSU)” provides the highest highest accuracy,

accuracy, followed by “Alt2 “Alt2 (MCET)”
(MCET)”
and “Alt1 (Avg)”
(Avg)” (see
(see Section
section 3.3
3.3 for
for more
more details
details on
on the
the ranking
ranking that
that is
is based
basedononoverall
overallKappa).
Kappa).
In the Doñana
Doñana complete area, “Alt3 (MCET)” provides the highest accuracy, followed by “Alt1
(MCET)”
(MCET)” andand “Alt2
“Alt2 (MCET)”
(MCET)” (see
(see Section
section 3.3 for more details on the ranking that is based on overall
Kappa). “Alt1 (avg)” provides the most consistent Kappa values values across
across the
the study
study period.
period.
Figure
Figure 9 shows an example of the Landsat derived inundation reference map onJanuary,
9 shows an example of the Landsat derived inundation reference map on 27 20182018
27 January, and
the
andinundation mapsmaps
the inundation estimated basedbased
estimated on theonalternative approaches
the alternative usingusing
approaches Alt3 input. GreenGreen
Alt3 input. and pink
and
colors indicate
pink colors the classification
indicate differences
the classification between
differences the reference
between map and
the reference map the alternative
and maps. maps.
the alternative “Alt3
(MCET)” (Figure(Figure
“Alt3 (MCET)” 9d) agrees
9d) well
agreeswith the with
well reference map (Figure
the reference map9a), while “Alt3
(Figure (OTSU)”
9a), while (Figure
“Alt3 9b)
(OTSU)”
largely
(Figureoverestimates the water presence
9b) largely overestimates the waterand “Alt3 (avg)”
presence (Figure
and “Alt3 9c) shows
(avg)” (Figureintermediate results.
9c) shows intermediate
results.
 Remote
Remote Sens.
Sens. 2019,
2019, 11,11,2251
x FOR PEER REVIEW 1212ofof20
20

Figure9.9.(a)(a)Landsat
Figure Landsat derived
derivedinundation
inundation map
mapon on
27 January,
27 January,2018. Inundation
2018. mapmap
Inundation on 27
onJanuary, 2018
27 January,
using
2018 (b) “Alt3
using b)(OTSU)”, (c) “Alt3(c)
“Alt3 (OTSU)”, (avg)”,
“Alt3and (d) “Atl3
(avg)”, and (MCET)”. Blue and gray
(d) “Atl3 (MCET)”. Bluecolors correspond
and gray colors
tocorrespond
inundated to and non-inundated
inundated areas, respectively.
and non-inundated Green color corresponds
areas, respectively. to areas classified
Green color corresponds as
to areas
inundated
classified inas (a) but non-inundated
inundated in (b–d), while
in (a) but non-inundated in the
(b),pink color
(c), (d), corresponds
while to areas
the pink color classified to
corresponds as
non-inundated
areas classifiedinas (a)non-inundated
but inundatedin in (a)
(b–d).
but The yellow in
inundated polygon
(b), (c),encloses
(d). Themarshland area. encloses
yellow polygon
marshland area.
3.3. Overall Kappa Per Approach and Examined Areas
3.3.Figure
Overall10
Kappa Per Approach
provides overalland Examined
Kappa Areas which are estimated when the number of true
coefficients,
positive pixels
Figure 10 of the inundated
provides overall class,
Kappa false positive pixels
coefficients, whichofaretheestimated
inundatedwhen
class,the
truenumber
positiveofpixels
true
ofpositive
the non-inundated class and false positive pixels of the non-inundated class are added
pixels of the inundated class, false positive pixels of the inundated class, true positive forpixels
dates
per
of approach and examined
the non-inundated class study.
and falseThepositive
combined Overall
pixels of theAccuracy of eachclass
non-inundated approach per study
are added case
for dates
isper
alsoapproach
given asand
a number close
examined the point
study. corresponding
The combined OveralltoAccuracy
its overall
of Kappa. Kappaper
each approach values
studycan be
case
classified according to Reference [53] as follows: “Moderate” when 0.40 < Kappa ≤
is also given as a number close the point corresponding to its overall Kappa. Kappa values can be0.60 (shortened
asclassified
“Mod”),according
“substantial” when 0.60
to Reference < as
[53] Kappa ≤ 0.80
follows: (shortened
“Moderate” as “Sub”)
when or “almost
0.40 < Kappa ≤ 0.60perfect” when
(shortened as
0.80 < Kappa ≤ 1 shortened as (“Alm Perf”).
 Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 20

Remote“Mod”), “substantial”
Sens. 2019, 11, 2251 when 0.60 < Kappa ≤ 0.80 (shortened as “Sub”) or “almost perfect” when 0.8013 of 20
< Kappa ≤ 1 shortened as (“Alm Perf”).

Figure
Figure 10. 10. Overall
Overall Kappaper
Kappa perapproach
approach and
and examined
examinedstudy
studyarea. TheThe
area. dotted line line
dotted delineates the the
delineates
threshold of the “Alm Perf” case. Yellow arrows indicate alternative approaches succeeding
threshold of the “Alm Perf” case. Yellow arrows indicate alternative approaches succeeding “Alm Perf” “Alm
Perf” results in Camargue and Doñana marshland at the same time, while blue arrows showcase the
results in Camargue and Doñana marshland at the same time, while blue arrows showcase the most
most successful overall alternative ones. The combined overall accuracy of each approach per study
successful overall alternative ones. The combined overall accuracy of each approach per study area is
area is given as a number close the point corresponding to its overall Kappa.
given as a number close the point corresponding to its overall Kappa.
“Alt2 (avg)” and “Alt3 (avg)” (see yellow arrows in Figure 10) have both “Alm Perf” overall
“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
Kappa whenarea,
complete considering Camargue
their Kappa and overall
is “Sub”, since Doñana marshland
Kappa for “Alt2study
(avg)” areas.
and “Alt3While,
(avg)”for theand
is 0.63 Doñana
complete area, their Kappa
0.73, respectively. is “Sub”,
On the other hand,since
“Alt2 overall
(MCET)”Kappa for (MCET)”
and “Alt3 “Alt2 (avg)” andarrows
(see blue “Alt3 in
(avg)”
Figureis 0.63
and 10)
0.73, respectively.
exhibit On the
more consistent other hand,
performance “Alt2
across (MCET)”
sites, since theirand “Alt3
Kappa (MCET)”
is “Alm (see
Perf” for blue
both arrows in
Doñana
marshland and Doñana complete area, while, for Camargue, their overall Kappa is very close
Figure 10) exhibit more consistent performance across sites, since their Kappa is “Alm Perf” for both to “Alm
Perf.”
Doñana marshland and Doñana complete area, while, for Camargue, their overall Kappa is very close
to “Alm Perf.”
4. Discussion
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
The aim of this paper is to examine the performance of various alternative automatic thresholding
Camargue and Doñana. These are considered as testbeds to evaluate the transferability of the
approaches for inundation mapping in two different wetland areas, which include Camargue and
different approaches and select the ones that favor the acquisition of consistently credible results for
Doñana.
both These are In
wetlands. considered
this study,asthe
testbeds to evaluate alternative
local thresholding the transferability
approaches of the different
utilize approaches
multispectral
and select
satellite data, while the vast majority of local thresholding approaches have been designed forIn this
the ones that favor the acquisition of consistently credible results for both wetlands.
study, the local
utilizing thresholding
radar alternativeMoreover,
data (e.g., [27,54–58]). approaches the utilize multispectral
estimation satelliterelies
of local thresholds data,onwhile the vast
patches
majority of local thresholding
of expanding size, contrary approaches
to most of thehave been designed
approaches for utilizing
using image subsets ofradar
fixed data (e.g., [27,54–58]).
size [55–58]. This
helps to
Moreover, increase
the robustness,
estimation of localwhen input data
thresholds have
relies ondifferent
patches spatial resolution
of expanding such
size, as when
contrary to inputs
most of the
from different satellites are used.
approaches using image subsets of fixed size [55–58]. This helps to increase robustness, when input
data haveThe overall spatial
different Kappa of the alternative
resolution suchapproaches varies from
as when inputs from0.73 to 0.84 satellites
different for Camargue, from 0.74
are used.
to 0.94 for Doñana marshland, and from 0.45 to 0.85 for the complete area of Doñana. The
The overall Kappa of the alternative approaches varies from 0.73 to 0.84 for Camargue, from 0.74
experimental results demonstrate that there is not a common approach across all study cases
to 0.94 for Doñana marshland, and from 0.45 to 0.85 for the complete area of Doñana. The experimental
achieving absolute top accuracy. “Alt3 (OTSU)” is best for Camargue, “Alt1 (OTSU)” is best for
results demonstrate
Doñana marshland, that
andthere
“Alt3is(MCET)”
not a common approach
is best for Doñana as across all The
a whole. studyusecases achieving
of “Alt2” absolute
and “Alt3”
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
Remote Sens. 2019, 11, 2251 14 of 20

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.
 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
RemoteRemote
water
Sens. 2019,
Sens. 2019, 11, x FOR PEER REVIEW
11, 2251
presence under Salicornia salt marshes classified as flooded in the reference map.
15 of 20
15 of 20

Figure 11. Zoomed part of the inundation map generated on 23 January, 2018 using: (a) “Alt3
(OTSU)”, and (b) reference map. In this region of the Camargue, the main habitat is halophilous
scrubs (Salicornia marshes). Areas of open water and bare ground are clearly visible on the Google
Earth image of the same area acquired on 22 April, 2018 shown in (c).

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
Figure 11.(avg)”
Zoomed(Figure
part of12).
the This misclassification
inundation (see on
map generated within the red2018
23 January, square of (a)
using: Figure
“Alt312b) was
(OTSU)”,
attributed
and to themap.
(b) reference presence
In thisof waves
region caused
of the by strong
Camargue, winds
the main during
habitat the S2 image
is halophilous scrubsacquisition.
(Salicornia
Presumably,
marshes). other
Areas ofsources of “false
open water andmisclassifications” could visible
bare ground are clearly arise from the Google
on the reference maps
Earth of Doñana
image of the
and Camargue because they were built from
same area acquired on 22 April, 2018 shown in (c). models that were not 100% accurate.

Figure
Figure 12. Inundation
12. Inundation mapmap generatedon
generated on1414 November,
November, 2017
2017using:
using:(a)(a)
“Alt2 (agv)”,
“Alt2 and and
(agv)”, (b) reference
(b) reference
map. Blue and gray colors correspond to inundated and non-inundated areas, respectively.
map. Blue and gray colors correspond to inundated and non-inundated areas, respectively.

Another
Another finding
finding is related
is related withwith
thethe rice
rice fields
fields withinthe
within thestudy
studyareas.
areas. The
The criteria
criteriafor
formapping
mapping the
the water-vegetation subclass (see paragraph 2.4.3) as defined for Doñana are not satisfied for
water-vegetation subclass (see paragraph 2.4.3) as defined for Doñana are not satisfied for Camargue
Camargue for any of the dates or any of the alternative approaches. On the contrary, for Doñana, the
for any of the dates or any of the alternative approaches. On the contrary, for Doñana, the criteria
criteria are satisfied for the dates falling within summer (e.g., 11 July, 2017 and 20 August, 2017) and
are satisfied for the dates
the water-vegetation falling
subclass is within summerin(e.g.,
mainly detected 11 July,
the upper east 2017 and 20
area where August,
rice 2017)
fields are and the
located.
water-vegetation
This result complies with the results in Reference [23] and is in agreement with the growing cycle of This
subclass is mainly detected in the upper east area where rice fields are located.
resultrice,
complies
since thewith the results
rice grows duringinthe
Reference [23] and
summer period is in agreement
in Doñana and, at thewith
samethe growing
time, cycleare
the paddies of rice,
sinceflooded
the rice[59].
grows during the
The difference summer
between period and
Camargue in Doñana
Doñana is and,
thatat the
the same time,
thresholds the paddies
examined in the are
criteria
flooded [59].forThe
detecting water-vegetation
difference between Camargue subclass andcan be detected
Doñana in thethe
is that Inp for Alt1 and
thresholds Alt2 andin the
examined
MNDVI histograms of Doñana, but not of Camargue (see Figure 4 and Figure
criteria for detecting water-vegetation subclass can be detected in the Inp for Alt1 and Alt2 and 5). This likely relates
MNDVI
histograms of Doñana, but not of Camargue (see Figures 4 and 5). This likely relates to the found
to the land cover synthesis of each area. In particular, for the rice paddies of Camargue, it is land cover
that there are other land cover types (e.g., reed marsh) with similar spectral behavior (i.e., similar Inp
synthesis of each area. In particular, for the rice paddies of Camargue, it is found that there are other
and MNDVI values) to the rice paddies. This is evident in a much smaller degree for Doñana.
land cover types (e.g., reed marsh) with similar spectral behavior (i.e., similar Inp and MNDVI values)
Moreover, the average MNDVI value of the Camargue rice paddies is lower than the average MNDVI
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.
 Remote Sens. 2019, 11, x FOR PEER REVIEW 16 of 20
Remote Sens. 2019, 11, x FOR PEER REVIEW 16 of 20
Remote Sens. 2019, 11, 2251 16 of 20
value of the Doñana rice paddies and, thus, the discrimination of Camargue rice paddies from other
value of the
vegetated Doñana
areas in the rice paddieshistogram
MNDVI and, thus,isthe discrimination
impeded. Due tooftheCamargue
previouslyricementioned
paddies from other
reasons,
vegetatedrice
emergent areas in the
cannot be MNDVI
detected histogram
in Camargue is impeded.
from Due toofthe
the middle previously
July to early mentioned reasons,
Furthermore, by comparing Alt1, Alt2, and Alt3, it is evident that Alt2September.
and Alt1 histograms
emergent rice cannot
Furthermore, by be detectedAlt1,
comparing in Camargue
Alt2, and from
Alt3, the
it is middle
evident of
thatJuly
Alt2to and
earlyAlt1
September.
histograms
allow the detection of Tupper (see Figures 4b and 13a), while, in the histogram of Alt3, Tupperallow
cannot
the Furthermore,
detection of T by comparing Alt1, Alt2, and Alt3, it is evident that Alt2 and Alt1 histograms
upper (see Figure 4b and Figure 13a), while, in the histogram of Alt3, Tupper cannot be
allow
be detected (see Figure 13b). As a consequence, when comparing the inundation maps presented
the detection of Tupper13b).
(see Figure 4b and Figurewhen 13a), while, in thethe histogram of Alt3, Tupper cannot inbe
in detected
Figure 14, (seeit Figure
is evident that As atheconsequence,
water-vegetationcomparing inundation
subclass, corresponding maps
mainly presented
to rice fields
detected (see Figure 13b). As a consequence, when comparing the inundation
Figure 14, it is evident that the water-vegetation subclass, corresponding mainly to rice fields located maps presented in
located
Figure in the upper
14, it isright
evident right part of the inundation map, cannot be detected when using the “Alt3
in the upper partthat theinundation
of the water-vegetation subclass,
map, cannot becorresponding
detected whenmainly using to
therice fields
“Alt3 located
(MCET)”
(MCET)”
in the approach.
upper right part of the inundation map, cannot be detected when using the “Alt3 (MCET)”
approach.
approach.

(a) (b)
(a) (b)
Figure
Figure 13.13.Histogram
Histogramofofthe the(a)
(a)Alt1
Alt1and
and (b)
(b) Alt3
Alt3 maps
maps forfor Doñana
Doñanaon on2020August,
August,2017.2017.Tinit is the
Tinit value
is the value
Figure
for which 13.the
Histogram
first deep ofvalley
the (a)inAlt1
the and
Alt1(b)
andAlt3
Alt3 maps for Doñana
histograms is on 20 August,
detected (black 2017. TThe
point). init is the value
red point
for which the first deep valley in the Alt1 and Alt3 histograms is detected (black point). The red point
forAlt1
on which and the firsthistograms
Alt3 deep valleycorresponds
in the Alt1 and
to TAlt3 histograms is detected (black point). The red point
final estimated relying on MCET. The deep valley right
on Alt1 and Alt3 histograms corresponds to Tfinal estimated relying on MCET. The deep valley right
on Alt1 and Alt3 histograms corresponds to T final estimated relying on MCET. The deep valley right
after the first valley in the Alt1 histogram corresponds to threshold Tupper. Tupper cannot be detected on
after the first valley in the Alt1 histogram corresponds to threshold Tupper . Tupper cannot be detected
after
the thehistogram.
Alt3 first valley (Histograms’
in the Alt1 histogram corresponds
values reach up to 255. to threshold upper. Tupper cannot be detected on
However, Tsince there are a few pixels with
on the
the Alt3
Alt3 histogram.
histogram.(Histograms’
(Histograms’valuesvaluesreach
reachupuptoto255.255.However,
However, since
since there
there areare a few
a few pixels with
value over 200 that do not cause fluctuations of the histograms’ curves, the upper limit ofpixels with
the x-axis
value over 200 that do not cause fluctuations
fluctuations of the histograms’ curves, the upper limit of the x-axis is
of the histograms’ curves, the upper limit of the x-axis
isvalue
set toover
200 200 that do not cause
for visualization reasons.).
setistoset
200 tofor
200visualization
for visualizationreasons.).
reasons.).

Figure 14. Inundation map generated on 20 August, 2017 using: (a) “Alt1 (MCET)”, (b) “Alt2
Figure
Figure 14. 14.
(MCET)”, andInundation
Inundation map
(c) “Alt3 map
(MCET)”generated
generated onAugust,
on 20
alternatives.20 August,
Blue 2017
and 2017
using:
gray using: (a) (MCET)”,
“Alt1
(a)correspond
colors “Alt1 to (MCET)”, (b)
(b) “Alt2
inundated and “Alt2
(MCET)”,
non-
(MCET)”,
inundated
and and
(c) “Alt3areas, (c)
(MCET)” “Alt3 (MCET)”
respectively. alternatives. Blue and gray colors correspond to inundated
alternatives. Blue and gray colors correspond to inundated and non-inundatedand non-
inundated areas,
areas, respectively. respectively.
This study proves that automatic thresholding can be applied to more than one study areas and
This
achieveThis study
study
high provesthat
proves
inundation that automatic
automatic
mapping thresholding
thresholding
accuracy, can
the be
withoutcan be applied
applied
need to
tomore
morethan
for simultaneousthan one
onestudy
ground study
truthareas
dataand
areas and
or
achieve
user’s high
achieve high inundation
intervention. Machine
inundation mapping
mapping accuracy,
learning without
approaches
accuracy, withoutthatthe need
theare
need for simultaneous
developed ground
based on ground
for simultaneous truth
ground truth data or or
data
user’s
derived
user’s intervention. Machine
from a specific
intervention. learning approaches
approaches
site [13,30,33,35,37]
Machine learning that
may perform are
aredeveloped
thatvery accurately based
developed well.onHowever,
asbased onground
ground truth data
a credible
truth data
derived from a specific site [13,30,33,35,37] may perform very accurately as well. However, a credible
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
Remote Sens. 2019, 11, 2251 17 of 20

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].

Author Contributions: G.A.K. and I.M. conceived and implemented the alternative approaches in this work,
which were fine-tuned based on suggestions from G.L. and B.P. G.A.K. performed an accuracy assessment of the
generated Sentinel-2 inundation maps. Reference maps for Camargue were provided by G.L. and B.P. Results
analysis was carried out by all authors. G.A.K. and I.M. mainly wrote the paper, with significant contribution
of G.L. and B.P. to the writing of the discussion and the description of the Camargue study area. All authors
reviewed, suggested improvements, and approved the manuscript.
Funding: The European Union’s Horizon 2020 research and innovation program under grant agreement No 641762,
ECOPOTENTIAL supported and funded this study.
Acknowledgments: We are grateful to Javier Bustamante, Ricardo Díaz-Delgado, and David Aragonés, who
are affiliated with the Remote Sensing and GIS Laboratory (LAST-EBD), Estación Biologica de Doñana (CSIC),
for kindly providing the reference Landsat inundation maps for Doñana and to Loïc Willm from Tour du Valat
Research Institute for providing data and maps from the Camargue.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. World Resources Institute. Millennium Ecosystem Assessment: Ecosystems and Human Well-Being: Wetlands and
Water Synthesis; World Resources Institute: Washington, DC, USA, 2005.
2. Park, E.; Lee, S.; Peters, D.J. Iowa wetlands outdoor recreation visitors’ decision-making process: An extended
model of goal-directed behavior. J. Outdoor Recreat. Tour. 2017, 17, 64–76. [CrossRef]
3. Davidson, N.C. How much wetland has the world lost? Long-term and recent trends in global wetland area.
Mar. Freshw. Res. 2014, 65, 934–941. [CrossRef]
4. Huang, W.; DeVries, B.; Huang, C.; Lang, M.W.; Jones, J.W.; Creed, I.F.; Carroll, M.L. Automated Extraction
of Surface Water Extent from Sentinel-1 Data. Remote Sens. 2018, 10, 797. [CrossRef]
5. Töyrä, J.; Pietroniro, A. Towards operational monitoring of a northern wetland using geomatics-based
techniques. Remote Sens. Environ. 2005, 97, 174–191. [CrossRef]
6. Marti-Cardona, B.; Dolz-Ripolles, J.; Lopez-Martinez, C. Wetland inundation monitoring by the synergistic use
of ENVISAT/ASAR imagery and ancilliary spatial data. Remote Sens. Environ. 2013, 139, 171–184. [CrossRef]
Remote Sens. 2019, 11, 2251 18 of 20

7. Martinis, S.; Plank, S.; Ćwik, K. The use of Sentinel-1 time-series data to improve flood monitoring in arid
areas. Remote Sens. 2018, 10, 583. [CrossRef]
8. Giustarini, L.; Hostache, R.; Matgen, P.; Schumann, G.J.P.; Bates, P.D.; Mason, D.C. A change detection
approach to flood mapping in urban areas using TerraSAR-X. IEEE Trans. Geosci. Remote Sens. 2012, 51,
2417–2430. [CrossRef]
9. Schumann, G.; Di Baldassarre, G.; Bates, P.D. The utility of spaceborne radar to render flood inundation
maps based on multialgorithm ensembles. IEEE Trans. Geosci. Remote Sens. 2009, 47, 2801–2807. [CrossRef]
10. Shen, X.; Wang, D.; Mao, K.; Anagnostou, E.; Hong, Y. Inundation Extent Mapping by Synthetic Aperture
Radar: A Review. Remote Sens. 2019, 11, 879. [CrossRef]
11. Markert, K.N.; Chishtie, F.; Anderson, E.R.; Saah, D.; Griffin, R.E. On the merging of optical and SAR satellite
imagery for surface water mapping applications. Results Phys. 2018, 9, 275–277. [CrossRef]
12. Manakos, I.; Kordelas, G.A.; Marini, K. Fusion of Sentinel-1 data with Sentinel-2 products to overcome
non-favourable atmospheric conditions for the delineation of inundation maps. Eur. J. Remote Sens. 2019,
1–14. [CrossRef]
13. Davranche, A.; Poulin, B.; Lefebvre, G. Mapping flooding regimes in Camargue wetlands using seasonal
multispectral data. Remote Sens. Environ. 2013, 138, 165–171. [CrossRef]
14. Donchyts, G.; Schellekens, J.; Winsemius, H.; Eisemann, E.; van de Giesen, N. A 30 m resolution surface water
mask including estimation of positional and thematic differences using landsat 8, srtm and openstreetmap:
A case study in the Murray-Darling Basin, Australia. Remote Sens. 2016, 8, 386. [CrossRef]
15. Du, Y.; Zhang, Y.; Ling, F.; Wang, Q.; Li, W.; Li, X. Water bodies’ mapping from Sentinel-2 imagery with
modified normalized difference water index at 10-m spatial resolution produced by sharpening the SWIR
band. Remote Sens. 2016, 8, 354. [CrossRef]
16. Kyriou, A.; Nikolakopoulos, K. Flood mapping from Sentinel-1 and Landsat-8 data: A case study from river
Evros, Greece. In Proc. SPIE 9644, Earth Resources and Environmental Remote Sensing/GIS Applications VI;
International Society for Optics and Photonics: Toulouse, France, 2015; p. 964405. [CrossRef]
17. Guo, Q.; Pu, R.; Li, J.; Cheng, J. A weighted normalized difference water index for water extraction using
Landsat imagery. Int. J. Remote Sens. 2017, 38, 5430–5445. [CrossRef]
18. Buma, W.; Lee, S.I.; Seo, J. Recent Surface Water Extent of Lake Chad from Multispectral Sensors and GRACE.
Sensors 2018, 18, 2082. [CrossRef]
19. Zhang, F.; Li, J.; Zhang, B.; Shen, Q.; Ye, H.; Wang, S.; Lu, Z. A simple automated dynamic threshold extraction
method for the classification of large water bodies from landsat-8 OLI water index images. Int. J. Remote
Sens. 2018, 39, 3429–3451. [CrossRef]
20. Lee, K.S.; Kim, T.H.; Yun, Y.S.; Shin, S.M. Spectral characteristics of shallow turbid water near the shoreline
on inter-tidal flat. Korean J. Remote Sens. 2001, 17, 131–139.
21. Feyisa, G.L.; Meilby, H.; Fensholt, R.; Proud, S.R. Automated Water Extraction Index: A new technique for
surface water mapping using Landsat imagery. Remote Sens. Environ. 2014, 140, 23–35. [CrossRef]
22. Lefebvre, G.; Davranche, A.; Willm, L.; Campagna, J.; Redmond, L.; Merle, C.; Guelmami, A.; Poulin, B.
Introducing WIW for detecting the presence of Water In Wetlands with Landsat and Sentinel satellites.
Remote Sens. 2019, 11, 2210. [CrossRef]
23. Kordelas, G.A.; Manakos, I.; Aragonés, D.; Díaz-Delgado, R.; Bustamante, J. Fast and Automatic Data-Driven
Thresholding for Inundation Mapping with Sentinel-2 Data. Remote Sens. 2018, 10, 910. [CrossRef]
24. Díaz-Delgado, R.; Aragonés, D.; Afán, I.; Bustamante, J. Long-term monitoring of the flooding regime and
hydroperiod of Doñana marshes with Landsat time series (1974–2014). Remote Sens. 2016, 8, 775. [CrossRef]
25. Bangira, T.; Alfieri, S.M.; Menenti, M.; van Niekerk, A. Comparing Thresholding with Machine Learning
Classifiers for Mapping Complex Water. Remote Sens. 2019, 11, 1351. [CrossRef]
26. Ban, H.J.; Kwon, Y.J.; Shin, H.; Ryu, H.S.; Hong, S. Flood monitoring using satellite-based RGB composite
imagery and refractive index retrieval in visible and near-infrared bands. Remote Sens. 2017, 9, 313. [CrossRef]
27. Chini, M.; Hostache, R.; Giustarini, L.; Matgen, P. A Hierarchical Split-Based Approach for Parametric
Thresholding of SAR Images: Flood Inundation as a Test Case. IEEE Trans. Geosci. Remote Sens. 2017, 55,
6975–6988. [CrossRef]
28. Otsu, N. A threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern. 1979, 9,
62–66. [CrossRef]
29. Kittler, J.; Illingworth, J. Minimum error thresholding. Pattern Recognit. 1986, 19, 41–47. [CrossRef]
Remote Sens. 2019, 11, 2251 19 of 20

30. Ko, B.C.; Kim, H.H.; Nam, J.Y. Classification of potential water bodies using Landsat 8 OLI and a combination
of two boosted random forest classifiers. Sensors 2015, 15, 13763–13777. [CrossRef]
31. DeVries, B.; Huang, C.; Lang, M.; Jones, J.; Huang, W.; Creed, I.; Carroll, M. Automated quantification of
surface water inundation in wetlands using optical satellite imagery. Remote Sens. 2017, 9, 807. [CrossRef]
32. Skakun, S. A neural network approach to flood mapping using satellite imagery. Comput. Inform. 2012, 29,
1013–1024.
33. Nandi, I.; Srivastava, P.K.; Shah, K. Floodplain mapping through support vector machine and optical/infrared
images from landsat 8 OLI/TIRS sensors: Case study from Varanasi. Water Resour. Manag. 2017, 31, 1157–1171.
[CrossRef]
34. Jakovljević, G.; Govedarica, M.; Álvarez-Taboada, F. Waterbody mapping: A comparison of remotely sensed
and GIS open data sources. Int. J. Remote Sens. 2019, 40, 2936–2964. [CrossRef]
35. Miao, Z.; Fu, K.; Sun, H.; Sun, X.; Yan, M. Automatic water-body segmentation from high-resolution satellite
images via deep networks. IEEE Geosci. Remote Sens. Lett. 2018, 15, 602–606. [CrossRef]
36. Isikdogan, F.; Bovik, A.C.; Passalacqua, P. Surface water mapping by deep learning. IEEE J. Sel. Top. Appl.
Earth Obs. Remote Sens. 2017, 10, 4909–4918. [CrossRef]
37. Acharya, T.; Lee, D.; Yang, I.; Lee, J. Identification of water bodies in a Landsat 8 OLI image using a J48
decision tree. Sensors 2016, 16, 1075. [CrossRef]
38. Yousefi, P.; Jalab, H.A.; Ibrahim, R.W.; Noor, N.F.M.; Ayub, M.N.; Gani, A. Water-body segmentation in
satellite imagery applying modified Kernel K-means. Malays. J. Comput. Sci. 2018, 31, 143–154. [CrossRef]
39. Sivanpillai, R.; Miller, S.N. Improvements in mapping water bodies using ASTER data. Ecol. Inform. 2010, 5,
73–78. [CrossRef]
40. Jiang, H.; Feng, M.; Zhu, Y.; Lu, N.; Huang, J.; Xiao, T. An automated method for extracting rivers and lakes
from Landsat imagery. Remote Sens. 2014, 6, 5067–5089. [CrossRef]
41. Zhou, Y.; Dong, J.; Xiao, X.; Xiao, T.; Yang, Z.; Zhao, G.; Zou, Z.; Qin, Y. Open Surface Water Mapping
Algorithms: A Comparison of Water-Related Spectral Indices and Sensors. Water 2017, 9, 256. [CrossRef]
42. Haibo, Y.; Zongmin, W.; Hongling, Z.; Yu, G. Water body extraction methods study based on RS and GIS.
Procedia Environ. Sci. 2011, 10, 2619–2624. [CrossRef]
43. Schwatke, C.; Scherer, D.; Dettmering, D. Automated Extraction of Consistent Time-Variable Water Surfaces
of Lakes and Reservoirs Based on Landsat and Sentinel-2. Remote Sens. 2019, 11, 1010. [CrossRef]
44. Chauvelon, P.; Tournoud, M.G.; Sandoz, A. Integrated hydrological modelling of a managed coastal
Mediterranean wetland (Rhone delta, France): Initial calibration. Hydrol. Earth Syst. Sci. 2003, 7, 123–132.
[CrossRef]
45. Grillas, P. The Camargue: Rhone River Delta (France). In The Wetland Book: II: Distribution, Description and
Conservation; Finlayson, C., Milton, G., Prentice, R., Davidson, N., Eds.; Springer: Dordrecht, The Netherlands, 2016.
46. Lefebvre, G.; Germain, C.; Poulin, B. Contribution of rainfall vs. water management to Mediterranean
wetland hydrology: Development of an interactive simulation tool to foster adaptation to climate variability.
Environ. Model. Softw. 2015, 74, 39–47. [CrossRef]
47. Lefebvre, G.; Redmond, L.; Germain, C.; Palazzi, E.; Terzago, S.; Willm, L.; Poulin, B. Predicting the
vulnerability of seasonally-flooded wetlands to climate change across the Mediterranean Basin. Sci. Total
Environ. 2019, 692, 546–555. [CrossRef] [PubMed]
48. Martí-cardona, B. Spaceborne SAR Imagery for Monitoring the Inundation in the Doñana Wetlands.
Ph.D. Dissertation, Polytechnic University of Catalonia, Barcelona, Catalonia, Spain, 2014.
49. Green, A.J.; Bustamante, J.; Janss, G.F.E.; Fernández-Zamudio, R.; Díaz-Paniagua, C. Doñana Wetlands
(Spain). In The Wetland Book: II: Distribution, Description and Conservation; Finlayson, C., Milton, G., Prentice, R.,
Davidson, N., Eds.; Springer: Dordrecht, The Netherlands, 2016.
50. Kloskowski, J.; Green, A.J.; Polak, M.; Bustamante, J.; Krogulec, J. Complementary use of natural and artificial
wetlands by waterbirds wintering in Doñana, south-west Spain. Aquat. Conserv. Mar. Freshw. Ecosyst. 2009,
19, 815–826. [CrossRef]
51. Kullback, S. Information Theory and Statistics; Dover Publications: New York, NY, USA, 1968.
52. Comaniciu, D.; Meer, P. Mean shift: A robust approach toward feature space analysis. IEEE Trans. Pattern
Anal. Mach. Intell. 2002, 24, 603–619. [CrossRef]
53. Landis, J.R.; Koch, G.G. The measurement of observer agreement for categorical data. Biometrics 1977, 33,
159–174. [CrossRef]
Remote Sens. 2019, 11, 2251 20 of 20

54. Nakmuenwai, P.; Yamazaki, F.; Liu, W. Automated Extraction of Inundated Areas from Multi-Temporal
Dual-Polarization RADARSAT-2 Images of the 2011 Central Thailand Flood. Remote Sens. 2017, 9, 78.
[CrossRef]
55. Martinis, S.; Kuenzer, C.; Wendleder, A.; Huth, J.; Twele, A.; Roth, A.; Dech, S. Comparing four operational
SAR-based water and flood detection approaches. Int. J. Remote Sens. 2015, 36, 3519–3543. [CrossRef]
56. Li, J.; Wang, S. An automatic method for mapping inland surface waterbodies with Radarsat-2 imagery. Int. J.
Remote Sens. 2015, 36, 1367–1384. [CrossRef]
57. Bioresita, F.; Puissant, A.; Stumpf, A.; Malet, J.P. A method for automatic and rapid mapping of water
surfaces from sentinel-1 imagery. Remote Sens. 2018, 10, 217. [CrossRef]
58. Bovolo, F.; Bruzzone, L. A split-based approach to unsupervised change detection in large-size multitemporal
images: Application to tsunami-damage assessment. IEEE Trans. Geosci. Remote Sens. 2007, 45, 1658–1670.
[CrossRef]
59. Toral, G.M.; Aragones, D.; Bustamante, J.; Figuerola, J. Using Landsat images to map habitat availability for
waterbirds in rice fields. Ibis 2011, 153, 684–694. [CrossRef]
60. Padró, J.-C.; Pons, X.; Aragonés, D.; Díaz-Delgado, R.; García, D.; Bustamante, J.; Pesquer, L.;
Domingo-Marimon, C.; González-Guerrero, Ò.; Cristóbal, J.; et al. Radiometric correction of simultaneously
acquired Landsat-7/Landsat-8 and Sentinel-2A imagery using pseudoinvariant areas (PIA): Contributing to
the Landsat time series legacy. Remote Sens. 2017, 9, 1319. [CrossRef]
61. Martins, V.; Barbosa, C.; de Carvalho, L.; Jorge, D.; Lobo, F.; Novo, E. Assessment of atmospheric correction
methods for Sentinel-2 MSI images applied to Amazon floodplain lakes. Remote Sens. 2017, 9, 322. [CrossRef]

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).