2.1.

Study Area

This study focuses on the significant flood event that occurred in Kerala, a southern state of India, during August 2018, as shown in Fig. 1. Kerala is situated between 8.291°N to 12.795°N latitudes and 74.863°E to 77.412°E longitudes. The spatial extent of the Kerala state is about $\mathrm{38,863}{\mathrm{km}}^2$. The state is majorly categorized into three geographical regions, i.e., low lands, mid lands, and high lands. Low lands are mainly characterized by backwaters, river deltas, and Arabian Sea shore.⁶⁹ Unlike other metropolitan cities, the urban areas in Kerala state are covered with a large amount of tree cover.

Fig. 1

Study area showing the Kerala state of India. Sentinel-1 swath coverage over the test area is shown in green rectangles. PWB is derived from GEE and the details are explained in Sec. 3.3.

Kerala state received two abnormal rainfall spells during August 7 to 9, 2018, and August 14 to 19, 2018. A total of 857.4 mm of rainfall, which is 18% above normal, occurred in July 2018. It is the first onset of flooding in Kerala state. With this heavy rainfall in July 2018, the rainwater started flowing into the dams. Later, the state received a total rainfall of 758.6 mm between August 1, 2018, and the end of second rainfall spell, which is 164% more than normal rainfall.⁷⁰ By the end of the first rainfall spell, the water level in several dams was up to full reservoir level. The water from the dams was released due to heavy rainfall at the time of second rainfall spell. Also, there were numerous landslides throughout the Kerala state. This unprecedented massive rainfall combined with dam water release,⁷¹ landslides, and improper urban planning contributed to severe flooding in 13 out of 14 districts in Kerala state. According to the National Disaster Management Authority, Kerala State Disaster Management Authority (KSDMA), and other secondary sources, a large part of Kerala’s population is affected by floods and its related vector-borne diseases. According to KSDMA, 504 people died, and 3.4 million people were accommodated in 12,300 relief camps. The other socio-economic losses include damage to houses (10,319 complete damaged and 100,000 partially damaged), damage to 10,000 km of transport network, and $600{\mathrm{km}}^2$ of agriculture crops.⁷² All these facts about Kerala flood prove that there is an immense requirement of flood maps, which can help in postdisaster management activities and resilient disaster planning.

2.2.

Datasets

For mapping the spatial extent of 2018 Kerala flood event, both SAR and optical images are used. The Sentinel-1 SAR images that are available before, during, and after the flood are used. The Sentinel-1 Ground Range Detected (GRD) images obtained from GEE platform are already preprocessed⁷³ to terrain-corrected ${\sigma }^0$ images with a pixel size of $10\times 10\mathrm{m}$. As the Sentinel-1 orbital information is accurate, the terrain-corrected SAR images also have subpixel geolocation accuracy. The specification of the SAR images used in the research is given in Table 1. To validate the results, WorldView-3 optical images that are freely available from DigitalGlobe under open data program,⁷⁴ Sentinel-2, and Landsat-8 optical images are used. The WorldView-3 image acquired before the flood, i.e., on March 9, 2018, was only used to understand the land use pattern of the study area. Also, the field data collected during and after the flood are used for validation. The spatial resolution of WorldView-3, Sentinel-2, and Landsat-8 optical images is 0.5 m, 10 m, and 30 m, respectively.

Table 1

Sentinel-1 images fetched from GEE for 2018 Kerala flood mapping.

Apart from these imaging satellite datasets, the Global Precipitation Measurement (GPM) precipitation data, i.e., Integrated Multisatellite Retrievals for GPM (IMERG), and World Wildlife Fund (WWF) HydroSHEDS flow accumulation data are fetched into the GEE platform for further understanding of flooding pattern in the state. The spatial resolutions of IMERG precipitation and flow accumulation layers are 0.1 deg and 15 arc sec, respectively, and the detailed description of these datasets are given in Earth Engine Data Catalog.⁷⁵

3.1.

Sentinel-1 Data Fetching

