An automatic method for mapping inland surface waterbodies with Radarsat-2 imagery

The use of synthetic aperture radar (SAR) imagery is generally considered to be an effective method for detecting surface water. Among various supervised/unsupervised classification methods, a SAR-intensity-based histogram thresholding method is widely used to distinguish waterbodies from land. A SAR texture-based automatic thresholding method is presented in this article. The use of texture images substantially enhances the contrast between water and land in intensity images. It also makes the method less sensitive to incidence angles than intensity-based methods. A modified Otsu thresholding algorithm is applied to selected sub-images to determine the optimal threshold value. The sub-images were selected using k-means results to ensure a sufficient number of pixels for both water and land classes. This is critical for the Otsu algorithm being able to detect an optimal threshold for a SAR image. The method is completely unsupervised and is suitable for large SAR image scenes. Tests of this method on a Radasat-2 image mosaicked from 8 QuadPol scenes covering the Spritiwood valley in Manitoba, Canada, show a substantial increase in land–water classification accuracy over the commonly used SAR intensity thresholding method (kappa indices are 0.89 vs. 0.79). The method is less computationally intensive and requires less user interaction. It is therefore well suited for detecting waterbodies and monitoring their dynamic changes from a large SAR image scene in a near-real time environment).


Introduction
Freshwater is crucial for much of life on Earth and is an essential part of the natural environment. Almost 9%, or 900,000 km 2 , of Canada's total landmass is covered by freshwater in the form of lakes, rivers, streams, etc. There are over 2 million lakes in Canada, of which 910,400 have an area greater than 0.1 km 2 , 84,516 have an area greater than 1 km 2 , and 564 have an area greater than 100 km 2 (Minns et al. 2008). Lakes and rivers are important for water supply, fishing, and recreation. Many large lakes in Canada support commercial and subsistence fisheries which provide vital sources of food and income for people, especially the First Nations (Minns et al. 2008). The inland waterbodies play an important role in the terrestrial water cycles and surface water budget (Wang, Huang, Rivera, et al. 2014;Wang, Huang,Yang, et al. 2014;, which strongly affect the atmosphere and surface/subsurface processes such as cloud development (Molders and Raabe 1996), surface albedo (Wang et al. 2006), evapotranspiration (Wang et al. 2013), stream flow (Koster and Milly 1997), and groundwater recharge (Sophocleous 2002). Surface waters are also integral parts of groundwater flow systems and hence are indicators of the status of the overall freshwater resource. For example, lakes gain water from groundwater systems and also are a source of groundwater recharge (Winter et al. 1998). Decreasing groundwater levels may result in reduced areal extent of lakes and less stable water temperatures. The available freshwater resources are under increasing pressure from a wide range of human activities. The greatest threat is the changing climate which is bringing large changes to surface water such as falling water levels, changing lake distribution, and even the complete disappearance of lakes, especially in permafrost regions of Canada (Riordan, Verbyla, and McGuire 2006;Smith et al. 2005;Yoshikawa and Hinzman 2003). To successfully manage this valuable resource requires an understanding that can only be achieved by improved monitoring. To effectively conserve and manage the fresh surface water resources, it is essential to have up-to-date information of their spatial and temporal variability. Due to the large extent of surface water in Canada, with a significant portion in remote areas, satellite remote sensing is the only practical approach that can map surface water cost-effectively and in a timely manner (Rundquist, Narumalani, and Narayanan 2001).
Various optical satellite sensors such as the Advanced Very High Resolution Radiometer (AVHRR), the Moderate-Resolution Imaging Spectroradiometer (MODIS), Landsat's Multispectral Scanner (MSS) and Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+), as well as the Satellite Pour l'Observation de la Terre (SPOT) High Resolution Visible (HRV) instrument have been employed for surface waterbody detection (Campos, Sillero, and Brito 2012;Du et al. 2012;Giardino et al. 2010;Huang, Li, and Xu 2012;Jain et al. 2006;Ma et al. 2007;Sheng, Gong, and Xiao 2001;Tulbure and Broich 2013). The multispectral nature of optical sensors provides some advantages for water detection; however, their applications in detecting surface water are constrained by several environmental factors, such as cloudy sky conditions, cloud shadows, smoke from wildfires, and haze, etc. Synthetic aperture radar (SAR) is able to penetrate cloud, haze, and smoke, and hence observe the Earth's surface in all weather conditions, day and night. SAR sensors thus have the ability to provide data for surface waterbody detection that can overcome the limitations of optical sensors. In fact, SAR is generally considered an effective tool for detecting surface water (Brisco et al. 2009) and has been used for flood detection (Giustarini et al. 2013;Kuenzer et al. 2013;Lu et al. 2014;Martinis, Twele, and Voigt 2009), monitoring open water dynamics (Bartsch et al. 2012), and delineating shorelines (Shu, Li, and Gomes 2010). Smooth water surfaces usually provide a specular reflection of microwave radiation, and hence very little energy is scattered back. In contrast, land surfaces scatter much more energy back to the radar due to, for example, surface roughness and volume scattering. The difference in the energy received back leads to a high contrast between water and land.
The contrast between waterbodies and their surrounding land is highly affected by the roughness of water surfaces and SAR parameters including wavelength, incidence angle, and polarization (Martinis 2010). The specular reflector is a simplified model for smooth surface water. The effects of wind and current can lead to rough water surfaces (e.g. waves), which cause a higher backscatter signal (Martinis 2010). This may lead to a lower contrast ratio between surface water and land (Solbo et al. 2003), which affects the ability to identify the respective signal return from waterbodies and land. Several studies have shown that the cross-polarizations HV and VH are less affected by water surface roughness induced by winds and currents than HH and VV polarization (Henry et al. 2006;Schumann et al. 2007). Over smooth water surfaces, like-polarization (e.g. HH) offers improved land/water separability compared to crosspolarization. An increase in surface roughness reduces the ability to discriminate between water and land in VV more than in HH polarization (Martinis 2010). The choice of polarization thus plays an important role in surface water detection using SAR data. A shorter SAR wavelength has greater backscattering from smooth water surfaces because of increased sensitivity to diffuse scattering. A longer SAR wavelength is less sensitive to the short amplitudes of landscape features and consequently has a reduced signal separation between land and water (Drake and Shuchman 1974). The strongest signal separation and highest contrast between water and land occurs in a SAR image with shorter wavelength (Martinis 2010). C-band is thus more suitable for delineating surface waterbodies from land than L-band SAR.
Among the various supervised/unsupervised classification methods, histogram thresholding is one of the most popular approaches used to delineate waterbodies from land in SAR intensity imagery (Brisco et al. 2009;Brivio et al. 2002;Martinis 2010). The threshold value approach to separating water from land in a SAR image depends on a number of factors that can vary temporally (e.g. wind, soil moisture), spatially (landscape roughness, vegetation zonation), and by sensor (e. g. frequency), and also on factors that are processing specific (SAR incidence angle and polarization choice) (Martinis 2010). The threshold value is determined by visual inspection of image histograms, quick checking of the classification results, and refinement of the threshold if necessary until the classification result is satisfactory (Brisco et al. 2009). Such an interactive thresholding method is analyst dependent and time consuming, and thus extensive research has been carried out into semi-automatic or automatic thresholding of SAR imagery (Brisco et al. 2009;Fan and Lei 2012;Hahmann and Wessel 2010;Kandus et al. 2001;Solbo et al. 2003). Most applications of histogram thresholding methods have been limited to specific areas, and the methods mostly require some manual interaction.
Generally, the grey contrast between surface water and land in a SAR image decreases with decreasing incidence angle (Drake and Patton 1980;Malnes, Guneriussen, and Høgda 2002;Martinis 2010). This means that higher incidence angle SAR data are preferable for accurately extracting surface waterbodies. This requirement reduces by a considerable amount the proportion of SAR data that is suitable for surface water detection (Solbo et al. 2003). This suggests that a method for surface water detection independent of the SAR incidence angle is needed.
An automatic method for detecting inland surface water using Radarsat-2 imagery is presented in this article. The method is based on histogram thresholding of the texture image, which is not significantly influenced by the incidence angle.

