Synergistic Approach of Remote Sensing and GIS Techniques for Flash-Flood Monitoring and Damage Assessment in Thessaly Plain Area, Greece

: This paper describes the synergetic use of earth observation satellites optical and radar data with a high-resolution digital elevation model (DEM) to detect ﬂooded areas and explore the impacts of a ﬂood event. A ﬂash ﬂood episode took place in May 2016, in the central-eastern part of West Thessaly (Central Greece). Landsat-7 ETM+ and a Sentinel-1 SAR images were acquired. For Landsat-7, several water indices were applied and for the Sentinel-1 a threshold method was implemented. Elevation data were also used to improve the delineation of the inundated areas, and to estimate ﬂood water depth. Furthermore, Sentinel-2 images were utilized so as to record the land use/cover of the ﬂooded area. The inundated areas and the affected cultivations were delineated with high precision, and the ﬁnancial effects were evaluated.


Introduction
Flash floods are considered to be among the most frequent and destructive types of natural disasters worldwide, with significant consequences including: (a) human and animal life losses, (b) agricultural crop destruction and soil loss, (c) damage to infrastructure and communication networks, and (d) transport of sediment loads and pollutants [1][2][3][4][5][6][7][8][9].Particularly, it is estimated that floods provoke about 40% of the damage caused by all natural disasters [10,11].Definitions and classifications vary for different types of floods, including catastrophic, flash, extreme, river, coastal, pluvial, and large floods [12][13][14].Flash floods are the most common type of flood in the Mediterranean areas, which occur due to the small size of the river basins (short flow concentration time), the geomorphology (high slopes or/and low soil and geological formations permeability), and intense rainfall [15].Therefore, mapping and monitoring flood events and their impacts are essential for any flood risk mitigation plan, disaster detection, and compensation services [16].This information is often difficult to produce using traditional techniques because inundated areas are transient, changing fast, cover large areas, and by their nature restrict the accessibility [17].
Earth observation (EO) data from space can provide valuable and timely information aiding the effective response and mitigation of natural hazards such as floods.Satellite observations enable the acquisition of data for large and hard-to-reach areas, they may provide continuous measurements, and facilitate the detection of the flood extent, which is crucial information for damage assessment and risk management [7,9,18,19].
Currently, the most common approach is to use both optical and radar data.Satellite images captured by optical sensors, such as those on board the Landsat-7 (L7), or Landsat -8 and Sentinel-2 (S2) satellites; these satellites acquire data from several spectral bands that cover the visible through shortwave infrared (SWIR) spectrum, and provide high spectral resolution and temporal data as well as a substantial archive to cover a range of hydrological conditions [20,21].Optical sensors detect energy naturally reflected or emitted by the earth's surface in visible and infrared spectral bands.Where clouds, trees, and floating vegetation do not obscure the water surface, high-resolution visible/infrared sensors provide good delineation of inundated areas [18,22,23].Several different water indices have been used up to now to delineate flooded areas in studies using Landsat imagery.Specifically, the Normalized Difference Water Index (NDWI), which can outline and enhance open water features by means of near infrared (NIR) and green channels of Landsat, was initially used [24,25]; it was later modified (Modified Normalized Difference Water Index) by substituting the middle infrared band for the NIR band [26][27][28].The MNDWI can enhance open water features while efficiently suppressing and even removing built-up land noise as well as vegetation and soil noise [26].Tasseled cap transformations (TCT) of L7 ETM+ data have also been used to map water [29][30][31].TCT summarizes Landsat data into fewer components, by rotating pixel vectors in multidimensional space, using all bands.The third component of the L7 tasseled cap transformation relates to wetness, and tasseled cap wetness (TCTW) has been used as an input to delineate the flooded area [32][33][34].Furthermore, several other differentiated and sophisticated water indices have been more frequently utilized, such as the Normalized Difference Water Index (NDWI) of McFeeters, the Difference Between Vegetation and Water Index (DVW), the Index of Free Water (IFW), the Water Impoundment Index (WII), the Modified Water Index (MWI), the Red and Short Wave Infrared index (RSWIR), and many others [24,26,27,[34][35][36][37][38][39].
Since the early 1970s, the Landsat program has provided us with the longest continuous global measurements of the Earth's surface, creating a historical collection unparalleled in quality, coverage and length.L7 is the seventh satellite of Landsat mission equipped with Enhanced Thematic Mapper Plus (ETM+); it operates on six multispectral bands with a spatial resolution of 30 meters, a thermal band at 60 meters and a panchromatic band at 15 meters.
The European Space Agency (ESA) through the Copernicus program offers accurate, timely, and easily accessible SAR data of Sentinel-1 satellite system.Sentinel-2 (S2) is a new land monitoring and classification mission that provides high spatial resolution optical imagery to perform terrestrial observations with global coverage of the Earth's land surface every five days.S2 is equipped with the MSI (Multispectral Imager) sensor operating on 13 different bands at a spatial resolution of 10 (4 bands), 20 (6 bands), and 60 (3 bands) meters.Its data have significant use in studies related to land use/cover mapping and the quantification of related changes [40,41].
Optical sensors are unable to procure clear images during cloudy weather conditions and therefore its application potential is limited in case of flood events accompanied by clouds.In such cases, SAR images are valuable, as several previous studies have been demonstrated by Kussul et al. [7]; Matgen et al. [42]; and Klemas [18], because they are not dependent on daytime and weather conditions, they provide data through clouds and hazy atmospheric circumstances and they can differentiate waterlogged and dry land [1,7,9,18,21,[42][43][44][45][46][47].Cloud penetration capacity turns out to be considerably significant for monitoring flood events as they customarily take place during periods of extended rainfall [18,43,48].Several commonly used SAR flood extent mapping techniques include image texture algorithms [47], histogram thresholding [9], and various multi-temporal change detection methods [45].Henry et al. [49] showed that the detection of inundated areas can be optimized by selecting the most suitable polarization of the radar waves.It has been observed that usually HH polarization (H=Horizontal) provides a more efficient distinction of regions flooded than HV or VV (V=Vertical).VV polarized data were highly influenced by surface roughness conditions, but in some cases, they provide good results in flooded area delineation [18,49,50].The main disadvantage of SAR data is that the resulting backscattered signal is a complex combination of effects which depend on incidence angle, vegetation density and orientation, relative water height, and wind effects.This could cause opposite backscatter responses for similar conditions, and the disentangling of such effects requires additional information [23].
Sentinel-1 (S1) is a synthetic aperture radar (SAR) mission of the Copernicus Program.The contribution of S1 to the application of flood mapping is due to the sensitivity of the backscatter signal to open water [40].Through radar intensity imagery, in the absence of wind, the specular reflection of C-band signals over open water means that the signal is significantly lower than average.In the case of flash floods though, wind could be a problem.Furthermore, vegetation that is not entirely covered by water can be also a problem (for optical also) [9, 21,40,51,52].
Having as a goal to overcome the inherent limitations of the various remote sensing (RS) techniques (cloudiness, trees or floating vegetation, long revisit cycles of satellites, etc.) [18,22,23,53], various researchers investigated the use of additional spatial data such as land cover and elevation [54,55].As an example, Matgen et al. [46] developed image processing methodologies that allow for the extraction of water levels at any point of the floodplain in order to further improve the exploitation of SAR images for flood mapping that is crucial in hydraulic modeling and near real-time crisis management.Chen et al. [56] investigated the use of a LiDAR survey during a major flood event and a digital terrain model to document flood water depths.Qian et al. [57] also used LiDAR-derived DEM for the extraction of uncultivable trails, ditches, and cultivated field parcels within flood plains in order to evaluate flood risk.Furthermore, Vesakoski et al. [58] investigated the effect of various sources of digital terrain models (DTM) on the delineation of topographic depressions and their hydrological impacts.
The main objective of this study is to further investigate the application of remote sensing techniques, using optical and SAR open data in conjunction with other data sources and GIS techniques, for monitoring and evaluating the impacts of a natural disaster.The case study was a flash flood event that took place from 20 to 22 May 2016, following an intense and heavy rainfall that occurred in the central-eastern part of West Thessaly, in Central Greece, including a part of the Thessalian plain, which is characterized by intensive agricultural activity.A L7 image, acquired on May 23, one day after the phenomenon and two Sentinel-1 SAR images (S1), acquired on May 21 (in the middle of the storm event), and were processed.For L7 the Tasseled Cap Transformation (TCT), the Modified Normalized Difference Water Index (MNDWI), the Difference between Vegetation and Water (DVW), and the Red and Short Wave Infrared (RSWIR) index, were applied and tested.Their results were compared, and from their analysis, the inundated area was delineated.As far as the S1 SAR data are concerned, a Level-1 Ground Range Detected (GRD) product has utilized.Also, a method setting a threshold value (that is established manually) of radar backscatter, which is followed by a binary algorithm to determine whether a given raster pixel is flooded or not, was implemented.A detailed DEM with 5 m spatial resolution was also used to assess the flood risk in the affected area, to improve the delineation of the inundated areas, and to estimate flood water depth.Finally, four Sentinel-2 (S2) images were utilized so as to map the land use/cover and classify the existing crop types in the flooded area [59].
Specifically, this study set out to: (a) examine the utility, the accuracy, and the limitations of high resolution optical and SAR data for the detection and delineation of floodplain water, (b) investigate the use of additional spatial data to improve the obtained results, (c) recommend a suitable methodology for mapping inundated area extent and flood water depth to better assess flood impacts, (d) determine the best data processing required technique for such cases, and (e) explore the advantages of the synergistic use of remote sensing (optical and radar) and GIS (detailed DEM) techniques to overcome limitations in flood monitoring and damage assessment.

Study Area
The study area lies between 39 • 11' N and 39 • 58' N latitude and 21 • 52' E and 22 • 45' E longitude (Figure 1), in the Pinios river basin (GR16), the largest catchment in the Thessaly Water District (GR08) with an area of 10,700 km 2 , mean annual rainfall 779 mm and mean annual runoff 3500 hm 3 (327 mm) [60].It is bounded by various mountainous regions, including the highest mountain of Greece, Olympus (altitude: 2918 m), encircling a low plain.The Thessaly plain lies in the center of the catchment and it is mostly a flat terrain which was subsequently filled with alluvial deposits that form the largest and one of the most productive agricultural lands in Greece [61,62].Specifically, 49.5% of the Thessaly Water District area is being used for agricultural purposes while it is estimated that approximately 14.2% of the primary products of Greece are being produced on the Thessalian plain [63].Pinios river originates from the Pindus sierra, the biggest mountain range of Greece, and empties into the Aegean Sea.The affected areas reside in Karditsa, Larissa, and Trikala prefectures (geographic administrative region of Thessaly) along with a slight segment in the northern part of Domokos municipality that belongs to Fthiotida prefecture (geographic administrative region of Central Greece).

Background
The geomorphological characteristics of Pinios river basin mentioned above, e.g., the circular shape of the catchment and flat to gentle slopes favor inundation conditions and thus render it prone to flood events.Particularly, according to the EU Floods Directive (2007/60/EC)-Preliminary Flood Risk Assessment that was implemented by the Ministry of Environment and Energy of Greece, several notable historical floods had been documented, such as those that occurred on 4 June 1907; 27 October 1980; 23 March 1987; 22 October 1994, while 31.7% of the total area of the Thessaly Water District is characterized as an area of potential significant flood risk [63].

Intensity of Rainfall from 20 to 22 May 2016
In the broader area, three meteorological stations are operated by the National Observatory of Athens [64], namely Kalabaka, Trikala, and Karditsa.In the 2013-2017 period (except for the study year 2016), the height of rain output for May varies between 13.2 and 66.8 mm.On the contrary, the recorded rainfall from 20 to 22 of May 2016 for Kalabaka, Trikala, and Karditsa station was 78.2, 87.8, and 91.8, respectively.Kalabaka station was the most representative station since it is located in the most crucial part of the flooded areas.The intense storm lasted from noon of 20 May until the dawn on the 22 May while its peak was reached on 21 May.

Materials
Several data sets were investigated to be acquired for the purpose of this study.A Level 1T L7 image (path/row: 184/33), free of charge was obtained from the USGS Earth Explorer (EE) tool.The image was geometrically corrected and registered to a UTM 34N WGS84 ellipsoid with standard terrain correction applied.The L7 image was acquired on 23 May, hence one day after the two-day flood event [65].It should be noted that Landsat-8 images were not available during the flood event.
As well, S1 and S2 datasets were provided free via ESA's Sentinels Scientific Data Hub [66].S1 data was an L1 GRD, high resolution (10 meters) Interferometric Wide (IW) image, of VV and HV polarization, using the Terrain Observation with Progressive Scans SAR (Table 1).The SAR image was obtained from S1 C-band during the flood, on 21 May (Time 07:05:10 pm).The available S2 products are in processing Level-1C, geometrically corrected and registered to the UTM 34N WGS84 ellipsoid with elevation correction applied.Four S2 images were acquired, from December 2015 till August 2016 (Table 1), so as to cover all the plant growth stages of the crops over a year and create the land cover classification map (Table 1).Topographic maps of the Hellenic Military Geographical Service (HMGS) were used to extract different types of info-layers: geomorphology, rivers, roads, railway tracks etc.A detailed DEM with 5 m spatial resolution and 0.5 m in geolocation, derived from topographical maps (scale 1:5000, contour interval 4m) of HMGS, was also used [67].
Detailed spatial data of the utilized agricultural area (UAA) and the specific crop cultivated in each farm of the affected area during the study period were obtained by the Integrated Administration and Control System (OPEKEPE) [68] in the form of vector files with scale 1:10,000.Finally, compensation data of the Greek Agricultural Insurance Organization (ELGA) for the flood event period of the 20-22 May 2016 were processed and compared with the results of this study.
The digital image processing and analysis of the satellite data were carried out using ENVI (v.5.3) and SNAP (v.5.0) software (Harris Geospatial Solutions, USA and European Space Agency, respectively), while the manipulation of the spatial information was made using ArcGIS (version 10.5.1, Environmental Systems Research Institute, Redlands, CA, USA).

Methodology
Optical remote sensing data collected by satellites are subjected to various deleterious effects related to sensor observation and solar illumination geometry, solar radiance/atmosphere conditions, detector calibration and settings, topography, as well as date and time of acquisition.It was necessary to achieve the comparison of images acquired from different sensors and at different times, to normalize the data sets to a standard set of conditions (normally reflectance values) [16].
The intervening atmosphere between the satellite and the Earth's surface can degrade the quality of the remote sensing data, by altering the intensity of the spectral response of ground features received at the sensor from their actual spectral characteristics.Several algorithms are used for removing the effects of the atmosphere on the data and by which sensor radiance values are thereby converted to reflectance values.Various levels of input data are used to drive these atmospheric models, ranging from detailed atmospheric profiles acquired by radiosondes to standardized models representative of general latitudinal areas [16,69,70].As atmospheric conditions can vary both spatially and temporally, standardised atmospheric models-such as that provided by QUAC atmospheric correction model in ENVI-is crucial for the atmospheric correction of the L7 image [71,72].Moreover, the treatment of different image data sets derived from different sensor systems in multiple combinations requires a high degree of geometric accuracy for proper feature information extraction and use in geographical information or environmental monitoring systems.Thus, an image to image registration was made between the two different satellite systems images, to achieve a well-built precision that will help the impending comparison of the data.
In the present study, three different processing methodologies were followed for the three-different kinds of satellite data used (Figure 2).The final pre-processing step of the L7 ETM+ image included the gap filling that scan line corrector (SLC) failure creates.For that reason, Landsat gap fill plugin of ENVI was used to retrieve approximately 22% of data loss (Figure 3).This process deals with the implementation of a gap mask for each band that points the missing data in the scan gap and classifies these regions as 0 and the existing data as 1.Then, the digital number values of the L7 image have been converted to top of atmosphere reflectance and then to surface reflectance.This procedure was fulfilled utilizing the digital elevation model (DEM) of the area, geographic longitude/latitude and the weather conditions that were present at the time of image acquisition.These latter parameters were automatically calculated by importing the images through the suitable import algorithm.The basic processing step concerned the calculation of the four water indices, the MNDWI, TCT, DVW, and the RSWIR (Table 2).Modified Normalized Difference Water Index enhances open water features while suppressing noise from built-up land, vegetation, and soil (Figure 4a).Xu [26] reported that the MNDWI produced better results than the Normalized Difference Water Index [35] in enhancing and extracting water from a background that is dominated by built-up land areas.The tasseled cap transformation [73] is a useful tool for compressing spectral data into a few bands associated with physical scene characteristics [30,31].In particular, it converts the readings in a set of The basic processing step concerned the calculation of the four water indices, the MNDWI, TCT, DVW, and the RSWIR (Table 2).Modified Normalized Difference Water Index enhances open water features while suppressing noise from built-up land, vegetation, and soil (Figure 4a).Xu [26] reported that the MNDWI produced better results than the Normalized Difference Water Index [35] in enhancing and extracting water from a background that is dominated by built-up land areas.The tasseled cap transformation [73] is a useful tool for compressing spectral data into a few bands associated with physical scene characteristics [30,31].In particular, it converts the readings in a set of bands into composite values, using the linear regression method-similar in concept to principal components analysis.The coefficients used to create the TCT are derived statistically and vary among the imaging sensors.The composite values define a new coordinate system with a new set of orthogonal axes in order to represent the brightness of each pixel, the degree of greenness/yellowness of vegetation and the wetness which is correlated to soil moisture, water, and other moist features.The other additional components contain image noise and atmospheric influences, such as clouds, haze, sun angle differences, etc.For L7 data, the tasselled cap transformation produces six output bands: brightness, greenness, wetness, fourth (haze), fifth, sixth.This type of transformation should be applied to atmospherically corrected and calibrated reflectance data rather than applying to raw digital number imagery.The third component relating to the tasseled cap wetness (TCW) has been used as a tool for flooded area delineation.Furthermore, the DVW and RSWIR indices were implemented.The first one functions by subtracting the simple Normalized Difference Water Index (NDWI) from the Normalized Difference Vegetation Index (NDVI), while the second one uses the short wave infrared (SWIR) band in its equation.
One of the major advantages of using SAR images, apart from its all-weather capability, lies in its ability to sharply distinguish water from other classes.Water bodies act as a mirror reflecting surface, their response is low (low backscatter coefficient in SAR images) and thus look like a dark area.The land, for its part, gives a much higher amount of radar energy due to the surface roughness and this generates the high contrast between surfaces: soil and water [75].The S1 SAR image obtained from Sentinel-1C-band comprises an image which was acquired at the early stages of the flood.The processing method sets a threshold value (which was established manually) of radar backscatter, which is followed by a binary algorithm (band math equation) to determine whether a given raster pixel is flooded or not.The data processing of the S1 VV polarization SAR image followed specific steps.The first step is related to the pre-processing of the data.Initially, a subset of the images was made, to reduce the size of the data and make them more easily managed.Then, a calibration of the images was carried out, comparing the instrument's measuring accuracy to a known standard.Finally, a speckle filter (Lee filter with window size 7 by 7) was applied [76].Speckle appears as a grainy 'salt and pepper' texture in an image.This is caused by random constructive and destructive interference from the multiple scattering returns that occur within each resolution cell.Speckle is essentially a form of noise which degrades the quality of an image and may make interpretation (visual or digital) more difficult.Thus, it is generally desirable to reduce speckle before interpretation and analysis.The second step deals with the main process of the data.Thus, a threshold was selected, analyzing first the histogram of the filtered backscatter coefficient.The histogram shows a peak of different magnitude.Low values of the backscatter correspond to the water, and high values correspond to the non-water class.To binarize the image band arithmetic was applied, putting as logical value (true) for values less than the chosen threshold and false for higher values, producing the final 'water' image.
As far as it concerns the four S2 images, the atmospheric correction was applied to convert the top-of-atmosphere reflectance values (TOA) to corrected bottom-of-atmosphere reflectance values (BOA).For the atmospheric correction, ESA's Sen2Cor plugin was used.Finally, bands 1, 9, and 10 were separated from the dataset.Then, the images were resampled using a bilinear method and were also subset to the study site extending [41,66].Then, a layer stack of the four images was made and a supervised maximum likelihood classification method was implemented, utilizing approximately 200 sampling points (Figure 4a-d).The resulted land cover classification was compared with the UAA data obtained by the Integrated Administration and Control System OPEKEPE [68] in order to cross check their accuracy.
Following the delineation of the inundated areas using remote sensing, the detailed DEM with 5 m spatial resolution was introduced to improve the inundation extraction results, to estimate the inundation depth as well as to assess the flood risk in the affected area.The 5 m digital elevation model (DEM) is shown in Figure 5.As it can be seen, most of the affected area, and especially the area which is located between the levees, is almost flat with small elevation differences.The altitude gradually decreases from upstream (west and south) to downstream (east).The use of other freely available DEMs like SRTM or ASTER GDEM was also considered; however, the low vertical resolution of these datasets (1 m) and especially the low vertical accuracy in flat areas [77] were critical limitations for this application.To this end, a flood surface elevation model was created, based on the ground elevation at the boundaries of the inundated areas obtained by RS.Specifically, the elevation along the boundary lines of the identified flooded areas was sampled from the DEM every 5 m (Figure 6a).Then, a low pass filter with a 5 × 5 pixel window was applied to smooth the obtained elevation data; the aim was to reduce the influence of outliers or extreme values which are due to possible errors in the inundated areas map, errors in the DEM, misalignment between the RS and the DEM pixels, etc.Then, the flood water elevation was spatially interpolated using a first-order local polynomial interpolation with barriers.An approximate spatial interpolation method was used instead of an exact method in order to further reduce the impact of outliers or extreme values in the obtained flood surface.Furthermore, the levees along the rivers were set as barriers in the spatial interpolation in order to consider the sharp changes of water surface elevation caused by their presence (Figure 5).It should be noted that initial tests assuming a constant water surface elevation in the flooded polygons indicated that neglecting the slope of the water surface along the river leads to large errors.The flood depth map is obtained by subtracting the DEM from the flood surface elevation model (Figure 6b).Based on the flood depth map a fine tuning of the inundated areas delineation was also attempted.To this end, pixels with flood depth lower than 0.05 m were considered as not flooded.This threshold was used in order to avoid the characterization of pixels having very small differences between the land surface and the water surface as flooded areas, which may be the result of very small inaccuracies.The chosen threshold (0.05 m) is a round number and at the same time a reasonable amount that could be normally applied as irrigation water, so it is assumed that may not cause problems to the crops.In any case, other depth thresholds can be easily considered for specific purposes (e.g., for specific crops) based on the flood depth map.Additionally, the obtained flood surface elevation model was extended outside the boundaries of flooded polygons in order to further improve the inundated areas delineation by considering areas near the flooded polygons with lower ground elevation than the water elevation.
As an additional step in this analysis, topographic depressions in the affected area were also identified based on the DEM.Topographic depressions are areas with lower ground elevation than their surrounding areas (local elevation minima) that have no direct connection to the downslope flow paths [78].Accordingly, topographic depressions are susceptible to flooding in the absence of artificial drainage systems or in case that drainage systems cannot outflow the excessive water to the natural streams (e.g., due to high water levels).Normally, depressions are filled prior to hydrologic analyses requiring hydrologically connected flow networks [79].However, in this analysis the topographic depressions are identified and analyzed as a means to assess areas prone to flooding.For this purpose, the method developed by Jenson and Domingue [78] was used to fill the topographic depressions in the affected area.In this approach, the presence of artificial drainage systems in the affected area was not considered.Then the original DEM was subtracted from the filled DEM in order to obtain the depth of the topographical depressions.The resulted output was classified with equal intervals to provide a measure of flood risk.
The obtained flood depth and flood risk raster datasets were clipped using the crop parcel polygons obtained by the Integrated Administration and Control System OPEKEPE [68].This dataset contains detailed information on the exact boundaries and the crops cultivated in each farm and it is updated on a yearly basis [80].It was thus made possible to assess in great detail the UAA affected as well the specific crops cultivated in the affected farms, neglecting for example uncultivable trails and ditches.Accordingly, the obtained results were able to provide a clearer picture of the flood impact and were directly comparable with the information provided by the Greek Agricultural Insurance Organization (ELGA) [81] on the compensations paid for this flood event.

Results and Discussion
During floods, timely and detailed situation reports are required by the disaster management authorities to locate and identify the affected areas and to implement the corresponding damage mitigation; this is the most delicate management category since it involves rescue operations and the safety of people and property [82,83].Many of the existing and future satellite missions and airborne platforms provide rich data with great potential for enhanced monitoring, measuring and mapping of floods [50,84].In this study, the flood dynamics and its impact were revealed, calculated, and analyzed using multisource data (L7, S1, and S2 EO data; DEM and public domain data).
S1 satellite mission data are valuable for flood monitoring and damage assessment.However, in this case study, the outputs of water extraction using S1 data (Figure 7) represent solely numerous dispersed spots given that the image was obtained at the early stages of the phenomenon.Specifically, the transformation of the precipitation to runoff and thence to flood routing, requires a certain period that is related to the hydrological and geomorphological characteristics of the watershed.Thus, the major part of the study area was not under water at the S1 acquisition time at 07:03:18 a.m. of 21 May, because the phenomenon was at its early stages.The four calculated water indices based on the L7 image (MNDWI, TCW, DVW, and RSWIR) are presented in Figure 8.The evaluation results of the inundated area which were extracted from the water indices, utilizing raster calculator of ArcGIS [equation: Con ("raster"> = 0, "raster",0)], reveal that MNDWI (Figure 8a), TCW (Figure 8b) and RSWIR indices (Figure 8d) resulted in a higher accuracy of surface water mapping as compared to DVW (Figure 8c).The three indices calculated a flooded area of 5381.3 ha while by the DVW an area of 5125.2 ha was estimated.The flood extent map obtained from the four different water indices derived from L7 data has been compared and evaluated.The flood map depicted that during 20-22 May 2016 an area of 5381.3 ha was inundated in the study site (Figure 8).Specifically, from the MNDWI, TCW and RSWIR indices the estimated total flooded area was approximately the same (about 54 km 2 ), while from the DVW there was an underestimation and it was assessed nearly at 49 km 2 .As already mentioned, the L7 image was acquired on 23 May 2016, i.e., one day after the end of the storm event.Hence, there was a reduction of the total surface area covered by water in comparison to the maximum inundated area during 20-22 May 2016.Another limitation is that flood water depth, which is a critical factor for the assessment of floods impact on the agricultural production, cannot be directly estimated using remote sensing techniques [56].Considering other possible limitations, such as trees or vegetation not entirely covered by water, two types of areas can be distinguished: (a) the flooded areas derived from the water indices and (b) the broader flooded areas (areas that were probably affected by the flood).
Given these considerations the detailed DEM with 5 m spatial resolution was used to improve the inundation extraction results (using the MNDWI), to estimate the inundation depth as well as to assess the flood risk in the affected area.At a first step, a flood surface elevation model was calculated as it is described in detail in the methodology section and then the flood depth map was created by subtracting the DEM from the flood surface elevation model.The obtained flood depth map is presented in Figure 9. Specifically, in Figure 9a the flood depth map is limited to the flooded polygons delimited by RS, while in Figure 9b the flood depth map is extended outside the boundaries of flooded polygons; thus, the inundated area delineation was further improved by considering areas near the flooded polygons with lower ground elevation than the water elevation.The flood depth maps presented in Figure 9 were clipped on the utilized agricultural area in order to provide a better representation of the flood impacts neglecting uncultivable trails, ditches, etc.The corresponding flooded areas statistics for both the original and the clipped flood depth maps are presented in Table 3.As it can be seen the fine-tuned total affected area obtained by the flood depth map (48,536,100 m 2 ) is a little smaller than the total flooded area obtained by RS (53,836,700 m 2 ); this is because parts of the inundated area (according to RS) had actually higher or equal ground elevation than the flood water elevation.
Furthermore, the total affected area clipped to the cultivated parcels (neglecting uncultivable trails, ditches, etc.) is about 18% smaller than the original flooded area delimited by RS, revealing that the use of detailed spatial information for the UAA is important for the accurate estimation of the impacts on the agricultural production.Another important finding is that the extended area, including areas near the flooded polygons with lower ground elevation than the water elevation, is considerably larger than the total flooded area obtained by RS.This is due to the fact that the available RS image acquired one day after the end of the storm event.During this one-day period, large parts of the initially flooded area naturally drained as the water level in the main streams dropped.The extended map, therefore, may provide a hint of the possible extent of the flooded areas at the end of the storm event.However, it should be noted that this assumption cannot be tested with the applied methodology, because the extent of the area drained in one day depends on many factors (e.g., floodplain morphology, drainage network characteristics, stream flow dynamics) and its estimation requires hydraulic modeling techniques.It is also interesting that the total area which received compensations (38,238,400 m 2 ) by the Greek Agricultural Insurance Organization (ELGA) is very close to the total affected area delimited by RS and clipped to the UAA (39,963,400 m 2 ).This observation indicates that acquiring the RS image one day after the end of the event provided a sufficient time window for the estimation of areas that were inundated long enough to significantly affect the cultivations (i.e., very short inundation may have a limited impact on production depending on the cultivation).It should be noted that compensations are only paid under certain circumstances, namely the total damaged production is greater than 20%, the total annual compensation of each parcel does not exceed 80%, etc. [81].Accordingly, considering the temporal dynamics of the flood as well as the flood depth instead of only the maximum extend of the flooded areas, may provide a better estimation of the inundated area.
Detailed statistics of the various crops affected by the flood based on ELGA compensations, on the flooded areas limited to the flooded polygons delimited by RS, and on the flooded area extended outside the boundaries of the delimited polygons are presented in Table 4.The comparison of the affected areas for each crop according to ELGA and to RS in Table 4 indicates that there is a direct correspondence between the flooded area that was delineated with the present methodology and the area calculated and compensated by ELGA based on field surveys.The observed differences for individual crops (e.g., fodder crops) could be a result of the different sensitivity of each crop to temporal inundation.For example, fodder crops could suffer less damage for the same inundation time and depth than cotton or vegetables (they are a significant part of the other crops).In order to examine the accuracy of the land use data provided by the Integrated Administration and Control System OPEKEPE [68] a land use/cover classification of four multi-seasonal S2 images (Figure 4a-d) using a supervised maximum likelihood classification was also made, utilizing approximately 200 sampling points (mainly declared parcels).The overall accuracy of the classification was 89%, with a kappa of 0.87, due to high homogeneity of the different land use/cover types of the area.The lowest accuracy was observed mainly in the peri-urban areas where there is an intense mixture of classes and secondly the deficient or not updated declarations of the farmers (Figure 10a,b).10a).a final step in this analysis, topographic depressions in the affected area that are more susceptible to flooding were also identified.Flood risk was due to the absence of artificial drainage systems or to their inability to outflow the excessive water to the natural streams (e.g., due to high water levels in the main streams).This process was based on the DEM as a means to assess areas prone to flooding and to provide a measure of flood risk.The obtained flood risk map is presented in Figure 11.As can be seen, the areas with higher flood risk are visually correlated to the areas with greater flood water depths in Figure 9. Furthermore, it is interesting that the total area of the identified depressions is 148,432,800 m 2 and it is similar to the maximum estimated flooded area (extended), which is 151,872,700 m 2 .

