Automated Extraction of Inundated Areas from Multi-Temporal Dual-Polarization RADARSAT-2 Images of the 2011 Central Thailand Flood

This study examines a novel extraction method for SAR imagery data of widespread flooding, particularly in the Chao Phraya river basin of central Thailand, where flooding occurs almost every year. Because the 2011 flood was among the largest events and of a long duration, a large number of satellites observed it, and imagery data are available. At that time, RADARSAT-2 data were mainly used to extract the affected areas by the Thai government, whereas ThaiChote-1 imagery data were also used as optical supporting data. In this study, the same data were also employed in a somewhat different and more detailed manner. Multi-temporal dual-polarized RADARSAT-2 images were used to classify water areas using a clustering-based thresholding technique, neighboring valley-emphasis, to establish an automated extraction system. The novel technique has been proposed to improve classification speed and efficiency. This technique selects specific water references throughout the study area to estimate local threshold values and then averages them by an area weight to obtain the threshold value for the entire area. The extracted results were validated using high-resolution optical images from the GeoEye-1 and ThaiChote-1 satellites and water elevation data from gaging stations.


Introduction
Floods occur almost every year in Thailand and cause unfavorable situations.The worst flooding in the last five decades occurred in 2011 [1].The World Bank has estimated the damage and the losses due to this flooding at approximately THB 1.43 trillion (USD 46.5 billion) in total, while the recovery and reconstruction needs were estimated to be THB 1.5 trillion (USD 50 billion) over the five-year period [2].This flood event spread throughout the northern, northeastern, and central provinces of the country.The flooding caused heavy economic impacts by disturbing industrial production in the affected areas and the supply chains of industries worldwide [1][2][3].In this study, satellite imagery data, which can effectively extract information in large-scale disasters, were introduced to evaluate the extent of the flood.Among other types of sensors, synthetic aperture radar (SAR) sensors can operate during the day and night and under all weather conditions [4].RADARSATs, which are Canadian SAR satellites with C-band radars operating at a wavelength of 5.6 cm, have been mainly used for flood monitoring in Thailand since 2000 (RADASAT-1, RS1) and 2008 (RADARSAT-2, RS2) [5,6].ThaiChote-1 (TH1), the first satellite of Thailand, has been used to provide optical support under clear sky conditions since 2004.
Earth terrain surfaces are considered to be rough at radar wavelengths and exhibit diffuse scattering with moderate backscatter.In contrast, water surfaces are generally smooth and can be regarded as specular reflectors that yield small backscatter at radar wavelengths.As a consequence, the surrounding terrain corresponds to brighter intensities in SAR images, whereas water is regarded as low intensity area.Therefore, SAR images are considered to be very effective and have been extensively used for water and flood mapping [7,8].
Single co-polarized HH (horizontal transmit and horizontal receive) SAR images are the most common and are useful for determining water and flood areas.Especially for C-band radars, this is the preferred polarization for mapping flooded vegetation because it maximizes canopy penetration and enhances the contrast between forests and flooded vegetation [9][10][11][12].Although it has been used in many cases, but with respect to water surfaces, obstacle cover, floating objects, and wind ripples, affect SAR backscatter, preventing C-band from returning good results.Dual and full polarizations have been employed to enhance capability [11][12][13].Dual-polarizations can potentially be used to detect and map vegetation water content (VWC) in forested areas and more reliably distinguish open water surfaces affected by wind.When an SAR image is acquired in lighter winds or under smooth water surface conditions, the HH co-polarization has been shown to be the most suitable for mapping surface waters.However, when wind or surface water roughness is present, the single cross-polarized HV (horizontal transmit and vertical receive) often yields better results for water extraction [14][15][16].Unfortunately, the Thai government only uses HH polarization in most cases to monitor flood events.Therefore, the intent of this study is to improve the effectiveness of mapping surface waters by combining the depolarization information in HV with HH as the total backscatter.
Deriving the extent of inundation from a single SAR image has been carried out using several methods, e.g., pixel-based classification [17][18][19][20][21][22][23][24], segment-based classification with region growth [25][26][27][28], and mixing between the two methods [19,29,30].The most common and efficient way is thresholding, which is a pixel-based operation.In computer vision and image processing, thresholding was introduced to reduce a gray level image to a binary image, foreground and background.The algorithm assumes that the image histogram is distributed in two classes or has a bimodal distribution.Flooded areas or foregrounds are separated from the background by a constant threshold value.Manual threshold-value selection may be faced with a problem; it is hard to judge the most suitable value objectively.Automatic thresholding methods have been introduced to overcome this issue and to improve the classification speed or efficiency.Several techniques have been proposed to determine threshold values for SAR images [17][18][19][20][21][22].The optimum global threshold value can be obtained from the minimum within-class variance [17,18] or the minimum error thresholding [19][20][21], whereas some techniques look for local threshold values using auxiliary data, e.g., elevation and slope [22][23][24].In this study, we selected a threshold value by a modification of the Otsu's method.The Otsu's method is a clustering based thresholding, one of the most referenced methods [31].This technique establishes an optimum threshold by minimizing the weighted sum of within-class variances of the foreground and background pixels [17,18,31].This technique is robust for noisy images with Gaussian noise and is the best for presenting the inter-region contrast of SAR images [17,18].
To detect floods over large-scale areas such as the Chao Phraya river basin, identifying the global threshold value from the histogram of the entire image is almost impossible because this histogram has a unimodal distribution.The global threshold value is estimated in an indirect way as the arithmetic average of local threshold values for small areas.First, the image is divided into small portions, and then only those portions that have a bimodal distribution histogram are selected as representative.This technique can be carried out by the fully automated division of an image by a certain pattern, the bi-level quad tree [19][20][21].However, the advantage of systematic hierarchical image division is location independent.It is time consuming, particularly for processing multi-temporal images in the same place, and tiled pieces at different times may not be in the same place.A permanent representative area is proposed in this study to decrease processing time and to monitor local water bodies in a time series.This is much more suitable for Thailand because floods occur almost every year.
Because the Chao Phraya river basin is located in a large flat plain, inundation always remains for an extended period of time.The damage caused by flooding depends not only on the water depth but also on the flood duration.In 2011, the inundation depths ranged between 0.6 m and 4.9 m, with an average of 2.2 m, and the mode (most frequent value) was 2.5 m.Fifty-seven (57) days was the average inundation duration, which ranged from 3 to 120 days and had a mode of 60 days [3].The flooded areas could be captured remotely by a single satellite image, and the duration could be obtained by monitoring the flooded areas in a sequence of time.However, water depth requires auxiliary data for processing, e.g., a digital elevation model (DEM).In this study, estimating water depth was set aside for future work.Thus, only the inundation area and duration are studied in this article.Working with multi-temporal SAR images, it is very difficult to obtain images in the exact same conditions, e.g., satellite position, look-side (right/left), and incidence angle.Several methods must be introduced to normalize multi-temporal images prior to analysis [32,33].In contrast, the proposed method, which obtains the threshold value from water references, is independent of the acquisition condition, and there is no need to normalize the images [19][20][21].

