High-resolution habitat and bathymetry maps for 65,000 sq. km of Earth’s remotest coral reefs

With compelling evidence that half the world’s coral reefs have been lost over the last four decades, there is urgent motivation to understand where reefs are located and their health. Without such basic baseline information, it is challenging to mount a response to the reef crisis on the global scale at which it is occurring. To combat this lack of baseline data, the Khaled bin Sultan Living Oceans Foundation embarked on a 10-yr survey of a broad selection of Earth’s remotest reef sites—the Global Reef Expedition. This paper focuses on one output of this expedition, which is meter-resolution seafloor habitat and bathymetry maps developed from DigitalGlobe satellite imagery and calibrated by field observations. Distributed on an equatorial transect across 11 countries, these maps cover 65,000 sq. km of shallow-water reef-dominated habitat. The study represents an order of magnitude greater area than has been mapped previously at high resolution. We present a workflow demonstrating that DigitalGlobe imagery can be processed to useful products for reef conservation at regional to global scale. We further emphasize that the performance of our mapping workflow does not deteriorate with increasing size of the site mapped. Whereas our workflow can produce regional-scale benthic habitat maps for the morphologically diverse reefs of the Pacific and Indian oceans, as well as the more depauperate reefs of the Atlantic, accuracies are substantially higher for the former than the latter. It is our hope that the map products delivered to the community by the Living Oceans Foundation will be utilized for conservation and act to catalyze new initiatives to chart the status of coral reefs globally.