GEE data catalog provides a wide range of remote sensing datasets, such as optical, thermal, infrared, and microwave images. In this work, the readily available level-2 processed high-resolution GRD Sentinel-1 images are utilized. A massive volume of 264 time-series Sentinel-1 images, which are available between January 1, 2015, to December 31, 2017, with a revisit time of 12 days interval are used for PWB mask creation. A total of seven Sentinel-1 images (Table 1) that are available at preflood, flood, and postflood time duration are used for Kerala flood mapping. All these images are directly fetched from GEE data catalog to GEE processing platform using Earth Engine operator (ee.ImageCollection) as an image collection. As the Sentinel-1 images available on GEE are already preprocessed to terrain-corrected radar backscatter ${\sigma }^0$, the time required in data downloading, SLC to GRD conversion, and data preprocessing time is saved in GEE cloud platform.

3.2.

Metadata Filtering

The Sentinel-1 image collection fetched from the previous step posses enormous information about each image. The metadata information, such as data acquisition mode, acquisition time, satellite pass, polarization, incident angle, and orbit number, is stored with image collection. Thus, the user can create a filter to restrict the data volume according to the objective of the work. In this work, all the Sentinel-1 images are filtered with descending pass, Interferometric Wide swath mode, and VV polarization parameters using .filterMetadata operator. Also, the .filterDate operator is used to fetch the images according to the image acquisition date. Subsequently, a spatial filter is also created using .filterBounds operator to limit the processing extent within the minimum bounding box of the study area. Thus metadata and spatial filtering play a crucial role in restricting the data volume and processing time.

3.3.

Permanent Water Body Mask Creation

The PWB mask used in the study was created using two datasets namely, (1) global surface water mapping (GSWM) layer with a spatial resolution of 30 m and (2) temporal mean SAR image with a pixel size of 10 m. The water occurrence band available with GSWM layer was prepared by Joint Research Centre using time-series Landsat optical images from 1984 to 2015. A threshold of $>80$ on water occurrence band was considered to obtained optical-based PWB. This indicates where surface water was present at 80% occurrence over 32 years⁶⁵ capturing interannual and intra-annual changes and variability. However, optical-based PWB has few limitations such as cloud cover and spatial resolution (30 m). Also, the water bodies may be expanded during monsoon season due to high inflow of rainwater and even change their course in due time after 2015. Therefore, the optical-based PWB mask needs to be updated for accounting seasonality change of water bodies after 2015, with the images acquired after 2015.

To update the optical-based PWB mask, all the 264 Sentinel-1 images are reduced to generate a single mean SAR image using .reduce operator. Several region of interests (ROIs) over water bodies are created in GEE platform that is spatially covering throughout the study area. An automatic Otsu’s thresholding algorithm is implemented on the mean SAR image in which water ROIs are used to find an optimum threshold. The SAR-based PWB mask is created in GEE platform by classifying the mean SAR image with all the pixels less than the optimum threshold.

Some water bodies that are not detected by optical remote sensing images are identified by SAR images due to its high sensitivity to dielectric constant.^76,77 Finally, the optical PWB mask is updated with SAR-based PWB mask to obtain the final PWB mask in GEE platform using .updateMask operator. The GEE .updateMask operator resamples the optical-based PWB to 10 m and merges with SAR-based PWB mask.

3.4.

Rapid Flood Mapping

The constellation of Sentinel 1A and 1B satellites can cover the globe with a revisit time of 3 to 12 days depending on the location and ESA program. For the southern part of India, the Sentinel-1 SAR images are acquired for every 12 days in most cases. During the disasters, the images may be acquired with even low revisit time (6 days) from both ascending and descending pass directions. There is a time difference between the SAR image acquisition time and ingestion time into GEE catalog. At present, this latency time varies from 30 to 70 h. The GEE developers have plans to reduce this latency time in coming future.⁷⁸ The SAR image ingestion latency time is not considered into GEE4FLOOD processing chain time. Therefore, the processing chain is named as rapid flood mapping rather than real-time or NRT flood mapping because of following points.

The workflow of rapid flood mapping is divided into four submodules, i.e., coregistration, speckle filtering, optimum threshold identification, and flood mapping. Based on the flood duration, the images are grouped as preflood, flood, and postflood images. Among the seven Sentinel-1 images used to understand the flood increasing and decreasing pattern, the images acquired on January 5, 2018, and September 2, 2018, are referred to as preflood and postflood images, respectively. The images acquired between July 15, 2018, and August 28, 2018, are designated as flood images.