Study Area and Imagery Data
This paper focuses on the Chao Phraya river basin in the central part of Thailand, which has an area of approximately 7400 km 2 .The basin was assigned as a regularly inundated area by the Thai government.The city of Ayutthaya, which is approximately 16.5 km in width and 21.0 km in length and includes the Ayutthaya Historical Park and the Rojana and Hi-Tech Industrial Parks, was selected as a validation area for extracting water bodies during the 2011 Thailand flood, as shown in Figure 1.The imagery data acquired by RADARSAT-2 (RS2) during that event have several beam types, i.e., Wide1 (W1), Wide2 (W2), Wide3 (W3), SCNA (W1 + W2), and SCNB (W2 + S5 + S6).The radar frequency, resolution, and incident angle varied depending on the beam type.Most of the images were taken in the HH and HV dual-polarization mode, and 30 images using that polarization were selected.The images were observed either from the ascending (ASC) or descending (DES) path.A summary of the image properties is presented in Table 1, and their color composite maps are shown in Figure 2. Because the swath width and revisit time of the satellite are limited, only some of the images could capture the whole study area.Two optical images were also used for validation.A GeoEye-1 (GE1) image was prioritized because it had a higher resolution and was taken on 22 November 2011, which was the nearest in time to one of the RS2 images.The pan-sharpened GE1 image had four multispectral bands with a 1.0-m resolution.ThaiChote-1 (TH1) images are also important because they are used as common data.A pan-sharped TH1 image has a multispectral band with a 2.5-m resolution.A TH1 image was taken on 25 November 2011.Another TH1 image was taken prior the flood on 12 December 2009.That image was used to assess the environment in the dry season but was not used in the analysis process.All the optical images will be shown in the accuracy assessment section.

Methodology and Results
The methodology used in this study is comprised of two parts.The main part is an automated system that extracts water bodies and produces flood duration maps.The second part is an accuracy assessment, which is a process to measure the efficiency of the main part using optical images as truth data.That part did not need to be implemented in the automated system.
The main part begins with preparing RS2 images by the radiometric calibration using the lookup tables provided in these products, the Refined Lee speckle filtering [34,35] with a window size of 5 × 5 pixels, and the terrain correction using an SRTM 90-m DEM.Radiometric calibration provides images in which pixels can be directly related to the radar backscatter of the scene by applying the sigma naught, beta naught, and gamma naught lookup tables provided in the product [32].The pre-processing was performed using the Sentinel Application Platform (SNAP) software program.Next, RS2 images were clipped by the study area and then temporarily cut again into specific segments throughout the study area.Those particular segments were defined as water references, as mentioned in the next section.The automatic thresholding was applied to each water reference one at a time to estimate a local threshold value.In that step, water references with unimodal distribution histograms were rejected.The local threshold values were weight-averaged by area to estimate the global threshold value.Finally, the global threshold value was applied to the whole image to extract the water bodies.The flood duration map is an auxiliary product for flood management and is calculated from the cumulative days of water bodies.All these processes were developed using a Python script and applicable modules, e.g., the Geospatial Data Abstraction Library (GDAL) for reading and writing image data, NumPy for numerical calculation, and Matplotlib for plotting histograms.NumPy is the fundamental package for scientific computing with powerful N-dimensional array objects and linear algebra.Matplotlib is 2D and 3D plotting library which produces quality figures under interactive environments.
The purpose of the accuracy assessment was to measure the efficiency and accuracy of the final results by comparing the water areas extracted from the SAR images (RS2) with those from the optical images (GE1 and TH1).The backscattering coefficient was used for the SAR images, and the Normalized Difference Water Index (NDWI) was used for the optical images.Although the RS2 images were processed using automatic thresholding, the GE1 image was processed using manual thresholding.The RS2 and TH1 results were compared with the GE1 image to determine the most accurate thresholding.The work flow diagram used in this study is shown in Figure 3.
results by comparing the water areas extracted from the SAR images (RS2) with those from the optical images (GE1 and TH1).The backscattering coefficient was used for the SAR images, and the Normalized Difference Water Index (NDWI) was used for the optical images.Although the RS2 images were processed using automatic thresholding, the GE1 image was processed using manual thresholding.The RS2 and TH1 results were compared with the GE1 image to determine the most accurate thresholding.The work flow diagram used in this study is shown in Figure 3.