Study area and data sets
The study area covers the Spiritwood buried valley (red polygon in Figure 1) which extends from southwest Manitoba, Canada, to North Dakota, USA. The dominant landcover class is agricultural crops, intermingled with grassland, surface waterbodies and forests. Surface water is fairly abundant in this area, in numerous small waterbodies of the prairie landscape.
Radarsat-2 QuadPol (fully polarimetric mode) wide fine-beam-mode images (with a swath of 50 km and an azimuth resolution of 7.6 m) were used in the analysis. To cover the entire study area, two neighbouring paths (4 scenes to each path) were acquired on 26 July 2013 (FQ1W with incidence angle~20.0°) and 2 August 2013 (FQ3W with incidence angle~21.5°) and processed into a single-look complex product type. All eight scenes of the images were first processed to generate calibrated multilook backscatter images and then auto-orthorectified using the GAMMA program. The eight orthorectified images were mosaicked (with radiometric normalization) to one image with a resolution of 30 m. A 5 × 5 enhanced Lee filter was applied to the mosaic image for speckle noise reduction. The background image in Figure 1 shows the composite of the HH (red), HV (green), and VV (blue) polarized images.
To evaluate the performance of the proposed method, a reference map was created for a small portion of the Radarsat-2 mosaic image from the 5 m-resolution SPOT-5 pansharpened multispectral image by visual interpretation and manual digitizing. The highlighted blue square in Figure 1 indicates the location and coverage of the reference map. 10% of the area covered by the reference map consists of surface water, shown in yellow, and 90% is land.