For any change detection study, coregistration of master and slave images must be implemented to ensure proper spatial alignment of images. The master image has to be carefully chosen from archived images based on the parameters that influence the radar backscatter (${\sigma }^0$), such as the same satellite track, and atmospheric conditions. Here, the preflood image is treated as the master image since it is acquired on the clear sky, which is free from clouds and with the same satellite track of slave images. The other six sentinel-1 SAR images are treated as slave images.

As the first step of rapid flood mapping procedure, the slave images are precisely coregistered to a preflood image by local rubber sheet deformations in GEE platform using .register operator. Then, the coregistered flood and postflood ${\sigma }^0$ images are converted to natural scale from decibel (dB) scale. A user-defined function is created for refined Lee speckle filter with a $3\times 3$ window size. Subsequently, to reduce the inherent speckle noise in the SAR data, a refined Lee speckle filter is applied on coregistered flood and postflood images in GEE platform. After speckle filtering, the ${\sigma }^0$ images are converted back to dB scale. An Otsu’s automatic thresholding algorithm is implemented on the speckle filtered images. The detailed explanation is given in Sec. 4.1. The criterion function required for optimal threshold identification is the minimization of within-class variance of the nonflood and flood class pixels in the SAR image. It is the same as maximizing the between-class variance as given in Eq. (1):

3.5.

Validation Approach

The Kerala flood event happened at a regional level, which covers a larger spatial area. Thus, it is extremely challenging to collect the validation data throughout the entire state as it requires enormous resources. Therefore, the validation is conducted by both direct (field data) and indirect (secondary data) methods.

Secondary data include optical remote sensing images (Landsat-8 and Sentinel-2), 3-h GPM v5 IMERG rainfall data, WWF HydroSHEDs flow accumulation, and hydrologically conditioned digital elevation model (DEM) that are available in GEE catalog have been used in the study. A function is created in GEE platform to generate the daily rainfall by summation of hourly data. All these secondary data are clipped to the required geographical extent in GEE.

The field data were collected at various locations that are highly affected by flood on different time scales. A narrative survey was conducted during and immediately after the Kerala flood that recorded information such as location, flood depth, and photographs. Also, some photographs were obtained through online sources.⁷⁹ Thus, with the limited field data, validation of the results is confined to region C of Fig. 3. With the help of photographs and survey data, a ground truth flood area polygon vector file was digitized in GEE platform to compare against the results. Various performance indices were defined as follows to find the overall accuracy. The detailed validation approach is given in Sec. 4.3.

Fig. 3

Spatio-temporal flood maps derived from GEE on July 16, 28, August 9, 21, and 27, 2018, for three regions A, B, and C as highlighted on the Kerala state.

4.1.

Threshold Identification

A function is developed in GEE platform for Otsu’s thresholding algorithm as described in Ref. 43. From the water ROIs, the histograms are obtained as shown in Fig. 4 for all the SAR images that were given as the inputs to the Otsu’s function. In each plot of Fig. 4, the speckle filtered backscattered values (orange color histogram) are plotted against the nonfiltered values (blue color histogram). From Fig. 4(a), it is inferred that a bimodal histogram with the values ranging from $-4$ to $-26\mathrm{dB}$ is seen for the preflood image. The lower values ($-18$ to $-26\mathrm{dB}$) and higher values ($-14$ to $-4\mathrm{dB}$) of the histogram correspond to backscattering from the stable/calm water bodies and dry riverbed/dry sand, respectively. Subsequently, after the first heavy rainfall occurred in July 2018, the backscatter values started shifting toward the lower values of the histogram, as shown in Fig. 4(b). Here, the histogram is a left-skewed histogram rather than bimodal. This is due to the change of land surface condition from dry to wet and the movement of floodwater to the low elevated areas. From the visual interpretation of three optical images in the first row of Fig. 5, it is evident that the land surface condition is dry at the time of preflood image acquisition and persisted as wet surface even after the July 16, 2018, SAR image acquisition.

