# Improved Consistency of an Automated Multisatellite Method for Extracting Temporal Changes in Flood Extent

Authors: Husniyah Binti Mahmud; Vaibhav Katiyar; Masahiko Nagai
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Malaysia is affected by floods almost every year. In this situation, high-frequency flood monitoring is crucial so that timely measures can be taken. However, the low revisit time of the satellites, as well as occlusion cast by clouds in optical images, limits the frequency of flood observation of the focused area. Therefore, this study proposes utilising multisatellite data from optical satellites such as Landsat 7, Landsat 8, and Moderate Resolution Imaging Spectroradiometer (MODIS), as well as Synthetic Aperture Radar (SAR) images from Advanced Land Observation Satellite (ALOS-2) and Sentinel-1, to increase observation of flood. The main objective was to utilize Otsu image segmentation over both optical and SAR satellite images to distinguish water and nonwater areas in each image separately. For this, modified normalized difference water index (MNDWI) for the optical satellite and total dual-polarization backscatter for SAR satellite images were estimated. The focused area has been divided into Universal Transverse Mercator (UTM) square-size grids of 30 pixels, and each satellite image was reprojected and resampled with a pixel size of 0.001° to standardize the flood map resolution. The second objective was to assess the potential of image fusion for increasing the consistency of water area extraction. Two pairs of satellite images with the same observation period covering a flood event in September 2017 in Perlis, Malaysia, were processed using 2D wavelet transform. Lastly, the temporal changes of the integrated surface water extent were evaluated by comparing the output from both multisatellite and fused images with the observed water level data from the Department of Drainage and Irrigation. The results showed that the proposed model can be used to estimate flood duration as well as to estimate the flood-related losses, especially in ungauged or data-poor regions.


---

## Body

## 1. Introduction

Floods are recurrent natural disasters in Malaysia as a result of heavy rainfall especially during the monsoon season. Massive efforts have been performed to mitigate flood disasters in the country [1]. In 2017, Perlis experienced two flood events, in May and September. The flood consequences caused damage to paddies, resulting in a 65% reduction in crop weight and causing farmers an average loss of about Malaysian Ringgit (RM) 1800/h, with almost 194 farmers affected [2]. Flood mapping is one of the sources of information used by decision makers, planners, and response teams to mitigate flood disasters. Multitemporal coverage of large areas using Earth observation satellites is widely used to observe differences in land cover before, during, and after flood events [3–6]. Large-scale information of flood situations can be retrieved from optical and radar satellite images. The Moderate Resolution Imaging Spectrometer (MODIS) onboard the Aqua and Terra satellites [7, 8], Land Remote-Sensing Satellite (Landsat) [9, 10], and Sentinel-2 [11] are example of optical sensors useful for mapping inundation areas. Sentinel-1 [12, 13] and ALOS-2 [14–16] are both carry radar sensors and valuable source of information for flood mapping. Optical images provide detailed information on land cover, while radar data can be used to provide prompt observations of flooded areas. Combining optical and radar images would improve the number of flood surface observations [17].During a flood event, the presence of a cloud is typically unavoidable, especially in the tropical region where cloud cover tends to vary with time [18, 19]. Synthetic aperture radar (SAR) is extremely useful for consistently monitoring flooded regions due to its wavelength that can penetrate clouds [17]. In general, the interaction of microwaves with the Earth’s surface depends on the SAR sensor’s wavelength, polarization, resolution, look angle, surface roughness, dielectric properties, and slope. The smooth water surface acts as a specular reflector, resulting in low backscatter and dark appearance on SAR images [20]. The selection of SAR polarization is important for obtaining high-accuracy surface water mapping. Single-polarization SAR is either horizontal transmission and horizontal reception (HH), horizontal transmission and vertical reception (HV), vertical transmission and vertical reception (VV), or vertical transmission and horizontal reception (VH). HH polarization has shown promising results for the delineation of open water bodies [21–23]. Manjusree et al. [24] assessed the backscatter response in multipolarizations using the C-band Radar Satellite (RADARSAT-2) SAR data and found a difference in the response of the flooded area discrimination depending on the polarization. The radiometric profile of HH is less backscattered than the profiles of HV and VV; however, the accuracy of HH is lower in heterogeneous regions than HV and VV [25]. Henry et al. [26], Horritt et al. [27], and Irwin et al. [28] addressed the use of dual polarization in flooded vegetation and their findings confirmed the advantages of dual-polarization observations in the detection of inundated area in vegetated regions. In dual-polarization SAR, there are two channels of intensity and phase information which allow the discrimination of scattering mechanisms [25]. The combined use of HH and HV (hereafter, HH + HV) has been shown to achieve the most reliable result for automatic classification mapping using a large dataset [28]. Furthermore, the Otsu thresholding [29] method shows good performance if the histogram has a bimodal distribution because a threshold will be chosen at the bottom of a sharp valley linking two tops that represent the object and the background. Moreover, HH + HV retains bimodality as opposed to using single-polarization (HH or HV) [30, 31].Optical images at visible and near-infrared (NIR) wavelengths can be used to extract water surfaces through visual interpretation. However, unlike SAR images, optical images are affected by weather conditions. The quality of the information in optical images is strongly affected by sunlight and inability of penetrating clouds. Additionally, compared with SAR images, optical images lack information about surface structure and roughness. The easiest and most reliable method to identify water areas from an optical image is using a water index calculated using two or more bands [32]. Several water indices have been developed in recent years [33–36]. The most widely used are the Normalized Difference Water Index (NDWI) and the Modified Normalized Difference Water Index (MNDWI) [37]. The NDWI uses green and NIR bands to enhance water features in remote-sensing images [33]. However, the MNDWI is a more reliable and stable water index since it uses short‐wave infrared (SWIR) instead of NIR to enhance open water features and effectively remove noise in images of urban areas, vegetation, and soil [34]. Husniyah and Nagai [38] compared NDWI, MNDWI, Automated Water Extraction Index with no shadow (AWEInsh), and Automated Water Extraction Index with shadow (AWEIsh) based on optical images from Sentinel-2, Landsat 7, Landsat 8, and MODIS MOD09A1; the study deduced that MNDWI was relatively stable and able to lessen the inconsistency of flood extent detection when combined with multiSAR images.Generally, the approaches to estimate flood extent based on remote-sensing images include simple visual interpretation, supervised classification, image texture algorithms, histogram thresholding, several multitemporal change detection methods, and active contour models. Thresholding is the preferable method to obtain an automatic binary classification of an image [12]. In most cases, thresholding is achieved by either global thresholding or local thresholding. Global thresholding uses the information of the gray-level histogram of an entire image dataset [39]. In contrast, local thresholding divides an entire region into subregions and uses a gray-level histogram of images in the subregions to derive threshold values [29]. Surface water extraction is mostly implemented using global thresholding. However, global thresholding often results in misclassification error since there is no difference in backscatter among different land cover types. Meanwhile, local thresholding can minimize the classification error between water and other land cover classes. The Otsu method, which uses uniformity and shape measures of an image, is preferable to determine threshold values than the other thresholding techniques. The selection of threshold values using the Otsu method was based on maximized separability of the resultant classes in the gray scale [29].In general, the study of flood monitoring comprises spatial and temporal dimensions [40]. This study proposes the use of the conceptual model to provide key information at different stages of a flood event: preflood, during flood, and postflood (Figure 1). The dynamic nature of flood in an inundated area is shown in Figure 1. The conceptual model is aligned with findings from the previous studies [32, 41], where the flood extent was maximum during flood and decreased drastically a few days after.Figure 1
Conceptual model of spatial and temporal variations for preflood, flood, and postflood.Presently, even though an abundance of near real-time data is available, the data and information still appear to be underutilized by decision makers in the disaster response phase due to a few constraints, such as the limitation on (1) time and capacity of a person to process new geospatial datasets, (2) accessibility of near real-time data, (3) compatibility of user platforms and geospatial data formats, (4) knowledge of the data availability and inadequacy of data latency, and (5) understanding of the end-user demand on the product and timing needs [42]. Despite various notable efforts, to this end, existing satellite technology or any single data product could not solve the current challenge from the decision-making standpoint. Hence, to fulfill decision makers’ needs, there is a demand to develop a process for improving the end product for an effective disaster response rather than developing a new system to improve flood operation. This study focuses on the standardization of multisource remote-sensing data by using the same processing platform. In this context, to obtain a single image that retains essential features of original images, the simple and robust fusion of images with the same observation period has been proposed in this study.The main objectives of this study are as follows: (i) classify the water extent based on the average Otsu threshold values calculated for individual grids with a bimodal histogram; (ii) then, the temporal changes of the flood extent are assessed by analyzing the grids that indicate flood-affected areas; (iii) the analysis proceeds to apply 2D wavelet transform to fuse two images with the same observation period. Finally, the temporal distribution of water classification was compared with the ground water level data.


