Automated global water mapping based on wide-swath orbital synthetic-aperture radar

This paper presents an automated technique which ingests orbital syntheticaperture radar (SAR) imagery and outputs surface water maps in near real time and on a global scale. The service anticipates future open data dissemination of water extent information using the European Space Agency’s Sentinel-1 data. The classification methods used are innovative and practical and automatically calibrated to local conditions per 1× 1 tile. For each tile, a probability distribution function in the range between being covered with water or being dry is established based on a long-term SAR training dataset. These probability distributions are conditional on the backscatter and the incidence angle. In classification mode, the probability of water coverage per pixel of 1 km× 1 km is calculated with the input of the current backscatter – incidence angle combination. The overlap between the probability distributions of a pixel being wet or dry is used as a proxy for the quality of our classification. The service has multiple uses, e.g. for water body dynamics in times of drought or for urgent inundation extent determination during floods. The service generates data systematically: it is not an on-demand service activated only for emergency response, but instead is always up-to-date and available. We validate its use in flood situations using Envisat ASAR information during the 2011 Thailand floods and the Pakistan 2010 floods and perform a first merge with a NASA near real time water product based on MODIS optical satellite imagery. This merge shows good agreement between these independent satellite-based water products.


Introduction
The consequences of inland and coastal flooding can be devastating and flooding needs to be detected and mapped as accurately and quickly as possible, so that appropriate measures can be taken by governments or disaster management agencies, pre-warnings may be issued, and downstream forecasts may be initiated (Carsell et al., 2004;Werner et al., 2005).In situ networks of hydrological gauges are increasingly being complemented by satellite imagery, which plays an important role in the European Global Monitoring of Environment and Security (GMES; Brachet, 2004) Emergency Response Core Service.That service is meant to provide "Rapid Mapping": fast retrieval of information from satellite imagery in order to map consequences related to hazards and civil protection.
Fast retrieval and systematic retrieval are different terms.Thus, a number of commercial and non-commercial agencies can respond to flood disasters within a short amount of time (fast retrieval).However, these agencies react on demand, when an emergency response has already begun.Also, due to the required manual expertise and labour requirements, such response cannot be accomplished on a daily basis and, commonly, not within a processing time comparable to the capture time of satellite images (Hostache et al., 2012)."Systematic" water mapping can instead be developed; wherein water extent information is routinely provided through the consistent and automated generation of maps and associated GIS (geographic information system) data.In surface water Published by Copernicus Publications on behalf of the European Geosciences Union.
mapping, these maps can then be used within a GMES Service for different purposes such as flood status, environmental monitoring of lake and reservoir extents or initializing hydrodynamic models.In this article we focus on the development and use of such an automated system specifically for use in flood response.

Relative advantages of SAR and optical imaging
The quality of C-band (e.g.Envisat) synthetic-aperture radar (SAR) images is independent of the time of the day and cloud cover.Water can often be visually distinguished due to the low backscattering exhibited by relatively flat water surfaces (with very low return to the side-looking sensor due to specular reflection, as compared to brighter "backscatter" from rougher surfaces).In contrast, and although not capable of observing through clouds, the MODIS optical sensor on NASA's Terra and Aqua satellites has some important advantages: MODIS visible and near IR (NIR) bands 1 and 2 provide global, twice daily coverage at 250 m spatial resolution, and optical multispectral classification methods may better distinguish land and water in some areas, including in deserts where SAR backscatter may be very low and highly variable (a.o.Ridley et al., 1996;Raghavswamy et al., 2008).The utility of MODIS for flood-related work has been repeatedly demonstrated by maps disseminated from the Dartmouth Flood Observatory (http://floodobservatory.colorado.edu/).These water area products are usefully compared to numerical 2-D model output in the case of catastrophic storm surges (Brakenridge et al., 2012).Improvements in wideswath SAR data processing can be undertaken to the same end as the addition of all-weather, day-night imaging capability to hydraulic models provides a powerful approach (Schumann et al., 2009).
The anticipated data output from ESA's Sentinel-1 satellites (Attema et al., 2009) will further open opportunities.However, at present, semi-automatic classical water extraction techniques, such as thresholding or change detection applied on SAR images, may fail due to windy conditions, or partially submerged vegetation, resulting in higher backscattering values (Yesou et al., 2007;Prathumchai and Samarakoon, 2005).According to Silander et al. (2006) misclassifications may also be caused by the dependency of backscattering on detection angle.Mason et al. (2010) mention the problem of misclassification due to topography, vegetation or canopy.O' Grady et al. (2011) conclude that misclassification due to low backscatter values from non-flooded areas can be reduced via image differencing approaches.Matgen et al. (2011) and Giustarini et al. (2012) present a method relying on the calibration of a statistical distribution of "open water" backscatter values inferred from SAR images of floods.Given the many circumstances that can affect classification results, it is difficult to derive a consistent classification technique that, ideally, also includes an error or accuracy assess-ment, and for all incidence angles.Up to the present, for example, some manual interpretation is still normally required to translate SAR data into water maps.However, Hostache et al. (2012) research an automated way of selecting the best reference image for change detection.