Fig. 4

Temporal histograms of ${\sigma }_{\mathrm{vv}}^0$ before and after speckle filtering of the images acquired on (a) January 5, (b) July 16, (c) July 28, (d) August 9, (e) August 21, and (f) August 27, 2018.

Fig. 5

Optical images acquired within the time duration (before/after) for every corresponding SAR image acquisition. B and C are the same regions given in Fig. 3.

In Figs. 4(c) and 4(d), the histograms are shifted toward left with the majority of the values in the range of $-24$ to $-18\mathrm{dB}$. At this phase, the land surface condition is saturated, and the floodwater receded through the drainage channels. Even though the Sentinel-2 optical images acquired on July 28, 2018, and August 7, 2018, as shown in Fig. 5 are partially cloudy, it is understood that agriculture fields in the low lying areas are affected by floodwater. In Fig. 4(e), all the values in the histogram are concentrated around $-22\mathrm{dB}$. This is due to a major rainfall spell that occurred during the August 14, 2018, to August 17, 2018, where the land surface got flooded and yielded low backscattering values. This is evident from the visual interpretation of optical images acquired on August 22, 2018, which is a day after the SAR image acquisition on August 21, 2018. In Fig. 4(f), the histogram is right skewed due to high backscattering values as the land surface is started to regain its original condition.

After generating the histograms of all the SAR images, optimum threshold values are obtained by implementing the Otsu’s function. The thresholds values obtained for temporal SAR images as shown in Table 1 are $-16.99$, $-14.94$, $-19.25$, $-19.11$, $-19.27$, $-21.12$, and $-21.12\mathrm{dB}$. However, there are some limitations with Otsu’s algorithm.^43,80 Those limitations include (1) the histogram needs to be bimodal; (2) size of flood and nonflood classes should be comparable; and (3) if the variances of flood and nonflood classes are larger than mean values of flood and nonflood classes, there is possibility of identifying an incorrect threshold. Therefore, due to these limitations, the threshold identified for the temporal images might yield underestimation of flood area. The overestimation of flood area is an extreme low possibility due to masking of water-like reflecting surfaces, such as layover, shadows, fore-shorting regions, and wet lands using preflood image.

4.2.

Spatio-Temporal Flood Mapping

The initial flood extent is obtained by considering the pixels whose values are lower than the optimum thresholds. The initial flood extent obtained from flood images and postflood images includes flood area along with water bodies and radar-like reflecting surfaces, such as highways and airport runways. The inclusion of water bodies in the flood area identification contributes to overestimation. Therefore, the PWB are masked out using PWB mask in GEE platform to reduce the overestimation. Apart from this, the radar-like reflecting surfaces and wetlands that are detected as flood area in the preflood image are masked out to obtain final flood extent. For a detailed understanding, three regions were selected based on the flood severity, location accessibility, and land use. The spatio-temporal flood extent of these regions A, B, and C are shown in Fig. 3, where red and blue color represent the flood area and PWB, respectively. The processed GPM daily rainfall data obtained in GEE cloud are clipped to the three regions A, B, and C in Fig. 3 separately and spatially averaged to get a single value. The quantitative estimation of flood area for the respective regions, as shown in Fig. 3, is plotted against the GPM daily rainfall obtained from GEE, as shown in Fig. 6.

Fig. 6

Flood area obtained from three regions A, B, and C, as shown in Fig. 3, is plotted against GPM IMERG daily rainfall from July 1, 2018, to September 5, 2018.