## 2. Materials and Methods

### 2.1. Study Area

The study area is located in the state of Perlis, Malaysia, which was affected by floods in 2017. The study area extents ranges between 6°43′19″ N latitude and 100°07′59″ E longitude, covering a land area of about 795 km2. There are two distinct seasons: a wet season with heavy rainfall and a hot season with prolonged drought. The wet season, which is associated with the southwestern monsoon, usually starts between April and May and lasts until September or October. Meanwhile, the dry season, which is influenced by the northeast monsoon, lasts from December to March. The mean annual rainfall is between 2000 and 3000 mm, and the average annual temperature is 27°C. Figure 2 shows an elevation map of the study area, including the locations of the affected villages as specified in a Department of Drainage and Irrigation (DID) report [43].Figure 2
Elevation map of the study area. Location and flood information of some affected residential areas represented by red dots are shown on the map. The water level data from the ground station are indicated by blue triangles.In two districts in Perlis, Kangar and Arau, the DID reported two occurrences of flooding in 2017 and 2018. In these districts, continuous rainfall was reported from 15 September 2017 until early October 2017. High atmospheric humidity caused the rainfall during this time as well as Typhoon Doksuri which passed over Northern Vietnam on 15 September. Malaysia was affected by the relatively weak southwestern monsoon and the Northern Hemisphere summer monsoon in this monsoon’s final phase. The tail effect of the Typhoon Doksuri caused wind gusts in Malaysia, especially in northern Peninsular Malaysia [44]. This situation contributed to the episodes of heavy rain that occurred throughout Perlis in September 2017 and caused some areas to experience floods that lasted for two weeks. Additionally, the floods were exacerbated by high tides. Figure 3 shows photographs of the flood event in Perlis on 22 September 2017 and 23 September 2017.Figure 3
Photographs of the flood disaster in Perlis on 22 September 2017 and 23 September 2017.


### 2.2. Data

For the mapping of flood extent and comparison with ground truth, the following data were used. The mapping of flood extend for before, during, and after the event was acquired from two sets of remotely sensed data, optical and SAR. The daily changes of the flood extent were compared with water level ground station data provided by DID.


#### 2.2.1. Satellite Data