Need for automated data processing and map generation
Automation is required for any systematic mapping approach to avoid subjective, time-consuming and expensive manual interpretation for each flood.Consider the case of a single ESA Sentinel-1 satellite: (1) the expected amount of data for only Level 0 data across all acquisition regions will reach 320 TB per annum, amounting to 2.3 PB (petabytes) in the course of 7.5 yr (Snoeij et al., 2009;Attema et al., 2008); and (2) when processed further, Hornacek et al. (2012) expects the matching Level 1 data volumes for baseline soil moisture products to be 4 to 5 times larger than those for Level 0. In order to cope with these amounts of data, the need for automation, and to fully utilise the very high information content of these new sensor data streams, new techniques are needed.

Methodology
We present a prototype automated technique, embedded in an online service, which classifies SAR imagery to probability of water for each image pixel, in near real time (NRT) and at global scale, which we call "Global Flood Observatory" (GFO).The service used Envisat ASAR data while that sensor was operating (it failed 8 April 2012, after 10 yr), which was made available in Level 1 format by the European Space Agency (ESA) in near real time from a 15-day rolling archive.The data were processed in near real time, i.e. within 3 h (but usually faster), after the data had been put on the ESA NRT Rolling Archive.Output results were subsequently placed on an open data server in widely-used data formats (i.e.NetCDF and Google Earth KML files).Previous employment of relatively high resolution (narrow swath) SAR for flood classification includes Kasischke et al. (1997) and McCandless and Jackson (2004).Three attributes of SAR that are of importance for the present algorithm are total backscatter, incidence angle, and signal polarisation.
In this paper we have chosen the axes in some of our plots as the amplitude of the backscatter measured by the Envisat-ASAR sensor, being the direct output of the ESA NEST (Next ESA SAR Toolbox) software used.Backscatter is the portion of the outgoing satellite radar signal -usually looking sideways in different incidence angles (as shown in Fig. 1) -that the target redirects back towards the radar receiver antenna.If the target is horizontal, the backscatter is a measure of the electromagnetic roughness of the first very thin layer of the subsurface (a.o.Verhoest et al., 2008).As already widely known from ground penetrating radar and other microwave techniques, the electromagnetic roughness, creating a "subsurface microtopography" is depending on the physical contrasts between the conductivity and permittivity within this layer, causing a reflection coefficient: a measure of the reflective strength of a radar target.Usually for the solid Earth, this contrast is caused by differences in soil moisture, whereas differences in soil type within this thin layer play a minor role (Beres and Haeni, 1991).
Because the beamed radar is to the side of the sensor, an incidence angle applies, and a lesser amount of total energy returns to be recorded by the sensor.Such radar backscatter is dependent on both the incidence angle α and on land cover, land topography, and soil moisture.Incidence angles in operating SAR sensors commonly range between between 15 • (closest to the satellite) and 45 • (furthest from the satellite), as shown in Fig. 1.
Our algorithm calculates the probability of a pixel within a satellite imaging swath being water, by matching its backscatter signal to a probability distribution of the pixel being dry, or being wet.These probability distributions are conditioned on geographic location, incidence angle and polarisation of the signal and were established using a training dataset of three years of Envisat ASAR data (Global Mode -GM, Wide Swath Mode -WSM, Image Mode -IMM, and Alternating Polarisation Mode -APM).The probability distribution is distributed over each 1 × 1 • latitudelongitude tiled dataset of the land covered globe.In general, for most incidence angles, the backscatter-incidence angle (σ − α) pair for land targets is different than that for water.In practice, an empirical distribution function is estimated per geographical area by building 2-D histograms of σ − α pairs for (a) pixels across the entire globe, which are permanently wet; and (b) pixels within a 1 × 1 • tile which are dry under average climatological circumstances.An example of a trained histogram is shown in Fig. 2, where for the Netherlands thousands of σ − α WSM pairs have been gathered, gaussian-smoothed and plotted for land and water in

