MODIS-Based Investigation of Flood Areas in Southern Cambodia from 2002–2013

: In Cambodia and the Vietnamese Mekong Delta, ﬂoods commonly occur during the rainy season, and a better understanding of their spatio-temporal distribution is important for both disaster prevention and the improvement of agricultural production. This study investigated spatio-temporal ﬂood inundation and land cover change from 2002 to 2013 in the southern part of Cambodia using Terra satellite on-board Moderate Resolution Imaging Spectroradiometer (MODIS) images. The algorithm for ﬂood inundation detection, WFFI (Wavelet-based Filter for detecting spatio-temporal changes in Flood Inundation) was used, and the parameters were modiﬁed to ﬁt the present study. The estimated inundation areas were validated using eight Landsat images. In a comparison between the original and modiﬁed WFFIs, the modiﬁed WFFI (70–96%) exhibited better accuracy than the original WFFI (30–70%). Overall, the temporal change in the ﬂood inundation area presented a decreasing trend, and a link to the in-situ observed water level showed a decreasing trend during the rainy season. Furthermore, the estimated ﬂood inundation exhibited a signiﬁcant delay since 2008. Based on the yearly land cover MODIS product, the permanent water body and wetland areas decreased, whereas the cropland areas increased. This was as a result of increased agricultural productivity. However, water shortage was the major obstacle to increasing agricultural productivity, and it also had a negative impact on aquatic ecology, such as ﬁsh spawning grounds.


Introduction
In Cambodia, as in the most of the countries along the Mekong River, which is the seventh longest river in Asia, floods commonly occur during the rainy season, and they have both beneficial and

MODIS Data
This study used the following two MODIS products: (1) The MODIS/Terra Surface Reflectance 8-Day L3 Global 500 m SIN Grid V005 (MOD09A1) and (2) The Land Cover Type Yearly L3 Global 500 m SIN Grid (MCD12Q1). The MOD09A1 products yield the best surface spectral-reflectance data over an 8-day period, with the least effects from atmospheric water vapor (EOS, 2014). This study involved an analysis of MOD09 8-day composite data acquired from 2002 (day of year (DOY) 1) to 2013 (DOY 365). The spatial distribution of the start dates varied year by year [16]. Although the level of cloud cover in the wet season means that it is difficult to obtain cloud-free images every day, 8day composite products can provide an image that is almost free of clouds and that is atmospherically corrected. Although there is a limitation that spatial resolution becomes coarser than Landsat, the time series MODIS data could be used to determine these dates to the nearest week and to map the spatial extent of a flood.
The MCD12Q1 product is supplied annually, and it covers land areas with a 500 m spatial resolution. There are five classification schemes in this product that display land cover properties from observations, which include one year's input of Terra and Aqua MODIS data. In this study, we used the primary land cover scheme (17 land cover classes) defined by the International Geosphere Biosphere Program (IGBP), which includes 11 natural vegetation classes, three developed and mosaicked land classes, and three non-vegetated land classes.
The default projection of original MODIS data is a MODIS sinusoidal tiling system. After mosaicking six tiled images to cover the study area, all the data were transferred into the Geographic Tagged Image File Format (GeoTIFF) with geographic projection (Latitude/Longitude, WGS84) using the MODIS reprojection tool [19].

MODIS Data
This study used the following two MODIS products: (1) The MODIS/Terra Surface Reflectance 8-Day L3 Global 500 m SIN Grid V005 (MOD09A1) and (2) The Land Cover Type Yearly L3 Global 500 m SIN Grid (MCD12Q1). The MOD09A1 products yield the best surface spectral-reflectance data over an 8-day period, with the least effects from atmospheric water vapor (EOS, 2014). This study involved an analysis of MOD09 8-day composite data acquired from 2002 (day of year (DOY) 1) to 2013 (DOY 365). The spatial distribution of the start dates varied year by year [16]. Although the level of cloud cover in the wet season means that it is difficult to obtain cloud-free images every day, 8-day composite products can provide an image that is almost free of clouds and that is atmospherically corrected. Although there is a limitation that spatial resolution becomes coarser than Landsat, the time series MODIS data could be used to determine these dates to the nearest week and to map the spatial extent of a flood.
The MCD12Q1 product is supplied annually, and it covers land areas with a 500 m spatial resolution. There are five classification schemes in this product that display land cover properties from observations, which include one year's input of Terra and Aqua MODIS data. In this study, we used the primary land cover scheme (17 land cover classes) defined by the International Geosphere Biosphere Program (IGBP), which includes 11 natural vegetation classes, three developed and mosaicked land classes, and three non-vegetated land classes.
The default projection of original MODIS data is a MODIS sinusoidal tiling system. After mosaicking six tiled images to cover the study area, all the data were transferred into the Geographic Tagged Image File Format (GeoTIFF) with geographic projection (Latitude/Longitude, WGS84) using the MODIS reprojection tool [19].

