1. Introduction
The National Hydrological Service (NHS) of Environment and Climate Change Canada is the national agency responsible for the collection, interpretation and dissemination of standardized data on water resources and information in Canada. Currently NHS’s primary mechanism to collect and process data is by means of approximately 2100 active gauging stations that are employed to measure lake level or river stage used to calculate water discharge. This system has been in place for over 100 years, and although valuable, provides only a limited picture of water resources within Canada. As a part of the Government of Canada’s strategy to leverage earth observation technology to further broaden the amount and value of environmental data and products [
1], the NHS is working with the Canadian Space Agency on the development of technology and methods based on satellite imagery for water monitoring [
2].
Water monitoring using Synthetic Aperture Radar (SAR) has been the object of study for many researchers [
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15,
16,
17,
18,
19]. Due to its all-weather capabilities, and its image acquisition capacity during day or night or in cloudy conditions, SAR imagery offers better alternatives for water mapping than optical imagery [
5]. As a result of the radar’s unique response to water, water mapping using intensity thresholding methods on SAR image had been extensively used [
4,
5,
6,
7,
20,
21]. In this method, water is separated from land in intensity images, but the accuracy of the results relies on the ability to differentiate land
vs. water pixels in the intensity domain, which becomes especially difficult when the water backscatter is affected by wind-induced roughness. Since the values of backscatter vary depending on the incidence angle, image quality and wind-induced roughness, the threshold needs to be modified on a scene per scene basis [
4,
5,
6]. In some other instances, image thresholding is combined with texture information [
7] or region-growing segmentation algorithms used to determine the water pixels [
14,
16,
21,
22]. Segmentation algorithms are quite successful but computationally expensive. In many cases, statistical information is combined with digital elevation models [
8,
14,
16], where pixels likely to contain water are first determined based on topographic data and the probability of water is based on histograms of water
vs. land pixels in the image, but this requires a coarse water mask to determine the statistics of water pixels. This dependence on a pre-determined water mask and high-resolution digital elevation models makes these methods unusable for mapping small ephemeral water bodies. More recently, active contour algorithms, also called snakes [
23], have been applied with some success for water mapping from radar images [
17,
22]. These algorithms, however, rely on ancillary data to determine candidate pixels for water as well as on morphological operators, which results in longer processing times.
This paper presents a novel threshold-based algorithm for automated water mapping based on Radarsat-2, which does not require ancillary data or
a priori knowledge of the response of water bodies in each scene. It has been developed in such way that it will be transferable to the Radarsat Constellation Mission (RCM) upon launch. The RCM will consist of three satellites that will provide daily worldwide coverage in lower resolution modes. The constellation will have a revisit period of 4 days, nominal resolution going from 3 m in very high resolution mode to 100 m in ScanSAR mode, and swaths ranging from 20 to 500 km. The system will provide imagery in single, dual and compact polarimetry modes and is expected to deliver 300,000 scenes per year [
24]. RCM is owned and will be operated by the Government of Canada and will provide Earth Observation radar data continuity for Canada to fulfill the government priorities in maritime surveillance, disaster and natural resources management [
25]. A combination of a thresholding technique with a texture indicator has been applied in a fully automated way on 82 Radarsat-2 images over the Canadian Prairies, to produce polygons delineating ephemeral water bodies (potholes) as well as permanent open water, which will be integrated into the Canadian Land Data Assimilation System (CaLDAS) [
26], to improve numerical weather prediction [
27]. The results were evaluated using nearly synchronous high resolution optical imagery for three areas of interest in the study site.
4. Results
Radarsat-2 derived polygons were evaluated against water polygons derived from cloud-free, temporally coincident high-resolution optical imagery over three areas of interest (AOIs) in Alberta, selected based on the availability of the imagery and landscape characteristics.
Due to the difficulty of accessing potholes’ areas and their ephemeral nature, concurrent optical imagery with sub-meter spatial resolution was used whenever available. If this type of imagery was not available, then optical imagery with a pixel size comparable to the SAR image was employed. Although the evaluation against imagery instead of ground surveys is recognized as being limited to some degree, this was considered an appropriate proxy for an expensive field campaign that did not align well with the priorities of the Radarsat-2 acquisition plan.
In order to assess the veracity of the method to map surface water with different conditions, the three pairs of concurrent optical-SAR image were chosen because they represent a variety of conditions in which non-forested surface water areas can be found: open water, water with some vegetation and areas of very small (<1 ha) closed drainage flows, commonly found in the agricultural landscape of the Canadian Prairie region. The radar/optical image pairs were:
Pair 1: WV2 and F0W3 (0.5 m and 8 m resolution respectively). Surface water polygons were derived from a cloud-free orthorectified and pan-sharpened WorldView2 (WV2) image taken on 12 May 2014 over the east side of the city of Red Deer and with a coverage 116 km long by 20 km wide. The polygons delineating open water were produced by thresholding of the near and short-wave infrared bands and manually edited. The resulting polygons were compared against the ones derived from a Wide Fine (F0W3) scene taken 1 day before (i.e., 11 May 2014) and fully containing the WorldView-2 scene. The area covered by the optical image is cropland, which hydrological features are rivers, potholes and many shallow drainage flows that can be perceived mainly in sub-meter optical imagery, but are still considered surface water.
Pair 2: WV2 and U76: (0.5 m and 2 m resolution respectively). Another WorldView-2 image taken on 10 May 2015 over an area 10 km west of the city of Schuler, Alberta was employed to derive water polygons using the same thresholding and manual editing procedure described above. The resulting polygons were compared against water polygons derived from an Ultra-Fine image (U76) taken 6 days after. The landscape in this area is mainly characterized by flooded vegetation.
Pair3: RapidEye and FQ19 (5 m and 7 m resolution respectively). A RapidEye (RE) Image from 8 September 2012 over Elk Island National Pak in Alberta was employed to derive water polygons using thresholding and manual editing. The resulting polygons were compared against the water polygons obtained from a fully overlapping Fine Quad image (FQ19) taken the same day. Surface water in this area is mainly composed of open water bodies larger than 1 ha.
The minimum area, maximum area, number of shapes and cumulative area of the resulting water polygons from the optical and SAR images were derived. To facilitate the comparison, polygons bigger than 25 m
2 were selected, their statistics were computed using 5 area intervals (represented by the rows in the tables below) and the number of polygons as well as the contribution (expressed as percentage) of their cumulative area to the total area was calculated for each area range. Polygons smaller than 25 m
2 were excluded as they are not perceived by Radarsat-2 operational beam modes with repeated pass and continuous coverage. The comparison of surface water delineation based on area intervals was motivated because missing a high number of small polygons would not be detrimental for water quantification, as long as, together, they do not hold the biggest percentage of surface water in a particular area. Having an area for which missed small water polygons hold the majority of surface water will render the method or the data inappropriate for mapping ephemeral water bodies.
Table 1,
Table 2 and
Table 3 show the accuracy in area quantification obtained from polygons, when surface water is mapped by manual vector editing on the optical image
vs. applying our operational procedure on the SAR images:
The accuracy achieved in quantifying the area of water bodies, when compared to high-resolution optical imagery for different beam modes, showed that:
Our procedure fails to map open water bodies smaller than 1 ha when applied to Wide Fine mode. For the first paired optical-SAR dataset evaluated (
Table 1), small water bodies were mostly missed by the algorithm. Also, the cumulative area contained in polygons smaller than 2 ha contributed to 52% of total surface water area in this particular AOI, which explains why the accuracy of polygons extracted from the SAR image drops significantly for this beam mode—see
Figure 6.
The quantification of the area on large water features were missed from Wide Fine mode due to fragmentation. This occurs due to discontinuity of polygons delineating water edges, when vegetation patches occur at the edges (e.g., riparian forest)—see
Figure 7. On the other hand, two big water polygons that are seen as separate entities in the RapidEye image can be joined together and form one in the SAR image, due to the missing separation by small low vegetation patches (which are visible on the RapidEye image)—this changed the distribution of the area contribution for the water polygons derived from Fine Quad between 2 ha and 0.5 km
2 (
Table 3).
As expected, the algorithm also fails to detect flooded vegetation: a closer look at the polygons missed from the Ultra-Fine image
vs. the polygons obtained from WolrdView 2 of the same time period showed that many polygons were larger than 1 ha and could have been seen in the Ultra-Fine image due its fine resolution, but were missed because of their high backscatter, which is characteristic of vegetation—see
Figure 8 and
Table 2.
Fine mode imagery seems to provide the best results, as it quantified 88% of the total surface water area and picked up 97% of the total number of polygons larger than 1 ha when compared to polygons obtained from RapidEye (
Table 3).
Some water polygons that are selected by the algorithm from the SAR image are not shown as discernible open water bodies in the optical image (especially when the optical image pixels size is 5 m or more). They could be seen as false positives, but considering that SAR is more sensitive than optical imagery to water content, these areas could also be areas of low vegetation (where the vegetation cover is not high enough to influence backscatter but its chlorophyll content does influence reflectance) with particularly high water content—higher than its surroundings.
In addition to this comparison of area estimates, randomly stratified points (50% of them on water pixels and 50% of them on non-water pixels) were generated and error matrices as well as kappa statistics [
44] were computed using the same SAR-optical image pairs described above (see
Table 4 and
Table 5, respectively). The accuracy of the reference optical images to map water was evaluated through visual inspection, and their kappa values varied between 91.6% and 98.8%.
Considering the uniform spatial distribution, number and representation of each class in the sample points, we can state that the confidence interval provides a precise approximation of the estimated overall accuracy. The kappa statistic, as a measure of how closely the instances classified by the algorithm matched the reference data (compared to a random classifier) confirms the results presented in
Table 1 through
Table 3—despite the likelihood of large water bodies attracting more random points. Our method has low accuracy to map water from Wide Fine beam mode over areas with a high number of small potholes (less than 1 ha) or over flooded vegetation (
i.e., the dominant landscape in UF76). However, the method succeeds in automatically mapping non-vegetated water bodies larger than 1 ha from Fine Quad beam modes when compared to optical imagery of equivalent spatial resolution.
5. Discussion
The areal extent of water bodies is important for a variety of reasons including a better understanding of the land–atmosphere boundary for meteorological and climate modeling [
27]. The size and extent of water bodies across Canada vary both geographically and temporally throughout the year. This variation depends significantly on local climate and topography: the prairie region has many shallow water depressions (prairie potholes) that are filled through snow redistribution, snow melt, infiltration, and precipitation, and can evaporate quickly and exhibit “spill and fill” movement [
45]. Due to local topography, small amounts of contributed water influence the areal extent of water bodies significantly. On the other hand, in the Canadian Shield, and mountainous regions, a large influx of water may significantly raise the stage, but areal extent does not significantly change.
The method presented in this paper includes procedures for geometric correction, calibration, filtering and thresholding of image values and derived texture. For surface water detection, a threshold-based algorithm that requires dual polarization is less restrictive than polarimetry-based algorithms that rely on fully polarimetric data, and therefore is easier to implement operationally, as more data can be obtained from scenes with larger coverage than fully polarimetric scenes.
Having accurate orbit and satellite ephemeris information as well as ancillary information (e.g., RPCs and GCPs) improves the efficiency of any operational system designed to retrieve data from satellite imagery: processing times were significantly lower when the system could skip the point collection routines, and eliminate manual quality control process on point collection. The nominal geometric model provided with Radarsat-2 was sufficient to perform orthorectification for the Prairie region. For the Radarsat Constellation Mission, we expect that the orbit accuracy as well as the rational polynomial coefficients are going to have accuracy that is at the same level or better than Radarsat-2, so that processing times remain as they are. Even for mountainous areas, where ground control point collection is required, the accuracy of the orbit and RPCs can significantly affect the processing time and quality of the results, as the image-matching algorithm uses this information to define the initial position of a searching window and to define a search radius [
35].
A threshold-based technique offers an efficient and non-expensive way to do automatic mapping, but it has limitations. Work done by Manjusree
et al. [
20] proved that a water mapping threshold needs to be decreased from –8 to –12 decibels in the near range and from –15 to –24 in the far range. In our experiments, using a single threshold resulted in missing water bodies in the far range, because there were no seeds captured: the threshold value excluded more pixels on open water in the far range, which usually has lower signal to noise ratios and higher calibration uncertainty. Work by O’Grady
et al. [
21,
46] also concluded that there are variations in the backscatter based on the incidence angle, given the varying noise levels of the scene. Another drawback of threshold-based methods for water mapping from SAR is that the threshold used for water mapping is severely affected by the noise (and noise floor) of the scene. Hardcoded threshold values are always “subject to change”—based on image quality and beam mode specifications. Wide beam modes with fine resolution (such as Wide Fine and Wide Ultra-Fine) are achieved by increasing the Noise Equivalent Sigma Zero (NESZ) which regulates the minimum signal that a SAR can measure, in order to maintain a constant resolution [
36,
40]. Also, to reduce data volumes when wide swaths modes with high resolution are acquired, a 2-bit block adaptive quantization (BAQ) compression technique is used [
40]. Predominantly, the 2-bit BAQ does not provide enough detail to capture the difference in signal strength (and therefore compensate for it) when there is a large change in the radar return as a function of time. This results in artificial backscatter variability from the near to the far edge, and increased noise levels, especially in darker areas (calm water) that are close to bright features (vegetation), hence the necessity to vary the threshold used for water based on the incidence angle.
The traditional method of accuracy assessment based on a confusion matrix and random points (or stratified random points) did not adequately evaluate the success (or failure) of surface water mapping for different beam modes since an indication of the accuracy on obtaining the complete area of surface water is required for numerical weather prediction models [
47]. Also, large water bodies that are easily identified in any SAR beam mode, have more pixels and therefore are more likely to get more random points, hence introducing some bias on the results of traditional accuracy assessment procedures. Furthermore, the elaboration of the confusion matrix requires a direct comparison of all points, without considering their contribution to the total surface water area. Quantifying water area based on polygons of a defined range (
i.e., 1 to 2 ha) better described the capacity of different SAR imaging modes for water mapping, as well the effectiveness of the method given the area distribution of ephemeral and permanent water bodies for the AOIs in our study region.
While omission errors in the SAR-derived polygons are relatively easy to evaluate, false positives (some of which could be wet surface or just homogeneous areas of low return, unrelated to water content) are much harder to assess. If we define the mapping unit as non-vegetated potholes that have a characteristic shape (normally high compactness), the commission errors can be assessed based on high-resolution optical image. However, if we were to define surface water as any standing water that might or might not contain any type of vegetation, the false positives are very difficult (if not impossible) to discriminate from high-resolution optical imagery, specially without a short-wave infrared band. Water areas with some vegetation are as important as open water for numerical weather prediction [
26]. Nonetheless, doing fieldwork in these areas might neither be possible (due to accessibility) nor relevant if the field work campaign is not executed at the exact same time of the image acquisition, as water can evaporate very quickly, especially during the summer season. Radar is more sensitive to water content than optical imagery, and therefore it can detect water that cannot be easily discriminated in the visible or NIR part of the spectrum [
48]. Furthermore, the chlorophyll activity and leaf structure can also affect the reflectance in optical imagery whereas for SAR these factors play a secondary role on the backscatter registered by the sensor. More work is underway to gather concurrent UAV and LiDAR data along with more SAR scenes, to obtain relative accurate ground truth that will serve to evaluate the commission errors. However, the validation of dynamic targets for multi-temporal, multi-source data remains a big challenge under a limited budget.
Further work is required to improve accuracy for small potholes and flooded vegetation. The upcoming Radarsat Constellation mission will offer compact polarimetry [
24], which means that polarimetric data can be easily obtained in wider areas with more resolution than its predecessor, Radarsat-2. The acquisition of polarimetric data in the standard coverage that will be available with RCM, means that imagery will be acquired on a periodical basis, without the overhead of image acquisition planning and consolidation of conflicting acquisition requests from different users. If this type of imagery can be provided in fine resolution modes for inland applications, such as forestry, agriculture, moisture analysis, biomass estimation, flooding and landslides, it will open many new opportunities not only for more accurate and frequent water monitoring systems but also for better natural resource mapping, in general.
6. Conclusions
A threshold-based procedure to automatically extract surface water polygons from SAR imagery has been presented here. For the AOI covered by the Wide Fine scene, our analysis based on WorldView-2 imagery at 50 cm pixel size showed that 52% of the total water area was contributed by polygons smaller than 2 ha, whereas in the AOI covered by the Fine Quad scene, most of the potholes and permanent water bodies (87%) were larger than 1 ha. For areas where small potholes are dominant, mapping with Wide Fine can only reach 21% accuracy when estimating total surface water area, but for areas where potholes larger than 1 ha are dominant, Radarsat-2 Fine mode captured 88% of the total water area extracted from RapidEye. This provides much better information than currently available for regions where seasonal and ephemeral changes are expected on surface water. More efforts should be made to improve the performance and meet the requirements for operational applications on smaller water bodies and low-resolution modes.
Despite its limitations, our procedure is a fully automated method to derive polygons of open water from Radarsat-2 Fine beam modes, based on pixel values, which is, in terms of processing time, much faster than running more complex segmentation procedures. Processing time becomes a key factor when the information on the land surface needs to be derived to feed the models used in weather forecasting, which run as frequently as every 6 h [
47]. The procedure has been implemented in Python as a processing chain (see WaterExtents.py, provided as
Supplementary Material to this paper), and it uses image processing software with multi-threat algorithms from PCI Geomatics, where only two steps for quality control and analysis are required: first, when a collection of points is required for geometric correction, and at the end, when water polygons are generated. The topological intersection of polygons to generate the final water mask is also being integrated into the same script using Python-based algorithms and open-source GIS software, which will make it easier to integrate into operational and/or web-based systems.