Training method details
As noted, a training period is first used to derive a spatially distributed probabilistic model for distinguishing land and water.Then the application of this model in near real-time is accomplished.
First it should be noted that the local incidence angle α as used here does not take into account local topographic features.Next, for building the histogram training set, an ancillary dataset is used, called the "water mask".The water mask is derived from the NASA SRTM Water Body map (SWBD), documented by USGS (2005) to provide water and For each 1 × 1 • tile (in latitude-longitude), a histogram is made for land for discrete values of the backscatter and α, for each polarisation.For the freshwater dataset, one global histogram file is made.The resulting multidimensional trained histograms consist of one or more hydrological years of SAR data.The trained land and water histograms are used as the reference set for classifying newly downloaded SAR data to land or water pixels.The training sets are built as separate entities for each SAR mode (e.g.ASAR-GM, ASAR-WSM, ASAR-IMM, ASAR-APM).It should be mentioned that by doing this, some "noise" is created, since flood events that occur while building the training set are not filtered out.

Classification and quality assessment details
Because of the difference found in the σ − α pair 2-D histograms (for dry land and water), a distinction can be made between dry land and water.This is shown in a visual example in Fig. 3.The probability that a pixel in a SAR dataset is wet or dry is established using Bayes' law, also used by e.g.Frigessi and Stande (1994) for satellite image classification.We consider two empirical distribution functions for wet and dry pixels as posterior distributions.The procedure to establish a single probability of a pixel being wet is given below.All equations are written as if continuous probability distributions are used and are applicable on a limited area within the earth's surface for which the set of empirical distributions for land/water apply.
Bayes' law in its general form can be written as where P [M|D] represents the probability of a model M given demonstrative data D. P [D|M] is the probability of data D occurring when model M applies and P [M] is the prior distribution.In our case, we may write this as where s = w means that a pixel s is classified as water (w) and b represents a certain σ − α pair.P [b|s = w] is the probability that a certain σ − α combination is experienced when a pixel is classified as water.This probability distribution is approximated empirically based on discrete slices of the trained 2-D histograms per discrete incidence angle value as described in Sects.4 and 5. P [s = w] represents prior knowledge that the pixel within the SAR scene is water.Since we have no prior knowledge about this, and a pixel can only have two states (land or water), this probability is set on 0.5.One may argue that the prior knowledge of a pixel being wet or dry should be equal to the amount of permanent water pixels within a limited area for which the training is established, divided by the total number of pixels within the area.However, during the recording of the analysed scene, we should not skew the probability towards such a prior, since the goal of this exercise is to distinguish wet events over pixels that are in fact normally dry.For example, let us assume that we have established a training probability distribution for a desert area that has no permanent wet pixels.If we would now analyse a scene during a catastrophic flash flood event, we would estimate the probability of a pixel being wet zero everywhere.Finally, P [b] is the normalisation constant.The same equation can be established for the probability that a pixel should be classified as dry land, being where s = d means that a pixel s is classified as dry/land.Equations ( 2) and (3) both share the same denominator.Furthermore, both priors have the same value being 0.5.Therefore, the following applies: and Furthermore, it is known that the sum of probabilities of a pixel being dry land or water is equal to unity, given that dry land or water are the only two states possible: Substituting Eq. ( 6) in Eqs. ( 4) and ( 5) gives Substituting Eq. ( 7) in Eq. ( 4) gives Equation ( 8) is used to determine the probability that a pixel is water.Finally, knowing the empirical probability distributions for dry land and water, we define a quality indicator q at a certain latitude, longitude and polarisation, as defined in Eq. ( 9): in which the shared area of the normal distributions of the land and water probabilities P (b|s = d) and P (b|s = w) are a measure for the quality of the probability calculation.For example, after reading in a NRT SAR pixel we already know what the two probability distributions of the training set are for the polarisation, angle and backscatter of the SAR data in this pixel, knowing the 2-D histograms of the 1 × 1 • tile.When these two probability distribitions overlap completely -which could for example happen at some low and intermediate local incidence angles and in very dry areas like deserts q will be close to 0 %.If the two distributions are separated completely q will be 100 %.The indicator is dependent on backscatter, incidence angle and geographical location (latitude-longitude).It can be used as a post-processing tool to filter out data that are already pre-defined as inferior for calculating water probability by defining threshold values (e.g. to only show probabilities where the quality indicator is higher than 70 %).

