Mapping Coral Reef Resilience Indicators Using Field and Remotely Sensed Data

In the face of increasing climate-related impacts on coral reefs, the integration of ecosystem resilience into marine conservation planning has become a priority. One strategy, including resilient areas in marine protected area (MPA) networks, relies on information on the spatial distribution of resilience. We assess the ability to model and map six indicators of coral reef resilience—stress-tolerant coral taxa, coral generic diversity, fish herbivore biomass, fish herbivore functional group richness, density of juvenile corals and the cover of live coral and crustose coralline algae. We use high spatial resolution satellite data to derive environmental predictors and use these in random forest models, with field observations, to predict resilience indicator values at unsampled locations. Predictions are compared with those obtained from universal kriging and from a baseline model. Prediction errors are estimated using cross-validation, and the ability to map each resilience indicator is quantified as the percentage reduction in prediction error compared to the baseline model. Results are most promising (percentage reduction = 18.3%) for mapping the cover of live coral and crustose coralline algae and least promising (percentage reduction = 0%) for coral diversity. Our study has demonstrated one approach to map indicators of coral reef resilience. In the context of MPA network planning, the potential to consider reef resilience in addition to habitat and feature representation in decision-support software now exists, allowing planners OPEN ACCESS Remote Sens. 2013, 5 1312 to integrate aspects of reef resilience in MPA network development.


Introduction
Throughout the past few decades, severe climate-driven threats to coral reefs have been identified, including increased frequency and severity of acute disturbances, such as storms, floods and mass coral bleaching from elevated water temperatures [1,2], as well as chronic environmental changes, such as ocean acidification [3,4].The effects of climate change combine with more direct human impacts, like overfishing and land-based runoff, to critically threaten coral reefs globally.Strategies for management of reef ecosystems often rely on the creation of marine protected areas (MPAs) that limit the direct human impacts, while strategies to also take climate-related threats into account have only recently been introduced [5][6][7].In coral reef management, the "resilience" of reef ecosystems has emerged as a key concept.The term "resilience" is often considered by reef scientists and managers to be a composite of "resistance", the ability to maintain structure and function while under stress, and "recovery", the ability to restore structure and function following a temporary disturbance [8,9], a definition we adopt in this study.This simplistic definition is useful for application of the resilience concept, though we recognize that it ignores complex and interactive ecosystem dynamics.Given the expected increase in climate-related effects on coral reefs, there is an urgent need to measure and map both components of resilience, as well as the level of coral reef exposure to climate-driven stresses [10].When mapped, data on reef exposure and resilience can be used in MPA network planning, for example, by including reef areas with low exposure and high resilience as targets for protection.
Based on global data sets, Maina et al. [10,11] developed a spatially explicit global coral reef exposure model that identifies regional hotspots where coral reefs are vulnerable to climate-related stressors and quantifies local stress-reinforcing and stress-reducing factors.The data can be used to assess cumulative impact, as well as change and variability of exposure, all of which may affect corals' abilities to cope with climate change.However, the coarse spatial scale of the model (gridded cells of approximately 21 km 2 ) presents real limitations, because its cell size is larger than the majority of reef features [12], so each model cell may contain several reef features, even several reefs, for which exposure values would be averaged.The large cell size also impacts this model's utility in designing MPA networks, particularly in areas like the Western Pacific, where many MPAs are small (median sizes ranging between 0.1 km 2 in Samoa to 11.8 km 2 in Tuvalu) [13,14], and the planning region or tenure boundary may only cover a handful of grid cells in the global model.
Other authors have sought to improve the spatial scale of reef exposure and resilience maps.Working in the Saudi Red Sea, Rowlands et al. [15] used a suite of remote sensing data to map indicators of both reef stress (fishing, industrial development and temperature stress) and resilience (coral abundance, framework abundance and water depth variability) and combined these into a single index of reef resilience.At a grid cell size of 1 km 2 and covering a large area, these maps supply important information on reef exposure and resilience for MPA network planning.However, further investigation is needed into the choice and relevance of the indicators themselves, the effects of the weightings used to combine individual indicators in the resilience index and the errors with which they can be mapped, to shed light on the sensitivity of the final product to those issues.
"Resilience", being a property of the reef ecosystem as a whole, is determined by different states and processes acting in concert [16,17].Various frameworks, e.g., [7], have emerged to assess what confers resistance or recovery potential for coral reefs, and a recent review suggests 11 key indicators, listed in Table 1, for which there is both expert consensus and considerable empirical evidence [18].Although some of these indicators can be mapped with current methods, several others have so far only been measurable in the field and at very small spatial extents.For example, coral recruitment and survival is typically measured using <1 m 2 quadrats [19].This makes them difficult to incorporate into the spatially explicit conservation planning software used to design MPAs and MPA networks, e.g., Marxan [20], which relies on spatially complete data layers.

