Automated inundation history mapping over large areas using Landsat and Google Earth Engine

: Accurate inundation maps for flooded wetlands and rivers are a critical resource for their management and conservation. In this paper we automate a method (thresholding of the short-wave infrared band) for classifying inundation, using Landsat imagery and Google Earth Engine. We demonstrate the method in the Okavango Delta, northern Botswana, a complex case study due to the spectral overlap between inundated areas covered with aquatic vegetation and dryland vegetation classes on satellite imagery. Inundation classifications in the Okavango Delta have predominately been implemented on broad spatial resolution imagery. We present the longest time series to date (1990-2019) of inundation maps at high spatial resolution (30m) for the Okavango Delta. We validated the maps using image-based and in situ data accuracy assessments, with accuracy ranging from 91.5 - 98.1%. Use of Landsat imagery resulted in consistently lower estimates of inundation extent than previous studies, likely due to the increased number of mixed pixels that occur when using broad spatial resolution imagery, which can lead to overestimations of the size of inundated areas. We provide the inundation maps and Google Earth Engine code for public use.


Introduction
Inundation maps can be used to study the past and present state of wetlands, to predict their future transformations, and to understand how natural processes, climate change, and human resource use affect them [1]. Further, inundation maps can be incorporated into management strategies and biodiversity studies. If the area of interest is large or difficult to access, if multiple maps are required, or human resources are limited, then inundation maps can be created using a historical suite of satellite imagery. A range of spatial and temporal resolution products are available, and higher spatial resolution information can increase confidence in associated decision-making [2].
The Okavango Delta (the delta) in northern Botswana is a wetland of international and domestic significance [3][4][5], yet pressures on its water resources from water abstraction (for agriculture and human consumption), damming for power generation, and climate change are growing [3,6,7]. This large wetland consists of a panhandle region, a channel system surrounded by permanent swamps, and a large, low gradient alluvial fan [6,8,9]. The delta is a flood pulsed system subject to an annual flood event which is asynchronous with the local rainy season; rainfall in the highlands of Angola flows into the Okavango River, entering the Botswana panhandle in April/May, and slowly moves down the fan, reaching maximum inundation extent in July-September [4,6,8,[10][11][12]. Intra and inter-annual variations in the frequency, duration, and extent of the inundation produces a complex mosaic of vegetation, supporting a vast number of ecological niches and a rich diversity of flora and fauna [3,13,14].
The hydrology of the delta, including temporal and spatial changes in its inundation history, have been investigated through inundation maps [3]. While differentiating open water (e.g. channels, lagoons) from dryland vegetation (e.g. shrublands, grasslands) on satellite imagery is relatively simple, there is spectral overlap between inundated areas covered in aquatic vegetation (e.g. floodplain) and some dryland vegetation classes, making separation of these classes difficult and traditional water classification methods unviable [4,[10][11][12]15,16]. Therefore a range of classification methods (unsupervised, supervised, band thresholding, band ratios, indices, and combinations of these methods) have been trialled [4,6,8,[10][11][12]16]. Recently, the relatively simple method of band thresholding has been successful [4,10,16], with thresholding of the short wave infrared (SWIR) band producing high accuracy results on MODIS imagery [11].
The majority of delta inundation studies have used imagery with broad spatial resolution (MODIS (250m, 500m, and 1km) [10,11,16] and NOAA AVHRR (1km) [8,12], but see [4,16]), taking advantage of the high temporal resolution of these sensors, which allows sub-monthly analysis of inundation [11]. However, dependent on the intended use of the inundation maps, such broad spatial resolution may result in unacceptable simplification of the complex mosaic of the delta [10]. In addition, broad spatial resolution increases the likelihood of mixed pixels (e.g. half inundated, half dry) which can confuse classification attempts [6,11].
Computational power, data procurement, management, and storage, and processing times have also traditionally been a motivation for using broad spatial resolution images, particularly when creating time series over large areas [6,10,11,16]. Recent advances in computing power and cloud-processing infrastructure (e.g. Google Earth Engine [17]) have enabled much wider access to satellite image time series, along with capacity to process and analyse these data.
In this paper we utilised the family of Landsat satellite sensors to create the longest time series of inundation maps for the delta at high spatial resolution (30m pixels) to date. We adapted a previously developed method based on thresholding of the SWIR band [11], and implemented an automated version in Google Earth Engine [17], a cloud-based geospatial analysis platform. We created a time series of peak inundation for the last 30 years, up to and including the flood event of 2019, thought to be the lowest flood on record [18]. Further, we provide validation results that confirm the accuracy of the SWIR thresholding method. The inundation maps and Google Earth Engine code are provided publicly for use and adjustment by stakeholders, land managers, and researchers.

Annual (July -September) Landsat composites
Unless otherwise stated, all processing was conducted using Google Earth Engine. Every tier 1 atmospherically corrected surface reflectance Landsat (Landsat 5 TM, Landsat 7 ETM+, Landsat 8 OLI) scene covering the Okavango Delta ( Figure 1) for the peak inundation period (1 July to 30 September) from 1990-2019 were used. These sensors capture scenes with 30m spatial resolution, containing eight (Landsat 5 and 7) and eleven (Landsat 8) bands, including the SWIR band (band 7) used in this study. For each scene, pixels classed as cloud or cloud shadow on the cfmask band were masked and then these pixels were filled using the median value for the pixel from a year before and after the scenes' date, using a pre-existing algorithm (www.stackoverflow.com/questions/55256739).
For each year, a composite was then created from the median value of all the scenes for that year (Figure 2b). Annual composites with large areas missing (e.g. whole scenes) were filtered out.

Annual (July -September) inundation maps
The composites were transformed into inundation maps using a SWIR thresholding technique. To apply this method, we assessed and digitised areas that were permanently inundated (e.g. permanent swamp or main channels) or permanently dry, based on Wolski et al.'s [11] designated areas but altered slightly to suit our use of higher spatial resolution imagery ( Figure 1). The median SWIR value for the inundated (SWIRwet) and dry (SWIRdry) areas were calculated for each individual composite, and a composite-specific SWIRthreshold value calculated using Eq. 1 (adapted from Wolski et al. [11]).
Pixels with a SWIR value below the threshold were classified as inundated and vice-versa for dry pixels (Figure 2c). Calculating the threshold value separately for each image accounts for the dynamic (seasonal and annual) nature of the inundation patterns of the delta [11]. Further details about this classification method and its development are in Wolski et al. [11].

Validating inundation maps (image-based accuracy assessment)
Gathering in situ data for historical time series is difficult, and in many cases these data simply do not exist. To enable validation for at least some of the time series, we used a visual interpretation approach to generate reference data, a generally accepted approach [19]. We validated the accuracy of a subset (2000 -2016) of the inundation maps (created without filling masked pixels) by visual assessment of true colour versions of the Landsat composites used to make the inundation maps, and high-resolution satellite images. High resolution imagery taken between July and October was accessed via Google Earth's historical imagery function and the Digital Globe collection (obtained from the DigitalGlobe Foundation). Given the study area ( Figure 1) included large tracts of permanently dry areas (i.e. the Kalahari Desert) which we predicted would rarely be misclassified [11], we used an amended area for the validation, removing some of the larger dry areas ( Figure S1). Using this amended area, fifty sample points ( Figure S1) were randomly generated using the sampleRandom function (raster package [20]) in R [21]. All 50 sample points were visually assigned as inundated or dry on the Landsat and high-resolution imagery for each year before progressing onto the next year to avoid bias of temporally proximate samples. The classification (inundated or dry) of each sample point was extracted from the relevant inundation map and an error matrix created.

Validating inundation maps (in situ data accuracy assessment)
To further validate the accuracy of the classification method we conducted an accuracy assessment using in situ data of an inundation map created from one Landsat scene (scene 175/73, 25 July 2018). Due to accessibility constraints during the peak flood period, we only sampled from a small section of the scene (within the Abu Concession, Figure S2) and due to safety constraints only sampled inundated areas that were within 100m of dry land (that could be accessed by wading). Sampling points were chosen using a random stratified sampling approach, where inundated and dry were the stratification levels. The sampling area was created in QGIS [22] by outlining islands that were accessible by vehicle, applying a 100m buffer, and clipping the inundation map to this shape. The raster was imported into R [21] and 55 points in each stratum were randomly selected ( Figure S2) using the sampleStratified function (raster package [20]), which were then exported and uploaded to a handheld Garmin GPSMAP® 64 GPS. This number of sample points was chosen to ensure that if some points were inaccessible (e.g. vegetation too thick to drive through, unsafe to wade into water, wildlife within close proximity) the recommended minimum number of 50 [19] could still be obtained. Data collection occurred within two days of the sensor's collection of the scene (25 -27 July 2018), with each point classified as either inundated (standing water) or dry based

Validating inundation maps
Using the true colour Landsat composites, a total of 691 points were visually classified; fifty points per year, except for years without inundation maps and 2010 and 2011, where 5 and 4 validations pixels were masked, respectively. There were fewer points (123) visually classified using the high-resolution imagery due to a lack of available data. Based on visual assessment of Landsat composites, the inundation maps had few misclassified pixels (1.9%), but this was slightly higher based on the visual classification of high-resolution imagery (4.1%)( Table 1).
Out of the 110 sample points that were generated for the in situ validation, we classified 106 (four were inaccessible). The inundation map had an accuracy of 91.5% (nine sample points were misclassified) (Table 1, Figure S2). Of these nine misclassified points, seven were noted as uncertain in the field as they had approximately equal areas of inundation and dry. There was an additional uncertain point that was correctly classified. When these uncertain points were removed, accuracy increased to 98.0% (Table 1).

Discussion
Our study details the longest ever time series of flooding extents for the Okavango Delta at high spatial resolution (30m), demonstrating the remarkable inter-annual variability of this system; the largest inundation extent recorded was almost three times that of the smallest (Figures 3 and 4). Based on available inundation maps, the flood event of 2019 represents the smallest inundation extent since 1985, being around 769km 2 smaller than the previous record in 1996 (Figure 4). Further, predictions of inundation extent (based on annual in-flow and rainfall) going back to 1934 calculated the lowest inundation to have occurred in 1996 [8], suggesting 2019 is the smallest flood in 85 years.
The annual inundation extent estimates in our study were systematically smaller than previous studies but showed a similar trend (Figure 4). The most likely cause for this effect is our use of higher resolution imagery, an effect that is also evident in other study systems when comparing estimates from different sensor resolutions [23][24][25][26]. Broad spatial resolution imagery (as used in existing studies) increases the number of mixed pixels, and can lead to overestimations of the size of the inundated areas [6], which we reduced by using Landsat images. Another potential cause for smaller estimates is our use of three-month composites rather than individual consecutive images, concealing the time of true maximum inundation, although our use of median values should be robust to this effect.
Accuracy assessments of delta inundation maps have generally not included in situ data (but see [16]), instead using comparisons to high resolution aerial orthomosaics [4], other inundation maps [10,11,16], or hydrological observations [8,11]. The accuracy of the method used in this study as determined by the image-based assessment (95.9 -98.1%) is comparable to these other studies, being predictably higher for the Landsat composites than the independent high-resolution images, as these are the true colour version of the images that the inundation maps were based on. Also contributing to the reduced accuracy, high-resolution images covered a slightly longer time period (July -October) than was used to make the inundation maps (done to increase available data), and were individual images, as opposed to composites, meaning they may have occurred before/after the full extent of the flood.
The accuracy from in situ data validation (91.5%) was slightly lower than from the image-based assessment, although not so (98.0%) when we removed points we had flagged as uncertain in the field. These were points that were approximately half inundated and half dry (i.e. a mixed pixel) and had remaining moisture in the soil (i.e. were muddy) where the flood had recently receded (but which we classified as dry as there was no standing water) ( Figure S2). In situ data collection in this area is difficult due to logistics, accessibility, and safety issues, particularly during the high flood period [4,6,12,27]; so, while our in situ validation was small in scale, it represents a rarely conducted true accuracy assessment of delta inundation mapping. In addition, given the sampling area was centered on small islands and edges of the floodplain, it fittingly represents the boundary between dry and inundated areas, the area where most classification errors are likely to occur [11]. Therefore, the high accuracy within this sampling area suggests the overall delta wide classification is likely to be reliable, confirmed by the delta-wide image-based validation.
In addition to inaccuracy caused by mixed pixels, we noted, as did Wolski et al. [11], the presence of some true misclassification when using this method. Visual inspection of our maps suggests riparian woodland vegetation is sometimes misclassified as inundated area, a known problem in the delta where these classes can have overlapping spectral signatures [4] and where riparian woodland can saturate Landsat pixels [15]. A potential solution to minimise this would be filtering out pixels that are discontinuous from the larger inundated area [4], although for the sake of simplicity we have not attempted to do this. Misclassifications can also occur where there is a small difference between the SWIR value of the dry and inundated areas used for the threshold calculation [11].
Choosing a sensor is a compromise between spatial and temporal resolution and the computational time and power required to process the images. Google Earth Engine allowed us to take advantage of high-resolution images with minimal effort; all images are called directly to the software without downloading and functions (e.g. cloud masking) can be automated. In an ecosystem as complex as the delta, broad spatial resolution maps may have restricted utility [8]. Using 500m spatial resolution imagery, Wolski et al. [11] noted that some important terminal rivers of the delta were not well represented on their inundation maps as they were narrower than the resolution of the imagery. However, by using Landsat imagery, the mosaic of floodplains and islands are well represented (Figure 2), and rivers that are important indicators of the hydrology of the delta, as well as essential to local communities (e.g. the Thamalakane River), are clearly mapped (Figure 3). Access to high spatial resolution maps is also important for spatial planning [28][29][30], assessing land use change [31][32][33], water resource planning [34,35], and wildlife habitat studies [36][37][38].

Conclusions
Wolski et al. [11] developed a simple method (thresholding of the short wave infrared band) to accurately classify inundation in the Okavango Delta using broad spatial resolution (500m) satellite imagery, noting the method was suitable for automation, but also cautioning that creating inundation maps using Landsat imagery was "laborious … making creation of a consistent, long time series of inundation maps difficult". In this paper, we have shown that periodic, accurate inundation maps can be created using relatively high-resolution imagery (Landsat) suitable to capture the complexity of this important ecosystem, by utilizing Google Earth Engine, a cloud-based platform. We provide the longest time series (1990-2019) of inundation maps at a relatively high spatial resolution (30m) to date.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Figure S1: Location of validation points for image-based accuracy assessment. Dashed line represents amended area for point sampling, Figure S2