Creating a topography mask
With proper correction for topography, SAR classification methods can be improved in mountainous areas (e.g.van Zyl et al., 1993).Resulting errors in this topography correction will depend on the spatial resolution and quality of this topography data.Instead of correcting for topography, however, and because our concern is surface water, we have chosen to improve efficiency of the automated NRT calculation methods by using a pre-processing filter or mask, prior to classification to water probability.By using threshold values of the height above nearest drainage (HAND) index (Rennó et al., 2008), areas that are unlikely for long-term flooding are filtered out.The HAND index is calculated by expressing the relative height of a location to its drainage outlet in an associated channel.It has now been calculated globally based on the HYDRO1k dataset, developed at the US Geological Survey ( 2008) and based on GTOPO30, a global Digital Elevation Model (DEM) at 30 arc second (approximately 1 km at the equator) resolution (Gesch et al., 1999).An example of the HAND index for Thailand is shown in Fig. 4.

Preliminary results
The downloaded SAR files are temporarily stored in digital (NetCDF) grid files.The trained histograms for each 1 × 1 • (in latitude-longitude) tile have been generated for land, for discrete values of the backscatter and the local incidence angle, and for each polarisation.Figure 5 shows a global compilation for ASAR GM data and three examples of discrete histogram training sets: rainforest, desert, and freshwater.Different areas in the world show different backscatter characteristics.For example, deserts (as well as savannah regions) have a low backscatter, rainforest generally has a rather constant backscatter value, whereas ice (not shown in the figure) exhibits a very high backscatter.This regional dependency is the reason that the probability distributions of land are stored per 1-by-1 • latitude-longitude tile.
Water probabilities and quality indicators are calculated and stored together in (daily) folders containing all calculated water probabilities, and designed to be made publicly available.The water probabilities are made available in gridded (NetCDF) files of 10 × 10 • tiles.An example of an output in a resolution of 0.009 × 0.009 latitude-longitude degrees (equivalent to app. 1 × 1 km on the equator) is shown in Fig. 6, where areas are shown in which the water probability is set to 0 and the quality indicator to 100, regardless of satellite data availability.This is the result of the incorporation of the HAND filter, thereby automatically setting water probabilities as a pre-processing step before classification.9 Validation of the method in two case studies