Stress-tolerant coral taxa
Resistance Limited ability to map coral taxa directly [21], but identification of some species possible locally [22].Fine

Coral diversity Resistance
Beta dissimilarity or diversity can be derived from habitat maps [23,24], and a texture measure [25] can describe habitat composition.

Resistance
Existing sea surface temperature (SST) data products can be used to evaluate historical SST variability [11,26].Coarse Nutrients (pollution)

Resistance and recovery
The concentrations of chlorophyll a (Chl-a) and coloured dissolved organic matter (CDOM) can be used to evaluate historical nutrient levels [27], but accuracy is relatively poor in the shallow case-2 waters typical of coral reefs areas [28].

Resistance and recovery
Progress has been made toward defining an appropriate algorithm for coral reef areas [29].Coarse

Herbivore biomass Recovery
The biomass of herbivorous fish can be mapped from LIDaR or multispectral data combined with field data [30,31].Fine

Physical impacts Recovery
Physical disturbances that substantially alter the reef structure, such as ship groundings [32], tsunamis [33], hurricanes [34,35] and exposure to waves [24], can be mapped with high resolution satellite data.

Coral Disease Resistance and Recovery
Direct mapping of coral disease is not currently possible.Disease outbreak risk, however, has been predicted spatially on the basis of SST [36].

Macroalgae Recovery
The location of dominant macroalgae cover can be mapped, albeit with moderate accuracy, using standard classification techniques [37,38].

Coral Recruitment Recovery
Although coral recruitment cannot be directly mapped, environmental correlates with the potential for mapping include storms [39], suspended sediment [40], distance from river mouths coupled with current speed [41], herbivores [42] and the availability of suitable substrate for settlement, such as dead coral and crustose coralline algae [43][44][45].Dead coral is mappable under ideal conditions [46].Crustose coralline algae have not been mapped directly with remote sensing.
Fine/Medium/Coarse Fishing Pressure Recovery Fishing pressure cannot be mapped directly, but may be quantified with proxy variables, such as distance from settlements and markets [15].Could potentially also be inferred from residuals of models predicting biomass of fished species [47].

Fine/Medium
Our study investigates the ability to model and map four of these resilience indicators using high resolution satellite data and field observations.Specifically, we use field data to test the ability to map the following indicators: stress-tolerant coral taxa, coral diversity, herbivore biomass and coral recruitment (Table 1).We also include two other resilience indicators that complement those listed in Table 1, for which field data were available: (a) herbivore functional group richness, used in a recent study by Cheal et al. [48], as resilience may depend on having the full complement of ecological processes derived from different forms of herbivory [49][50][51], and (b) the cover of live coral and crustose coralline algae, promoted as an indicator within the Coral Health Index [52], both as an indicator of present-state algal-coral dynamics when coupled with macroalgal cover [53,54] and a proxy for the amount of substrate available for coral settlement following a mass bleaching event and, hence, related to coral recruitment potential [43][44][45][46][47][48][49][50][51][52][53][54][55].We do not investigate those indicators for which data products are already available (historical temperature variability and nutrients), those that can be mapped with existing methods (sedimentation, physical impacts and macroalgae) or those for which we do not have spatially distributed field data (coral disease and fishing pressure).We adopt and expand upon methodology that we previously employed to map aspects of reef fish communities in Fiji [47], in which we used high-spatial resolution satellite imagery to create spatial data layers of environmental variables that inform models predicting response variables.The fine spatial resolution of the satellite data (≤4 m pixels covering more than 260 km 2 ) enables mapping of the resilience indicators with a level of spatial detail meaningful for customary management systems in the Western Pacific and comparable to other maps of the area [47].

Study Area
The study area (Figure 1) comprises the traditional fisheries management area of Kubulau, Vanua Levu, Fiji, that covers a complex reef system with fringing and barrier reefs, lagoons and very deep water (>500 m depth) off the outer reef edge.The villages in Kubulau District are among the poorest in Fiji and heavily reliant on marine resources for both subsistence and cash income [56].Since 2005, a network of MPAs, including three large no-take reserves and 17 smaller periodically harvested fishing closures (tabu areas), has provided moderate protection for fish stocks in the area [57].To date, the marine ecosystem remains relatively intact, and outer reef areas, in particular, support high fish biomass and catch rates [57,58].