Water References and Histogram Analysis
Applying automatic clustering-based thresholding to a large area may return unsatisfying results because its histogram is not bimodally distributed [31].This problem occurs when the ratios of the water and non-water areas are very different.A novel technique was implemented in this study to overcome this problem by applying automatic thresholding to specific smaller areas located throughout the study area.Those small areas were defined as water references, which can be used for any satellite images acquired on a different date and time.
All the water references were selected using the following criteria: having an area larger than 320,000 m 2 (512 pixels for a resolution of 25.0 m, and 2048 pixels for a resolution of 12.5 m), containing water throughout the year, not facing a flood situation, located on flat ground as much as possible, and having a water and non-water cover ratio of nearly 1:1.An irregular shape is allowed.In this study, an elliptical shape was preferred because it was easy to maintain the ratio of the water and non-water proportions.The 87 water references were selected from different types of water bodies, natural and man-made.Their locations are shown in Figure 1 (right).
Because not all of the RS2 images covered the entire study area and included all the water references, only those water references whose entire areas were covered by each RS2 image were taken into account.The number of water references and their covering areas for each image are

Water References and Histogram Analysis
Applying automatic clustering-based thresholding to a large area may return unsatisfying results because its histogram is not bimodally distributed [31].This problem occurs when the ratios of the water and non-water areas are very different.A novel technique was implemented in this study to overcome this problem by applying automatic thresholding to specific smaller areas located throughout the study area.Those small areas were defined as water references, which can be used for any satellite images acquired on a different date and time.
All the water references were selected using the following criteria: having an area larger than 320,000 m 2 (512 pixels for a resolution of 25.0 m, and 2048 pixels for a resolution of 12.5 m), containing water throughout the year, not facing a flood situation, located on flat ground as much as possible, and having a water and non-water cover ratio of nearly 1:1.An irregular shape is allowed.In this study, an elliptical shape was preferred because it was easy to maintain the ratio of the water and non-water proportions.The 87 water references were selected from different types of water bodies, natural and man-made.Their locations are shown in Figure 1 (right).
Because not all of the RS2 images covered the entire study area and included all the water references, only those water references whose entire areas were covered by each RS2 image were taken into account.The number of water references and their covering areas for each image are shown in Table 2.It was impossible to show all of them in this article, and thus only 10, taken on 28 February 2011, are shown in Figure 4 (see supplementary materials for the full version).Smooth water surfaces are shown in black for the HH, HV, and HH + HV polarizations and in deep blue for the color composites of HH, HV, and HH/HV."HH + HV" denotes the sum of the backscattering coefficients as the total backscattering, and "HH/HV" denotes the ratio of the backscattering coefficients as the relative backscattering.Both values were calculated on a linear scale but are presented on a logarithmic scale (dB).
The histogram plots shown in Figure 4 were obtained from the RS2 images and were clipped by the selected water reference within the red elliptical boundary.The blue curve shows the HH+HV backscattering coefficient, the red curve shows the HH backscattering coefficient, and the green curve shows the HV backscattering coefficient.Theoretically, these histograms should display a bimodal distribution, but some of them were found to have unimodal distributions.Examples include water surfaces covered by floating or emerged plants or surface waves caused by strong winds.These situations cause a larger than normal amount of SAR energy to be reflected back to the sensor [11][12][13].This effect can be observed in the sample images for water references 024, 028, 063, and 076.
By considering the peaks, valleys, and curvatures of the smoothed histograms, the water references with unimodal distributions were rejected, and only those with bimodal distributions were taken into account.The bimodal occurrences for the HH and HH + HV polarizations are shown in Table 2.This number indicates the occurrence probability of the bimodal distribution for each image.In that sense, HH + HV is more likely to have a bimodal distribution and is more suitable for automatic classification.In other words, the HV polarization can improve the efficiency of water surface extraction.Thus, the extracted water areas presented in this article were derived from HH + HV.