Introduction
Humans have been damaging reefs since they first started to interact with them (e.g., Pandolfi et al. 2003;McClenachan et al. 2017), but it is only in the last 40 yrs, or so, that impacts such as overfishing, pollution, and climate change have precipitated their global collapse (Jackson et al. 2001;Bellwood et al. 2004). Targeted intervention can reverse this demise using tools ranging from marine protected areas to reef restoration, but to be effective, it is necessary to understand the location of Earth's reefs and their status. largest coral reef studies in history, visiting a global transect of remote shallow-water reef sites (Fig. 1). The primary goals of the GRE were to map and characterize coral reef ecosystems, identify their status and major threats, examine factors affecting their resilience, and to promote local and regional conservation efforts through data sharing, outreach, and education. A key stipulation of the endeavor was that the Living Oceans Foundation was invited by each host nation into their territorial waters. By operating under invitation, it was deemed that the chances of nourishing local conservation efforts would be maximized. Further to this aim, every effort was made to include local scientists, managers, educators, as well as representatives from not-for-profit organizations, as part of the shipboard party. Many of these efforts have been described in high-level post-cruise reports (e.g., Bruckner et al. 2016;Purkis et al. 2017Purkis et al. , 2018, documentary films (e.g., Barrat 2013Barrat , 2014Barrat , 2015Barrat , 2016, and scientific papers (see bibliography at www.livingoceansfoundation.org/pub lications/scientific-articles-accessed 11/27/2018). With the exception of the Red Sea surveys (Bruckner et al. 2011), however, little has been published on the methods used or accuracy of the KSLOF-GRE remote sensing products. Such is the overall purpose of this paper. The KSLOF-GRE employed remotely sensed imagery along with contemporaneous field data to produce both Fig. 1 a-j Location of the sites visited on the Khaled bin Sultan Living Oceans Foundation Global Reef Expedition where habitat and bathymetric maps were produced. Red polygons emphasize extent of mapping and encompass a total area of 65,000 sq. km of habitat situated shallower than 25 m water depth. Accompanying site names in red also. GBB in e abbreviates Great Bahama Bank. North is top in all maps; scales as noted habitat maps and bathymetry at 2 m spatial resolution over an area of 65,000 sq. km. This area corresponds to coral reefs found in water shallower than 25 m. The completed work reflects the remarkable increase in accuracy of satellite-derived reef maps over the past 20 yrs and represents an important milestone toward mapping all of Earth's reefs at meter-scale spatial resolution.
Early remote sensing studies of coral reefs used government-operated sensors such as Landsat (e.g., Ahmad and Neil 1994;Andréfouët et al. 2001;Purkis et al. 2002;Naseer and Hatcher 2004) or SPOT (Satellite Pour l'Observation de la Terre; Loubersac et al. 1991;Capolsini et al. 2003). Although the 20-30 m pixel sizes of those sensors were considered high resolution for the time and were adequate to give the gist of seabed character, results using those data were incapable of capturing the heterogeneity of a typical reef environment as it would be experienced in situ. This shortcoming was largely overcome with the launch of IKONOS with 4 9 4 m pixels in the visible spectrum in 2000, followed by QuickBird with 2.4 9 2.4 m pixels in 2001, both of which were swiftly assigned to mapping reefs, albeit over areas of just a few hundred square kilometers Purkis, 2005;Hernández-Cruz et al. 2006;Purkis et al. 2006;Rowlands et al. 2008). More recent progress has been incremental and led by the WorldView series of satellites which offer visible-spectrum spatial resolutions between 1 and 2 m, two orders of magnitude smaller pixels than Landsat or SPOT of a generation ago. Beyond enhanced spatial resolution, the WorldView program delivers data in eight spectral bands, of which five are water penetrating, facilitating improved separation of seabed types and more accurate bathymetry derivation (Collin and Hench 2012;Goodman et al. 2013;Roelfsema et al. 2014Roelfsema et al. , 2018Glynn et al. 2015;Hedley et al. 2016;Warren et al. 2016;Kerr and Purkis 2018;Purkis 2018).
Although orbital sensors are now able to image Earth with spatial and spectral resolutions that could only be achieved using aircraft a decade ago, it is only recently that reef-mapping projects have started to tackle regional scales, as opposed to individual or small collections of reefs at specific locations. Computational limitations somewhat explain this local focus; regional image datasets are large and laborious to classify, not to mention complicated by broad variations in environmental conditions, such as tides, waves, and water clarity. The lack of research funds allocated to global mapping projects is an equal culprit, however. A small number of reef-mapping programs have, nevertheless, made the important jump to regional audits suitable for countrywide marine spatial planning initiatives. These programs can be categorized as to whether they deliver maps at Landsat-type resolutionthat is, with minimum mapping units (MMUs) measured in hundreds of square meters-versus WorldView-type resolution, which is at least one order of magnitude finer.
There have been two regional-scale reef-mapping programs in the first category (Landsat-scale spatial resolution). First was the Biogeography Reef-Mapping Program of the US National Oceanic and Atmospheric Administration (NOAA) which was limited to US territorial waters and tendered at the minimum mapping unit of 1000 m 2 (Monaco et al. 2012). Second, and of similar resolution because it was developed from Landsat imagery, was the near-global reef database compiled by the United Nations Environmental Programme-World Conservation Monitoring Center (UNEP-WCMC). The Millennium Coral Reef-Mapping Project (Andréfouët et al. 2006) was the dominant constituent here.
There have also been two regional-scale reef-mapping programs in the second category (meter-scale spatial resolution). The first program to deliver meter-resolution reef maps at regional scales was the KSLOF-GRE and is the focus of this paper. The second program in this category is a recently launched initiative called the Allen Coral Atlas (www.allencoralatlas.com-accessed 03/07/2019). This endeavor aims for global coverage, is funded by Microsoft co-founder and philanthropist Paul Allen, and is utilizing satellite imagery provided by Planet Labs. The methods employed are based on work by Roelfsema et al. (2018) in the Great Barrier Reef. To date, the Allen Coral Atlas has completed various reef globally, including Heron Island (Australia), Karimunjawa (Indonesia), Kayankerni Reef (Sri Lanka), Moorea (French Polynesia), Lighthouse Reef (Belize), and West Hawai'i (USA.). All these data are accessible through the project's Web portal-www.allen coralatlas.com/atlas (accessed 03/11/2019). The KSLOF-GRE builds forward from earlier programs by providing coverage of much of Earth's major reef provinces, as accomplished by the Millennium Mapping Program, but at meter scale. Furthermore, the work conducted by KSLOF through the GRE might be considered as complimentary to new initiatives, such as the Allen Coral Atlas, by providing insight as to how field assessments and benthic habitat mapping can be scaled regionally.

Global Reef Expedition
The KSLOF-GRE simultaneously examined reef geomorphology, habitat, and satellite-derived bathymetry. Bathymetry is traditionally partnered with habitat maps because it has been demonstrated to hold predictive power over several ecologically important aspects of the reef system, such as the use of rugosity to forecast the diversity and biomass of reef fish Mellin et al. 2009;Knudby et al. 2011). Bathymetry maps are also the primary Coral Reefs (2019) 38:467-488 469 input to the calculation of local hydrodynamic exposure (e.g., Purkis et al. 2012a, b;Callaghan et al. 2015). The underlying philosophy of the KSLOF-GRE was to marry meter-resolution remote sensing and ground verification with traditional field surveys of reefs, such as those which have been continuously developed since 1998 by the Atlantic and Gulf Rapid Reef Assessment (AGRRA) Program (Ginsburg et al. 2000). AGGRA, and its expansive partner network, provides a standardized assessment of key structural and functional indicators that can be applied to reveal spatial and temporal patterns of reef condition. Based on a modified version of the AGRRA protocols, the KSLOF-GRE used SCUBA-diver surveys to systematically collect data at multiple depths for all visited sites to quantify, at a minimum, live coral and algal cover, as well as reef fish biomass and diversity.
From 2006 through 2009, the KSLOF-GRE operated along the Red Sea coastline of the Kingdom of Saudi Arabia during which four cruises were accomplished for the purpose of developing and refining the field and remote sensing protocols which would later be deployed globally. In this initial phase, reef geomorphology as it pertains to reef resilience was examined (Hamylton 2011;Riegl et al. 2012;Rowlands et al. 2014Rowlands et al. , 2016Rowlands and Purkis 2015) alongside the sedimentology and Pleistocene development of the Red Sea (Purkis et al. 2010. Approximately 32,000 sq. km of shallow-water (\ 25 m water depth) reef habitat and bathymetry were mapped from satellite and aircraft data (Fig. 1f), work which was summarized in a marine atlas of the Red Sea (Bruckner et al. 2011), a format which would later also be used to disseminate geographic products for the KSLOF-GRE. All of these data can be viewed on the interactive Living Oceans Foundation GIS Data Portal (https://maps.lof.org/ lof-accessed 11/7/2018). Since it has already been published on extensively, the Red Sea component of the GRE is not the focus of this study and will not be considered further.
KSLOF-GRE surveys from 2011 to 2015 used a standardized survey protocol to collect baseline data on reef extent, habitat distribution, and health using a combination of diver, satellite, and other observations. The aims of this paper are fourfold: 1. To emphasize the economies of scale that can be achieved by object-based interpretation of Digi-talGlobe satellite data. 2. To highlight trends and patterns in error for the delineation of benthic habitats from orbit across diverse reef geomorphologies, seafloor types, water depths, and environmental settings.
3. To initiate a public repository of coral reef maps generated at appropriate scales to support regionalscale marine spatial planning initiatives. 4. To promote awareness of the KSLOF-GRE map products and initiate their widespread usage by the community.

Methods
Diver surveys for training data The field component of the KSLOF-GRE was conducted between 2006 and 2015, and Table 1 provides an overview of the quantity of data acquired by country visited. For each of the 1000 individual reefs visited, the benthic cover of major functional groups and substrate type were assessed along 10 m transects using both diver-recorded observations, point-intercept counts, and photographic assessments. A minimum of four transects were completed at each dive site, and surveys were completed at 25, 20, 15, 10, and 5 m water depths. Via these methods, the following parameters were quantified: corals identified to genus, other sessile invertebrates identified to phylum or class, and six functional groups of algae. Reef fish surveys were also conducted at each dive site at depths stratified between 5 and 20 m via visual census as described by English et al. (1997). For more detail, the reader is directed to the Foundation's Field and Final Country Reports which are available online (www.livingoceansfoundation.org/publica tions/final-reports/-accessed 03/11/2019). These reports contain exhaustive lists of all sites and survey protocols. Diver-collected data were used to aid in the definition of map classes, as described in ''Definition of habitat map classes'' section. In addition, the dominant habitat type was extracted from each dive site to serve as labeling (training) data for satellite mapping, as described in ''Level 4 biological cover and Level 5 habitat maps'' section. A total of 1240 dive sites across the KSLOF-GRE were treated in this way. This number of dives equates to approximately 15,000 hours of underwater data collection achieved by the [ 200 scientists involved in the expedition.

Small-vessel surveys for training and validation data
A small vessel was used to collect several datasets at each of the * 1000 visited reef sites. A total of 30-million tidecorrected single-beam sonar soundings were acquired throughout the KSLOF-GRE. These measurements were used to create the bathymetry maps (''Satellite-derived bathymetry maps'' section). Additional ground-truth data were collected in the form of 2000 surficial sediment samples and 150 linear km of low-fold subbottom geophysical profiles obtained with a 5 kHz SyQwest Stratabox subbottom profiler-protocols for each detailed by Purkis et al. (2014). These datasets were used in conjunction with the diver data just described to help define habitat classes and segment geomorphological structures (''Definition of habitat map classes'' and ''Development of habitat maps'' sections).
A total of 11,000 seabed videos were captured across all sites via a tethered SeaViewer 'drop' camera integrated with a differential global positioning system (dGPS). This video system allowed seabed observations to be obtained from the intertidal to approximately 50 m water depth at a frequency far exceeding that achievable via SCUBA. The drop camera videos were analyzed in the laboratory and used for map validation (''Accuracy assessment'' section).

WorldView-2 satellite imagery
The KSLOF-GRE employed the DigitalGlobe Inc. WorldView-2 (WV2) satellite to image each visited reef site. The instrument images in eight multispectral bands with pixel widths of 1.85 m for images acquired with look angles \ 20°off-nadir coarsen to 2.07 m for look angles exceeding 20°. Pixel brightness values are digitally encoded with 11-bit radiometric resolution. WV2 is particularly adept at imaging the shallow seabed since five of the eight spectral bands are of sufficiently short wavelength to have meaningful penetration in water-these five are the coastal blue band (400-450 nm), blue (450-510 nm), green (510-580 nm), yellow (585-625 nm), and red (630 -690 nm). Experience across the KSLOF-GRE suggested that under ideal conditions, the seabed could routinely be imaged for habitat mapping down to water depths of 25 m.
The tropics are often cloudy and therefore challenging to image. To address this difficulty, at least eight months prior to each of the 15 field missions, the WV2 was tasked to acquire imagery at look angles\ 15°off-nadir to minimize sun glint. At 1 month prior, all acquired imagery was purchased from DigitalGlobe Inc. and assembled to support mission planning and subsequent fieldwork. If insufficient cloud-free data had been obtained for mapping a given country, the sensor was tasked for an additional two-month post-cruise, to fill areas that remained stubbornly cloud contaminated. In this way, the majority of imagery was acquired within four months of fieldwork, but with a maximum differential of eight months. For large sites, such as the 6000 sq. km Cay Sal Bank (Bahamas- Fig. 1e), up to 50 individual WV2 acquisitions were assembled to deliver an image mosaic with \ 3% cloud cover, which was the threshold deemed as the maximum tolerable for mapping. In many cases, cloud cover was further reduced by replacing individual cloud-contaminated areas with a portion of a cloud-free acquisition from an alternative date and equivalent tidal state, a process termed 'cloud patching.' Adjacent image scenes were selected to have a similar tidal state and equivalent water clarity.
Prior to mosaicking the individual scenes, each was processed to units of above-water remote sensing reflectance, which encompasses radiometric, solar geometry, and atmospheric correction, as described in detail by Kerr and Purkis (2018) and corrected for sun glint following Hedley et al. (2005). At this point, the processed satellite scenes were stitched into a mosaic using the image-processing software ENVI (v. 5.4, Harris Geospatial Inc.), emergent areas identified using a threshold in the 860-1040 nm spectral band, and areas of deepwater identified also, defined as having \ 5% reflectance in the 450-510 nm band (Fig. 2a). The remainder of the imagery was considered as potentially containing shallow-water habitat, defined as \ 25 m water depth, and was passed forward to the mapping workflow ( Fig. 2b-d).

Satellite-derived bathymetry maps
Bathymetry maps were derived for all the KSLOF-GRE sites via spectral derivation of water depth from WV2 satellite imagery (workflow detailed in Fig. 2b). These products served as stand-alone data layers, but were also utilized in the habitat-mapping workflow to partition each reef site into zones, which in turn were populated with a zone-specific suite of habitat classes. Stumpf et al. (2003) offer the most widely adopted empirical algorithm for extracting bathymetry from multispectral imagery. This solution uses a ratio of reflectance from two spectral bands which is tuned against known water depths to yield a bathymetry map. Motivated by the fact that this method does not exploit all five water-penetrating bands of WV2 and its successors, Kerr and Purkis (2018) evolved the algorithm via multi-linear regression of five bands, a solution which provided enhanced estimates of water depth. Their algorithm allowed viable bathymetric models Fig. 2 Workflow for the production of bathymetry and benthic habitat maps. a Image preparation encompassed correction for solar, radiometric and atmospheric effects, and, if required, correction for sun glint also. Sites imaged by multiple satellite scenes were stitched into a seamless mosaic once these corrections had been implemented. The resulting mosaic was processed to yield a bathymetry map which was calibrated by sonar depth soundings (b). Via manipulation in eCognition software, the bathymetry map was used to create a map of reef zonation which was then combined with the multispectral image mosaic via hierarchical classification to yield a map of seabed habitat (c)-text for details. Finally, the accuracy of the habitat map was computed with reference to drop camera videos acquired in the field (d) to be derived even in cases where ground truth via sonar was limited, and, under ideal conditions, even absent. Mapping of water depth for the KSLOF-GRE sites followed the Kerr and Purkis (2018) methodology and was calibrated by the sonar soundings described in ''Smallvessel surveys for training and validation data'' section. Bathymetry maps were masked below the 25-m-depth contour, as derived from sonar soundings.

Definition of habitat map classes
The KSLOF mapping endeavor built forward from two noteworthy regional-scale programs; the Millennium Coral Reef-Mapping Project (Andréfouët et al. 2006) and the NOAA Biogeography Reef-Mapping Program (Monaco et al. 2012). Although our habitat map classes differed from these predecessors, we adopted a hierarchical scheme which allows for cross-comparison (as also done by Roelfsema et al. 2018). The Landsat-derived maps of Andréfouët et al. (2006) delineated reef geomorphology, not habitat, though it was implied in many cases. For instance, class 'fore reef' in the Andréfouët et al. (2006) scheme describes location within the benthic system but, importantly, does not address substrate or cover type at that location. A fore reef environment can reasonably be anticipated to be coral-dominated, however. The NOAA effort (Monaco et al. 2012) also captured geomorphology, termed 'structure' in their nomenclature, but developed two additional map layers, 'biological cover' and 'geographic zone.' The former described dominant biota (e.g., live coral, seagrass, etc.), whereas the latter referred to the location of the benthic community within the system (e.g., reef crest, back reef, etc.). Unlike NOAA, the KSLOF-GRE products do not provide three map layers for each area, but the classification scheme was hierarchically arranged such that geomorphological structure, geographic zone, and biological cover can be separated if required. As described in ''Development of habitat maps'' section, this cross-compatibility is implicit to the way that the maps are created; a bathymetric map was initially interpreted into geographic zones (termed the 'Level 1' output), which was subsequently populated with increasing detail of geomorphological structure (Levels 2 and 3), before addition of biological cover recorded in situ (Level 4), to produce a final homologated Level 5 'habitat' map in which zone, structure, and cover are aggregated.  Table 2. North is top Table 2 Hierarchy of habitat classes utilized for satellite mapping For each site, Level 1 classes were developed from a bathymetric map, Levels 2 and 3 were generated via hierarchical classification within eCognition software and Level 4 constituted in situ seabed observations of biological cover. Homologating Levels 1 through 4 delivered the KSLOF 36-class scheme (Level 5). Colors for the Level 5 classes correspond to those used in the example habitat maps (Figs. 3,4,5,6). Map accuracy was assessed against 16 consolidated classes for all non-Atlantic sites and seven classes for Atlantic sites (final column). Consolidation of the Level 5 classes was required for accuracy assessment because of the reduced level of detail that can be extracted from underwater video-text for details. Each consolidated class was given a class ID which spans 1 through 16 for the non-Atlantic sites (emphasized in red, final column) and 1 through 7 for the Atlantic classes (black). These IDs were used in the error matrices ( The combination of reef zone, geomorphological structure, and biotic cover resulted in 36 habitat classes used across the Red Sea, Pacific, and Indian Oceans (Table 2). In the Atlantic, the same scheme was used, but not all combinations of zone, geomorphology, and cover were found in this ocean basin; only 25 of the classes were represented in the Atlantic maps. For example, there was no difference defined between 'lagoon' and 'back reef' in the Atlantic sites visited by the GRE. The description of these classes should make intuitive sense based on their zone, structure, and cover (Table 2), but there are also lengthy descriptions and example photographs for each class in the field reports previously published by KSLOF (see, for example, Bruckner et al. 2016).

Development of habitat maps
The KSLOF-GRE used eCognition software (v. 5.2, Trimble Inc.) to segment the WV2 imagery into polygons that were then labeled by zone, structure, and ultimately habitat class. In contrast to pixel-based classifiers, which assign image pixels to map classes based on their spectral content (Purkis and Klemas 2011), eCognition follows an object-based approach (Knudby et al. 2011;Phinn et al. 2012;Purkis et al. 2012aPurkis et al. , b, 2014Roelfsema et al. 2013Roelfsema et al. , 2014Roelfsema et al. , 2018Zhang et al. 2013;Warren et al. 2016). In a workflow termed 'hierarchical classification,' edgedetection routines are used to segment imagery into eCognition 'objects,' which are precincts of the image set with similar spectral and/or textural attributes. These objects are subsequently assigned into one of several map classes based on rules which consider spectral/textural signatures, shape, and contextual relationships with surrounding classes.
Whereas recent progress has been made to automate the assignment of objects to map classes, such as by Saul and Purkis (2015) using multinomial logistic discrete choice models, we found the accuracy of the automated assignments to be consistently lower than that delivered manually by an expert user. For this reason, we elected to use manual assignment of eCognition objects to map classes in our workflow for the production of the KSLOF-GRE habitat maps (Fig. 2). The workflow required four steps to handle preprocessing of the satellite imagery, derivation of a bathymetry map, development of a habitat map, and accuracy assessment (Fig. 2a-d, respectively). This section deals solely with developing the habitat map (Fig. 2c); the other three steps are described in their corresponding sections.

Level 1 zone map
A Level 1 zone map for each site was created using eCognition by applying a multi-resolution segmentation algorithm to the bathymetry map. This algorithm, because it was described in detail by Baatz and Schäpe (2000), will only be treated briefly here. The general concept of multiscale image segmentation is to subdivide an image set into objects with spectral and/or textural homogeneity. The solution proposed by Baatz and Schäpe (2000) considers this task an optimization problem. In the first step, every image pixel is considered a separate image object. Each object is then visited iteratively and merged with its neighbors to form larger (multi-pixel) objects. With each iteration, the merging decision is based on local homogeneity criteria describing the similarity of adjacent image objects. In a process similar to the annealing function described by Purkis et al. (2012b), a cost function is tracked as each merge is conducted and objects cease to be further amalgamated at the point that the function ceases to reduce.
To create the Level 1 map, the multi-resolution segmentation was deployed on the bathymetry map, which has pixel values enumerating water depth. Once segmented, an expert user manually grouped the resulting image objects that correspond to five reef zones (lagoon, back reef, fore reef, reef crest, and shelf), plus two zones encompassing terrestrial areas (land and intertidal), and deep ocean ( Table 2). The upshot of this process was a Level 1 zone map.

Levels 2 and 3 geomorphological structure maps
The next step toward the final habitat map was the delineation of geomorphological zones which were first crudely defined (Level 2) and then refined in more detail (Level 3). For the Level 2 map, the inputs were (a) the Level 1 reef zone map produced from bathymetry and (b) the multispectral WV2 image mosaic. First, for each Level 1 zone, the multispectral imagery was segmented via the multiresolution method of Baatz and Schäpe (2000). Second, in a process termed 'labelling' and with reference to the surficial sediment samples and geophysical profiles acquired in the field, the expert user manually selected image objects and attributed them as belonging to one of the three Level 2 geomorphology classes (unconsolidated sediment, coral reef and hardbottom, or other; see Table 2). b Fig. 4 Satellite-derived map products for Gizo Island, Solomon Islands. a Location of Gizo Island in the New Georgia Group. b Enhanced true-color WorldView-2 (WV2) image of Gizo and surrounding reef systems. c Bathymetry map created via spectral derivation from the WV2 imagery calibrated by in situ sonar soundings. d Corresponding habitat map developed via object-based mapping in eCognition software. Colors in d correspond to those for the Level 5 map classes in Table 2. North as indicated in a Table 3 Error matrices for the consolidated classes used to assess the classification accuracy for Atlantic (a) and non-Atlantic sites (b) Columns capture the ground-truth data and rows the classified data. The integers in parentheses after the consolidated class names (black for Atlantic, red Third, based on these user-defined training sets for each Level 2 class in each Level 1 zone, eCognition was used to classify all of the objects in the image set into geomorphological structures based on spectral, textural, and neighborhood parameters. The upshot of this process was a map of major geomorphological structure primarily split into unconsolidated sediment-dominated areas (spectrally bright and texturally homogeneous) and coral reef and hardbottom-dominated areas (spectrally dark and texturally heterogeneous). Note that by conducting this process independently within each Level 1 zone, the bias introduced by varying bathymetry across the satellite imagery was mitigated by the fact that each zone occupies a limited range of water depths. This was important because the rapid attenuation of light by water tends to override the subtle spectral differences between reef habitats (e.g., Purkis 2005).
The detailed geomorphological structure maps (Level 3) were produced in the same way as in the preceding step, but the imagery was re-segmented on the basis of the Level 2 classes and for each, the expert user applied labels for the 11 Level 3 classes defining seabed character (mud, sand, rock, etc.; see Table 2) and in the case of reefs, their morphological type (pinnacle versus aggregate, etc.; see Table 2). The advantage of conducting this segmentation based on the Level 2 classes was a radical reduction in computational overhead since subsets of the overall image mosaic were segmented separately. As before, the user manually developed these labels with reference to known points on the ground visited during fieldwork, and, again, eCognition was used to classify the unlabeled image objects based on their similarity to the training set.  Table 2. North is top Fig. 6 Satellite-derived map products for a series of isolated reef platforms offshore Île des Pins. a Location of Île des Pins, part of the Île Loyauté, New Caledonia. b Enhanced true-color WorldView-2 (WV2) image of the reef complex. c Bathymetry map created via spectral derivation from the WV2 imagery calibrated by in situ sonar soundings. d Corresponding habitat map developed via object-based mapping in eCognition software. Colors in d correspond to those for the Level 5 map classes in Table 2. North is top   Fig. 9 Relationship between habitat map accuracy and map area (a) and map complexity (b). Map accuracy was quantified by the Tau coefficient and complexity via the proportion of edge pixels (text for details). In both plots, the (linear) correlation was computed for Atlantic (blue broken line) and non-Atlantic sites (brown). No correlation was observed between map accuracy and area for either grouping of sites. Only low correlation (R 2 = 0.35) exists between accuracy and map complexity for non-Atlantic sites. These variables are highly correlated for the Atlantic maps (R 2 = 0.99), meanwhile, but with only three sites in this ocean basin, the relationship should be not be overly emphasized. Site abbreviations in plots as follows: Galap. = Galápagos, Inag. and HS = Inagua and Hogsty, N. Cal = New Caledonia, Solo. = Solomon Isl., Aust. = Austral Isl., Gam. = Gambier

Level 4 biological cover and Level 5 habitat maps
In the final step in the mapping workflow, field observations of biological cover (termed 'Level 4' data) were convolved with the geomorphology map to yield a map of habitat (e.g., Figs. 3d, 4d, 5d, 6d). This step was again achieved via application of the multi-resolution segmentation algorithm (Baatz and Schäpe 2000), but this time, the Level 3 classes were individually segmented and, again with reference to field data, the user manually selected labels for objects characterized by the 12 Level 4 classes of biological cover. Again, eCognition was used to classify the remaining unattributed image objects on the basis of similarity to the training set. As laid out in Table 2, each object, now classified according to benthic cover, was attributed with the addition of its zone and geomorphological structure, which varied by location within the image set, as defined by the previously created Level 1 and 2 maps, respectively. As an example, an image object describing a patch reef in the lagoon would be reattributed as 'Lagoon-Patch Reef,' and so on. This reattribution process delivered the 36 'aggregate classes' of the final habitat map. To complete the map, boundaries existing between image objects of the same class were dissolved such that areas of a single habitat type were encompassed by a single polygon. At this stage, and again with reference to the diver observations, the evolving map was examined by an expert user and any obvious errors corrected in a process termed 'contextual editing' (as originally proposed by Mumby et al. 1998). To complete the process, the finished habitat map was exported as an ESRI shapefile for further analysis in a geographic information system (GIS).

Accuracy assessment
Accuracy assessment of the habitat maps was conducted using error matrices (Story and Congalton 1986;Congalton 1991) with reference to the 5106 dGPS-positioned seabed videos captured using a tethered 'drop' camera that remained independent from the map-making workflow. These drop camera videos had three advantages for the purposes of accuracy assessment: a large sample size, wide, consistent coverage across the entirety of the GRE sites, and independence from the training/labeling process of map creation. The drop camera dataset suffered a few limitations as well. First, some habitat types were undersampled due to physical constraints navigating the vessel. Second, the limited field-of-view of the camera created difficulties discriminating certain habitat classes. Third, there was some geographic uncertainty in camera location due to the tether length and the horizontal field-of-view. The first two of these limitations were addressed by eliminating or consolidating certain map classes for the purposes of accuracy assessment. The third was addressed by considering the neighborhood around each drop camera point using a technique we call 'lagged accuracy.' It is important to emphasize that field-operation logistical planning helped reduce these uncertainties by accounting for wind, as well as current magnitude and direction, when deploying the camera and capitalizing on precise boat handling techniques by the highly skilled skipper. This allowed us to accurately position and 'fly' the tethered camera over each habitat sampled. Terrestrial habitat classes were impossible to sample with the drop camera, for obvious reasons. Thus, the accuracy of terrestrial habitat classes was not quantitatively assessed for these maps. Nevertheless, we assume that the maps are very accurate for a consolidated 'terrestrial' class (i.e., consolidated map Class #1; Table 2), since segmenting land versus marine habitats is straightforward with the infrared channels of satellite imagery. Intertidal and reef crest classes also proved difficult to sample, due to their extremely shallow depths at the islands surveyed. Only two reef crest videos and no intertidal videos were captured. Thus, the accuracy for intertidal classes was not assessed and fore reef crest was insufficiently sampled to draw strong conclusions. To put this limitation in perspective, however, intertidal and reef crest classes were each found to have \ 1% of the total number of classified pixels (Fig. 7). Therefore, their omission from the accuracy assessment is unlikely to change overall conclusions about classifier performance.
The limited field-of-view of the drop camera prevented the discrimination of many of the fine details between Level 5 classes ( Table 2). For instance, the videos were adequate to classify the seabed in general as a 'Lagoonal Reef,' but the field-of-view was inadequate to resolve whether a given lagoonal reef was only 10 m in diameter, or smaller, which would correspond to the Level 5 map class 'Lagoon-Coral Bommies,' versus a much larger patch, which would be a Level 5 'Lagoon-Pinnacle Reef.' To compensate for this discrepancy in scale between the satellite data and the ground-truth data, we grouped the 36 Level 5 classes into a smaller number of 'consolidated classes' (Table 2). For most sites around the world, 16 consolidated classes were used, reflecting different combinations of geographic zone and substrate. In the Atlantic, however, geographic zone was not as easy to define, so additional classes were consolidated in the Atlantic, reducing the total to seven for those sites.
Overall, producer's and user's map accuracies were computed for each site using the consolidated classes (Table 2) via the error matrix approach (Story and Congalton 1986). In addition, the Kappa (Congalton 1991) and Tau (Ma and Redmond 1995) coefficients were computed to quantify the degree to which the accuracy of each map was better than random chance. Equal prior probability was used for calculating Tau because no a priori information on class probability was used in the hierarchical segmentation. It should be noted that the accuracies quoted in the error matrices (Table 3) are for the consolidated classes and cannot be extrapolated to speak to the accuracy of the individual classes prior to their consolidation.
Map accuracy as assessed via standard error matrices does not allow for geographic offsets between the habitat map and reference data. Such offsets are often unavoidable, however, and stem from the many vagaries of setting an exact position on the ocean during fieldwork. Sources of positional error include GPS inaccuracies, diver observations not made exactly beneath the position recorded when entering the water, and the drift of the tethered 'drop' camera away from the boat. These offsets might reasonably be expected to routinely exceed the 1.85 m pixel width of the WorldView satellite, with the result that the groundtruth data are not perfectly registered with the habitat map. Furthermore, with a horizontal field-of-view, as was the case with the drop camera used for this study, the video data are directional, which can have just as great an impact as positional uncertainty. Imagine the camera positioned on an edge between two classes; the class assigned to that ground-truth point would depend on the direction in which the camera was orientated, even if the camera position did not change. Whereas such offsets might legitimately be considered as inaccuracies for habitat maps produced on a local scale, we consider them to be acceptable when mapping across hundreds of thousands of sq. km of Earth's remotest reef systems. Thus, we wanted a way to assess accuracy that would account for uncertainty in the relative position of ground truth to satellite data.
To sensibly address geographic uncertainty in camera location, we used a metric called 'lagged accuracy' which, for each ground-truth point, collates the cumulative probability of encountering pixels mapped as the same class attributed to the assessment point, for lag distances between 0 and 300 m offset from that point, in all Cartesian directions. Providing that map pixels of the class sought exist within the specified lag distance around the accuracy assessment point, the cumulative probability of encounter will rise as a function of increasing lag, with the rate of that rise dictated by the density of pixels of that class in the queried portion of the habitat map. Of course, it is unreasonable to take the existence of a map pixel within the search radius with the same class assignment of that of the accuracy assessment point to justify scoring the map as accurate. Indeed, for large lag distances, the correct assignment will be recognized even if the queried habitat map is random. Hence, it is necessary to set sensible thresholds in both lag distance and cumulative probability that might rationally indicate that a positioning error has precluded an exact match at the location of the accuracy assessment point. While there are no precise answers, we felt 25 m was a sensible threshold for the accuracy of a regional-scale map used to support marine spatial planning.

Patterns of habitat classification accuracy
To explore possible causes of habitat map error, the persite Tau coefficients (a measure of map accuracy) were cross-plotted against mapped area and the complexity of the habitat maps (Fig. 9). Doubtless, the KSLOF-GRE dataset allows for all sorts of analysis of the spatial patterns among reef systems around the globe, and it is our hope that it will be used for such in the future. The present goal, however, was simply to check for broad and systematic patterns related to habitat classification accuracy.
There are many ways to quantify scene complexity. One of the simplest is to count the proportion of edge pixels, i.e., those which border a different class. Edge pixels are good metrics for assessing habitat classification because they are affected by a combination of class variety, spatial arrangement, and pixel mixing (Heydari and Mountrakis 2018). Furthermore, accuracy has sometimes been shown to decrease with increasing proportion of edge pixels (Heydari and Mountrakis 2018). Checking whether this pattern held for the KSLOF-GRE was valuable because datasets with a sufficient number of scenes with varying complexity to test this are rare.

Results
The KSLOF generated benthic habitat maps and bathymetry over a total of 65,000 sq. km during the 10 yrs of the Global Reef Expedition. Examples of these products are reproduced here for the Pacific and Indian Oceans (Figs. 3,4,5,6). The entire digital dataset can be explored on the KSLOF GIS data portal (https://maps.lof.org/lof-accessed 11/7/2018). As demonstrated by Figs. 3, 4, 5, and 6, the object-based workflow used for mapping seabed character scaled to tens of thousands of sq. km, while maintaining high spatial fidelity.
Sites in the Atlantic were comprised of \ 5% reef habitat, nearly 25% hardbottom, and 29% as seagrass (Fig. 7). Non-Atlantic sites, by contrast, contained \ 1% seagrass and [ 15% reef habitat. Even though biogeographic differences in benthic character were not the subject of this paper, these statistics underline the diversity of sites mapped using a common workflow across the KSLOF-GRE.
A quantitative assessment of classification error (Table 3) revealed the overall accuracy of the maps developed for the three Atlantic sites to be approximately 10% lower than that for the ten non-Atlantic sites (81% vs. 90%, respectively). The Kappa and Tau coefficients differed by nearly 20% between the Atlantic and non-Atlantic sites, however. Kappa and Tau both penalize the Atlantic results more than the other sites because of the fewer number of consolidated classes used to conduct the accuracy assessment in the Atlantic (4) versus elsewhere (9).
The producer's accuracies for the habitat classes in the Atlantic maps were approximately 70%, save for 'Hardbottom,' which was 94% (Table 3). Much of this discrepancy arose from confusion between macroalgal stands and seagrass meadows which occupied large swaths of the 6000 sq. km Cay Sal Bank (Fig. 1e). The user's accuracies for the Atlantic sites were generally higher than the producer's accuracies, ranging from 76% (Sediment with Macroalgae) to 89% (Reef).
Both producer's and user's accuracies for the 13 non-Atlantic sites were considerably higher than for the Atlantic. An outlier here was the class 'Reef Crest Hardbottom' (Producer's Accuracy = 50%, User's = 100%) which should be ignored as it was validated by a single ground-truth point, a function of the difficulty of safely navigating a small vessel across the shallow reef crest. Withstanding this class, the median producer's accuracy was 92% and user's accuracy was 89%, both reassuringly high values.
Cumulative distribution functions describing the lagged accuracy of the three Atlantic (Fig. 8a) and ten non-Atlantic sites (Fig. 8b) echoed the trend obvious in the error matrices (Table 3). The habitat maps for the Atlantic had lower accuracies than those developed for non-Atlantic sites. Maps created for Colombia and the Cay Sal Bank (Bahamas) were the worst performers here, with only a 75% probability of encountering a correctly classified pixel within 25 m of a ground-truth point (Fig. 8a). The maps for Great and Little Inagua, and Hogsty, were * 10% better, but still underperformed the non-Atlantic sites. For these, the probability of encountering a correctly classified pixel within 25 m of the ground-truth data exceeded 90% for the Solomon Islands, the four mapped archipelagos in French Polynesia (Australs, Gambier, Society, and Tuamotu), and the Cook Islands. The probability at the same lag distance for the remaining sites (Tonga, Fiji, New Caledonia, and Galápagos) exceeded 85%. By a lag distance of 50 m, which is likely on the upper limit of what is useful for a regional-scale map product, the probability of encountering a correctly identified pixel exceeded 90% for all sites except Colombia and the Cay Sal Bank.
Map accuracy and area were uncorrelated for both the Atlantic and non-Atlantic sites (Fig. 9a). For the latter grouping, however, a low level of correlation (R 2 = 0.35) was observed between accuracy and complexity (Fig. 9b). The correlation between these parameters was strong for the three Atlantic sites (R 2 = 0.99), but the result is unreliable because of the small number of sites in this grouping.

Discussion
Regional-scale coral reef mapping from remote sensing has a role to play in the widening portfolio of intervention measures that are being mobilized against the reef crisis. Of these measures, the establishment of large-scale marine protected areas (MPAs) has been particularly effective (Sheppard et al. 2012;Toonen et al. 2013;Wilhelm et al. 2014). The Big Ocean Network (https://bigoceanmanagers. org-accessed 10/25/2018) provides a rule of thumb as to what constitutes 'large scale,' with their 17 member sites ranging in size from approximately 150,000 sq. km to nearly 2000,000 sq. km. Beyond large size, the success of an MPA rises if it is nested within a network of ecologically connected protected areas, which as a collective encompass a full range of critical habitats Gleason et al. 2010). Designing ecological connectivity into MPA networks, however, requires careful consideration of available information on habitat distribution, larval dispersal patterns, adult movement ranges, and oceanography Gaines et al. 2003;Palumbi 2004). Marine spatial planning-MSP- (Douvere 2008) is a central component to balancing these and other considerations in the design of MPAs. Habitat and bathymetry maps lie at the base of the MSP workflow and are therefore critical to its success. With MPAs becoming ever larger, it is imperative that the mapping keeps pace.
Covering 65,000 sq. km and taking nearly a decade to complete, the KSLOF-GRE is the largest coherent reefmapping program accomplished to date. The results from the GRE are presently being used to train the next generation of image classifiers which hold the potential to deliver truly global audits of reef status through time. As developed by Chirayath and Earle (2016) and Chirayath and Li (2019), the delineation of reef habitat via machine learning holds particular promise. Further, the goal to 'map once, use many ways' underpins and justifies the KSLOF-GRE and the archiving of its outputs in the public domain. As the expedition was planned, the Foundation was agnostic to the degree to which each host nation was conducting MSP. Instead, the philosophy was to deploy the KSLOF-GRE seafloor mapping program to stimulate the creation of national and regional databases and information systems containing essential coral reef environmental data, else contribute to these if they already existed.
Accuracy of the KSLOF-GRE habitat maps, as computed with the traditional error matrix approach (Story and Congalton 1986), varied from * 70 to 90%. The technique used here of computing lagged accuracy, in addition to a traditional error matrix, resulted in cumulative probability functions for each map class (Fig. 8) which grant the map user an alternative means of judging map quality. In some ways, lagged accuracy conveyed the same information as overall accuracy computed using a traditional error matrix approach. For example, the error matrices showed that Atlantic sites had lower overall accuracy, in general, than non-Atlantic sites (Table 3). This same result was also clear from plots of lagged accuracy (Fig. 8). In other ways, however, lagged accuracy complemented the error matrix approach with new information about the spatial distribution of errors. For example, the sites with the steepest slopes of the lagged accuracy curve in the 0-25 m spatial scale had the highest fraction of edge pixels (Tuamotu, Gambier, and Society). Conversely, those with the shallowest slope over that range had the fewest edge pixels (Colombia, Cay Sal, Galápagos). Thus, information about the relationships among patchiness, scale, and accuracy was contained in the lagged accuracy cumulative distributions. This is a topic to be examined further in the future.
Lagged accuracy also provides a thematic analog to familiar specifications for spatial accuracy. Horizontal spatial data accuracy is typically reported in the following form: 'X meters (feet) horizontal accuracy at the 95% confidence level' (FGDC 1998). One could use lagged accuracy to report a thematic accuracy in an analogous form, for example: 'X% thematic overall accuracy within Y meters horizontal lag.' Users could specify an acceptable spatial lag, Y, according to their needs. A manager planning marine protected areas across one million square km of ocean might acceptably tolerate a larger Y than, say, an engineer planning a dredging operation.
All KSLOF-GRE sites except Cay Sal and Colombia were found to have at least 85% thematic overall accuracy within 25 meters horizontal lag; Cay Sal and Colombia had 75% thematic overall accuracy within 25 m horizontal lag. If faced with developing a network of MPAs across hundreds of thousands of sq. km of tropical ocean, as has already been accomplished by the 17 member sites of the Big Ocean Network, access to maps which correctly position the occurrence of critical habitats such as coral reefs, seagrass meadows, and mangrove stands to within 25 m would be of huge benefit, especially considering that most of these protected areas were defined without any precise knowledge of the locations, size, and architecture of benthic habitat throughout their range.
The fact that, regardless of which accuracy metric was used, the Atlantic sites were * 10% less accurate than those in the Pacific and Indian oceans was surprising because the diversity of seabed character is considerably lower in the Atlantic, and therefore fewer classes were used to map these sites. In contrast, several previous studies have shown that increasing the number of habitat classes decreases map accuracy (e.g., Andréfouët et al. 2003). One potential explanation for this inconsistency can be linked to the workflow used to produce the maps. Our workflow relied on a bathymetric map derived from satellite imagery to guide development of a Level 1 map of reef zones which was subsequently evolved to a Level 2 map of major geomorphic structure (Fig. 2). The sites considered outside the Atlantic were well poised for this approach as they were predominantly atolls with well-defined zones (fore reef, reef crest, lagoon, etc.). Once the seascape has been split into zones, the burden of mapping the habitats contained within them was eased because a limited number of benthic cover types were prescribed in advance for each zone (as detailed in Table 2). The Atlantic sites mapped by KSLOF-GRE were challenging to partition into reef zones, except for the diminutive Hogsty Reef, which is only 100 sq. km in area but atollic in morphology. The Level 1 zone map for the 5500 sq. km Cay Sal Bank, however, was more poorly defined because Cay Sal lacks platform-margin reefs and takes the form of a sediment-dominated flattopped carbonate bank ). Thus, the Atlantic maps placed greater emphasis on correctly ascribing benthic cover across a large area, from a large quantity of potential biotic classes. We anticipate that this disjoint in our workflow explains the reduced accuracy of the Atlantic habitat maps. This said, the 81% overall accuracy for Atlantic maps falls in line with, or exceeds comparable studies conducted at much smaller spatial scales (e.g., Phinn et al. 2012;Zhang et al. 2013;Collin et al. 2014;Hedley et al. 2016;Roelfsema et al. 2018) and is therefore not deemed to be a limiting factor.
A regression of map accuracy and area yielded no (linear) correlation in the Atlantic and at most a slight inverse relationship for non-Atlantic sites (Fig. 9a). The lack of a relationship was reassuring and indicated that the workflow was not confounded by very large sites. The positive correlation observed between map accuracy and map complexity (Fig. 9b) was counterintuitive. Rather than accuracy decreasing with increasing map complexity, as previously documented for pixel-based approaches to classification Heydari and Mountrakis 2018), our object-based approach yielded increasing accuracy for more complex benthic systems. The explanation for this result might be due to pixel-based versus object-based classifiers, or due to the zone-based classification scheme used. Note, however, that the correlation between accuracy and complexity was high for the Atlantic sites but rather modest for the non-Atlantic sites. Thus, we feel this result can be most easily interpreted as another manifestation of the lowered performance of our workflow in the Atlantic sites where reef morphology is less well developed than outside the Atlantic. Throughout our mapping endeavors, we found that some habitat types tended to be misclassified more often than others. For example, in the Atlantic, sediment with macroalgae was occasionally misclassified as hardbottom. Confusion existed too between seagrass and either sediment with macroalgae or hardbottom. Such errors are to be anticipated because of the spectral similarity of these classes (Hochberg and Atkinson 2003;Purkis 2005). For non-Atlantic sites, lagoonal sediment with macroalgae was sometimes misclassified as lagoonal reef, confusion which might variably be attributed to their spectral and textural similarity and the fact that turbidity in restricted atoll lagoons can be elevated (Kjerfve 1986). Forereef sediment with macroalgae was occasionally wrongly mapped as forereef hardbottom. In this case, the rapid downslope increase in water depth can likely be implicated as nearvertical morphology is challenging to image because of light attenuation and shadowing (Jay et al. 2017). Ways to more routinely separate live coral from macroalgae in multispectral imagery are of heightening importance given the large-scale regime shift of reefs to algal-dominated states (Graham et al. 2015;Hughes et al. 2017;Hempson et al. 2018).
The products discussed in this manuscript, and the spatial breath of coverage the KSLOF-GRE achieved, provide a major contribution to the science and management of coral reefs, and the ecosystem services they provide. Studies of this magnitude provide critical baseline data to benchmark the condition of reefs now, thereby enabling; (1) quantification of the rate and direction of future change at seascape scales, (2) enhanced understanding of how reefs should be managed to ensure their sustainability, and (3) documentation of how they change once management interventions are in place. Many countries visited during the GRE are Small Island Developing States, whose economies are wholly dependent on the submerged marine resources located within their EEZ. As such, resource management at the national level is critical to the local economies of these nations. In addition, most of the KSLOF-GRE sites are home to isolated villages of people (i.e., monthly ferry service and no airport) who are truly dependent on their adjacent reef environment for food security (Béné et al. 2016). Dissemination of the KSLOF-GRE mapping products and survey data to the host countries and communities provides the greatest opportunity for their use at both the national and local levels.
Coral reefs are icons of environmentalism because they have degraded so rapidly with causes easily linked to climate change and other human pressures. Despite iconic status, though, Earth's reefs have not been systematically mapped with the intensity of, for instance, tropical deforestation. This deficit means that even fundamental questions such as area covered by reefs globally are unknown. This includes the inability to formally assess coral reef health and status at country level. Though by no means covering every reef worldwide, the KSLOF-GRE covers a meaningful proportion of key reef provinces around the world and provides a baseline of their health prior to the 2017 mass bleaching event. According to Spalding et al. (2001), Earth's reefs cover nearly 285,000 sq. km, a figure which would suggest that the KSLOF-GRE mapping, which covers 65,000 sq. km, has characterized one-fifth of them. This proportion is tenuous, however, as the true global reef area is poorly constrained, not least because of the rarity of large-scale maps-a deficit which motivated this study. We hope that our large-scale maps will open new vistas of potential enquiry and motivate others to work toward a global reef audit. Many aspects of the reef crisis are presently intractable. We show that accurately mapping bathymetry and habitat at regional scale is not one of them.