Data
The data for this study consist of georeferenced fish and substrate survey data from the Wildlife Conservation Society (WCS)'s reef monitoring program, as well as IKONOS and QuickBird satellite images for Kubulau District.The field data were processed to derive resilience indicators for each field site (Sections 3.1.1and 3.1.2),while a series of environmental data layers were derived from the satellite images (Section 3.1.3).These data sets were then used in random forest and universal kriging predictive models to produce maps of each resilience indicator for the study area and associated error estimates (Section 3.2).An overview of the data processing chain is shown in Figure 2.

Fish Data for Resilience Indicators
Field data were available from a total of 73 sites; after removal of sites that had missing data or were not covered by satellite data, 66 sites remained with fish data and 72 with substrate data.At each site, between three and five 5 m × 50 m belt transects were used to collect data on the fish community.For each belt transect, WCS staff recorded the number and total length (TL) of observed fish in each species, using 5 cm length classes for fish less than 40 cm and exact length for each fish over 40 cm, as per Jupiter and Egli [57].Fish species were classified into herbivores and others, and herbivores were then subdivided into four functional groups: large excavators, browsers, grazers/detritivores and scrapers/small excavators, as per Green and Bellwood [59].For each site, "herbivore functional group richness" was calculated as the number of functional groups present.To calculate the biomass for each fish, we applied the mean length of the fish size class or the exact length if available, in the length-weight (L-W) expression W = a × L b , with a and b parameter values for each species preferentially selected from sites closest to Fiji in Fishbase [60].For L-W conversions requiring fork length (FL), a TL-FL conversion factor, also obtained from Fishbase [60], was applied before the weight calculation.For each site, "herbivore biomass" was then calculated as the total estimated weight of all herbivorous fishes, divided by the area covered by the transects.No invertebrates were considered in the calculations of "herbivore biomass" and "herbivore functional group richness".A positive outlier was found for "herbivore biomass", caused by roaming schools of large-bodied herbivores (Scarus prasiognathos, Acanthurus xanthopterus) observed in one of the three transects at that site, and subsequent modelling was done both with and without this data point.Outliers were not present for the other indicators.