Automatic Thresholding
Otsu's method (OT) is one of the best threshold selection methods for general gray-level images.This technique chooses the threshold value of the minimum within-class variance (σ 2 W ) or the maximum between-class variance (σ 2 B ) in Equation (1).Although this method can obtain satisfactory segmentation results in many cases, it is basically limited to images with background and foreground Gaussian distributions of equal variance.Therefore, images that do not meet this criterion may return unsatisfactory results, especially when the gray level histogram is unimodal or close to a unimodal distribution [33].
To address this weakness, many modifications of the Otsu method have been proposed.For example, the valley-emphasis method (VE), modified by weight σ 2 B with p(t), the complement of a probability at a threshold value t, causes the valley in the histogram to be more likely to be better determined.The neighborhood valley-emphasis (NE) is an improvement of the valley-emphasis method by weighting σ 2 B with the neighborhood information in n = 2m + 1 intervals at the threshold value t as in Equation ( 2).The result is closer to the valley of the histogram because it considers the neighborhood around the threshold point in addition to the threshold point.The optimal threshold is chosen by maximizing the between-class variance function [17] as in Equation (3). where is the between-class variance at threshold value t, µ is the mean value of all the intervals, µ w (t) is the mean value of the water portion at threshold value t, µ n (t) is the mean value of the non-water portion at threshold value t, p w (t) is the probability of the water portion at threshold value t, p n (t) is the probability of the non-water portion at threshold value t, m is the number of neighborhood intervals for threshold value t, p(t) is the probability of the interval at threshold value t, p(t) is the sum of the probabilities of the neighborhood interval around threshold value t, Arg max is the argument of the maxima for threshold value t in the function, and t * is the optimum threshold value.
From the previous study [17] and our empirical test for all the images using three thresholding methods (OT, VE, and NE), the NE method yielded the most suitable result for water extraction from the RS2 images.The threshold value of that method was close to the valley of the histogram and was able to extract water surfaces quite well.Consequently, the NE method was used for automatic thresholding in the following discussion.
All the histograms were calculated by dividing the entire range of values into 256 intervals and then setting the number of neighborhood intervals to 11 (m = 5).The results of the threshold values are displayed as vertical lines in the histogram plots in Figure 4.In the figure, the solid lines represent bimodal distributions, whereas the dashed lines represent unimodal distributions, which were excluded from the calculation of the global threshold.For example, water reference 076 was rejected when calculating the global threshold for HH + HV.The threshold values obtained from the bimodal distributions are the area-weighted average, resulting in the global threshold value.This threshold is expected to be effective for the entire image.All the global threshold values are listed in columns 10 and 11 of Table 2.
Among the 30 RS2 images, the global threshold values were slightly different.This may be because the RS2 images were acquired under different conditions.Different paths result in different azimuth angles, whereas different beam modes result in different incident angles, wave frequencies, and image resolutions.Those conditions produce different backscattering coefficients.The proposed technique has less bias for selecting the threshold value because it does not depend on the satellite mode and seasonal environment.Thus, the extracted water bodies for all the acquisition dates can be combined.

Water Area Extraction and Flood Duration Map
The water areas were simply extracted by applying a global threshold value to each RS2 image.Close-ups of the extracted result from the image acquired on 28 November 2011 are shown in Figure 5.The water boundaries obtained using this technique appear to be reasonable in comparison to the visually classified results from the original RS2 and TH1 images acquired in a dry season.All the RS2 images from the 30 dates were processed separately.Because some RS2 images did not cover the entire study area, the missing parts were estimated from the previous image.The results are shown in Figure 6

Accuracy Assessment
The accuracy of a classification must be assessed by comparing the results with truth data.In this study, optical images that captured ground surface activity and gaging stations that recorded water heights were introduced as truth data sources.During this flood event, almost all the country was covered by clouds, and optical satellite clear-sky images therefore were very rare.GeoEye-1 (GE1) and ThaiChote-1 (TH1) clear-sky images taken on 22 November 2011 over Ayutthaya city were selected as the truth data.Those data covered the Ayutthaya Historical Park and the Rojana and Hi-Tech industrial parks, which are very important economically and for tourism.
To extract the flooded areas from the optical images, using the Normalize Different Water Index (NDWI) calculated from the Green (G) and near-infrared (NIR) band values is the most popular and effective method [36].McFeeters proposed the NDWI in 1996 to detect surface waters in wetland environments and to allow for the measurement of surface water extent, and asserted that values of NDWI greater than zero are assumed to represent water surfaces, while values less than, or equal to, zero are assumed to be non-water surfaces [37].In this image interpretation, rivers and ponds were also extracted as flooded areas because nearly the entire study area was covered by water.The GE1 NDWI threshold for the water areas during the flood period was determined by visual Therefore, the period of inundation or a flood duration map is very important to the Thai government for controlling floods and developing remedial plans.The flood duration map was produced by stacking the interval of the extracted water surfaces.The final flood duration map is shown in Figure 7.   Similar to GE1, TH1 has an optical sensor with four bands with quite similar spectral ranges, but its spatial resolution is much lower.Before calculating the NDWI values, a TH1 image in Figure 8 (5A) was up-sampled and co-registered to the GE1 image.We then compared the NDWI in Figure 8 (5B) with that from the GE1 image.The most accurate TH1 NDWI threshold for water was found to be 0.28, which corresponds to 64.6% of the estimated water areas (W-W and N-W) in Figure 9(2A-2C).Among those extracted areas, 85.2% were similar to those from the GE1 (W-W and N-N), 9.6% were false negatives (W-N), and 5.1% were false positives (N-W).By applying this NDWI threshold to the TH1 pre-flood image, the water-covered ratio of this area was 11.9%.
For thresholding SAR sensors, the backscattering coefficient (the sigma naught value, σ °) is used more often than the NDWI for optical sensors.The best threshold for the RS2 HH+HV acquired on November 21, 2011 was determined to be 3.18 dB, which corresponds to 64.9% of the estimated water areas (W-W and N-W).Among those extracted water areas, 79.9% were similar to the results from GE1 (W-W and N-N), 12.2% were false negatives (W-N), and 7.9% were false positives (N-W).