Method
The proposed automated thresholding method is based on the Otsu algorithm, which is one of the best threshold selection methods for image binarization (Fan and Lei 2012). The Otsu algorithm selects the threshold value that maximizes the between-class variances of the histogram. This is optimal for thresholding large objects from the background, which means that it is ideal for thresholding a histogram with bimodal or multimodal distribution (Ng 2006). Generally speaking, SAR images do not have such characteristics. Water often occupies a small portion in a large SAR image scene and thus a SAR image histogram may be unimodal. The Otsu algorithm may also fail to select the optimal threshold value in a large SAR image scene. To resolve the issue of uneven distribution of water and land classes in SAR images, a sub-image selection process can be applied to divide the original image into smaller sub-image tiles with a user-defined size. This provides a data set that is likely to have a higher probability of containing a sufficient proportion of both classes. Only the selected sub-images are used for the Otsu thresholding. To fulfil the sub-image selecting process, it is necessary to obtain an initial classification map of water and land. The histogram thresholding method only works on single polarization images. Thus, the selection of optimal polarization as the method's input is necessary if multiple polarizations are available. Multiple criteria can be employed in this process, for example, weather on the SAR image acquisition date can be assessed and visual analysis of the wind effects on the images can be reviewed. An overview of the automated thresholding method for surface water detection in a SAR image is shown in Figure 2. The main steps of this method are summarized below: Step 1: k-means cluster analysis. An initial waterbody mask, as well as a lowbackscatter image mask, is created by applying k-means clustering to the SAR intensity image.
Step 2: Sub-image selection. This step selects sub-images containing a sufficient proportion of land and water classes.
Step 3: Entropy texture image histogram analysis. This step generates the entropy texture of the SAR intensity image. The histogram of the sub-images selected at step 2 is then computed.
Step 4: Otsu algorithm for automatically determining the threshold value. A mask image with low entropy is generated by applying a modified Otsu algorithm to the histogram computed in step 3.
Step 5: Generation of the final waterbody image. The low-backscatter image mask created in step 1 is used to refine the low-entropy mask image to obtain the final waterbody image.
The purpose of the first three steps is to optimize the input SAR image for an optimal auto-image thresholding, which is implemented with the Otsu algorithm. The specific details of this method are explained in the following sub-sections.