In this study, 22 SAR and optical images collected over the study region between 4 September 2017 and 30 October 2017 were acquired (Table1). Four ALOS-2 SAR images in the HH and HV polarization modes were acquired from the Japan Aerospace Exploration Agency (JAXA) (https://auig2.jaxa.jp/ips/home). The acquired ALOS-2 images were in the L band with a wavelength of 22.9 cm and spatial resolutions of 6.25 m and 25 m. The incidence angle of the ALOS-2 images varied between 8° and 70°. A total of 10 Sentinel-1 SAR images with VV and VH polarization with ascending and descending modes were used (https://scihub.copernicus.eu/dhus/#/home); these were Level-1 processed ground range detected (GRD) products collected in interferometric wide swath (IW) mode. The Sentinel-1 images were in the C-band with incidence angles between 29° and 46° and wavelength 5.6 cm.Table 1
Details of the SAR and optical satellite images used in this study.
SatelliteSensorAcquisition dateResolution (m)Pass (polarization)SourceProjectionCentral wavelength (nm)ALOS-2PALSAR-24/9/201725Descending (HH & HV)JAXAUTML bandALOS-2PALSAR-29/9/201725Descending (HH & HV)JAXAUTML bandALOS-2PALSAR-228/9/20176.25Ascending (HH & HV)JAXAUTML bandALOS-2PALSAR-221/10/201725Descending (HH & HV)JAXAUTML bandSentinel-1C-SAR6/9/201710Descending (VV & VH)ESAEllipsoidL bandSentinel-1C-SAR12/9/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR18/9/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR24/9/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR30/9/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR6/10/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR12/10/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR18/10/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR24/10/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR30/10/201710Ascending (VV & VH)ESAEllipsoidC bandLandsat 7ETM+7/9/201730USGSUTMGreen: 560SWIR1: 1648Landsat 8OLI1/10/201730USGSUTMGreen: 562SWIR1: 1610Landsat 8OLI17/10/201730USGSUTMGreen: 562SWIR1: 1610MODIS terraTerra7/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua12/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua30/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS terraTerra18/10/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua19/10/2017500USGSSinusoidalGreen: 555SWIR1: 1640The cloud-covered images less than 20% [22] were obtained from Landsat 7, Landsat 8, daily MODIS Terra MOD09GA, and MODIS Aqua MYD09GA. The Landsat 7 Enhanced Thematic Mapper (ETM+) and the Landsat 8 Operational Land Imager (OLI) have approximately 710 km sun-synchronous circular 98.2° inclined orbit with a revisit time of 16 days, but these overpass each other by eight days [45]. MOD09GA and MYD09GA version 6 provide bands 1–7 in a daily gridded L2G product in a sinusoidal projection, including reflectance values with a resolution of 500 m and observation and geolocation statistics with a resolution of 1 km.Three sets of SAR and optical images were used for the prefood, flood, and postflood stages (Figure4). The preflood stage covers from 4 to 21 September 2017, the flood stage covers from 21 September to 4 October 2017, and the postflood stage covers from 4 to 30 October 2017. One SAR image and one optical image are available for four dates, namely, 7 September (Landsat 7 and MODIS Terra), 12 September (Sentinel-1 and MODIS Aqua), 30 September (Sentinel-1 and MODIS Aqua), and 18 October (Sentinel-1 and MODIS Terra). The images from each of these dates were used for image fusion analysis.Figure 4
Temporal coverage of the satellite images used in this study.


#### 2.2.2. Ground Station Data

Water level data from ground station were available at DID website (https://publicinfobanjir.water.gov.my/aras-air/data-paras-air/?lang=en) and opened for public access. The only provided information is 15 minutes’ resolution water level, monthly and daily rainfall data. Unfortunately, information on flood extent and discharge was not available for free, and only two ground-measured water level stations are located in the flood-affected grid (Figure 2). However, there were no missing data during flood.


### 2.3. Image Preprocessing

The Sentinel-1 and ALOS-2 satellite image preprocessing was implemented using the Sentinel Application Platform (SNAP) version 7, while the images from Landsat 7, Landsat 8, and MODIS were preprocessed using Google Earth Engine (GEE). The image analysis was made using Python 3.6, and mapping was done in QGIS software. Image preprocessing was carried out to obtain high accuracy when implementing the Otsu method. From the SAR images, namely, the Sentinel-1 (VV/VH) and ALOS-2 (HH/HV) images, the standard preprocessing procedure was followed using the SNAP software [46]. The preprocessing method of Sentinel-1 comprises six steps to obtain a reliable and stable image for later use: (i) automatically downloads and applies an orbit file to each Sentinel-1 image to obtain precise satellite position and velocity information, (ii) thermal noise is removed from the SAR images by eliminating the noise from an intensity image using the SNAP thermal noise removal tool, (iii) the SAR data are calibrated for quantitative use so that the pixel values of the SAR images show the radar backscatter of the reflecting surface, (iv) the Lee filter (3 × 3 window) is applied to reduce speckle [47], (v) terrain correction is performed to correct distortion related to the side-looking geometry using the range-Doppler terrain correction in SNAP, and (vi) the last preprocessing step is to convert the unitless backscatter coefficient to decibel (dB) using a logarithmic transformation. The preprocessing of ALOS-2 images was performed through radiometric calibration, speckle filtering, and summation of the HH and HV polarization. Radiometric calibration was used to convert backscatter intensity (recorded as digital numbers by the satellite) into (dB) format, speckle filtering was performed using the Lee filter (3 × 3 windows), and the backscatter images for each polarization were summed (VV + VH for Sentinel-1 and HH + HV for ALOS-2).For the quantitative analyses of the Landsat, MODIS Terra, and Aqua images, atmospheric correction was performed [45]. The preprocessing of optical images was fully implemented in Google Earth Engine (GEE). Landsat 7 ETM+ contains data gaps due to the failure of the Scan Line Corrector (SLC), which were filled in GEE. For Landsat 7 masked for “cloud shadows,” it is bit 3, while “clouds” are bit 5, respectively. Cloud mask for Landsat 8 is bit 8-9, and for MODIS state_1 km, it was masked and assigned according to the following: cloudBitMask = 1, cloudshadowBitMask = 1 << 2, cirruslBitMask = 1 << 9, cloudInternalBitMask = 1 << 10, and adjacentBitMask = 1 << 13. The equation of MNDWI (equation (1)) was used to extract water features using the green and SWIR bands:(1)MNDWI=green−SWIRgreen+SWIR.


### 2.4. Grid Resolution and Image Resampling

This study uses multisatellite images with different orbital platforms. These satellites have different properties such as daily revisit time but lower spatial resolutions (500 m), such as MODIS Terra and Aqua. Other satellites with a better spatial resolution have a lower temporal resolution, such as Landsat 7, Landsat 8, Sentinel-1, and ALOS-2. Hence, image alignment and resampling is an essential step for this multitemporal and multisatellite study. In this study, the input remote-sensing images were resampled using bicubic spline interpolation with a pixel size of 0.001°. The size of the pixel was decided considering the spread of the ground data and convenience for the user since high spatial resolution might cause inconvenience during access [48]. A vector grid (Figure 5) was built as a reference for image projection and alignment. The vector grid was reprojected as a Universal Transverse Mercator (UTM) zone 47 N coordinate system. The proposed method required overlapping images, and to do this, the pixel of each image is aligned with the corresponding pixel of the referred vector grid. Then, each resampled and reprojected image was cropped to a 30 × 30 pixel (w×h) following the reference vector grid for further processing.Figure 5
A reference vector grid of 30 × 30 pixel (w×h) size.


### 2.5. Otsu Thresholding

The thresholding value is the minimum class variance of the image pixel frequency, and it is decided by separating the image pixels into two classes: foreground and background [20]. The first step of the Otsu algorithm is to draw the image pixel histogram and then calculate the various intensity levels. Next, the possible threshold to divide pixels into two classes, background and foreground, are determined. Wb, Wf, μb,μf, σb, and σfindicate background and foreground weights, mean, and variance, respectively. The equation for within-class variance, σW2, is simply a sum of variance with multiplication of its associated weight and is shown as [29](2)σw2=Wbσb2+Wfσf2.The equation for between-class variance is shown in equation (3). The between-class variance formula includes the multiplication of background and foreground weight with a square difference of the mean background pixels, μb, with the mean of the foreground pixels, μf [29]:(3)σB2=WbWfμb−μf2.The calculation of equations (2) and (3) was done automatically using a Python script with GDAL and Numpy-based Python library.


### 2.6. 2D Wavelet Transform for Image Fusion

The fusion process begins by decomposing the image volume to frequency bands. Wavelet transform decomposes images into four quadrants of high- and low-frequency bands. The input images were transformed into four wavelet coefficients called HH (containing high-frequency diagonal information), LH (containing high-frequency horizontal information), HL (containing high-frequency vertical information), and LL (containing low-frequency approximation information). Then, the frequency bands of each image are analyzed by the set fusion rules to determine which one can be combined and which one has to be removed from the final volume of coefficients. To produce the fused and combined coefficient map of the input images, “the average” was used, and then, an inverse wavelet transform was implemented to reconstruct the fused image [49]. Figure 6 shows the schematic diagram of 2D wavelet transform fusion process. The 2D wavelet transform was performed using PyWavelets which is wavelet transform software for Python.Figure 6
2D wavelet transform fusion process [50].


### 2.7. Flood Extent Detection Using SAR Dual Polarization

The process of flood extent detection using SAR dual polarization is shown in the flowchart in Figure7. The following three main steps were implemented to achieve surface water extraction using dual-polarization SAR: (1) preprocessing the image, (2) dividing the study area into an aligned grid, and (3) classifying surface water.Figure 7
Flowchart of the automatic flood mapping process.Preprocessing of Sentinel-1 was conducted as illustrated in Section2.3, and the image was overlaid and cropped per grid. Then, Otsu thresholding was implemented on each grid image. Grayscale digital images were converted to binary images by automatic clustering based on image thresholding using Python 3.6. The selected grid image with a bimodal histogram shape was averaged to represent the extent of flooding over the whole study area (Figure 8).Figure 8
Extent of the September 2017 flood extracted from the Sentinel-1 satellite images.


### 2.8. Grid-Based Flood Temporal Changes

The second analysis involved the assessment of the temporal changes of floodwater area in the affected grid. The following four main steps were taken to achieve this objective: (1) preprocessing the images, (2) resampling and dividing the study area into an aligned grid, (3) identifying flood water extent time changes, and (4) fusion by a grid.The same process, as mentioned in the previous section, was applied to classify water surfaces (Figure8).Referring to Figure4, the image with the same observation date was used for the image fusion process. There are two kinds of fusion based on the availability of data: optical with another optical image (for example, Landsat 7 and MODIS Terra on 7 September 2017) and optical image with SAR image (for example, Sentinel-1 and MODIS Terra on 12 September 2017). The overall fusion processing goes through image registration and preprocessing. MNDWI was calculated for optical images; however, the backscatter of SAR is not directly comparable with the MNDWI. Hence, the MNDWI information was rescaled by inverting the minimum and maximum index values so that the pixel close to the value of −1 represents water while that close to 1 was classified as nonwater [51]. Then, 2D wavelet transform was applied to the rescaled MNDWI optical image and SAR. This study used wavelet transform to fuse two images using Python 3.6. Finally, Otsu thresholding was applied on the fused images.


### 2.9. Performance Evaluation for Data-Scarce Areas

To assess the accuracy of the surface water area derived from the multisatellite time series, a correlation coefficient (equation (4)) was used to measure the correlation between the satellite-derived water area and the ground-measured water area:(4)Correlation coefficient=∑i=1NPSi−P¯SPOi−PO¯∑i=1NPSi−PS¯2∑i=1NPOi−PO¯2,where N is the number of samples, PSi is the estimated surface water extent value, P¯Si is the mean of estimated surface water extent, POi is the observed water level value, and P¯Oi is the mean of observed water level value. Additionally, the root mean square error (RMSE) was used to measure the magnitude of the average error in the extracted surface water area. The RMSE (equation (5)) is important for detecting undesirable large errors [52]:(5)Root Mean Square ErrorRMSE=∑i=1NPSi−POi2N1/2.Finally, the mean bias error (MBE) (equation (6)) was used to measure the mean absolute value of the error. If the predicted value is larger than the observed value, the mean bias error will be positive, and vice versa:(6)Mean Bias Error=∑i≡1NPSI−POIN.


## 2.1. Study Area

The study area is located in the state of Perlis, Malaysia, which was affected by floods in 2017. The study area extents ranges between 6°43′19″ N latitude and 100°07′59″ E longitude, covering a land area of about 795 km2. There are two distinct seasons: a wet season with heavy rainfall and a hot season with prolonged drought. The wet season, which is associated with the southwestern monsoon, usually starts between April and May and lasts until September or October. Meanwhile, the dry season, which is influenced by the northeast monsoon, lasts from December to March. The mean annual rainfall is between 2000 and 3000 mm, and the average annual temperature is 27°C. Figure 2 shows an elevation map of the study area, including the locations of the affected villages as specified in a Department of Drainage and Irrigation (DID) report [43].Figure 2
Elevation map of the study area. Location and flood information of some affected residential areas represented by red dots are shown on the map. The water level data from the ground station are indicated by blue triangles.In two districts in Perlis, Kangar and Arau, the DID reported two occurrences of flooding in 2017 and 2018. In these districts, continuous rainfall was reported from 15 September 2017 until early October 2017. High atmospheric humidity caused the rainfall during this time as well as Typhoon Doksuri which passed over Northern Vietnam on 15 September. Malaysia was affected by the relatively weak southwestern monsoon and the Northern Hemisphere summer monsoon in this monsoon’s final phase. The tail effect of the Typhoon Doksuri caused wind gusts in Malaysia, especially in northern Peninsular Malaysia [44]. This situation contributed to the episodes of heavy rain that occurred throughout Perlis in September 2017 and caused some areas to experience floods that lasted for two weeks. Additionally, the floods were exacerbated by high tides. Figure 3 shows photographs of the flood event in Perlis on 22 September 2017 and 23 September 2017.Figure 3
Photographs of the flood disaster in Perlis on 22 September 2017 and 23 September 2017.


## 2.2. Data

For the mapping of flood extent and comparison with ground truth, the following data were used. The mapping of flood extend for before, during, and after the event was acquired from two sets of remotely sensed data, optical and SAR. The daily changes of the flood extent were compared with water level ground station data provided by DID.


### 2.2.1. Satellite Data

In this study, 22 SAR and optical images collected over the study region between 4 September 2017 and 30 October 2017 were acquired (Table1). Four ALOS-2 SAR images in the HH and HV polarization modes were acquired from the Japan Aerospace Exploration Agency (JAXA) (https://auig2.jaxa.jp/ips/home). The acquired ALOS-2 images were in the L band with a wavelength of 22.9 cm and spatial resolutions of 6.25 m and 25 m. The incidence angle of the ALOS-2 images varied between 8° and 70°. A total of 10 Sentinel-1 SAR images with VV and VH polarization with ascending and descending modes were used (https://scihub.copernicus.eu/dhus/#/home); these were Level-1 processed ground range detected (GRD) products collected in interferometric wide swath (IW) mode. The Sentinel-1 images were in the C-band with incidence angles between 29° and 46° and wavelength 5.6 cm.Table 1
Details of the SAR and optical satellite images used in this study.
SatelliteSensorAcquisition dateResolution (m)Pass (polarization)SourceProjectionCentral wavelength (nm)ALOS-2PALSAR-24/9/201725Descending (HH & HV)JAXAUTML bandALOS-2PALSAR-29/9/201725Descending (HH & HV)JAXAUTML bandALOS-2PALSAR-228/9/20176.25Ascending (HH & HV)JAXAUTML bandALOS-2PALSAR-221/10/201725Descending (HH & HV)JAXAUTML bandSentinel-1C-SAR6/9/201710Descending (VV & VH)ESAEllipsoidL bandSentinel-1C-SAR12/9/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR18/9/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR24/9/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR30/9/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR6/10/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR12/10/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR18/10/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR24/10/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR30/10/201710Ascending (VV & VH)ESAEllipsoidC bandLandsat 7ETM+7/9/201730USGSUTMGreen: 560SWIR1: 1648Landsat 8OLI1/10/201730USGSUTMGreen: 562SWIR1: 1610Landsat 8OLI17/10/201730USGSUTMGreen: 562SWIR1: 1610MODIS terraTerra7/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua12/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua30/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS terraTerra18/10/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua19/10/2017500USGSSinusoidalGreen: 555SWIR1: 1640The cloud-covered images less than 20% [22] were obtained from Landsat 7, Landsat 8, daily MODIS Terra MOD09GA, and MODIS Aqua MYD09GA. The Landsat 7 Enhanced Thematic Mapper (ETM+) and the Landsat 8 Operational Land Imager (OLI) have approximately 710 km sun-synchronous circular 98.2° inclined orbit with a revisit time of 16 days, but these overpass each other by eight days [45]. MOD09GA and MYD09GA version 6 provide bands 1–7 in a daily gridded L2G product in a sinusoidal projection, including reflectance values with a resolution of 500 m and observation and geolocation statistics with a resolution of 1 km.Three sets of SAR and optical images were used for the prefood, flood, and postflood stages (Figure4). The preflood stage covers from 4 to 21 September 2017, the flood stage covers from 21 September to 4 October 2017, and the postflood stage covers from 4 to 30 October 2017. One SAR image and one optical image are available for four dates, namely, 7 September (Landsat 7 and MODIS Terra), 12 September (Sentinel-1 and MODIS Aqua), 30 September (Sentinel-1 and MODIS Aqua), and 18 October (Sentinel-1 and MODIS Terra). The images from each of these dates were used for image fusion analysis.Figure 4
Temporal coverage of the satellite images used in this study.


### 2.2.2. Ground Station Data

Water level data from ground station were available at DID website (https://publicinfobanjir.water.gov.my/aras-air/data-paras-air/?lang=en) and opened for public access. The only provided information is 15 minutes’ resolution water level, monthly and daily rainfall data. Unfortunately, information on flood extent and discharge was not available for free, and only two ground-measured water level stations are located in the flood-affected grid (Figure 2). However, there were no missing data during flood.


## 2.2.1. Satellite Data

In this study, 22 SAR and optical images collected over the study region between 4 September 2017 and 30 October 2017 were acquired (Table1). Four ALOS-2 SAR images in the HH and HV polarization modes were acquired from the Japan Aerospace Exploration Agency (JAXA) (https://auig2.jaxa.jp/ips/home). The acquired ALOS-2 images were in the L band with a wavelength of 22.9 cm and spatial resolutions of 6.25 m and 25 m. The incidence angle of the ALOS-2 images varied between 8° and 70°. A total of 10 Sentinel-1 SAR images with VV and VH polarization with ascending and descending modes were used (https://scihub.copernicus.eu/dhus/#/home); these were Level-1 processed ground range detected (GRD) products collected in interferometric wide swath (IW) mode. The Sentinel-1 images were in the C-band with incidence angles between 29° and 46° and wavelength 5.6 cm.Table 1
Details of the SAR and optical satellite images used in this study.
SatelliteSensorAcquisition dateResolution (m)Pass (polarization)SourceProjectionCentral wavelength (nm)ALOS-2PALSAR-24/9/201725Descending (HH & HV)JAXAUTML bandALOS-2PALSAR-29/9/201725Descending (HH & HV)JAXAUTML bandALOS-2PALSAR-228/9/20176.25Ascending (HH & HV)JAXAUTML bandALOS-2PALSAR-221/10/201725Descending (HH & HV)JAXAUTML bandSentinel-1C-SAR6/9/201710Descending (VV & VH)ESAEllipsoidL bandSentinel-1C-SAR12/9/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR18/9/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR24/9/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR30/9/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR6/10/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR12/10/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR18/10/201710Ascending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR24/10/201710Descending (VV & VH)ESAEllipsoidC bandSentinel-1C-SAR30/10/201710Ascending (VV & VH)ESAEllipsoidC bandLandsat 7ETM+7/9/201730USGSUTMGreen: 560SWIR1: 1648Landsat 8OLI1/10/201730USGSUTMGreen: 562SWIR1: 1610Landsat 8OLI17/10/201730USGSUTMGreen: 562SWIR1: 1610MODIS terraTerra7/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua12/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua30/9/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS terraTerra18/10/2017500USGSSinusoidalGreen: 555SWIR1: 1640MODIS aquaAqua19/10/2017500USGSSinusoidalGreen: 555SWIR1: 1640The cloud-covered images less than 20% [22] were obtained from Landsat 7, Landsat 8, daily MODIS Terra MOD09GA, and MODIS Aqua MYD09GA. The Landsat 7 Enhanced Thematic Mapper (ETM+) and the Landsat 8 Operational Land Imager (OLI) have approximately 710 km sun-synchronous circular 98.2° inclined orbit with a revisit time of 16 days, but these overpass each other by eight days [45]. MOD09GA and MYD09GA version 6 provide bands 1–7 in a daily gridded L2G product in a sinusoidal projection, including reflectance values with a resolution of 500 m and observation and geolocation statistics with a resolution of 1 km.Three sets of SAR and optical images were used for the prefood, flood, and postflood stages (Figure4). The preflood stage covers from 4 to 21 September 2017, the flood stage covers from 21 September to 4 October 2017, and the postflood stage covers from 4 to 30 October 2017. One SAR image and one optical image are available for four dates, namely, 7 September (Landsat 7 and MODIS Terra), 12 September (Sentinel-1 and MODIS Aqua), 30 September (Sentinel-1 and MODIS Aqua), and 18 October (Sentinel-1 and MODIS Terra). The images from each of these dates were used for image fusion analysis.Figure 4
Temporal coverage of the satellite images used in this study.


## 2.2.2. Ground Station Data

Water level data from ground station were available at DID website (https://publicinfobanjir.water.gov.my/aras-air/data-paras-air/?lang=en) and opened for public access. The only provided information is 15 minutes’ resolution water level, monthly and daily rainfall data. Unfortunately, information on flood extent and discharge was not available for free, and only two ground-measured water level stations are located in the flood-affected grid (Figure 2). However, there were no missing data during flood.


## 2.3. Image Preprocessing

The Sentinel-1 and ALOS-2 satellite image preprocessing was implemented using the Sentinel Application Platform (SNAP) version 7, while the images from Landsat 7, Landsat 8, and MODIS were preprocessed using Google Earth Engine (GEE). The image analysis was made using Python 3.6, and mapping was done in QGIS software. Image preprocessing was carried out to obtain high accuracy when implementing the Otsu method. From the SAR images, namely, the Sentinel-1 (VV/VH) and ALOS-2 (HH/HV) images, the standard preprocessing procedure was followed using the SNAP software [46]. The preprocessing method of Sentinel-1 comprises six steps to obtain a reliable and stable image for later use: (i) automatically downloads and applies an orbit file to each Sentinel-1 image to obtain precise satellite position and velocity information, (ii) thermal noise is removed from the SAR images by eliminating the noise from an intensity image using the SNAP thermal noise removal tool, (iii) the SAR data are calibrated for quantitative use so that the pixel values of the SAR images show the radar backscatter of the reflecting surface, (iv) the Lee filter (3 × 3 window) is applied to reduce speckle [47], (v) terrain correction is performed to correct distortion related to the side-looking geometry using the range-Doppler terrain correction in SNAP, and (vi) the last preprocessing step is to convert the unitless backscatter coefficient to decibel (dB) using a logarithmic transformation. The preprocessing of ALOS-2 images was performed through radiometric calibration, speckle filtering, and summation of the HH and HV polarization. Radiometric calibration was used to convert backscatter intensity (recorded as digital numbers by the satellite) into (dB) format, speckle filtering was performed using the Lee filter (3 × 3 windows), and the backscatter images for each polarization were summed (VV + VH for Sentinel-1 and HH + HV for ALOS-2).For the quantitative analyses of the Landsat, MODIS Terra, and Aqua images, atmospheric correction was performed [45]. The preprocessing of optical images was fully implemented in Google Earth Engine (GEE). Landsat 7 ETM+ contains data gaps due to the failure of the Scan Line Corrector (SLC), which were filled in GEE. For Landsat 7 masked for “cloud shadows,” it is bit 3, while “clouds” are bit 5, respectively. Cloud mask for Landsat 8 is bit 8-9, and for MODIS state_1 km, it was masked and assigned according to the following: cloudBitMask = 1, cloudshadowBitMask = 1 << 2, cirruslBitMask = 1 << 9, cloudInternalBitMask = 1 << 10, and adjacentBitMask = 1 << 13. The equation of MNDWI (equation (1)) was used to extract water features using the green and SWIR bands:(1)MNDWI=green−SWIRgreen+SWIR.


## 2.4. Grid Resolution and Image Resampling

This study uses multisatellite images with different orbital platforms. These satellites have different properties such as daily revisit time but lower spatial resolutions (500 m), such as MODIS Terra and Aqua. Other satellites with a better spatial resolution have a lower temporal resolution, such as Landsat 7, Landsat 8, Sentinel-1, and ALOS-2. Hence, image alignment and resampling is an essential step for this multitemporal and multisatellite study. In this study, the input remote-sensing images were resampled using bicubic spline interpolation with a pixel size of 0.001°. The size of the pixel was decided considering the spread of the ground data and convenience for the user since high spatial resolution might cause inconvenience during access [48]. A vector grid (Figure 5) was built as a reference for image projection and alignment. The vector grid was reprojected as a Universal Transverse Mercator (UTM) zone 47 N coordinate system. The proposed method required overlapping images, and to do this, the pixel of each image is aligned with the corresponding pixel of the referred vector grid. Then, each resampled and reprojected image was cropped to a 30 × 30 pixel (w×h) following the reference vector grid for further processing.Figure 5
A reference vector grid of 30 × 30 pixel (w×h) size.


## 2.5. Otsu Thresholding

The thresholding value is the minimum class variance of the image pixel frequency, and it is decided by separating the image pixels into two classes: foreground and background [20]. The first step of the Otsu algorithm is to draw the image pixel histogram and then calculate the various intensity levels. Next, the possible threshold to divide pixels into two classes, background and foreground, are determined. Wb, Wf, μb,μf, σb, and σfindicate background and foreground weights, mean, and variance, respectively. The equation for within-class variance, σW2, is simply a sum of variance with multiplication of its associated weight and is shown as [29](2)σw2=Wbσb2+Wfσf2.The equation for between-class variance is shown in equation (3). The between-class variance formula includes the multiplication of background and foreground weight with a square difference of the mean background pixels, μb, with the mean of the foreground pixels, μf [29]:(3)σB2=WbWfμb−μf2.The calculation of equations (2) and (3) was done automatically using a Python script with GDAL and Numpy-based Python library.


## 2.6. 2D Wavelet Transform for Image Fusion

The fusion process begins by decomposing the image volume to frequency bands. Wavelet transform decomposes images into four quadrants of high- and low-frequency bands. The input images were transformed into four wavelet coefficients called HH (containing high-frequency diagonal information), LH (containing high-frequency horizontal information), HL (containing high-frequency vertical information), and LL (containing low-frequency approximation information). Then, the frequency bands of each image are analyzed by the set fusion rules to determine which one can be combined and which one has to be removed from the final volume of coefficients. To produce the fused and combined coefficient map of the input images, “the average” was used, and then, an inverse wavelet transform was implemented to reconstruct the fused image [49]. Figure 6 shows the schematic diagram of 2D wavelet transform fusion process. The 2D wavelet transform was performed using PyWavelets which is wavelet transform software for Python.Figure 6
2D wavelet transform fusion process [50].


## 2.7. Flood Extent Detection Using SAR Dual Polarization

The process of flood extent detection using SAR dual polarization is shown in the flowchart in Figure7. The following three main steps were implemented to achieve surface water extraction using dual-polarization SAR: (1) preprocessing the image, (2) dividing the study area into an aligned grid, and (3) classifying surface water.Figure 7
Flowchart of the automatic flood mapping process.Preprocessing of Sentinel-1 was conducted as illustrated in Section2.3, and the image was overlaid and cropped per grid. Then, Otsu thresholding was implemented on each grid image. Grayscale digital images were converted to binary images by automatic clustering based on image thresholding using Python 3.6. The selected grid image with a bimodal histogram shape was averaged to represent the extent of flooding over the whole study area (Figure 8).Figure 8
Extent of the September 2017 flood extracted from the Sentinel-1 satellite images.


## 2.8. Grid-Based Flood Temporal Changes

The second analysis involved the assessment of the temporal changes of floodwater area in the affected grid. The following four main steps were taken to achieve this objective: (1) preprocessing the images, (2) resampling and dividing the study area into an aligned grid, (3) identifying flood water extent time changes, and (4) fusion by a grid.The same process, as mentioned in the previous section, was applied to classify water surfaces (Figure8).Referring to Figure4, the image with the same observation date was used for the image fusion process. There are two kinds of fusion based on the availability of data: optical with another optical image (for example, Landsat 7 and MODIS Terra on 7 September 2017) and optical image with SAR image (for example, Sentinel-1 and MODIS Terra on 12 September 2017). The overall fusion processing goes through image registration and preprocessing. MNDWI was calculated for optical images; however, the backscatter of SAR is not directly comparable with the MNDWI. Hence, the MNDWI information was rescaled by inverting the minimum and maximum index values so that the pixel close to the value of −1 represents water while that close to 1 was classified as nonwater [51]. Then, 2D wavelet transform was applied to the rescaled MNDWI optical image and SAR. This study used wavelet transform to fuse two images using Python 3.6. Finally, Otsu thresholding was applied on the fused images.


## 2.9. Performance Evaluation for Data-Scarce Areas

To assess the accuracy of the surface water area derived from the multisatellite time series, a correlation coefficient (equation (4)) was used to measure the correlation between the satellite-derived water area and the ground-measured water area:(4)Correlation coefficient=∑i=1NPSi−P¯SPOi−PO¯∑i=1NPSi−PS¯2∑i=1NPOi−PO¯2,where N is the number of samples, PSi is the estimated surface water extent value, P¯Si is the mean of estimated surface water extent, POi is the observed water level value, and P¯Oi is the mean of observed water level value. Additionally, the root mean square error (RMSE) was used to measure the magnitude of the average error in the extracted surface water area. The RMSE (equation (5)) is important for detecting undesirable large errors [52]:(5)Root Mean Square ErrorRMSE=∑i=1NPSi−POi2N1/2.Finally, the mean bias error (MBE) (equation (6)) was used to measure the mean absolute value of the error. If the predicted value is larger than the observed value, the mean bias error will be positive, and vice versa:(6)Mean Bias Error=∑i≡1NPSI−POIN.


## 3. Results and Discussion

### 3.1. Extraction of Water Area

Sentinel-1 images of the Arau and Kangar areas on 24 September 2017 were used as flood stage images. The methods described in Section2.7 were used to extract the surface water area, and an automated thresholding method was developed to differentiate water and nonwater areas. Figure 9 shows an example of the water area extracted from the Sentinel-1 images.Figure 9
Examples of water area detection in (a) lake and (b) river (bottom). The water area is shown in blue and the background image shown is the google earth image.
(a)(b)Figure7 shows the flood extent in September 2017 extracted from the Sentinel-1 satellite images.The result shows that the flood extent information is consistent with the classification results of flood susceptibility in previous studies, in which major floods were found to be located in the south of Perlis and paddy field areas [38, 53]. As shown in Figure 7, the extracted flood extent (blue color) covered 41% of the total combined area (150.95 km2) of the Kangar and Arau districts. The high flood susceptibility in the south of Perlis is due to the presence of alluvial deposits, paddy fields, low-humidity clay soil, and flat terrain. The flood event in 2017, due to Typhoon Doksuri, severely damaged rice crop. Most paddy field areas were submerged in water for over three days, preventing seedlings from growing. Based on the result, this study shows agreement with the other previous study which indicated that flood mostly occurs due to its terrain structure, where low-lying areas like paddy fields are located [50, 54]. In general, it is challenging to perform flood classification using radar sensors when the land pixels have similar characteristics to water pixels. The results of this study show that the cause of overestimation of flooded areas is shadow due to terrain characteristics or the incidence angle of satellite images; the backscatter value of the shadow which is similar to that of water causes the misclassification of shadow as floodwater [55].


### 3.2. Grid-Based Image Fusion

Subsequently, the method explained in Section2.7 was applied to the MNDWI images, since it is less complex than applying the fusion method to the green and SWIR bands, which are used to calculate the MNDWI. Figure 10 shows the extracted water area for grid 21 (Figure 5) based on two optical images acquired on 7 September 2017: one from Landsat 7 and one from MODIS Terra.Figure 10
Landsat 7 optical image (a), MODIS terra optical image (b), both acquired on 7 September 2017, and a fused image created using these two images (c). The water area is shown in blue.
(a)(b)(c)Figure11 shows the scatterplots between the MNDWI of Landsat 7 optical images and that of MODIS optical images, the MNDWI of Landsat 7 optical images and that of fused images, and the MNDWI of MODIS optical images and that of fused images. The correlation coefficients were 0.31, 0.56, and 0.66, respectively. The amount of spectral content is preserved.Figure 11
Scatterplots between (a) MNDWI of Landsat 7 optical images and that of MODIS optical image, (b) MNDWI of Landsat 7 optical image and that of fused image, and (c) MNDWI of MODIS optical image and that of the fused image.
(a)(b)(c)Figure12 shows the results of the fusion of a MODIS Terra optical image and a Sentinel-1 SAR image acquired on 12 September 2017 for grid 21. The results show that the fused image can preserve the information of the multispectral image while maintaining the structural and geometrical information of the SAR image. Although the Sentinel-1 image has been filtered to remove noise, a few pixels of noise still appear, and the boundaries of objects are more apparent in the fused image than in the Sentinel-1 or MODIS Terra images.Figure 12
Sentinel-1 (a) and MODIS Terra (b) images, both acquired on 12 September 2017, and a fused image created using these two images (c). The water area is shown in blue.
(a)(b)(c)Figure13 shows the correlation between the backscatter of Sentinel-1 and that of MODIS optical images, the backscatter of Sentinel-1 and that of fused images, and the MNDWI of MODIS optical images and that of fused images. The correlation coefficients were 0.06, 0.33, and 0.18, respectively. Compared with the fused images created using two optical images, the fused images created using optical and SAR images have a lower correlation coefficient owing to the differences in data types, data structures, spatial resolution, and geometric characteristics.Figure 13
Scatterplots between (a) MNDWI of Sentinel-1 backscatter images and that of MODIS optical images, (b) MNDWI of Sentinel-1 backscatter images and that of fused images, and (c) MNDWI of MODIS optical images and that of fused images.
(a)(b)(c)The classification of the surface water area using the integration of optical, SAR, and fused images from 4 September to 30 October 2017 for some affected grids is shown in Figure14. In general, the water area in all of the grids is consistent with the ideal flood variation modelling curve (Figure 1). In grid 5, the maximum surface water extent occurred on 1 October, with a water area of 6.89 km2, and abruptly declined after 11 days. The maximum surface water extent in grids 20, 21, 29, and 39 was observed on 28 September 2017, with water areas of 7.19, 8.07, 7.47, and 10.02 km2, respectively. The surface water area abruptly declined in the next period, 30 September 2017, and subsequently decreased slowly until the end of the observation period.Figure 14
Time series of the surface water extent extracted from multiple satellite images for some of the flood-affected grids. (a) Grid 5, (b) grid 20, (c) grid 21, (d) grid 29, and (e) grid 39.
(a)(b)(c)(d)(e)The fused image on 7 September 2017 for grid 5 showed a zero amount of surface water recorded. This might be due to the fused optical images seriously affected by cloud occlusion, causing loss of data during the cloud masking process. In such circumstances, alternatively a single image with lower contaminated data is a better option to represent data on that certain date. Discrepancies of each image can cause fluctuations in a certain period. Since this study aimed to use all the images regardless of image condition, standardization is required. For improvement, polynomial modelling could be utilized to achieve low bias and low variance time-series data.


### 3.3. Evaluation Results

Because of the scarcity of rain gauge data for the September 2017 flood event, ground observations of the water level were only available for two grids. The correlation coefficient (r), RMSE, and MBE were calculated to evaluate the extraction results. Figure 15 shows the temporal variations of the extracted surface water area and the ground-measured water level for grids 21 and 29.Figure 15
Time series of (a) the extracted surface water area in grid 21 and the water level measured at the Kg Repoh station, and (b) the extracted surface water area in grid 29 and the water level measured at the Kg Bakau station.
(a)(b)Based on the observed water level in grid 21, on 24 September, the water level was at a dangerous level (>3.5 m), with a recorded level of 4.91 m. Additionally, in grid 21, the peak of the ground-measured water level (24 September) occurred before the peak of the extracted water area (28 September). The RMSE between the extracted water area and the ground-observed water level is 1.55, and the MBE is −0.57 (underestimation). Similarly, for grid 29, the peak of the ground-measured water level also occurred before the peak of the extracted water area, and the values of RMSE and MBE for this grid are 2.56 and −2.15 (underestimation), respectively. Overall the two data did not have a good relationship, but there is an assumption that there is a delay of flooding event associated with the increase of river water level [56]. Therefore, the results are acceptable for further analysis.


## 3.1. Extraction of Water Area

Sentinel-1 images of the Arau and Kangar areas on 24 September 2017 were used as flood stage images. The methods described in Section2.7 were used to extract the surface water area, and an automated thresholding method was developed to differentiate water and nonwater areas. Figure 9 shows an example of the water area extracted from the Sentinel-1 images.Figure 9
Examples of water area detection in (a) lake and (b) river (bottom). The water area is shown in blue and the background image shown is the google earth image.
(a)(b)Figure7 shows the flood extent in September 2017 extracted from the Sentinel-1 satellite images.The result shows that the flood extent information is consistent with the classification results of flood susceptibility in previous studies, in which major floods were found to be located in the south of Perlis and paddy field areas [38, 53]. As shown in Figure 7, the extracted flood extent (blue color) covered 41% of the total combined area (150.95 km2) of the Kangar and Arau districts. The high flood susceptibility in the south of Perlis is due to the presence of alluvial deposits, paddy fields, low-humidity clay soil, and flat terrain. The flood event in 2017, due to Typhoon Doksuri, severely damaged rice crop. Most paddy field areas were submerged in water for over three days, preventing seedlings from growing. Based on the result, this study shows agreement with the other previous study which indicated that flood mostly occurs due to its terrain structure, where low-lying areas like paddy fields are located [50, 54]. In general, it is challenging to perform flood classification using radar sensors when the land pixels have similar characteristics to water pixels. The results of this study show that the cause of overestimation of flooded areas is shadow due to terrain characteristics or the incidence angle of satellite images; the backscatter value of the shadow which is similar to that of water causes the misclassification of shadow as floodwater [55].


## 3.2. Grid-Based Image Fusion

Subsequently, the method explained in Section2.7 was applied to the MNDWI images, since it is less complex than applying the fusion method to the green and SWIR bands, which are used to calculate the MNDWI. Figure 10 shows the extracted water area for grid 21 (Figure 5) based on two optical images acquired on 7 September 2017: one from Landsat 7 and one from MODIS Terra.Figure 10
Landsat 7 optical image (a), MODIS terra optical image (b), both acquired on 7 September 2017, and a fused image created using these two images (c). The water area is shown in blue.
(a)(b)(c)Figure11 shows the scatterplots between the MNDWI of Landsat 7 optical images and that of MODIS optical images, the MNDWI of Landsat 7 optical images and that of fused images, and the MNDWI of MODIS optical images and that of fused images. The correlation coefficients were 0.31, 0.56, and 0.66, respectively. The amount of spectral content is preserved.Figure 11
Scatterplots between (a) MNDWI of Landsat 7 optical images and that of MODIS optical image, (b) MNDWI of Landsat 7 optical image and that of fused image, and (c) MNDWI of MODIS optical image and that of the fused image.
(a)(b)(c)Figure12 shows the results of the fusion of a MODIS Terra optical image and a Sentinel-1 SAR image acquired on 12 September 2017 for grid 21. The results show that the fused image can preserve the information of the multispectral image while maintaining the structural and geometrical information of the SAR image. Although the Sentinel-1 image has been filtered to remove noise, a few pixels of noise still appear, and the boundaries of objects are more apparent in the fused image than in the Sentinel-1 or MODIS Terra images.Figure 12
Sentinel-1 (a) and MODIS Terra (b) images, both acquired on 12 September 2017, and a fused image created using these two images (c). The water area is shown in blue.
(a)(b)(c)Figure13 shows the correlation between the backscatter of Sentinel-1 and that of MODIS optical images, the backscatter of Sentinel-1 and that of fused images, and the MNDWI of MODIS optical images and that of fused images. The correlation coefficients were 0.06, 0.33, and 0.18, respectively. Compared with the fused images created using two optical images, the fused images created using optical and SAR images have a lower correlation coefficient owing to the differences in data types, data structures, spatial resolution, and geometric characteristics.Figure 13
Scatterplots between (a) MNDWI of Sentinel-1 backscatter images and that of MODIS optical images, (b) MNDWI of Sentinel-1 backscatter images and that of fused images, and (c) MNDWI of MODIS optical images and that of fused images.
(a)(b)(c)The classification of the surface water area using the integration of optical, SAR, and fused images from 4 September to 30 October 2017 for some affected grids is shown in Figure14. In general, the water area in all of the grids is consistent with the ideal flood variation modelling curve (Figure 1). In grid 5, the maximum surface water extent occurred on 1 October, with a water area of 6.89 km2, and abruptly declined after 11 days. The maximum surface water extent in grids 20, 21, 29, and 39 was observed on 28 September 2017, with water areas of 7.19, 8.07, 7.47, and 10.02 km2, respectively. The surface water area abruptly declined in the next period, 30 September 2017, and subsequently decreased slowly until the end of the observation period.Figure 14
Time series of the surface water extent extracted from multiple satellite images for some of the flood-affected grids. (a) Grid 5, (b) grid 20, (c) grid 21, (d) grid 29, and (e) grid 39.
(a)(b)(c)(d)(e)The fused image on 7 September 2017 for grid 5 showed a zero amount of surface water recorded. This might be due to the fused optical images seriously affected by cloud occlusion, causing loss of data during the cloud masking process. In such circumstances, alternatively a single image with lower contaminated data is a better option to represent data on that certain date. Discrepancies of each image can cause fluctuations in a certain period. Since this study aimed to use all the images regardless of image condition, standardization is required. For improvement, polynomial modelling could be utilized to achieve low bias and low variance time-series data.


## 3.3. Evaluation Results

Because of the scarcity of rain gauge data for the September 2017 flood event, ground observations of the water level were only available for two grids. The correlation coefficient (r), RMSE, and MBE were calculated to evaluate the extraction results. Figure 15 shows the temporal variations of the extracted surface water area and the ground-measured water level for grids 21 and 29.Figure 15
Time series of (a) the extracted surface water area in grid 21 and the water level measured at the Kg Repoh station, and (b) the extracted surface water area in grid 29 and the water level measured at the Kg Bakau station.
(a)(b)Based on the observed water level in grid 21, on 24 September, the water level was at a dangerous level (>3.5 m), with a recorded level of 4.91 m. Additionally, in grid 21, the peak of the ground-measured water level (24 September) occurred before the peak of the extracted water area (28 September). The RMSE between the extracted water area and the ground-observed water level is 1.55, and the MBE is −0.57 (underestimation). Similarly, for grid 29, the peak of the ground-measured water level also occurred before the peak of the extracted water area, and the values of RMSE and MBE for this grid are 2.56 and −2.15 (underestimation), respectively. Overall the two data did not have a good relationship, but there is an assumption that there is a delay of flooding event associated with the increase of river water level [56]. Therefore, the results are acceptable for further analysis.


## 4. Conclusions

In mid-September 2017, Typhoon Doksuri caused heavy rain and high winds in northern Peninsular Malaysia, which resulted in a flood event in Kangar and Arau districts. In this study, Landsat 7, Landsat 8, MODIS, ALOS-2, and Sentinel-1 images were integrated to estimate changes in the floodwater extent during the September 2017 event. A simple and operational flood monitoring methodology is presented. The integration of data from different sensors enables a higher observation frequency, which is favorable for flood response and increases the reliability of the estimation of flood duration. The degree of economic loss and structural damage is highly dependent on the flood duration.This study found that the proposed automatic flood-detection method is efficient for mapping the water surface area during flood events. The water area extraction results are acceptable for lakes that are large and homogenous. However, the classification showed a lower performance for rivers, especially when the rivers are narrow, close to urban structures, or follow a meandering path. Therefore, mapping flood occurrence in urban areas using SAR images remains challenging. It was impossible to extract the flood area in urban areas using Sentinel-1 data due to the high radar backscatter. Furthermore, the results showed that incorporating a flood variation spatial and temporal model gives a useful overview of the flood inundation process, especially in data-scarce areas. By using a grid system, unflooded areas can be ignored, thus reducing the processing time.In the future, we plan to introduce a more reliable method to detect floods in urban areas and enhance the flood-detection ability in areas with small rivers. The implementation of a well‐defined method that considers spatial and temporal data about flood depth using digital elevation models (DEMs) and topographic information should be considered in future research.



---

*Source: 6164161-2021-11-01.xml*