Conclusions
Remotely sensed data are being increasingly used for flood extent mapping and damage assessment.Optical sensor data are readily available free of charge and they can be used with available well-defined, robust data processing techniques.However, as it was shown in this study, long revisit cycles of satellites can be a serious limitation for the investigation of flooding temporal dynamics using solely remote sensing data.Temporal evolution of the inundated areas can be critical for the assessment of the impact of floods in agricultural production, especially in the case of flash floods.Another limitation is that flood water depth, which is a critical factor for the assessment of floods impact on the agricultural production, cannot be directly estimated using remote sensing techniques.Accordingly, in this study case it was attempted to overcome these limitations by using multisource data (L7, S1, and S2 EO data; DEM and public domain data) to calculate and analyze flooding dynamics and their impact on agricultural production.
Radar data provided by the S1 satellite mission are quite valuable for flood monitoring and damage assessment.However, at this case study, inundated areas delineation using S1 data resulted solely in numerous dispersed spots because the available image was acquired at the early/middle stages of the phenomenon, highlighting the significant limitation of six-day revisit cycle.Revisit time is expected to decrease with the launch of new missions and the integration of different satellites constellations.This will likely make the water techniques for flood monitoring more reliable in the near future.
The most suitable RS data in terms of acquisition time was an L7 image acquired on 23 May 2016, one day after the end of the storm event, which was a very fortunate situation considering that usually cloudy conditions prevail and flooded areas are not visible using optical data.Accordingly, in this case, a reduction of the total surface area covered by water in comparison to the maximum inundated area during 20-22 May 2016 is expected.
The thorough analysis of the four different water indices calculated in this study based on the L7 image (MNDWI, TCW, DVW, and RSWIR), showed that all indexes were efficient at classifying pure water pixels; however, MNDWI and RSWIR indices appeared to be more suitable for delineating water bodies during floods.
A detailed DEM with 5 m spatial resolution was used to fine tune the inundation extraction results, to estimate the inundation depth and to assess the flood risk in the affected area.The fine-tuned total affected area obtained by the flood depth map was a little smaller than the total flooded area obtained by RS because parts of the inundated area (according to RS) had actually higher or equal ground elevation than the flood water surface elevation.Interestingly, the total area having received compensations by ELGA is very close to the fine-tuned total affected area.This observation indicates that acquiring the RS image one day after the end of the event provided a sufficient time window for the estimation of areas that were inundated long enough to significantly affect the cultivations.Accordingly, considering temporal dynamics and flood depth instead of only the maximum extent of the flooded areas, may provide a better insight on the actual impacts.Furthermore, a direct correlation between the flooded area that was delineated with the present methodology and the area calculated and compensated by ELGA, based on field surveys, was observed.
Finally, concerning the assessment of flood risk based on the analysis of topographic depressions it was found that the areas with higher estimated flood risk were visually correlated with the areas having greater flood water depths.
The approach used in this study may provide a suitable alternative for the complementary use of multisource spatial data to assess the impacts of floods.However, further investigation in more case studies with various characteristics and using more alternative data sources are needed to validate the feasibility and the accuracy of this approach.
Author Contributions: E.P. and M.Z.contributed to the acquisition of the appropriate data and to their pre-process and process stages, as well as to the analysis of the results.K.X.S. contributed to the analysis of the DEM and the calculation and analysis of the flood surface elevation model, the flood depth and the flood risk maps.E.P., K.X.S., N.D. and M.Z.wrote the paper.
Funding: This research received no external funding.