From Fig. 6, it is inferred that the flood area obtained from temporal SAR images follows the rainfall pattern obtained from GPM IMERG data. For example, after the major rainfall spell from August 14 to 19, 2018, the massive runoff generated from rainfall started receding. The same decreasing pattern is also observed in the flood area obtained from temporal SAR images acquired after the major rainfall spell. For the small geographical extent, as shown in Fig. 3, the total time consumed for spatio-temporal flood maps generation is $\sim 1\min$. The geographical area of region A in Fig. 3 is $33.62{\mathrm{km}}^2$, which is characterized by water bodies, semiurban areas mixed up with tree cover, and agricultural land uses. The geographical area of region B in Fig. 3 is $33.63{\mathrm{km}}^2$, which is mainly characterized by agricultural land use and a small portion of semiurban land use. The geographical area of region C in Fig. 3 is $11.78{\mathrm{km}}^2$, which is characterized by water bodies, urban areas, open spaces, and agricultural land uses, as shown in Fig. 7. These land features provide unique scattering phenomena when radar electromagnetic wave interactions with the targets are considered. Water has a high dielectric constant and acts as a specular reflector, which yield very low backscatter power return to the sensor. It makes the water bodies to appear as a dark feature in a SAR image (irrespective of polarization channel). For majority of other natural land features, particularly over agricultural areas, the radar wave produces diffused scattering. However, the dominant scattering of radar wave changes with respect to different crop growth stages. In the case of semiurban areas, double-bounce scattering becomes dominant, which is likely due to corner reflection from buildings and the structure of surrounding coconut trees in the study area.

Fig. 7

Optical remote sensing images acquired before the flood event for the regions, as shown in Fig. 3.

In the study area, SAR images acquired in July and August are used, where the crops are in a growing stage, which led to a high amount of volume scattering. From region A of Fig. 3, it is observed that the flood is predominantly detected among all the land uses except for urban areas. Many studies^50,55,81 stated that the major challenges for flood mapping in urban areas are layover, shadows, and the presence of radar-like reflecting surfaces. In this study, the test site is characterized by presence of radar-like reflecting surfaces and heterogeneous built forms, which results in layover and shadows. Therefore, the automatic thresholding algorithm is applied on preflood images to extract these land surfaces that lead overestimation of flood area. After identifying the flooded areas from flood images and post-flood images, the land surfaces obtained from preflood image are masked out to avoid overestimation. Also, the increased backscatter from flooded urban streets or high-density built-up areas results in underdetection, which needs to be tackled by high-resolution SAR imagery and interferometric SAR coherence.

After the first rainfall spell occurred in the second and third weeks of July 2018 as shown in Fig. 6, the land surface got flooded in all the regions of Fig. 3, which is evident from Fig. 5. Even though the geographical area and land use of region A and B of Fig. 3 are the same, the flood area obtained from region B is higher than region A. This is due to the low elevation of region B, high rainfall, and flow of drainage water to Vembanadu Lake. As the floodwater receded, there is a decrease in flood area until the major rainfall spell occurred during the third week of August 2018. After the major rainfall spell, there is a steep increase in the flood area and a decreasing trend afterward. Among all the selected locations, the region B is profoundly affected by the flood. For the better insight of water flow dynamics at a larger spatial area, the regional level flood mapping is carried out and described in Sec. 4.2.1.

4.2.1.

District-level flood mapping

To carry out district-level postdisaster strategic planning activities, a district-level flood damage assessment is required. Therefore, both quantitative and qualitative flood assessment are carried out in the GEE cloud platform at the district level. The reducer and quantification functions are developed in GEE to generate the spatio-temporal flood maps at the district level and to quantify the flood area, respectively. The district-wise temporal flood area obtained from GEE is shown in Fig. 8, and it is observed that Alappuzha is the worst affected, followed by Kottayam and Thrissur districts. Therefore, the flooding pattern (increasing and receding) in these three districts is studied in detail. The spatio-temporal distribution and quantitative estimation of flood area with the rainfall for these three districts are shown in Figs. 9 and 10, respectively. From Fig. 10, it is inferred that the rainfall pattern for all three districts is almost similar with an offset. The flood area derived from SAR images also follows the district-level rainfall pattern.

Fig. 8

District wise temporal flood area obtained from GEE for Kerala state.

Fig. 9

Spatio-temporal distribution of flood area for Thrissur, Alappuzha, and Kottayam districts. Red and blue color represent flood and perennial water bodies, respectively. Preflood ${\sigma }_0$ image derived from Sentinel-1 is used as a background.

Fig. 10

Flood area obtained from three districts as shown in Fig. 9 is plotted against their respective district wise GPM IMERG daily rainfall from July 1, 2018, to September 5, 2018.