Comparison between the Best Accuracy Method and the Proposed Method
The threshold value for RS2 from the most accurate thresholding method was greater, covered a larger area and was more accurate than the proposed weight averaged neighborhood valley-emphasis thresholding method.The RS2 threshold value acquired for 21 November 2011, (image 24) from the most accurate thresholding method was 3.18 dB (64.9% area coverage), whereas that by the proposed method was 2.08 dB (44.0%area coverage).For the proposed method, 71.9% of the results were similar to the results of GE1 (W-W and N-N), 26.6% were false negatives (W-N), and 1.5% were false positives (N-W).From Figure 9(3B,3C), the blue color shows the water areas extracted by the proposed method, and the green color shows the difference from the best accuracy method.
Although the threshold value of the best accuracy method is apparently more eligible, applying that local threshold value to the entire study area returned an inappropriate result and was difficult to implement as a practical procedure.The extraction of water surfaces using that threshold value for the entire area resulted in overestimation.For example, the water areas in Figures 5 and 6 became noisier.Moreover, it was almost impossible to find suitable optical images to be used as references for the SAR images throughout the event.

Accuracy Assessment
The accuracy of a classification must be assessed by comparing the results with truth data.In this study, optical images that captured ground surface activity and gaging stations that recorded water heights were introduced as truth data sources.During this flood event, almost all the country was covered by clouds, and optical satellite clear-sky images therefore were very rare.GeoEye-1 (GE1) and ThaiChote-1 (TH1) clear-sky images taken on 22 November 2011 over Ayutthaya city were selected as the truth data.Those data covered the Ayutthaya Historical Park and the Rojana and Hi-Tech industrial parks, which are very important economically and for tourism.
To extract the flooded areas from the optical images, using the Normalize Different Water Index (NDWI) calculated from the Green (G) and near-infrared (NIR) band values is the most popular and effective method [36].McFeeters proposed the NDWI in 1996 to detect surface waters in wetland environments and to allow for the measurement of surface water extent, and asserted that values of NDWI greater than zero are assumed to represent water surfaces, while values less than, or equal to, zero are assumed to be non-water surfaces [37].In this image interpretation, rivers and ponds were also extracted as flooded areas because nearly the entire study area was covered by water.The GE1 NDWI threshold for the water areas during the flood period was determined by visual interpretation of Figure 8

Accuracy Assessment
The accuracy of a classification must be assessed by comparing the results with truth data.In this study, optical images that captured ground surface activity and gaging stations that recorded water heights were introduced as truth data sources.During this flood event, almost all the country was covered by clouds, and optical satellite clear-sky images therefore were very rare.GeoEye-1 (GE1) and ThaiChote-1 (TH1) clear-sky images taken on 22 November 2011 over Ayutthaya city were selected as the truth data.Those data covered the Ayutthaya Historical Park and the Rojana and Hi-Tech industrial parks, which are very important economically and for tourism.
To extract the flooded areas from the optical images, using the Normalize Different Water Index (NDWI) calculated from the Green (G) and near-infrared (NIR) band values is the most popular and effective method [36].McFeeters proposed the NDWI in 1996 to detect surface waters in wetland environments and to allow for the measurement of surface water extent, and asserted that values of NDWI greater than zero are assumed to represent water surfaces, while values less than, or equal to, zero are assumed to be non-water surfaces [37].In this image interpretation, rivers and ponds were also extracted as flooded areas because nearly the entire study area was covered by water.The GE1 NDWI threshold for the water areas during the flood period was determined by visual interpretation of Figure 8(4B) as NDWI ≧ 0.02 or 69.1% of the image area shown in Figure 9(1A-1C).This NDWI threshold value was used as the truth data when obtaining the most accurate threshold values for NDWI from TH1 and HH + HV from RS2.The green pixels are the determined best accuracy results from RADARSAT-2, which was larger than that determined using the neighborhood valley-emphasis method.

Finding the Best Accuracy Thresholds for the ThaiChote-1 and RADARSAT-2 Images
The comparison of two (the truth data and estimation) two-class spatial images, water areas (W) or non-water areas (N) results in four combinations: W-W, N-N, W-N, and N-W.When the threshold for a client's image is set to the minimum value, all of the results will be non-water areas (N).Some of them are N-N, which represents the same N values as those from the master image (GE1), whereas the others are W-N, which represents a water extraction omission error (false negative).When the threshold of the client image moves to higher values, some values will be considered to be water (W).Some of these values are W-W, which represent the same W values, whereas the others are N-W, which represent overestimation (commission error, false positive) in the water extraction.The best threshold value is the point where the summation of the W-W and N-N areas becomes largest.This point corresponds to the result most similar to that from the GE1 image.The results of this approach are shown in Table 3.Similar to GE1, TH1 has an optical sensor with four bands with quite similar spectral ranges, but its spatial resolution is much lower.Before calculating the NDWI values, a TH1 image in Figure 8(5A) was up-sampled and co-registered to the GE1 image.We then compared the NDWI in Figure 8(5B) with that from the GE1 image.The most accurate TH1 NDWI threshold for water was found to be 0.28, which corresponds to 64.6% of the estimated water areas (W-W and N-W) in Figure 9(2A-2C).Among those extracted areas, 85.2% were similar to those from the GE1 (W-W N-N), 9.6% were false negatives (W-N), and 5.1% were false positives (N-W).By applying this NDWI threshold to the TH1 pre-flood image, the water-covered ratio of this area was 11.9%.
For thresholding SAR sensors, the backscattering coefficient (the sigma naught value, σ • ) is used more often than the NDWI for optical sensors.The best threshold for the RS2 HH+HV acquired on November 21, 2011 was determined to be 3.18 dB, which corresponds to 64.9% of the estimated water areas (W-W and N-W).Among those extracted water areas, 79.9% were similar to the results from GE1 (W-W and N-N), 12.2% were false negatives (W-N), and 7.9% were false positives (N-W).