k-means clustering analysis
The initial step in the analysis is the application of the k-means algorithm to a single polarization SAR image to generate an initial waterbody mask for further analysis. The kmeans algorithm is first applied to the SAR intensity image to obtain a map with an initial k clusters. The cluster map encodes each cluster with a unique grey-level value. The cluster number is represented by the grey level. For example, cluster 1 is assigned a grey-level value of 1 corresponding to dark appearance in the SAR intensity image, and cluster 2 is assigned a grey-level value of 2, which may correspond to a mixture of water and land. The initial waterbody mask is generated by assigning cluster 1 as water and other clusters as land, in which the water class may include some land pixels. A major problem with the kmeans algorithm is that the number of clusters (k) should be decided a priori. A limitation exists in finding the number of homogeneous objects (k) for large natural scenes, but this is not the case in this study because only an initial (not accurate) waterbody map is needed at this step. Tests of different values of k ranging from 10 to 20 on several different SAR scenes over different areas show that all values produce similar results and capture most (>85%) of the areal extent of the water (cluster 1), which is sufficient for the initial water mask. In this study, k is set to 15. In addition, a low-backscatter image mask is obtained by regrouping clusters 1-7 as low backscatter and others as high backscatter.

Sub-image selecting
To complete the sub-image selection, a SAR image is divided into M non-overlapping sub-images with an analyst-defined size of w × w. The selection of w, which determines M, depends on the extent of the two classes, water and land, within the SAR image (Martinis, Twele, and Voigt 2009). Due to the fact that the Ostu algorithm is optimal for thresholding a histogram with a bimodal distribution, only sub-images that contain a sufficient number of pixels from both water and land classes are selected for the threshold computation. Martinis, Twele, and Voigt (2009) introduced two statistical measures for the contrast within an image, on which the selection of appropriate sub-images was based.
However, for SAR images the empirical threshold values of the measures vary from image to image. In this study, since the initial water/land classification is available, appropriate sub-images can be easily selected without defining the threshold values for the measures of the contrast in the SAR image. According to Bazi, Bruzzone, and Melgani (2007), for accurate detection of threshold values, it is sufficient if each class has at least 10% of the pixels in an image belonging to it. Very highly reflecting pixels (cluster number > 7), for example those corresponding to urban areas, are excluded from the computation of the proportions of the water and land classes. This process prohibits the mis-selection of subimages that have tri-modal histograms containing more than two classes. The criterion for the water proportion in the sub-images to be selected is set to 10%-90%. If no sub-image meets this criterion, the image partitioning process is repeated by decreasing w by 10 each time until sub-images with adequate water content are successfully extracted.

Entropy texture analysis
Texture is defined as the tonal variation within a neighbourhood and thus reflects the spatial relations between pixels in an image. The texture information derived from a SAR image is a valuable feature for discriminating between different land-cover types, and thus texture analysis has been widely used for image segmentation and land-cover classification. Some studies have indicated that texture information is more useful than the intensity image for SAR image classification and segmentation (Song, Sohn, and Park 2007). Solbo et al. (2003) demonstrated that a texture-based surface water detection method performed better than the SAR intensity image thresholding method. Among the commonly used texture measures (e.g. variance, entropy, contrast, etc.), computed using grey-level cooccurrence matrix (GLCM), the entropy texture is often used in SAR image segmentation (Kekre and Gharge 2010;Samanta and Paul 2011;Samanta and Sanyal 2012). In this study, the surface water auto-detection strategy is proposed by using the entropy texture derived from the GLCM of SAR imagery The GLCM is a two-dimensional array, indicating the frequency of a pair of pixels within a local window. Each element p(l 1 ,l 2 ) represents the probability of the occurrence of the pair of grey levels (l 1 ,l 2 ) separated by a given distance d at angle θ. The entropy H of a SAR intensity image is computed by: where L is the number of grey levels of the SAR image. The entropy is a measure of nonuniformity in the SAR image. The histogram of the sub-images selected in step 2 is then computed.