Bangkok floods, Thailand, 2011
The region along the Chao Phraya River, north of Bangkok, Thailand, suffered from severe flooding in the fall of 2011 caused by heavy rains in upstream areas of the catchment.The flooded area is rather flat, but surrounded by hills.The Envisat satellite covered the flood propagation over de Chao Phraya on 13 October 2011 with the WSM mode.The flooded area on the WSM image is covering the incidence angles α from 24 to 49 • .The resulting water probability map generated by our algorithms is shown in a Google Earth flood map in Fig. 7 (left panel).Examination of the flood map in more detail leads to some interesting observations: -The detected water extent on 13 October 2011 by our algorithm, when the threshold is applied with a conservative 70 % probability and 70 % quality and superimposed on the HAND image in Fig. 7 (right panel), shows a strong visual correlation between flooded areas and low HAND index values between 0 and 1.
-Elevated features, such as roads, embankments and railways can be distinguished in our image.These objects constrain the flood water and are marked by a sharp boundary between pixels detected as dry land, and pixels detected as water.
-Just East of Bangkok, there are many flooded rice paddies.These paddies and their borders are visible in our image as square like features.
When applying a conservative low-pass HAND-filter from 0 to 5, i.e. all values higher than 5 are filtered out, we continue the validation by looking at the use of different threshold values of our probabilities and qualities for computing a binary water extent map. Figure 8 shows the raw backscatter map of 13 October image (top left panel).The city of Bangkok is recognised in the white zone showing high-amplitude direct reflections and double bounced backscatter and reflections.The validation map (Fig. 8, top right panel) has been created using an algorithm that follows a similar procedure to that of Rémi and Hervé (2007) based on backscatter ratios of the flooded and non-flooded image but was modified to include an increased weighting on low backscatter areas where the ratio was high but the backscatter low because of surface variability.This was applied on a pixel by pixel basis and areas in both the reference stack and the flooded image were classified as permanent water.The GFO and validation datasets are then compared using a simple algorithm to compare a binarised version of the GFO dataset at a series of thresholds for quality and probability (as required) based on a pixel by pixel comparison to the validation result.The output of this assessment is a series of statistics which identify how much of the flood area is identified (%) on the GFO as well as providing an agreement factor (0-1) to account for total area in agreement of the result (flooded and non flooded area).
The ideal output would be a flood area identified at 100 % and an agreement factor of 1.The results of our GFO probabilities thresholds are applied at 60 % (Fig. 8, middle left panel) and 75 % (Fig. 8, middle right panel).The areas with high quality indicator values are indicated in red (Fig. 8, bottom left panel: 50 % and bottom right panel: 60 %).Table 1 shows the proportion of the flood identified with the validation map and the proportion of the flooded area in agreement between the GFO and the validation map.Although areas have been flooded on the eastern side of Bangkok, the lower q values indicate that the GFO results are less reliable in that area.Places where the HAND index are higher than 15 are set to q = 100, water probability = 0 in a pre-processing phase.
The Dartmouth flood observatory (DFO) heavily utilises the two MODIS sensors aboard the NASA Terra and Aqua satellites.Currently, a team at NASA is also assisting this effort by performing the classification procedure in an automated way (their NRT Flood product, as described in Brakenridge et al., 2012).In this automated process, the NRT processor collects and combines 4 images over each 10 × 10 • latitude-longitude subset, and over a forward running period of 2 days (thus, two images/day worldwide; four images/two days).The resulting GIS file shows surface water as boundary polygons: each such "daily" file actually includes two days of imagery, using a MODIS band 1/band 2 threshold approach to detect water and requiring at least two detections per pixel in order to exclude cloud shadows (which have similar spectral characteristics to water, but which change location over time).Because the DFO approach suffers from cloud cover and the present SAR-based approach provides less frequent temporal coverage, it is desirable to merge the two independent approaches for flood mapping.A first attempt to merge the two products was accomplished using the Envisat-ASAR WSM data collected on 13 October 2011, and the DFO MODIS based data from 13 October 2011.Because DFO provides a binary map (flooded, or non flooded) and our GFO produces a probability map of flooding, along with a quality indicator for this probability, it was decided to establish a binary map from the GFO product as well.This was done by thresholding the GFO probabilities in the area north of Bangkok, where probabilities and q values are high, to high thresholds of 70 % probability and 70 % q to compare the results to optical data.The results are shown in Fig. 9.The image shows the visual correlation between the two independent mapping methods.In particular the far right image shows that there are flood areas located by the SAR processor, clearly following elevated features in the landscape such as roads and railways, which were cloud-obscured in the DFO product.It is hydrologically plausible that the flood  water is localised by these elevated features and we therefore conclude also that our algorithm is useful to detect floods in cloud-covered areas when used in tandem with MODIS.

Pakistan floods, 2010
The This algorithm can be modified for different areas based on the most appropriate thresholds for that location.The GFO flood probabilities that are most comparable to the validation image are applied at a threshold as higher than 40 % (Fig. 10, middle left panel) and higher than 50 % (Fig. 10, middle right panel).Especially in the east, corresponding to lower incidence angles, misclassifications can be seen, which are also recognised by the low q values.Indeed, looking at the quality values of higher than 50 (bottom left panel) and higher than 60 (bottom right panel), almost the entire area would have been removed if the threshold is applied above 50 % quality.The GFO algorithm marks almost the entire flooded area as too low quality for reliable flood classifications.Apparently, the combination between the land type and the incidence angles causes the large low quality zone in this area.