Comparison between the Best Accuracy Method and the Proposed Method
The threshold value for RS2 from the most accurate thresholding method was greater, covered a larger area and was more accurate than the proposed weight averaged neighborhood valley-emphasis thresholding method.The RS2 threshold value acquired for 21 November 2011, (image 24) from the most accurate thresholding method was 3.18 dB (64.9% area coverage), whereas that by the proposed method was 2.08 dB (44.0%area coverage).For the proposed method, 71.9% of the results were similar to the results of GE1 (W-W and N-N), 26.6% were false negatives (W-N), and 1.5% were false positives (N-W).From Figure 9(3B,3C), the blue color shows the water areas extracted by the proposed method, and the green color shows the difference from the best accuracy method.
Although the threshold value of the best accuracy method is apparently more eligible, applying that local threshold value to the entire study area returned an inappropriate result and was difficult to implement as a practical procedure.The extraction of water surfaces using that threshold value for the entire area resulted in overestimation.For example, the water areas in Figures 5 and 6 became noisier.Moreover, it was almost impossible to find suitable optical images to be used as references for the SAR images throughout the event.

Comparison with the Gaging Station Data
Because the satellite images were not taken from the same sensors and were not acquired on the same date, it is difficult to explain the cause of the different extracted water extents.The difference in sensor types and spatial resolutions and changes in water height over time may have contributed to the discrepancy.The heights of the flood waters were difficult to project because water flows downward without stopping, although water flows rather slowly in this area.To understand the flood event, the daily average water heights above the mean sea level (MSL) that were collected from the three nearest telemetry gaging stations (C35, C37, and S5) and recorded by the Royal Irrigation Department (RID) [38] and the water depths from the surface of the road at four checkpoints reported by Rojana Industrial Park Public Co., Ltd.[39] were introduced as truth data, as shown in Figure 10.The four water checkpoints were located in the Rojana Industrial Park at the power plant (RJ0), Rojana phase-1 gate-A in front of the head office (RJ1), Rojana phase-2 gate-B (RJ2), and Rojana-3 around the flyover (RJ3).Among them, RJ0 was closest to the Honda automobile factory shown in Figure 9.

Comparison with the Gaging Station Data
Because the satellite images were not taken from the same sensors and were not acquired on the same date, it is difficult to explain the cause of the different extracted water extents.The difference in sensor types and spatial resolutions and changes in water height over time may have contributed to the discrepancy.The heights of the flood waters were difficult to project because water flows downward without stopping, although water flows rather slowly in this area.To understand the flood event, the daily average water heights above the mean sea level (MSL) that were collected from the three nearest telemetry gaging stations (C35, C37, and S5) and recorded by the Royal Irrigation Department (RID) [38] and the water depths from the surface of the road at four checkpoints reported by Rojana Industrial Park Public Co., Ltd.[39] were introduced as truth data, as shown in Figure 10.The four water checkpoints were located in the Rojana Industrial Park at the power plant (RJ0), Rojana phase-1 gate-A in front of the head office (RJ1), Rojana phase-2 gate-B (RJ2), and Rojana-3 around the flyover (RJ3).Among them, RJ0 was closest to the Honda automobile factory shown in Figure 9.
The solid lines in the left graph show the water heights at the three telemetry gaging stations over a one-year period (April 2011 to March 2012).The period over which the satellite images were acquired was at the end of the flood event, when the water level had dramatically decreased.Although the water heights were slightly different on 22 November 2011, they had significantly dropped more than 30 cm in three days by 25  MSL (m) The solid lines in the left graph show the water heights at the three telemetry gaging stations over a one-year period (April 2011 to March 2012).The period over which the satellite images were acquired was at the end of the flood event, when the water level had dramatically decreased.Although the water heights were slightly different on 22 November 2011, they had significantly dropped more than 30 cm in three days by 25 November 2011.
Because Ayutthaya is situated on a relatively flat plain, a small increase in the water height may cause the water to spread over a wide area.At that time, only the water at station C37 was higher than the levee, and the water at this station spread beyond the river.Although the water heights at the other stations (S5 and C35) were below the height of the levee and no water spread outside, the water remained on the ground.At water checkpoint RJ0, in front of the power plant and close to the Honda automobile factory, the water depth above the road surface was 0.80 m on November 21 and deceased to 0.69 m on 22 November and the water had completely drained by November 25.The extracted areas determined in this study may have been slightly different because the satellite images were acquired on different dates.