Otsu algorithm
The Otsu algorithm is an automatic threshold selection method for the reduction of a greyscale image to a binary image containing two classes. It selects an optimal threshold value separating those two classes so that the between-class variance is maximized. The Otsu algorithm does not depend on modelling the probability density function. It assumes a bimodal distribution of grey-level values or two classes of pixels in an image. The advantage of this algorithm is that only the grey-level histogram is needed to derive an image threshold without any other a priori knowledge. For a given image with L grey levels, the number of pixels with grey-level i is denoted by n i . If N is the total number of pixels, the probability p i of the occurrence of grey-level i is then calculated as: The total mean (µ T ) grey-level of the entire image is calculated as: In the case of a single threshold, k, the pixels in the image are divided into two classes C 1 = [0, 1, ……, k] and C 2 = [k + 1, k + 2, ……, L − 1] based on grey levels. The probabilities of two classes are calculated as: The mean values of the two classes are: The between-class variance σ B 2 (k) of the two classes, C 1 and C 2 , is defined as: The optimal threshold k* can be determined by maximizing the between-class variance σ 2 B k ð Þ as: The Otsu algorithm works well when the image histogram is close to a bimodal distribution with equal variances (Fan and Lei 2012). In this study, although the sub-image selection process is employed to choose sub-images with bimodal distribution, the two classes may not have equal variances, which may result in an incorrect threshold value when applying the Otsu method. To improve the threshold selections in these cases, we use a modified Otsu method named valley-emphasis method proposed by Ng (2006). The valley-emphasis method makes the threshold closer to the actual valley of the histogram. It selects a threshold value that has small probability of occurrence (valley in the greylevel histogram) and also maximizes the between-class variance (Ng 2006). The optimal threshold k* is selected as: The multiplication of the weight 1 − p k by the between-class variance ensures that the selected threshold will always be a value that resides in the valley or on the bottom rim of the grey-level distribution. The histogram computed in step 3 first has to be normalized to a certain grey-level, L (e.g. 255). In this case, the entropy texture image will be transformed to contain values lying in the range 0-L (e.g. 255). An optimal threshold k* for the normalized entropy texture image can then be determined by applying the valley-emphasis Otsu method to the normalized histogram.

Generation of final surface waterbody image
A low-entropy mask image is obtained by applying the selected threshold k* to the normalized entropy texture image. In addition to the waterbodies, the low-entropy mask image may also contain some smooth wetland and agriculture areas, where standing water exists. These areas have a bright tone in the SAR intensity image due to double-bounce return of standing water and plants. Therefore, with the help of the low-backscatter mask image created in step 1, we can generate the final waterbody map by masking out the wetland and agriculture areas in the low-entropy mask image.