Landsat Data
Eight Landsat images (see Table 1) were downloaded from the Global Land Cover Facility GLCF [20]. In the study, five images were obtained from the Landsat 7 satellite on-board Enhanced Thematic Mapper plus (ETM+), and another three images were obtained by the Landsat 8 satellite on-board Operational Land Imager (OLI) within only one path/row (126/52). Atmospheric correction was performed using a Fast Line-of-sight Atmospheric Analysis of a Spectral Hypercube (FLAASH) module on ENVI software ver. 5.1 (Excelis Visual Information Solutions, Colorado, USA).

Daily Water Level
Daily water level data between 2002 and 2013 from five stations were provided by the Mekong River Commission (MRC). The stations are located at Bassac Chaktomuk (11 •

Detecting Water Surface Using MODIS Data with WFFI
The WFFI algorithm used in this study is a modified approach, as shown in Figure 2, based on the techniques originally developed by Sakamoto et al. [16] to detect spatio-temporal flood distributions in Cambodia and Vietnam. The WFFI is based on smoothed indices of MODIS images, including the Enhanced Vegetation Index (EVI), Land Surface Water Index (LSWI) and their difference (DVEL): where Blue, Red, NIR and SWIR are the surface reflectance value in blue (band 3; 459-479 nm), red (band 1; 621-670 nm), near infrared (NIR) (band 2; 841−875 nm) and short-wave infrared (SWIR) (band 6; 1628-1652 nm), respectively. In the original WFFI algorithm of MODIS sensor, a wavelet-based filter was used to interpolate any missing information. The production of the inundation maps involved a decision tree to classify each pixel into one of the following categories: flood, mixture, non-flood and a long-term water body. In the present study, the threshold values of the decision tree were modified based on the smoothed EVI, LSWI and DVEL values ( Figure 3). The modified WFFI was evaluated by comparing it with the original WFFI.

Validation with Landsat Data
Several studies have been conducted on detecting water bodies using Landsat data [21][22][23][24][25]. To detect water bodies from Landsat image, several band ratio indices have been developed. In this study, the water surface pixels were identified using the Modified Normalized Difference Water Index (MNDWI) [25] as a simple but more effective band ratio method than others. First, the nonanalysis areas including cloud cover and cloud shadow were determined from the blue reflectance value. The cloud-cover areas were derived from the majority analysis (window size: 3 × 3 pixels) for pixels whose blue reflectance (band 1 in Landsat 7 ETM+ or band 2 in Landsat 8 OLI) was greater than 0.2. The cloud-shadow areas were then masked by the cloud-cover areas, including 10-pixel buffer zones, where they were shifted in the opposite direction from the sun orientation. Next, water surface pixels were identified where the MNDWI was greater than or equal to 0. The equation for the MNDWI used in this study is as follows [25], where Green and MIR are the surface reflectance values in green (band 2; 530-600 nm for Landsat 5, 530-610 nm for Landsat 7) and middle-infrared (MIR) (band 5; 1550-1750 nm) bands

Validation with Landsat Data
Several studies have been conducted on detecting water bodies using Landsat data [21][22][23][24][25]. To detect water bodies from Landsat image, several band ratio indices have been developed. In this study, the water surface pixels were identified using the Modified Normalized Difference Water Index (MNDWI) [25] as a simple but more effective band ratio method than others. First, the non-analysis areas including cloud cover and cloud shadow were determined from the blue reflectance value. The cloud-cover areas were derived from the majority analysis (window size: 3 × 3 pixels) for pixels whose blue reflectance (band 1 in Landsat 7 ETM+ or band 2 in Landsat 8 OLI) was greater than 0.2. The cloud-shadow areas were then masked by the cloud-cover areas, including 10-pixel buffer zones, where they were shifted in the opposite direction from the sun orientation. Next, water surface pixels were identified where the MNDWI was greater than or equal to 0. The equation for the MNDWI used in this study is as follows [25], where Green and MIR are the surface reflectance values in green (band 2; 530-600 nm for Landsat 5, 530-610 nm for Landsat 7) and middle-infrared (MIR) (band 5; 1550-1750 nm) bands Lastly, 'not water-surface' pixels were identified where MNDWI was less than 0. The threshold value of 0 was adopted from previous literature by Xue [25]. The land-surface map was resampled to a grid with a resolution of 25 m using the nearest neighbor method, and then the areas of these different land-surface types were aggregated within each grid with a resolution of 500 m to compare the MODIS-derived results.

Start and End Dates of Flooding
The continuous flood period was determined by the following procedures: (1) The continuous inundation period was recorded by water-related pixels (both flood and mixture pixels) and computed with water-related or non-flood pixels over 12 years; (2) if the periods of continuous inundation were identified twice or more in a year, the longer period was taken as the period of annual flood inundation; and (3) the duration of annual flood inundation was detected using the start and end dates of annual flood inundation.

Water Level Trends
The annual changes in water levels in the rainy season (June-November) were examined. According to the slope factor of the linear regression function, the trend in water level time-series data could be incorporated into the calculation; the slope of the linear regression model fit resembles a straight line, with the water level data plotted on the y-axis as the response variable and the value of time (t = 1,2,3,...) plotted on the x-axis as an explanatory variable. Therefore, if the slope factor ≥0, the water has increased or not changed, and if the slope factor <0, the water has decreased. Figure 4 shows the surface water area derived by the modified WFFI and original WFFI on 1 January 2002. In the inundation map from the original WFFI, approximately 78% of the area is covered by water, suggesting overestimation. In a comparison between the original and modified WFFIs (see Table 1), the modified WFFI (70-96%) exhibited better accuracy than the original WFFI (30-70%). This might be considered to be differences in the sensitivity of the bands used for WFFI due to the change in MODIS product. The original WFFI was developed in version 4 of the MODIS product, while version 5 was used in this study, so it was necessary to adjust the threshold value for extracting water-related pixels in WFFI. This was in agreement with a previous study by Islam et al. [17] with MODIS product in version 5. They adopted and modified WFFI in Bangladesh to produce a flood inundation map and obtained similar accuracy (R 2 = 0.96) in validation using the RADARSAT derived surface water area.

Detecting Flood Areas
Spatial resolutions of 500 m could be another limitation in original WFFI, which are not always suited to small river catchments [26]. The threshold value in NDWI varies depending on the proportions of subpixel water or non-water components [27]. In this study, inundation maps derived from Landsat images, which was resampled in a resolution of 500 m to conform the MODIS-derived results were used for the validation. Better spatial resolution data, such as Landsat (30 m) or Sentinel-2 (10 m), could be reduced by the mixture issues, and provide more appropriate spatial information to water mapping applications [28][29][30]. However, their temporal frequency (revisit every 16 days in Landsat and 5 days in Sentinel-2) are not suited for monitoring spatial dynamics in flood events. Therefore, MODIS data is still used for mapping flood events due to the near-global spatial (250-1000 m) and high temporal (daily) characteristics.
Due to the limitation of field data, in this study, the validation was made by Landsat-derived results using MNDWI as the most suitable indices for detecting surface water [31]. To date, several indices have been developed to determine water-related pixels in optimal satellite data, such as normalized difference water index (NDWI) [32]. Future work needs to examine this by comparing other indices with field surveys data.

Temporal Characteristics of Annual Floods and Water Levels
Using the estimated inundation maps from the modified WFFI, the spatial distributions of the start and end dates of annual floods were computed ( Figure 5). The spatial distribution of the start and end flood dates varied year by year. Furthermore, seasonal changes in the extent of the estimated flood area were calculated for each year ( Figure 6). The estimated start date of the flood indicated that flood occurrence was late from 2010 to 2013. Moreover, the estimated flood inundation area showed 27% decrease over the last few years.
Water infrastructure developments, such as flood control, dam hydropower, and irrigation, were the main causes of reduced inundation area [33]. At the same time, the gradient of the water level at five stations in the rainy season exhibited a negative trend. Figure 7 shows the water levels at the Bassac Chamtomuk station in the rainy season. The water level showed decreased trends over the past 12 years (2002-2013). There are many water infrastructure plans in the upstream river, especially hydropower dams (Osborne 2008; Pearse-Smith 2012). Similarly, Arias et al. [34] found that water infrastructure development was expected to shrink the magnitude of floods by raising the water level during the dry season and reducing it during the rainy season. These developments resulted from the amount of water used and the demand increase due to economic development, especially hydropower dam construction upstream. Nik Hashim [35] noted that dam construction

Temporal Characteristics of Annual Floods and Water Levels
Using the estimated inundation maps from the modified WFFI, the spatial distributions of the start and end dates of annual floods were computed ( Figure 5). The spatial distribution of the start and end flood dates varied year by year. Furthermore, seasonal changes in the extent of the estimated flood area were calculated for each year ( Figure 6). The estimated start date of the flood indicated that flood occurrence was late from 2010 to 2013. Moreover, the estimated flood inundation area showed 27% decrease over the last few years.
Water infrastructure developments, such as flood control, dam hydropower, and irrigation, were the main causes of reduced inundation area [33]. At the same time, the gradient of the water level at five stations in the rainy season exhibited a negative trend. Figure 7 shows the water levels at the Bassac Chamtomuk station in the rainy season. The water level showed decreased trends over the past 12 years (2002-2013). There are many water infrastructure plans in the upstream river, especially hydropower dams (Osborne 2008; Pearse-Smith 2012). Similarly, Arias et al. [34] found that water infrastructure development was expected to shrink the magnitude of floods by raising the water level during the dry season and reducing it during the rainy season. These developments resulted from the amount of water used and the demand increase due to economic development, especially hydropower dam construction upstream. Nik Hashim [35] noted that dam construction and water demands led to a 4% decrease in the maximum water level in the rainy season and a catchment area decline of 5% in the Tonle Sap region. This negative impact would lead to decreased fish production and thereby threats to food security, restricted agriculture productivity from reduced water supply for agriculture land, water pollution that is harmful to the livelihoods of local communities, and ecological changes, especially in aquatic ecosystems, and environmental changes that will also affect economic development [36]. and water demands led to a 4% decrease in the maximum water level in the rainy season and a catchment area decline of 5% in the Tonle Sap region. This negative impact would lead to decreased fish production and thereby threats to food security, restricted agriculture productivity from reduced water supply for agriculture land, water pollution that is harmful to the livelihoods of local communities, and ecological changes, especially in aquatic ecosystems, and environmental changes that will also affect economic development [36].    catchment area decline of 5% in the Tonle Sap region. This negative impact would lead to decreased fish production and thereby threats to food security, restricted agriculture productivity from reduced water supply for agriculture land, water pollution that is harmful to the livelihoods of local communities, and ecological changes, especially in aquatic ecosystems, and environmental changes that will also affect economic development [36].    water supply for agriculture land, water pollution that is harmful to the livelihoods of local communities, and ecological changes, especially in aquatic ecosystems, and environmental changes that will also affect economic development [36].

Land Cover and Water Body Changes between 2002 and 2012
Based on the MCD12Q1 product, Table 2   The expansion of agricultural land is one of the reasons for decreasing water area. Senevirathne et al. [37] also found that water bodies significantly decreased due to increase croplands and built up areas in Cambodia. However, the increasing trend of cropland area may be closely linked to economic growth in Cambodia. Agriculture is one of the main sectors contributing to economic growth, and the expansion of agricultural land would lead to an increase in productivity, which would contribute to national food security and the alleviation of poverty [38]. The consistent trends in the estimated flood area and the decrease in the permanent water body land cover might be caused by hydropower dams, which take water for electricity generation, and other water uses fueling economic growth, such as industry and irrigation for agriculture [39,40].

Conclusions
This study investigated spatio-temporal flood inundation and land cover change in the southern part of Cambodia using the Terra MODIS products. The WFFI algorithm for flood inundation detection was used in this research, and the parameters were modified to fit with the current study. The flood inundation maps from MODIS 500 m resolution were evaluated by comparing them with flood inundation maps based on resampled Landsat data. In a comparison between the original and modified WFFIs, the modified WFFI (70-96%) showed better accuracy than the original WFFI (30-70%). This demonstrates that the modified WFFI methodology is very useful for detecting flood inundation and that the flood inundation maps from MODIS images can be used to characterize flooding, which is very useful for the integration of water resources, flood management and the maintenance of wetland ecosystems. Overall, the WFFI products are in good agreement with the water surface area derived from Landsat images and hydrological data in terms of temporal change and spatial distribution. Furthermore, the estimated flood inundation data exhibited a significant delay since 2008 caused by climate change. Based on the yearly land cover data, the permanent water body and wetland areas decreased, while the cropland areas increased. This was a result of increased agricultural productivity.
However, water shortage is one of the major obstacles to increasing agricultural productivity, and it also has a negative impact on aquatic ecology such as fish spawning grounds.
Funding: This work was partially supported by the Hiroshima Earth Environment Information Center for the paper submission and the Japanese Grant Aid for Human Resource Development Scholarship (JDS) and the Japan International Cooperation Agency (JICA) for financial support to the first author (N.V.) for the field survey in Cambodia.