Discussion and Future Work
The proposed method can improve the speed of automated thresholding because to determine the global threshold value from static local areas is much more efficient than to do from systematic hierarchical local areas.This technique is suitable for multi-temporal images in a specific place, like the Chao Phraya river basin.Although the processing for permanent water references is rapid, the selection of proper water references takes time at the first step.An improper water reference without a bimodal distribution histogram should not be considered because it does not contribute to the estimation of the global thresholding.
Open water from SAR images can be extracted quite well, but it is difficult to extract water in inundated urban areas and under trees [25,27,40].The most important factor that limits the use of satellite data are resolution.When a satellite records data in a high-resolution mode, the imaging swath becomes smaller than that in the normal mode.This fact is caused by limitations of the data recorder and transmitter.Observing an inundation over a wide area with the same sensor or the same satellite is almost impossible.Although RS2 can observe in the spotlight mode at a 1-m resolution, the image on November 22 was acquired in the SCNB mode at a 25-m resolution to cover a larger area.On the other hand, the resolutions of pan-sharpened GE1 and TH1 images are 1 m and 2 m, respectively.Therefore, the GE1 and RS2 results are quite different (approximately 23-33%), according to the accuracy assessment.For detecting the water surface, a large number of pixels were missing because the overall backscattering of the pixels surrounded by obstacles such as buildings was higher than for open spaces.
Sensor type can also cause differences.Although the NDWI obtained from optical sensors is suitable for the detection of water surfaces, objects having low infrared radiation are also classified as water surfaces.This misclassification can be observed in Figure 9(1C) and (2C).The roof-tops of buildings were sometimes identified as flood water by GE1, and there were more false-positive classifications from TH1. Water surfaces may also be misclassified by the backscattering coefficient from a SAR sensor.For example, non-water objects with smooth surfaces such as roads and runways are classified as water surfaces.Furthermore, side-looking transmission hinders observational ability.Water surfaces next to buildings produce double bounce backscattering, resulting in a total backscattering coefficient that is higher than usual and misclassifications as non-water.This effect can be observed in Figure 9(3C), where the west sides of the buildings were misclassified because the layover was projected onto the water surface.Similarly, the water surfaces between the buildings were misclassified.
Observing inundations using satellites has another drawback; water under roofs or inside buildings cannot be observed.Auxiliary data must be prepared to provide the missing information.The average ground elevation data were prepared from an aerial survey conducted by the Japan International Cooperation Agency (JICA) for the zone along the Chao Phraya river and by ESRI Thailand for the outer area in 2012.The combined data can cover the entire Chao Phraya river basin.Future research needs to focus on improving accuracy through the utilization of LiDAR DEMs as topographic data.

Conclusions
Multi-temporal RADARSAT-2 images with different acquisition conditions were used to extract water areas from the 2011 central Thailand flood along the Chao Phraya river.Considering the use of satellite SAR data in emergency situations, where validation data are scarce and optical images are hindered by cloud cover, an automated thresholding approach for water extraction was attempted.By introducing the global threshold value of the entire study area for each SAR image, the weight-averaged neighborhood valley-emphasis method was able to extract flooded areas automatically from the backscattering coefficient.The proposed method, which obtains the threshold value from the water references located throughout the Chao Phraya river basin, is more suitable for Thailand.Moreover, this technique is independent of the satellite acquisition condition, and there is no need to normalize the images before combining all the extracted water-areas as a flood duration map.In this case, the HH+HV dual-polarization achieved a higher accuracy than the HH single-polarization for open water extraction, which is affected by winds and floating/submerged plants.The extraction result for Ayutthaya city was similar to the visual inspection result from a GeoEye-1 image (approximately 67%), whereas the result obtained by the best accuracy thresholding method was approximately 77%.The results based on the ThaiChote-1 image were more similar to that from the GeoEye-1 image (approximately 82%) because they were both obtained by high-resolution optical satellites.On the other hand, the results from RADARSAT-2 had lower accuracies than those from GeoEye-1 and ThaiChote-1 because of its lower spatial-resolution and side-looking observational scheme.SAR images also have limitations for observing water areas covered by trees or adjacent to buildings.Despite these obstacles, the extraction of flooded areas from SAR intensity data can be improved by introducing pre-event topographic data such as LiDAR DEMs and building footprints.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/9/1/78/s1.For Figure 4: Automated threshold values determined using the neighborhood valley method for all the water references and for all the RADARSAT-2 images.

19 Figure 1 .
Figure 1.Map showing the study area (left), where the black border delimits the Chao Phraya river basin in central Thailand, and the red square delimits Ayutthaya City.Topographic map (middle), with main streams on the DEM.Water reference map showing 87 water references throughout the study area on a RADARSAT-2 image (right).

Figure 1 .
Figure 1.Map showing the study area (left), where the black border delimits the Chao Phraya river basin in central Thailand, and the red square delimits Ayutthaya City.Topographic map (middle), with main streams on the DEM.Water reference map showing 87 water references throughout the study area on a RADARSAT-2 image (right).

19 Figure 2 .
Figure 2. RADARSAT-2 color composite (HH as red, HV as green, and HH/HV as blue) images taken during the 2011 Thai flood event from the beginning to end.

Figure 2 .
Figure 2. RADARSAT-2 color composite (HH as red, HV as green, and HH/HV as blue) images taken during the 2011 Thai flood event from the beginning to end.