Results and discussion
The proposed histogram thresholding method (the 'texture method') only works on single polarization images, thus the selection of optimal polarization as the input is necessary. According to previous analysis, like-polarization (e.g. HH) images rather than cross-polarization images are selected if there is no wind or if the effects of the wind are small. If the water surfaces do not appear as smooth surfaces, a cross-polarization image (e.g. HV) is chosen. Figure 3 shows the effects of the wind on the QuadPol mosaics in different polarizations: (a) HH, (b) HV, (c) VH, and (d) VV. The HV-and VH-polarized images are less affected by the wind than the HH-and VV-polarized images. After checking the weather data for the acquisition dates of the SAR images and visually examining the wind effects, we selected the HV-polarized image for the Radarsat-2 QuadPol mosaic data as the input to the proposed method. For the purposes of performance assessment, we also tested a 'traditional' method (the 'intensity method'), in which the valley-enhanced Otsu algorithm is directly applied to the image intensity rather than the entropy texture after step 2 of the proposed method. The results obtained from both methods were then compared against the water polygons extracted from the 5 m-resolution SPOT-5 pan-sharpened multispectral image.
The initial waterbody mask and the low-backscatter mask were generated by applying the k-means algorithm to the QuadPol SAR HV image. Figures 4(a)-(c) show the HV polarization intensity image, initial waterbody mask, and the low-backscatter mask, respectively. By assigning cluster 1 to the water class, it was found that the water bodies defined using the initial water body mask made up about 85% of the total area of the water bodies. In the second step, the window size for the image splitting was set to 100. A total of eight 100 pixel × 100 pixel sub-images were selected. The cyan squares shown in Figure 4(a) indicate the location and coverage of these selected sub-images. Figures 4(d) shows the entropy texture image derived from the QuadPol SAR HV intensity image. The waterbodies are areas of low entropy, whereas most of the land areas correspond to areas of high entropy in the entropy texture image. Figures 4(a) and (d) indicate that the contrast between water and land in the entropy texture image is substantially higher than it is in the intensity image. It can also be seen that, as well as the waterbodies, a few wetland (marsh) areas (e.g. the area highlighted by the red circle) are also areas of low entropy. However, these areas have a bright tone in the intensity image since their main scattering mechanism is the double-bounce return from standing water and grass. Therefore, these wetland areas can be excluded from the histogram computation of the sub-images by using the lowbackscatter mask image. Figure 5 shows the histograms of the selected sub-images: (a) for the QuadPol HV intensity image and (b) for the QuadPol HV entropy image. The histograms have been normalized to 256 grey levels. The intensity and entropy texture images are then transformed so that all values are in the range 0-255. Threshold values of 50 for the HV intensity image and 105 for the entropy image were determined by the valley-enhanced Otsu algorithm. The dashed arrows in Figure 5 indicate the positions of these thresholds in the normalized histograms. The waterbody maps were generated by   setting pixels with grey levels less than the determined threshold values as water. Some wetlands were, however, misclassified as water due to their low-entropy values when the texture method was used. The waterbody map was further refined by using the lowbackscatter maskthis eliminates wetlands, which have high backscatter and a bright tone in the intensity image. Figure 6 shows the final waterbody maps generated by (a) the texture method and (b) the intensity method. It can be seen that the intensity method produces a greater area covered by surface water than the texture method. Error matrices were generated to assess the performance of the two methods using around 8000 samples selected from the reference waterbody map that had been derived from the 5 m-resolution SPOT-5 image.  Table 1 shows the error matrices for the classification results obtained using (a) the texture method and (b) the intensity method. From the error matrices, the kappa indices can be calculated in order to evaluate the accuracy of the distinguishing between water and land classifications. The kappa index measures the agreement between the actual land-cover classes and the classified classes (Congalton 1991). A value of 0 indicates no agreement between the two variables; a value of 1 indicates perfect agreement, with all the values falling on the diagonals. The values of the kappa indices obtained indicate that the texture method (kappa index = 0.89) has a higher accuracy than the intensity method (kappa index = 0.79) for the water/land classifications. For comparison, of 423 land pixels misclassified as water pixels by the intensity method, the texture method reduces the land/water misclassification error to 146 pixels. This misclassification is attributable to the miscalculation of the threshold value in the intensity method. From the histogram in Figure 5(a), it can be seen that the land class in the intensity image presents a flat histogram with a long tail on the right-hand side, which may result in a larger threshold value when the modified Otsu algorithm is used. Out of a total of 2686 water reference pixels, the texture method produces 409 (15.2%) error pixels (water misclassified as land and vice versa), in comparison with the 752 error pixels (28.0%) produced by the intensity method.
The texture method produces a lower classification error for the water pixels than the intensity method. However, about 10% of the error for water pixels that are misclassified as land pixels is still measured. It might be explained by the different acquisition dates of the SPOT-5 and Radarsat-2 data and temporal changes in water level. Another error source may come from the water and land edge pixels. As a result of the 3 × 3 window operation that is part of the computation of the entropy texture, the edge pixels of a waterbody may have higher entropy than the pixels inside, which may result in the water being misclassified as land. Therefore, the classification accuracy could possibly be improved by expanding one pixel of the detected waterbodies, but this would need to be verified by applying a morphological operator in future studies. Table 1 indicates that about 3% of the land pixels are misclassified as water. A probable reason for this is the confusion of water and roads. Roads present the same backscattering mechanism as water, and consequently roads also appear as a dark colour in a SAR image. This error could be further reduced by making use of the existing road network data.
One advantage of the proposed texture method is the usage of the sub-image selection process and the entropy texture. These are critical for the valley-enhanced Otsu algorithm to be able to detect an optimal threshold. Because of the small portion of water in a large SAR image scene, and the unimodal or close to unimodal histogram, the Otsu algorithm may not detect an optimal threshold value for a large SAR image scene. A sub-image selection process is applied to divide the original image into a tiling of smaller sub-images with a user-defined size and then to select sub-images that have a high probability of containing a sufficient proportion of water and land classes. To fulfil the task, an ancillary GIS data set containing the water-land boundary and statistical measures (Martinis, Twele, and Voigt 2009) is often used in the selection of the proper sub-images. However, the water/land boundary GIS data are not always available and the optimal thresholds for statistical measures must be predetermined by analysts. The k-means algorithm can be easily implemented to automatically generate an initial waterbody mask with sufficient accuracy. The initial waterbody mask ensures that the proper sub-images containing only water and land are selected. For simplicity, the kmeans algorithm is applied directly to the original intensity image in this study. When required, the initial waterbody map could be enhanced by a histogram clipping of the intensity image (Amitrano et al. 2014) because the clipping could relax the dynamics of the lowest value scatters and push the highest values towards the right-hand part of the distribution. It is noted that the Otsu algorithm still fails to detect an optimal threshold value for the SAR intensity image since the land class there is a lot of variability in the grey levels for this class ( Figure 5(a)). There is no such problem with the entropy texture, which substantially enhances the contrast between water and land (Figures 4(d) and 5(b)). Therefore, the use of the entropy texture image instead of the intensity image as an input to the Otsu algorithm can improve the results. In addition, higher incidence angle SAR data are preferred for accurately extracting surface waterbodies since the tonal contrast between surface water and land in a SAR image decreases with decreasing incidence angle. Since texture is characterized by the spatial distribution of grey levels in a neighbourhood (Kekre and Gharge 2010) and refers to the local variation in land-cover types, the texture is less sensitive to the incidence angle. Champion et al. (2008) performed a variance analysis to test the influence of the SAR incidence angle on texture values and showed that the SAR image incidence angle does not significantly influence the texture features (Champion et al. 2014), and that only signal intensity and signal-to-noise vary significantly with incidence angle. Thus, the texture method is less sensitive to the incidence angle and overcomes the limitations of intensity thresholding which requires the use of SAR imagery with a high incidence angle. It also overcomes the water detecting limitations of ScanSAR, which achieves a large swath coverage but has large variations in the incidence angles from the near range (~20°) to the far range (~60°). Therefore, the texture method increases the amount of SAR data suitable for surface waterbody detection. It may also provide an automatic method for flood detection and shoreline delineation using SAR data, since all these water-related applications are facilitated by the grey contrast between water and land in SAR imagery (Giustarini et al. 2013;Kuenzer et al. 2013;Lu et al. 2014;Martinis, Twele, and Voigt 2009;Shu, Li, and Gomes 2010). This study presents a fully unsupervised method that generated a satisfactory result. It should be noted that the human interpreters' high-level cognitive function and landscape recognition training should not be overlooked (Amitrano et al. 2014;Datcu and Seidel 2005;Madhok and Landgrebe 2002). The human-machine interaction could help identify and correct misclassifications made with an automatic method. A system implementing this automatic method may provide an interface for a user to adjust (if necessary) the threshold based on the output histogram. Such a system would optimize the strength of machine algorithms in solving very complicated mathematical tasks with the humans' capacity for discriminating and interpreting images to facilitate recursive operations with a desired accuracy.

