Automated Inundation Mapping Over Large Areas Using Landsat Data and Google Earth Engine

: Accurate inundation maps for ﬂooded 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 peak inundation in the Okavango Delta, northern Botswana, using Landsat imagery and Google Earth Engine. Inundation classiﬁcation in the Okavango Delta is complex owing to the spectral overlap between inundated areas covered with aquatic vegetation and dryland vegetation classes on satellite imagery, and classiﬁcations have predominately been implemented on broad spatial resolution imagery. We present the longest time series to date (1990–2019) of inundation maps for the peak ﬂood season at a high spatial resolution (30 m) for the Okavango Delta. We validated the maps using image-based and in situ data accuracy assessments, with overall accuracy ranging from 91.5% to 98.1%. Use of Landsat imagery resulted in consistently lower (on average, 692 km 2 ) estimates of inundation extent than previous studies that used Moderate Resolution Imaging Spectroradiometer (MODIS) and National Oceanic and Atmospheric Administration Advanced Very-High-Resolution Radiometer (NOAA AVHRR) imagery, likely owing 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. This classiﬁcation method can likely be adapted for inundation mapping in other regions.


Introduction
The Okavango Delta (the delta) in northern Botswana is a wetland of international and domestic significance [1][2][3], yet pressures on its water resources from water abstraction (for agriculture and human consumption), damming for power generation, and climate change are growing [1,4,5]. This large wetland consists of a panhandle region; a channel system surrounded by permanent swamps; and a large, low gradient alluvial fan [4,6,7]. The delta is subject to an annual flood event asynchronous with the local rainy season; rainfall in the highlands of Angola flows into the Okavango River, entering the Botswana panhandle, and slowly moves down the fan, reaching maximum inundation extent in July-September [2,4,6,[8][9][10]. Intra and inter-annual variations in the frequency, duration, and extent of the inundation produce a complex mosaic of vegetation, supporting a vast number of ecological niches and a rich diversity of flora and fauna [1,11,12].
The hydrology of the delta, including temporal and spatial changes in its inundation history, has been investigated through inundation maps [1]. These maps can be used to study the past and present state of the delta; to predict its future transformations; and to understand how it is affected by natural processes, climate change, and human resource use [13]. They may also be incorporated into management strategies and biodiversity studies. Inundation maps can be created using satellite imagery, which are available from a range of spatial and temporal resolution products. While differentiating open water (e.g., channels, lagoons) from dryland vegetation (e.g., shrublands, grasslands) is relatively simple, there is substantial overlap in the spectral values of inundated areas covered in aquatic vegetation (e.g., floodplain) and some dryland vegetation classes, making the separation of these classes difficult and traditional water classification methods unviable [2,[8][9][10]14,15]. For example, the normalized and modified normalized difference water index (NDWI and MNDWI), which were specifically developed to map waterbodies, had the least ability to classify inundation compared with six other methods [9]. Therefore, a range of classification methods (unsupervised, supervised, band thresholding, band ratios, indices, and combinations of these methods) have been implemented [2,4,6,[8][9][10]15]. Recently, band thresholding has been successful [2,8,15], with thresholding of the short wave infrared (SWIR) band producing high accuracy results on Moderate Resolution Imaging Spectroradiometer (MODIS) imagery [9]. The SWIR band is highly sensitive to moisture content [16], and can differentiate densely vegetated inundated areas from non-inundated vegetation [9]. As well as its accuracy, the advantage of this method is its relative simplicity, meaning it is easily automated, which reduces the time (and thus cost) of implementation compared with more complex methods.
The majority of delta inundation studies have used imagery with broad spatial resolution (MODIS (250 m, 500 m, and 1 km) [8,9,15] and National Oceanic and Atmospheric Administration Advanced Very-High-Resolution Radiometer (NOAA AVHRR) (1 km) [6,10]), taking advantage of the high temporal resolution of these sensors, which allows daily and sub-monthly analysis of inundation [9]. 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 [8]. Further, high spatial resolution information can increase confidence in associated decision-making [17][18][19][20][21][22][23][24][25]. Broad spatial resolution can be downscaled to achieved finer spatial resolution, if access to high temporal resolution data is a pertinent factor [8,[26][27][28]. In addition, broad spatial resolution increases the likelihood of mixed pixels (e.g., pixels containing both inundated and dry areas), which can confuse classification attempts [4,9,29], although methods exist to reduce this issue [26,30,31].
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 [4,8,9,15]. Recent advances in computing power and cloud-processing infrastructure (e.g., Google Earth Engine [32]) have enabled much wider access to satellite image time series, along with the 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 peak flood season for the delta at a high spatial resolution (30 m pixels) to date. We adapted a previously developed method based on thresholding of the SWIR band [9], and implemented an automated version in Google Earth Engine [32], 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 season on record [33]. 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 to 2019 was used (Step 1 in Figure 2). Scenes occurred within six Landsat path/row footprints (174/073, 174/074, 175/073, 175/074, 176/072, 176/073). These sensors capture scenes with 30 m 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 classified as cloud or cloud shadow on the Landsat cloud mask band were masked (Step 2 in Figure 2). These pixels were then filled using the median value for the pixel from a year before and after the scenes' date, using a gap-filling algorithm [34] (step 3 in Figure 2). The SWIR band was selected for each scene (step 4 in Figure 2) and then, for each year, a composite was created from the median value of all the scenes for that year (step 5 in Figure 2). Annual composites with large areas missing (e.g., there were no scenes available for a path/row) were filtered out (step 6 in Figure 2). This occurred five times (1993, 2000, 2009, 2010, and 2012 cloud shadow on the Landsat cloud mask band were masked (Step 2 in Figure 2). These pixels were then filled using the median value for the pixel from a year before and after the scenes' date, using a gap-filling algorithm [34] (step 3 in Figure 2). The SWIR band was selected for each scene (step 4 in Figure 2) and then, for each year, a composite was created from the median value of all the scenes for that year (step 5 in Figure 2). Annual composites with large areas missing (e.g., there were no scenes available for a path/row) were filtered out (step 6 in Figure 2). This occurred five times (1993,2000,2009,2010, and 2012).

Annual (July-September) Inundation Maps
The composites were transformed into inundation maps using an SWIR thresholding technique. To apply this method, we assessed and digitised areas that had permanent water (e.g., permanent swamp or main channels) or were permanently dry, based on Wolski et al.'s [9] designated areas, but altered slightly to suit our use of higher spatial resolution imagery (Figure 1, step 7 in Figure 2). The median SWIR value for the inundated (SWIRwet) and dry (SWIRdry) areas was calculated for each individual composite, and a composite-specific SWIRthreshold value was calculated using Equation 1 (taken from Wolski et al. [9]) (step 8 in Figure 2).
The relative frequency of SWIR values for the wet and dry areas is shown in Figure S1, with the SWIRwet, SWIRdry, and SWIRthreshold values marked. Pixels with an SWIR value below the threshold were classified as inundated, and vice-versa for dry pixels (step 9 in Figure 2). Calculating the threshold value separately for each image accounts for the dynamic (seasonal and annual) nature of the inundation patterns of the delta [9]. The multiplier of 0.3 represents the value needed to calculate the correct threshold to classify a pixel with an inundation fraction of 50% as inundated, as recommended by Wolski et al. [9]. Owing to the different imagery used in this study, we confirmed the value of 0.3 was appropriate, by also assessing inundation maps developed using values of 0.25 Figure 1. Okavango Delta study area (black outline) in northern Botswana. White and black hatched areas are the permanent water (permanent swamps and channels) and dry areas, respectively, used in calculating the threshold value, and blue lines are the major channels of the delta.

Annual (July-September) Inundation Maps
The composites were transformed into inundation maps using an SWIR thresholding technique. To apply this method, we assessed and digitised areas that had permanent water (e.g., permanent swamp or main channels) or were permanently dry, based on Wolski et al.'s [9] designated areas, but altered slightly to suit our use of higher spatial resolution imagery (Figure 1, step 7 in Figure 2). The median SWIR value for the inundated (SWIR wet ) and dry (SWIR dry ) areas was calculated for each individual composite, and a composite-specific SWIR threshold value was calculated using Equation (1) (taken from Wolski et al. [9]) (step 8 in Figure 2).
The relative frequency of SWIR values for the wet and dry areas is shown in Figure S1, with the SWIR wet , SWIR dry , and SWIR threshold values marked. Pixels with an SWIR value below the threshold were classified as inundated, and vice-versa for dry pixels (step 9 in Figure 2). Calculating the threshold value separately for each image accounts for the dynamic (seasonal and annual) nature of the inundation patterns of the delta [9]. The multiplier of 0.3 represents the value needed to calculate the correct threshold to classify a pixel with an inundation fraction of 50% as inundated, as recommended by Wolski et al. [9]. Owing to the different imagery used in this study, we confirmed the value of 0.3 was appropriate, by also assessing inundation maps developed using values of 0.25 and 0.35 (see next section). Further details about this classification method and its development are provided in Wolski et al. [9]. and 0.35 (see next section). Further details about this classification method and its development are provided in Wolski et al. [9].

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. Although this method can be subjective, it is a generally accepted approach, particularly when historical data are limited [2,[35][36][37][38][39][40][41]. We validated the accuracy of three sets of inundation maps (produced using 0.25, 0.3, and 0.35 in Equation 1) by visual assessment of true colour versions of the Landsat composites used to make the inundation maps, and high-resolution satellite images. The inundation maps were created without filling masked pixels, and only comprised a subset of the years (2000-2016). 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 that the study area ( Figure 1) included large tracts of permanently dry areas (i.e., the Kalahari Desert), which we predicted would rarely be misclassified [9], we used an amended area for the validation, removing some of the larger dry areas ( Figure S2). Using this amended area, fifty sample points ( Figure S2) were randomly generated using the sampleRandom function (raster package [42]) in R [43]. For each year, the same 50 sample points were visually assigned as inundated or dry on the Landsat and high-resolution imagery before progressing onto the next year. This prevented the assessors from making classifications based on a sample point's previous inundation history. The classification (inundated or dry) of each sample point was extracted from the relevant inundation map and an error matrix was created. Overall accuracy (the sum of the diagonal entries (correctly classified points) divided by the total sampled points), producer's accuracy (the diagonal entry divided by its respective column total), and user's accuracy (the diagonal entry of each row divided by its respective row total) were calculated.

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. Although this method can be subjective, it is a generally accepted approach, particularly when historical data are limited [2,[35][36][37][38][39][40][41]. We validated the accuracy of three sets of inundation maps (produced using 0.25, 0.3, and 0.35 in Equation (1)) by visual assessment of true colour versions of the Landsat composites used to make the inundation maps, and high-resolution satellite images. The inundation maps were created without filling masked pixels, and only comprised a subset of the years (2000-2016). 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 that the study area ( Figure 1) included large tracts of permanently dry areas (i.e., the Kalahari Desert), which we predicted would rarely be misclassified [9], we used an amended area for the validation, removing some of the larger dry areas ( Figure S2). Using this amended area, fifty sample points ( Figure S2) were randomly generated using the sampleRandom function (raster package [42]) in R [43]. For each year, the same 50 sample points were visually assigned as inundated or dry on the Landsat and high-resolution imagery before progressing onto the next year. This prevented the assessors from making classifications based on a sample point's previous inundation history. The classification (inundated or dry) of each sample point was extracted from the relevant inundation map and an error matrix was created. Overall accuracy (the sum of the diagonal entries (correctly classified points) divided by the total sampled points), producer's accuracy (the diagonal entry divided by its respective column total), and user's accuracy (the diagonal entry of each row divided by its respective row total) were calculated.

Validating Inundation Maps (In Situ Data Accuracy Assessment)
To validate the accuracy of the classification method, we carried out a field examination of inundated and dry regions within one Landsat scene (scene 175/73, 25 July 2018). Owing to accessibility and safety constraints, we only sampled from the Abu Concession ( Figure S3), where inundated areas could be accessed by field personnel by wading (within 100 m of dry land). Sampling points were chosen using a random stratified sampling approach, where inundated and dry were the stratification Remote Sens. 2020, 12, 1348 5 of 12 levels. The sampling area was created in QGIS [44] by outlining islands that were accessible by vehicle, applying a 100 m buffer, and clipping the inundation map to this shape. The raster was imported into R [43] and 55 points in each stratum were randomly selected ( Figure S3) using the sampleStratified function (raster package [42]), 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 [35] 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 on which class occurred over the majority of the 30 m 2 area centered on each point. Where the proportion of each class was approaching equality, the point was classified, but was also noted as an uncertain classification. The classification for each sample point was extracted from the inundation map created from the Landsat scene and an error matrix was created. Overall accuracy, user's accuracy, and producer's accuracy were calculated as above.

Results
The extent and distribution of the peak inundation varied annually (Figures 3 and 4

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 five and four validation pixels were masked, respectively. There were fewer points (123) visually classified using the high-resolution imagery owing to a lack of available data. For inundation maps created using a multiplier of 0.3, based on the 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-

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 five and four validation pixels were masked, respectively. There were fewer points (123) visually classified using the high-resolution imagery owing to a lack of available data. For inundation maps created using a multiplier of 0.3, based on the 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). Inundation maps created using a value of 0.25 had slightly lower overall accuracy, and using a value of 0.35 had almost identical overall accuracy (Table S1). Misclassified points were located predominately on the boundary of inundated and dry areas and where there was high inter-annual variation in inundation ( Figure S2). Table 1. Error matrices and overall accuracy of inundation maps using image-based accuracy assessment (Landsat and high-resolution imagery) and in situ points. Note: values in square brackets are based on points noted as 'uncertain' being removed.

Landsat
Hi 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 overall accuracy of 91.5% (nine sample points were misclassified) (Table 1, Figure S3). 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, overall accuracy increased to 98.0% (Table 1).

Discussion
Our study details the longest ever time series of peak flooding extents for the Okavango Delta at a high spatial resolution (30 m), demonstrating the remarkable inter-annual variability of this system; the largest inundation extent recorded was almost three times that of the smallest (Figures 4 and 5). There are also inter-annual variations in the spatial distribution of inundation (Figures 3 and 5, Data S1 (annual rasters)), driven by the volume of water discharged into the system, but also factors such as sedimentation, channel blockage from vegetation, and avulsion [6,8,45,46]. On the basis of the maps produced in this study, the 2019 flood event represents the smallest inundation since 1985, being around 769 km 2 smaller than the previous record in 1996 (Figure 4). Estimations of inundation extent going back to 1934 calculated the lowest inundation to be approximately 5100 km 2 [6]. Our estimate of inundation extent in 2019 was 3483 km 2 , making it the smallest flood in the last 85 years. This exceptionally low flood is likely driven by a multi-decadal (16-20 years) rainfall cycle in Southern Africa [6,47,48], being 14 years since the previous dry year. The delta is also at risk of drying owing to increases in temperature and evaporation and decreases in rainfall and river flow due to climate change and water abstraction and damming [4,49].
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 Remote Sens. 2020, 12, 1348 8 of 12 resolution imagery, an effect that is also evident in other study systems when comparing estimates from different sensor resolutions [50][51][52][53]. 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 [4,29], 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 (Murray-Hudson et al. [15] is an exception), instead using comparisons to high resolution aerial orthomosaics [2], other inundation maps [8,9,15], or hydrological observations [6,9]. The overall 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 on which the inundation maps were based. Also contributing to the reduced accuracy, high-resolution images covered a slightly longer time period (July-October) than that 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. Although only a subset of the maps (2000-2016) was validated, Figure S1 provides further evidence that the thresholding method can accurately detect the boundary between inundated and dry pixels, with non-validated and validated years following a similar pattern.
The overall 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 S3). In situ data collection in this area is difficult owing to logistics, accessibility, and safety issues, particularly during the high flood period [2,4,10,54]; 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 that 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 [9] (Figure S2), and which was under-represented in the sample points of the image-based accuracy assessment. 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. [9], 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 [2] and where riparian woodland can saturate Landsat pixels [14]. A potential solution to minimise this would be filtering out pixels that are discontinuous from the larger inundated area [2], 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, typically during the wet season [9], which this study did not measure. The anomalously low flood level in 2019 meant it was difficult to get suitable training data that were consistent with the other maps. Therefore, dry pixels within the permanent water polygon increased the range of SWIR wet values ( Figure S1). This may have led to a higher level of misclassification ( Figure 5).
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 [6]. Using 500 m spatial resolution imagery, Wolski et al. [9] noted that some important terminal rivers of the delta were not Remote Sens. 2020, 12, 1348 9 of 12 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, 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).

Conclusions
Wolski et al. [9] developed a simple method (thresholding of the short wave infrared band) to accurately classify inundation in the Okavango Delta using broad spatial resolution (500 m) 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 for the peak flood season at a relatively high spatial resolution (30 m) to date. The inaccessibly of remote sensing methods and processing capability has prevented wide-spread adoption of its use by non-experts. We anticipate that the methods/code and the data produced in this paper can be used and adapted by land managers, researchers, and other stakeholders, who require access to accurate high resolution inundation maps. Further, the classification method is likely to be suitable for mapping inundation in other regions, with only minimal adaptation of the methods and code presented here.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/8/1348/s1, Figure S1: Density of SWIR values for permanent water and dry areas, with median values (dashed line) and threshold value (solid line) for each year; Figure S2: Location of validation points for image-based accuracy assessment. Black points were always correctly classified, red points were incorrectly classified at least once (label displays the number of times incorrectly classified out of total times classified). Dashed line represents amended area for point sampling and background map is the summed annual inundation map; Figure Table S1: Error matrices and overall accuracy of inundation maps using image-based accuracy assessment (Landsat and high-resolution imagery) using alternative values of f in the threshold equation SWIRthreshold = SWIRwet + f * (SWIRdry − SWIRwet); Data S1: Rasters of individual inundation maps (1990-2019); Data S2: Raster of sum of all inundation maps; Data S3: Raster of variance of all inundation maps; Data S4: Annual inundation extents; Code S1: Google Earth Engine code and inundation rasters are available on Github (https://github.com/VictoriaInman/OkavangoDelta_flooding) and archived on Zenodo (DOI: 10.5281/zenodo.3693153); Code S2: Link to Google Earth Engine repository https://code.earthengine.google. com/?accept_repo=users/victoriainman/OkavangoDelta_TechnicalNote.