Using a global freshwater dataset as a training histogram
The use of a global freshwater histogram implies that local features may be neglected.An improved algorithm could use more local information, as some areas are more exposed to wind than others, and some water bodies may contain more sediment or salt.However, a separate training set per 1 × 1 • latitude-longitude tile is also not feasible, as not all of these tiles include freshwater bodies large enough to build a statistically robust training set.Using higher resolution will reduce this problem, as also smaller rivers can be taken into account in the training histograms.Even then, 1 × 1 • latitudelongitude areas in the world exist where there are no freshwater bodies.A recommendation after this research is that a validation of the algorithms in different climatological zones could improve the classification.

Misclassifications at low to intermediate incidence angles
In the Pakistan case study, we have analysed that the GFO algorithm marks almost the entire flooded area as too low a quality for reliable flood classifications.Ignoring the q value, we see that low to intermediate incidence angles α (lower than app.23-24 • ) are misclassified in this way.In the Thailand case study (Figs.7 and 8) the flooded area is in between α of 24 and 49 • .Figure 7 shows that for lower α (from 30 to 24 • ), the flood probability is still higher than 70 %, but decreasing as α decreases.If we interpolate this finding on the misclassification result of the Pakistan flood, it is reconfirmed that the probability threshold is related to the incidence angle.Depending on the 1 × 1 latitude-longitude degree tile on the Earth, flood classification is not reliable in an area where our land and water training histograms overlap and we cannot distinguish between land and water.This typically occurs at low to intermediate incidence angles.However, as also shown in Fig. 3, for even lower α values distinction could be possible again.Using the q value as defined in our algorithm masks out the lower quality areas, sometimes leading to a result where it shows that the flood map is not reliable at all.To consistently threshold the probability and quality values into a final water/land classification, while maintaining the automated character of the algorithm, more research is recommended to look at the behaviour of the q values in different climatological zones.

Improvements and present application of the topographic (HAND) index
The process to reach a topologically sound and accurate drainage network introduces occasional canyon-like artifacts into any DEM, as a result of aberrant height differences adjacent to the drainage network.These artifacts are transferred to the HAND grid during computation.An example of these artifacts is shown in detail in Fig. 11 in the area indicated by the ellipse.We processed the currently used HAND indices from the HYDRO1k dataset, based on GTOPO30 data.
In the currently used HAND, we use the empirically based data threshold of 15 (i.e.all data higher than 15 will not be processed and set to 0 % water probability and 100 % quality indicator).This is mainly done to be "on the safe side": to prevent flood-prone areas with artifacts to be wrongly filtered out.Rennó et al. (2008) state that when using original SRTM data for the HAND grid computation, these artifacts associated with the corrected DEM can be avoided.To filter more efficiently an improved version is recommended, using the HAND-data based on the HydroSHEDS 30 arc sec DEM (Lehner et al., 2006), in which the data are upgraded to streams that are "burned" less deeply in the DEM.
Also in regard to topographic effects, the simplified explanation shown in Fig. 12 (left panel) shows that terrain slopes, when assuming a small swath width and thus neglecting the ellipsoid of the geoid, cause most of the difference between incidence angle and local incidence angle.Note that a wrongly assumed incidence angle causes a different backscatter returned to the satellite: the σ − α pair shifts in our 2-D histograms.Figure 12 (right panel) shows a slice of the 2-D histogram at a certain incidence angle for a certain place on the globe.Correcting for this slope will shift the position of σ − α pairs and cause a noticeable shift in the 2-D histograms.For the present algorithm, however, this improvement is not deemed efficient, as we already filter out non-flood prone areas.Removing all HAND values higher than 15 in fact means that the pixels we do use are never higher than 15 m above the nearest drainage point.The largest shift of incidence angle is thus expected directly near the drainage point.
Looking again at Fig. 11, and assuming a GFO pixel size of 1 × 1 km, the shift in incidence angle is where β is the terrain slope, α is the incidence angle, θ the local incidence angle, z the elevation and x the length for which the slope is calculated.The shift is less than 1 • when z = 15 m and x = 1 km.Lastly, in regard to topographic corrections: when working with 1 × 1 km pixel scales, we can consider the errors in global topography models with roughly the same resolution, such as GTOPO30 or the newer 30-arc second global mosaic of the Shuttle Radar Topography Mission (SRTM, USGS, 2004).According to Rodríguez at al. (2006), SRTM can give average absolute height errors per continent up to almost 10 m and locally even higher.Harding et al. (1999) indicate that GTOPO30 can reach even higher errors of 30 m.When we consider this error in z in Eq. ( 11), it is clear that a global topography model can also cause shifts in incidence angles of the same order and higher than the maximum shift we expect in our 2-D histograms after applying the HAND-index based filter.Furthermore, when taking into account computer processing efficiency, it is practical to avoid the correction to local incidence angle, at least until higher quality global digital elevation models are available.