Figure 1 .
Figure 1.The study area located in Central Greece, with its principal geomorphological and hydrological characteristics.(a) Greece; (b) Central Greece.

Figure 2 .
Figure 2. Flow-chart of the implemented methodology.

Figure 3 .
Figure 3. Filling the wedge-shaped gaps of the L7 image, (a) image with gaps, (b) corrected image.

Figure 4 .
Figure 4.The four S2 images (False Color Composite: NIR, red, green band as RGB) used for the classification of the land use/cover of the area (a) 6 December 2015, (b) 04 April 2016, (c) 13 July 2016, and (d) 12 August 2016.

Figure 5 .
Figure 5. Detailed DEM of the affected area with 5 m spatial resolution.

Figure 6 .
Figure 6.Flood water depth calculation example.(a) Sampling points at the RS polygons boundaries; (b) flood water depth model.

Figure 7 .
Figure 7.The 'water' image derived from the processing of S1 SAR image.

Figure 9 .Table 3 .
Figure 9.The obtained flood depth map.(a) Clipped to the flooded polygons delineated with the RS methodology; and (b) extended in the entire affected area.Table 3. Flooded area statistics for the original and the clipped flood depth maps and for the cases that the flooded areas are limited to the flooded polygons delimited by RS or are extended outside the boundaries of the delimited polygons (extended).

Figure 10 .
Figure 10.(a) The land use/cover classified image; (b) Zoom in image showing the cross-checking of Sentinel-2 classified image and the land use/cover file of OPEKEPE for the accuracy assessment process (red polygon in Figure 10a).

Figure 11 .
Figure 11.Topographic depressions classified according to their depth altitude to provide a measure of flood risk.

Table 1 .
Satellite images data sets.

Table 4 .
Affected areas statistics according to the Greek Agricultural Insurance Organization (ELGA) compensations, to the flooded areas limited to the flooded polygons delimited by RS, and to the flooded area extended outside the boundaries of the delimited polygons (extended).