Conclusions
This study presents an automatic texture thresholding method for waterbody detection using large SAR image scenes. It combines the k-means clustering algorithm, sub-image selection process, and a modified Otsu thresholding algorithm. Using k-means clustering, sub-images that contain sufficient proportions of water and land classes are selected. It is critical for the Otsu algorithm to detect an optimal threshold for a SAR image. In contrast to the intensity image, the entropy texture image can enhance the contrast between water and land and reduce the variability in the grey level for the land class. It also makes the method less sensitive to incidence angles and thus increases the amount of SAR data that is useful (e.g. low incidence angle SAR images and ScanSAR images) for detecting surface waterbodies. An improved result is obtained when the Otsu algorithm is applied to the entropy texture image. The proposed approach was tested on a Radasat-2 QuadPol mosaic scene covering the Spritiwood valley in Manitoba, Canada. The results showed higher classification accuracy and lower classification error than the commonly used intensity method for overall land/water classification. The results presented in the article indicate that the waterbodies in a large SAR image can be delineated from land automatically by the proposed texture method. The method can be applied to a single geocoded SAR image, with no additional data (other than DEM data used for geoorthorectification) required. In addition, all algorithms used in the texture method provide fast execution and require minimal computational effort. These advantages mean that it can be easily implemented and used for extracting the waterbodies from a SAR image in near-real time.