Substrate Data for Resilience Indicators
The transects used for fish data collection were also used to collect point intercept data (0.5 m point spacing) on the substrate composition, as per [57].For the purposes of this study, point observations from these transects were classified as "coral", "crustose coralline algae" or "other", the latter category including all substrate not identified as coral or crustose coralline algae.The percentage cover of "live coral and crustose coralline algae" was then calculated as the combined cover of those two classes.The number of "Juvenile corals", defined as all visible corals with a diameter less than 5 cm, was counted using six 1 × 1 m quadrats placed every 5 m along separate 1 × 25 m belt transects [61].Coral genera dominance data were collected from measures of colony size structure along the opposing sides of the same 25 m belt transects that were used for "Juvenile corals", followed by timed 10 minute swims to search for rarer taxa in the general reef area [61].Dominance was ranked qualitatively as: 4-dominant, 3-abundant, 2-common, 1-uncommon and 0-absent."Coral diversity" was calculated as the Shannon Diversity for each site, using generic dominance as abundance.Dominance was then multiplied by a thermal resistance index for each genus (0-2 scale, 0 being the least resistant to thermally-induced coral bleaching) with relative values obtained through literature searches on taxa tolerances to bleaching stress [5,[62][63][64].A "Stress-tolerant coral taxa" indicator was calculated as the mean (dominance × thermal resistance) score for each site, producing a low value for sites dominated by coral types vulnerable to temperature-induced bleaching and a high value for sites dominated by more resistant coral types.To calculate the density of "Juvenile corals", the number of juvenile corals was summed across the six 1 m 2 plots per transect and scaled to 100 m 2 .All calculated resilience indicators and their quantification are shown in Table 2.

Table 2.
Resilience indicators derived from field data and their quantification."Herbivore biomass" was modelled both with and without the outlier.

Resilience Indicator Quantification
Stress-tolerant coral taxa Mean (dominance × thermal resistance) The processing of the satellite data to produce data layers for a range of environmental variables has been described in detail elsewhere [47]; the present description is, therefore, brief.The satellite data consist of three QuickBird images and two IKONOS images.Pre-processing of each image included calculation of at-surface reflectance with ENVI's FLAASH ® routine, masking of clouds, resampling to 4 m pixels and creation of mosaics from overlapping IKONOS and QuickBird images.Maps of substrate (33 classes) and geomorphology (9 classes) were produced and validated from 9,646 georeferenced photos of the substrate with object-based image analysis [65], and bathymetry was mapped from 24 depth points with Lyzenga's [66] empirical approach.
An index of coral cover was generated by assigning pixels with coral as the dominant substrate class a value of 2, pixels with coral as a secondary substrate class a value of 1 and remaining pixels a value of 0. Spatially averaged values of the coral index were calculated for circular areas around each pixel by treating the index value as a continuous variable.10 different radii were used for these calculations-4, 10, 30, 50, 75, 100, 200, 300, 500 and 1,000 m-resulting in 10 environmental data layers.Using the same range of radii, habitat richness was quantified for each pixel as the number of different substrate classes present, and structural complexity was quantified as the standard deviation of depth estimates, resulting in 10 environmental data layers for habitat richness and 10 for structural complexity.The radii used in these calculations were selected to range from the smallest possible radius given the pixel size of the satellite data to radii deemed at least large enough to quantify the environmental variables at a scale where they may influence the resilience indicators.Although a more targeted selection of radii and neighbourhood shape for each resilience indicator is desirable, sufficient knowledge of the spatial scales at which the physical and ecological processes that influence the resilience indicators operate does not exist, so the exploratory approach used here is necessary and commonly used in similar studies [30,31,47,67].
Three additional layers were calculated containing the distance of each pixel from land (excluding Namenalala Island), seagrass and mangroves.Land and seagrass were identified from the geomorphic and substrate classifications, while mangroves were visually identified by Fiji Department of Lands for all of Vanua Levu using a Landsat 5 TM scene from 2001 [68].A map layer describing the conservation status of each pixel (tabu area, no-take reserve, no protection) as of 2011 was derived from existing spatial data.All 38 environmental data layers are listed in Table 3.

Predictive Model Types
Two types of predictive model were used to produce maps of each resilience indicator: random forest [69] and universal kriging [70] (Figure 2).Random forest is a non-parametric ensemble classifier that predicts the value of a single response variable from the values of multiple predictors; this model type has previously been shown to outperform other regression model types for modelling ecosystem variables in coral reef environments [47,71].A random forest model is composed of multiple regression trees [72], where splits at each node in each tree are based on a random subset of the available predictors [69].For each resilience indicator, a model was trained on a data set consisting of the observed values of the indicator at the field sites (the response) and the values of all environmental data layers at the field site locations (the predictors).Model development was done using R [73], using the "randomForest" package [74].Default values were used for random forest parameters, except for the number of trees in each forest and the number of predictors used to partition the data at each node.These two parameters were optimized with internal 5-fold cross-validation for each training partition, where the combination producing the lowest error value, quantified as the root mean square error (RMSE), within the training set was used.
Universal kriging is a geostatistical method that does not rely on predictors from the environmental data layers, but rather fits a deterministic trend to the large-scale spatial variation in the response variable and then uses locally optimized spatial interpolation of the residuals to account for small-scale spatial variation [75]; the large-scale trend that can be modelled in universal kriging is suitable for our study area, where it may act as a proxy for variations in, e.g., the distance from land, proximity to deep water, fishing intensity, etc. Universal kriging models were developed for each resilience indicator, all using a linear trend model.Spherical, exponential, Gaussian and Matern models were automatically fitted to the empirical variograms, and the model type with the lowest residual sum of squares was used for kriging.Model development was done in R, using the "automap" and "gstat" packages [76,77].
As a baseline against which to compare the performance of the random forest and universal kriging models, an "average model" was also developed.As for the random forest model, the "average model" was also trained on a data set consisting of the observed values of the indicator at the field sites, but without any predictors included, because the "average model" always predicts the average value of the training set at all locations.

Evaluation of Model Performance
For each resilience indicator, the performance of the predictive models was quantified using the root mean square error (RMSE) of predictions for the field sites, estimated with 5-fold cross-validation.In cross-validation, the full data set, containing both observed and model-predicted values for all field sites, is split randomly into a number of groups (here 5) of equal size.A training set is then formed by combining all except one of these groups, with the last group forming the test set to derive performance measures [78,79].This was repeated 100 times with training and validation sets randomly sampled from the full data set without replacement [71].To provide an additional measure of how well each resilience indicator could be mapped, the median RMSE from the 100 cross-validation repetitions was also compared for all three model types.The ability to map each indicator was then quantified as the percentage reduction in RMSE of the best performing models, compared to the average model (henceforth: percentage reduction).

Resilience Indicator Maps
The best performing model types, i.e., those with lowest RMSE values, were then used to produce maps of each resilience indicator.For the map production, the full data set was used to develop each predictive model.In recognition of the decreasing accuracy of all environmental data layers with increasing water depth, a maximum depth threshold was used to remove predictions from pixels in deep water.The threshold value was determined at 14 m based on visual interpretation of each satellite image, where no discernible bottom reflectance could be identified beyond this depth.Although it could be expected that differences in turbidity would result in different thresholds for different areas, this was not the case in our data.

Results
Table 4 lists, for each indicator, the RMSE values for all model types, the percentage reduction and the range of values found in the field data.The random forest models had lowest RMSE values for the indicators "stress-tolerant coral taxa", "herbivore biomass" and "herbivore functional group richness", while universal kriging had the lowest values for "live coral and crustose coralline algae" and "juvenile corals", and neither improved on RMSE values from the average model for "Coral diversity".RMSE and percentage reduction for "Herbivore biomass" both change substantially when the outlier is excluded from the analysis.Apart from different error values, the map products resulting from each model type are very different.Maps based on the random forest and universal kriging models are illustrated in Figures 3-5, showing mapped predictions, using both models, for those indicators with percentage reduction >10%.Maps of the remaining resilience indicators are provided as supplementary material.Results from the average model, being spatially uniform, are not shown.In Figure 3, showing "herbivore functional group richness", the map produced from the universal kriging model shows only minor local variation, but captures a gradual spatial trend with lower values in the north and higher values in the south.This trend is in accordance with the field data.Other trends in the field data, such as higher values near reef slopes and lower values on reef flats, in lagoons and on reefs fringing Vanua Levu, are only captured by the random forest model.In Figure 4, showing "live coral and crustose coralline algae", the universal kriging model successfully captures two trends seen in the field data: a general increase in values from north to south and a local hotspot in the central western part of the study area.The random forest model captures the north-south trend and partly captures the local hotspot, while also predicting high values near the reef crests.In Figure 5, the universal kriging model captures local variation well for the "stress tolerant coral taxa" indicator, in addition to a weak north-south trend, while the random forest model, as for the other indicators, predicts higher values near reef slopes and lower values on reef flats and in lagoons.

Discussion
The produced maps do not quantify coral reef resilience per se, as we do not have time series data spanning reef ecosystem response to acute or chronic disturbance (e.g., [80]) from which to calibrate and validate the models.Instead, the maps quantify selected aspects of the local reef community that are likely to confer resilience.It would be desirable to combine these indicators (and others, such as hurricane/cyclone damage, which may strongly impact reef recovery [34,35]) to produce composite indices of resistance or recovery, as done for exposure [10], or to produce a single resilience index, as in [15].However, the knowledge to quantitatively determine how these indicators combine or interact does not currently exist, and the limited predictability of the individual indicators would combine to give little confidence in the usefulness of a combined resilience map.The characteristics of the resilience indicator maps, along with the potential to improve their accuracy, are discussed below, along with issues around their application for coral reef management.

Model Performance by Indicator
The random forest and universal kriging models produced map products of varying quality for the six indicators, with percentage reduction values ranging from 0% (coral diversity) to 18.3% (live coral and crustose coralline algae).The best performing model type varied between indicators.Random forest, an inherently aspatial model that relies solely on the environmental data for predictions, performed best for "stress-tolerant coral taxa" and "herbivore functional group richness", suggesting that the spatial distributions of these indicators are significantly influenced by their environment, as quantified in the data layers or by correlates.
For the random forest model of "stress-tolerant coral taxa", the most important satellite-derived predictor variables were the coral cover index at 200 m and 500 m radii and structural complexity at a 30 m radius.Structural complexity at this scale appears to capture the large amount of macro-structure noted at field sites with the highest values of "Stress-tolerant coral taxa", dominated by massive Porites, Pavona, favids and other more thermally-tolerant coral taxa.At each of these sites, field observers noted considerable reef macro-structure in both forereef and backreef habitats (WCS, unpublished data), often with abundant crevices providing habitat for species with shade preferences, which often fall within robust genera, such as Mycedium, Echinophyllia and Psammocora [81].
The highest values of "herbivore functional group richness" were mostly seen on steep, shallow forereef slopes far away from villages, where observers would have had greater probabilities of encountering large excavators, such as large individuals of Chlorurus microrhinos and Cetoscarus bicolor, that are highly vulnerable to fishing pressure [50,52].We note, however, that the absolute amount of herbivory performed by each functional group will be influenced by the numbers and biomass of fish within each group, as well as their community composition.For example, after removal of large-scale herbivore-exclusion cages on the Great Barrier Reef (GBR), herbivory of the dominant macroalgae (Sargassum) was primarily carried out by a single species, Platax pinnatus [51]; had it been missing from the assemblage, phase-shift reversal and reef recovery may not have been possible.Similarly, in an experimental herbivory assay study on the Coral Coast of Fiji, only four species (Naso lituratus, N. unicornis, Chlorurus sordidus and Siganus argenteus) accounted for 97% of macroalgal consumption.These fish exhibited strong feeding complementarity between functional groups, meaning that all feeding groups need to be present to effectively remove the full suite of macroalgae from reefs [82].Furthermore, although maximal herbivore functional richness has been documented from sites along the GBR that successfully resisted phase shifts, other sites with low functional richness did not succumb to phase shifts, suggesting that measures of herbivore functional richness and diversity should be evaluated in concert with other environmental measures that may impact resilience thresholds [48].
Universal kriging, an inherently spatial model that does not incorporate information from the environmental data layers, performed best for "live coral and crustose coralline algae".This suggests that the spatial distribution of this indicator is significantly influenced by a physical/biological variable or process that is poorly described by the environmental data layers, such as the availability of coral larvae for settlement [83,84], factors related to post-settlement survival of corals [44], the impact of previous disturbances and competition with other life forms for substrate space.The low percentage reduction of this indicator with the random forest model (2.5%) is surprising, because of the ability of remote sensing to produce reasonable estimates of live coral cover [46,85,86] and since live coral accounts for ~88% of the combined areal coverage of live coral and crustose coralline algae in our field data.
The ability to map aspects of a reef fish community from high spatial resolution optical remote sensing [30,67] suggests that we should be able to model and map "Herbivore biomass".However, even after exclusion of the outlier, "herbivore biomass" only had a percentage reduction of 8.2%.This result may be partially explained by the fact that field observations of this indicator are highly influenced by large individuals and roving schools of medium-sized fish and, therefore, sensitive to instantaneous variation in the fish community [87].The identified outlier is a dramatic example of this.In addition, local fishing practices that target high-biomass individuals have the potential to greatly distort the statistical relationships that may exist between herbivore biomass and the environmental variables in a natural setting.At least one intensive, prolonged harvesting event occurred in the area inside and adjacent to the Namena MPA around Namenalala Island between 2009 and 2010 [56].Because our data set includes observations before, during and after this event, the sustained fishing would have introduced variability in the field data themselves.In future studies, incorporation of information on both historic and current fishing activities has the potential to greatly reduce prediction errors and shed light on the relative importance of climate and direct human impacts on reef resilience.
The density of "Juvenile corals" was not mapped well with our approach, but this was expected, given the many factors influencing spatial patterns of coral larval dispersal, recruitment and post-settlement survival.Although moderate-scale (~150 m), integrated models of oceanographic and tidal currents derived from satellite measurements can be used to identify potential hot spots of larval settlement [88], recruitment and post-settlement survival are determined by a number of processes not easily mapped and variable in both time and space, including tropical cyclone and bleaching recurrence interval [89], presence of settlement-inducing species of crustose coralline algae [44] and post-settlement mortality from grazing [90].
The poor ability to map "coral diversity" was slightly surprising, given known relationships between coral species and generic richness with depth, exposure to waves and reef habitat type [91][92][93], all of which we are able to map for our study area.The low predictability could possibly be attributed to the mismatch between large field sampling areas that cover a broad range of structural complexity and fine-scale pixel-based measures of structural complexity, which did not contribute substantially to the model results.The model outputs are also likely to be sensitive to imperfect georeferencing of field data collection sites, as well as the absence of field data covering the full range of habitat features, both of which are discussed in more detail below.

Sources of Error and Uncertainty
Unless directly detectable in the remote sensing data, the ability to map a resilience indicator depends on it being statistically related to a set of (mappable) environmental variables and/or sufficiently spatially autocorrelated such that it can be predicted from the field observations with kriging or other spatial interpolation approaches.The statistical relationship with environmental variables relies on those variables being both related to the indicator in question and accurately mapped at the appropriate spatial scale(s).In addition, the field data must reflect the true mean of the resilience indicators at each site [94], and the georeferencing of both field sites and satellite data must be as accurate as possible.Temporal changes between acquisitions of different parts of the data set may also contribute to errors, as may the mismatch in spatial scale between field and satellite data [95].These multiple factors all contribute to limiting the predictability of the resilience indicators.Apart from statistics on the accuracy of the maps of depth, substrate and geomorphology, reported in [47], the errors introduced by each factor remain unquantified, but several should be considered carefully in future studies.If quantified, the sources of the largest errors can be targeted for cost-effective reduction in RMSE values, e.g., through additional data collection, improved georeferencing, purchase of new and improved satellite data or modelling refinements.Information on errors can potentially also be incorporated into decision-support software for MPA network design to give planners more confidence that sites containing features of interest are prioritized for protection [94].

Spatial Extent of Field Sites
For each field site, the data used to calculate the resilience indicators were collected using different methods (belt transects, point intercept transects, quadrats) covering different areas.Each site is, thus, variable in size, depending on the indicator, and additionally, does not correspond to the native pixel size of the satellite data (4 × 4 m for IKONOS, 2.4 × 2.4 m for QuickBird) for any of the indicators.Previous studies mapping aspects of the fish community have used a single 4 m × 25 m belt transect [31], a single 5 m × 25 m belt transect [96] or a 5 m radius point count [30,47] to collect the field data at each site.In comparison, our field sites are 6-15-times larger in their areal coverage, depending on the indicator.This may mean that our field data are more representative of the true mean of these indicators for the general area around each site, but also that each site incorporates local variation in both resilience indicators and environmental variables.It should be noted that the results of this study were obtained with data that were not specifically collected for the purpose of spatial predictive modelling and, therefore, not optimized for factors, such as the number, distribution, georeferencing and spatial extent of field sites.The optimal spatial extent of the field sites requires further investigation, for example, through the use of nested field sites of different sizes, to evaluate the influence of site size both on prediction errors and in terms of the subsequent utility of the maps for MPA network planning.

Georeferencing of Field and Satellite Data
Georeferencing for each field site ideally refers exactly to the geographic centre of the site, i.e., the middle of the central 50 m transect, but errors are introduced for each site from imperfect placement of the GPS receiver over this point.In addition, the handheld GPS receivers used for georeferencing in our study has positional errors of ~5 m, and the positional accuracy of the satellite data is ~15 m for IKONOS [97] and ~23 m for QuickBird [98].These errors occasionally combined to visibly displace field sites significantly from their actual location, e.g., a field site known to be located near the reef crest was shown in deep water tens of meters off the reef when superimposed on the satellite data.In such cases, the coordinates of the field site were manually adjusted to the nearest suitable pixel based on site environment notes taken during data collection.This was necessary to avoid gross errors in derivation of environmental variables with our data, but cannot fully compensate for errors introduced by imperfect georeferencing.Improved absolute positional accuracy from newer satellites (e.g.GeoEye-1 [97] and Worldview-2 [98]) and/or improved georeferencing of field sites (e.g. with differential GPS receivers or use of site-specific knowledge [95]) should reduce this problem in the future.

Number and Distribution of Field Sites
The number of field sites and their distribution determine the data on which the predictive models are trained and, hence, their ability to make predictions for unsampled locations.Without sufficient training data, the models are unable to accurately quantify the sometimes complex relationships between the numerous environmental variables and each resilience indicator, thus increasing prediction error.In addition to the amount of data, the distribution of values for each environmental variable should cover the full range of values found in the study area and, ideally, their combinations, to avoid extrapolations when predicting for unsampled locations.Although such considerations were used to guide some of the data collection for this study, the distribution of values for some environmental variables could have been improved with additional field sampling.One example is structural complexity at spatial scales of 4 and 10 m, where values are positively skewed (Figure 6(a)).Any influence of high structural complexity (at these spatial scales) on the resilience indicators is, therefore, poorly documented in the training data and, hence, in the models and their predictions.Another example is distance to land, where no values exist between 7 and 12 km, nor between 17 and 20 km (Figure 6(b)).These plots can be used, along with the spatial data layers, to selectively locate sites for future data collection that will improve the coverage of each variable.

Derivation of Environmental Predictors
The environmental data layers used as predictors were derived from a previous study [47] and broadly match those used in studies predicting aspects of the fish community from environmental variables.However, the use of somewhat generic predictors easy to derive from high spatial resolution satellite data ignores known relationships, such as the influence of past disturbance history on live coral cover, coral diversity and recruitment rate [92,[99][100][101], the influence of ocean and tidal currents and source populations on coral recruitment and, hence, juvenile coral density [41,88] and the influence of exposure to waves on coral diversity, morphology and thermal stress tolerance [91,102].
The incorporation of such targeted predictors in future studies and the appropriate data sets, and temporal and spatial scales used to quantify them, will require significant effort, but is also likely to reduce prediction errors and strengthen confidence in the ability of the models to make predictions for unsampled locations.

Conclusion
We assessed two possible modelling approaches (Random Forest and Kriging models) for mapping six indicators of coral reef resilience-stress-tolerant coral taxa, coral generic diversity, fish herbivore biomass, fish herbivore functional group richness, density of juvenile corals and the cover of live coral and crustose coralline algae-using field observations and high spatial resolution satellite imagery.All the maps, except the one showing coral diversity, constitute a better-than-random description of the spatial distribution of the resilience indicators and, as such, constitute a resource that enables incorporation of reef resilience into MPA network planning in the region.
There was a marked difference between both prediction errors and maps produced with the random forest and universal kriging models (Figures 3-5).The relative strengths of the two approaches depend on factors, including the density of field observations and the accuracy of the environmental data layers used in the random forest model, and the relative performance of the two model types should not be considered of general validity.However, their differences do suggest a potential for further error reduction with hybrid models, such as regression kriging [103].
The two model types are also different in terms of the cost of their implementation.In our study, the random forest models rely on IKONOS and QuickBird satellite imagery for spatial data, as well as on Definiens Developer 7.0 and ArcGIS 10.0 for data analysis, together representing a substantial cost.The universal kriging models, on the other hand, did not require remote sensing data and were implemented entirely in freeware (R and contributed packages), putting this model type within reach of any organization with field data and computer access.The utility of free satellite data (e.g., Landsat TM/ETM + data) or pre-classified maps from the Millennium Coral Reef Landsat Archive [104], in combination with freeware for spatial data analysis, for mapping resilience indicators has yet to be explored, but could potentially eliminate the current additional costs for random forest or other regression-type models relying on remote sensing data.The different data and software requirements should be taken into account when selecting a model for mapping of resilience indicators, in addition to the predictions errors each model produces.The extent to which the number, distribution and georeferencing of field sites, as well as the development of more targeted environmental predictors influence how well a resilience indicator can be mapped should also be investigated, and additional modelling approaches, including hybrid combinations of aspatial and spatial models should be explored.Given the budgetary constraint of both government and non-government organizations that may wish to produce maps of coral reef resilience, such analyses will provide important information on the cost-effectiveness and trade-off between additional field data collection, satellite image acquisition and increasing sophistication of analysis, both in terms of cost and the expected improvement in the accuracy of the resulting maps [105].
In the context of MPA planning, planners now have the potential to consider reef resilience in addition to habitat and feature representation in decision-support software, allowing them to capture properties that maximize the overall resilience of an MPA network.This would represent a considerable improvement over the current practice of designing or adapting MPAs based on resilience scores compiled from field data collected from relatively few sites across the planning region [18,106].

Figure 1 .
Figure 1.Kubulau District and fisheries management area, Vanua Levu, Fiji.The coverage of high-resolution satellite data inside the red boundary forms the study area, with central coordinates 16°51′S, 179°0′E.Red circles indicate field data collection sites.

Figure 2 .
Figure 2. Flow chart of the data processing chain used to create the environmental data layers, the predictive models, error estimates and the maps of each resilience indicator.

Figure 3 .
Figure 3. Observed and predicted values of "herbivore functional group richness" for the study area.(a) Predictions based on the Universal Kriging model.(b) Predictions based on the Random Forest model.

Figure 4 .
Figure 4. Observed and predicted values of "live coral and crustose coralline algae" for the study area.(a) Predictions based on the Universal Kriging model.(b) Predictions based on the Random Forest model.

Figure 5 .
Figure 5. Observed and predicted values of "stress-tolerant coral taxa" for the study area.(a) Predictions based on the Universal Kriging model.(b) Predictions based on the Random Forest model.

Figure 6 .
Figure 6.Histograms showing the distribution of values for (a) "structural complexity, 10 m radius", and (b) "distance to land" at the field sites.

Table 3 .
Environmental data layers derived from satellite data and their data type.
* Pixels with depth values > 14 m were masked out based on visual interpretation of the imagery (details in Section 3.2).

Table 4 .
For each resilience indicator, the table shows (a) the median cross-validated root mean square error (RMSE) values for each model (units in Table2), (b) percentage reduction and (c) the range of values found in the field data.The best performing model for each indicator is outlined in bold.