Conclusions
In preparation for the Sentinel-1 SAR satellite, and in order to address the urgent need for fast flood water detection and mapping, systematic and automated processing algorithms are needed.A binary product (e.g.water/not water or flooded/not flooded) is not optimum, as SAR-based classification products include noise and therefore an uncertainty indication is desirable.In this article an automated method to calculate probability of water, including a quality indicator, from Level 1 Envisat ASAR data is presented.The method is a new contribution, primarily because, in our approach, (a) data over a certain time span are stored in 2-D histogram training sets in the incidence angle -backscatter domain and (b) because the algorithms automatically calculate a water probability and a quality indicator for each image pixel.Also, the application of the HAND index as a pre-processing filter improves the final result.The results of our algorithm are evaluated using two different algorithms, using ASAR and MODIS data.A first merge with MODIS imagery in a case study in Thailand shows strong resemblance between the ASAR and MODIS derived results.At locations where MODIS suffers from clouds, ASAR shows hydrologically correct results, as observed through the clouds and as verified by other knowledge.We recommend more research in merging SAR and MODIS derived water imagery, to combine the strengths of both methods and improve the desired operational global surface water product.A second case study in Pakistan shows the added value of the quality indicator q to remove unreliable classifications, unfortunately also removing most of the data in the flooded area.In order to combine the product, one should be able to consistently apply a threshold with the probability from our algorithm into a final water/land classification.By looking at our two case studies, we see that the threshold values to make binary flood maps of our water probabilities varies, depending on the area in the world and the incidence angle.More research on the dependency of our algorithm to different climatological zones is recommended.

Fig. 1 .
Fig. 1.The backscatter is the portion of the outgoing satellite radar signal, usually looking sideways (left panel) in different incidence angles α (right panel) is highly depending on the backscatter characteristics of the subsurface (middle panel).Adapted from ESA (2012).

Fig. 2 .
Fig. 2. Trained and gaussian smoothed 2-D histograms for land and water (left panel) as derived from two years of ASAR WSM backscatter data in the Netherlands (right panel).

Fig. 3 .
Fig. 3. Explanation of the ability to distinguish between water and land for different incidence angles.

Fig. 4 .
Fig. 4. The height above nearest drainage (HAND) Index for an example location in Thailand.

Fig. 5 . 2 Figure 6 . 6 Fig. 6 .
Fig. 5. Different histograms as shown on a global backscatter map.Different areas in the world show different backscatter characteristics.In this figure, the backscatter characteristics for rain forest (left panel), desert (middle panel) and freshwater (right panel) are shown.1

Figure 7 .Fig. 7 .
Figure 7. Results of the GFO water probabilities over the Chao Phraya basin on 13 October 3 2011 as shown in Google Earth (left).HAND values with flooded areas (light blue) are shown 4 on the right.5 6 7

Fig. 11 .
Fig. 11.Detail of artifact along the river in the HAND index.

Fig. 12 .
Fig. 12. Left panel: the difference between incidence angle α and local incidence angle θ is mainly caused by the terrain slope β.Right panel: a slice of the 2-D histogram at α = 32.5 • (location latitude-longitude = 27.5-87.5).Corrections for the slope will shift the position of σ − α pairs.