During the first spell of rainfall in July 2018, the land surface got flooded among these three districts, and the floodwater started receding. An increase in the flood area is observed for Alappuzha and Kottayam districts, whereas decreasing trend for Thrissur district is obtained for the images acquired on July 28, 2018. Subsequently, a decreasing trend of flood area is observed in all three districts for images acquired on August 9, 2018. By the time of image acquisition on August 21, 2018, the land surface got saturated. After the heavy rainfall occurred during the third week of August 2018, there is a steep increase in the flood area among these three districts and subsequently decreasing trend. Even though the rainfall received by Thrissur during the third week of August 2018 is higher than the other two districts, the flood area obtained is lesser than the other two districts. This is due to the location of Thrissur in higher elevation region as shown in Appendix 2 (Sec. 6.2), and receding of the floodwater through the drainage channels is shown in Fig. 11.

Fig. 11

Flow accumulation maps generated from GEE cloud platform for Kerala state, Thrissur, Alappuzha, and Kottayam districts, respectively.

Apparently, Alappuzha has most flood-affected area due to its low elevation, backwater flow, dam water release, and the majority of the west-flowing rivers draining into Vembanadu Lake. Similarly, the state-level flood maps are generated and given in Appendix 1 (Sec. 6.1). At a district and a state level, the total time consumed for spatio-temporal flood maps generation is $\sim 2.5$ and 4 min, respectively.

4.3.

Validation

The results obtained from GEE for all image acquisition dates cannot be validated due to lack of ground truth data at that particular time. Most of the ground truth data are collected after the major rainfall spell in the third week of August 2018. So, validation is carried out on the flooded area obtained from the SAR image acquired at high flood stage on August 21, 2018. The ground truth data are generated from Sentinel-2 optical images acquired on August 22, 2018, along with field survey data. The ground truth vector file generated in GEE is rasterized with the same pixel size as Sentinel-1 ${\sigma }_0$ images and the same spatial extent as region C of Fig. 3. The PWB pixels are excluded from the performance metrics calculation. As mentioned in Sec. 3.5, various performance metrics are generated, as shown in Table 2.

Table 2

Flood mapping performance metrics for region C.

The TPR and TNR values are 66% and 98%, respectively, which indicate that for nonflood class, the ground truth is almost matching with algorithm output, whereas some discrepancies are seen for flood class. A 100% $F$-score value represents the perfect match between predicted flood area and ground truth flood area, i.e., perfect recall and precision for flood class, and a 0% represents a perfect mismatch. The accuracy and $F$-score values for the region C are 82.04% and 78.55%, respectively, which show the close match with the reference data for both classes and flood class, respectively. The value of $F$-score is promising considering the usage of Otsu’s automatic thresholding algorithm without any manual tuning or modification of algorithm. The difference in the amount of algorithm performance for both classes (flood and nonflood) and flood class alone is quite low. Further, high-resolution SAR images and timely acquisition of validation data may boost accuracy.

6.1.

State-Level Flood Mapping

The spatio-temporal flood maps generated in GEE cloud platform are shown in Fig. 12. From the figure it is observed that on September 2, 2018, the flood is receded and the land surface regained its original state.

Fig. 12

State-level spatio-temporal flood maps for 2018 flood event of Kerala state, India.

Fig. 13

Elevation map for Thrissur, Alappuzha, and Kottayam districts of Kerala state, India.

6.2.

Elevation

The HydroSHEDS hydrologically conditioned DEM available in GEE cloud platform with 3 arc sec spatial resolution is fetched and clipped to three districts as shown in Fig. 2. The elevation values for Thrissur, Alappuzha, and Kottayam are ranging from $-18$ to 1361, $-27$ to 107, and $-29$ to 1155 m, respectively (Fig. 13). While generating the flood maps at a regional scale, i.e., district or a state level, it is imperative that masking out higher elevation regions may optimize the computational time. Here, in this study, the higher slopes were masked out in GEE while computing the flood area at state and district level. This is achieved by fetching the HydroDEM into GEE from terrain information obtained using .Terrain() operator. The slope image was fetched from terrain dataset using .select() operator.