Figure 3 .
Figure 3. Work flow diagram for automatic flood extraction for RADARSAT-2.The optical image procedure shown by the dashed lines was used for a one-time accuracy assessment.

Figure 3 .
Figure 3. Work flow diagram for automatic flood extraction for RADARSAT-2.The optical image procedure shown by the dashed lines was used for a one-time accuracy assessment.

Figure 4 .
Figure 4. Automated threshold values using the neighborhood valley method of water references in the red elliptical areas (10 of 87) from the HH, HV and HH + HV sigma-naught values taken on the evening of 28 November 2011 PM (1 of 30 in Table2).The dashed lines represent unimodal distributions, and the solid lines represent bimodal distributions.The ThaiChote-1 images in the first column were acquired on different dates in the dry season.

Figure 4 .
Figure 4. Automated threshold values using the neighborhood valley method of water references in the red elliptical areas (10 of 87) from the HH, HV and HH + HV sigma-naught values taken on the evening of 28 November 2011 PM (1 of 30 in Table2).The dashed lines represent unimodal distributions, and the solid lines represent bimodal distributions.The ThaiChote-1 images in the first column were acquired on different dates in the dry season.
. The results from the current date are shown in deep blue, and the estimated results from the previous date are shown in light blue.The estimations coincided with the actual flood situation; flooding began in the north region on 2 October 2011 (image 1), the flooded areas spread to the south on 23 October 2011 (image 4), the flood in the north had abated by 21 November 2011 (image 28), and was completely finished on 14 February 2012 (image 30).

Figure 5 .
Figure 5. Close-ups of the extracted water extents from the RADARSAT-2 image acquired on 28 November 2011 using the HH+HV global threshold value (1.92 dB).The water areas are displayed in deep blue, and the red ellipses are the water references.

Figure 5 .
Figure 5. Close-ups of the extracted water extents from the RADARSAT-2 image acquired on 28 November 2011 using the HH + HV global threshold value (1.92 dB).The water areas are displayed in deep blue, and the red ellipses are the water references.

19 Figure 6 .
Figure 6.Extracted water areas after applying the global threshold value for each Radarsat-2 image.The extracted results from the actual date are shown in deep blue, and the results from the previous date are shown in light blue.

Figure 6 .
Figure 6.Extracted water areas after applying the global threshold value for each Radarsat-2 image.The extracted results from the actual date are shown in deep blue, and the results from the previous date are shown in light blue.

Figure 7 .
Figure 7. Water duration map calculated by stacking the extracted water surfaces.The red border shows Ayutthaya City, an enlargement of which is shown on the right side of the figure.The yellow border shows the Rojana industrial park.
Gaging Station  Ayutthaya Historical Park  Rojana Industrial Estate  Hi-Tech Industrial Estate

Figure 7 .
Figure 7. Water duration map calculated by stacking the extracted water surfaces.The red border shows Ayutthaya City, an enlargement of which is shown on the right side of the figure.The yellow border shows the Rojana industrial park.

19 Figure 7 .
Figure 7. Water duration map calculated by stacking the extracted water surfaces.The red border shows Ayutthaya City, an enlargement of which is shown on the right side of the figure.The yellow border shows the Rojana industrial park.

Figure 9 .Figure 8 .Figure 8 .
Figure 9. Histograms and cumulative probability plots (top) and extracted flooded areas (blue color) in Ayutthaya city (middle) and over part of the Rojana Industrial Park (bottom), from the visualization of the GeoEye-1 image (1A-1C), determining the best accuracy of the ThaiChote-1 image (2A-2C), and the emphasis of the neighborhood valley on the RADARSAT-2 image (3A-3C), plotted on a GeoEye-1 false color composite.The green pixels are the determined best accuracy results from RADARSAT-2, which was larger than that determined using the neighborhood valley-emphasis method.

Figure 9 .Figure 9 .
Figure 9. Histograms and cumulative probability plots (top) and extracted flooded areas (blue color) in Ayutthaya city (middle) and over part of the Rojana Industrial Park (bottom), from the visualization of the GeoEye-1 image (1A-1C), determining the best accuracy of the ThaiChote-1 image (2A-2C), and the emphasis of the neighborhood valley on the RADARSAT-2 image (3A-3C), plotted on a GeoEye-1 false color composite.The green pixels are the determined best accuracy results from RADARSAT-2, which was larger than that determined using the neighborhood valley-emphasis method.

Figure 10 .
Figure 10.Comparison of the water heights recorded at three telemetry gaging stations and the levee heights (left); enlargement of the left plot for 14-27 October and the water depths above the road surface observed at four checkpoints in the Rojana Industrial Park (right).

Figure 10 .
Figure 10.Comparison of the water heights recorded at three telemetry gaging stations and the levee heights (left); enlargement of the left plot for 14-27 October and the water depths above the road surface observed at four checkpoints in the Rojana Industrial Park (right).

Table 1 .
Summary of the RADARSAT-2 image data properties used in this study.

Table 1 .
Summary of the RADARSAT-2 image data properties used in this study.

Table 2 .
List of RADARSAT-2 images and threshold values for water extraction.

Table 3 .
Results of the accuracy assessment of the water extraction for Ayutthaya city.

Table 3 .
Results of the accuracy assessment of the water extraction for Ayutthaya city.