Multi-Temporal UAV Data and Object-Based Image Analysis (OBIA) for Estimation of Substrate Changes in a Post-Bleaching Scenario on a Maldivian Reef

Coral reefs are declining worldwide as a result of the effects of multiple natural and anthropogenic stressors, including regional-scale temperature-induced coral bleaching. Such events have caused significant coral mortality, leading to an evident structural collapse of reefs and shifts in associated benthic communities. In this scenario, reasonable mapping techniques and best practices are critical to improving data collection to describe spatial and temporal patterns of coral reefs after a significant bleaching impact. Our study employed the potential of a consumer-grade drone, coupled with structure from motion and object-based image analysis to investigate for the first time a tool to monitor changes in substrate composition and the associated deterioration in reef environments in a Maldivian shallow-water coral reef. Three key substrate types (hard coral, coral rubble and sand) were detected with high accuracy on high-resolution orthomosaics collected from four sub-areas. Multi-temporal acquisition of UAV data allowed us to compare the classified maps over time (February 2017, November 2018) and obtain evidence of the relevant deterioration in structural complexity of flat reef environments that occurred after the 2016 mass bleaching event. We believe that our proposed methodology offers a cost-effective procedure that is well suited to generate maps for the long-term monitoring of changes in substrate type and reef complexity in shallow water.


Introduction
Coral reefs are declining worldwide as a result of the effects of multiple natural and anthropogenic stressors, such as increased water temperature, ocean acidification, overfishing and land reclamation [1][2][3]. Among these factors, temperature-induced bleaching events have been reported to cause significant coral mortality in the Indo-Pacific Ocean, with two main destructive events in 1998 and 2016 [4][5][6][7][8].
One of the most critical consequences of a bleaching event is the loss of structural complexity, which seriously affects the ecological functionality of the ecosystem [9]. Biogenic reefs are formed by complex structures built by corals, that support a high diversity and abundance of related organisms, such as fish, mollusks, echinoderms and other invertebrates [10][11][12]. The resulting structural complexity of the reef Our work monitored four selected portions of the entire Magoodhoo reef system, hereafter named Areas 1, 2, 3 and 4 (Figure 1c), from a few months after the bleaching event (February and March 2017) to the end of 2018 (October/November), more than two years later. Logistics reasons bounded the period of field data collection. Area 1 and Area 2 are both parts of the same reef flat area, which includes the reef crest that faces the open ocean to the south. The surveyed reef system is here exposed to the waves and the currents that come from the oceanic channel that separates Faafu and Dhaalu atolls. Area 3 and Area 4 are instead located toward the lagoon and lack a direct connection to the crest. Area 3 is a small subcircular patch reef confined in a minor and sheltered small-scale lagoon system, formed within the reef rim and bounded to the north by a tiny subelongated reef representing Area 4. Areas 1 and 2 are subject to very similar environmental conditions because they are included in the same wide geomorphic sector of the reef system, while Areas 3 and 4 represent unique small-scale reef sectors located in an area that is less exposed to ocean currents.

UAV Data Collection and Field Survey
DJI Phantom 4 was chosen for data acquisition. This consumer-grade UAV is a quadcopter with high sensing characteristics. It is equipped with a 1/2.3″ CMOS camera sensor (12.4 MP) that can collect RGB images with a resolution of 4000 × 3000 pixels and an integrated GPS/GLONASS system. DJI Phantom 4 is lightweight and easy to carry, and it can efficiently fly at low altitudes to obtain good-quality ground-resolution images. Additionally, easy take-off and landing procedures make this UAV a remarkable solution for low-height and short-range surveys. Metadata of the acquired

UAV Data Collection and Field Survey
DJI Phantom 4 was chosen for data acquisition. This consumer-grade UAV is a quadcopter with high sensing characteristics. It is equipped with a 1/2.3" CMOS camera sensor (12.4 MP) that can collect RGB images with a resolution of 4000 × 3000 pixels and an integrated GPS/GLONASS system. DJI Phantom 4 is lightweight and easy to carry, and it can efficiently fly at low altitudes to obtain good-quality ground-resolution images. Additionally, easy take-off and landing procedures make this UAV a remarkable solution for low-height and short-range surveys. Metadata of the acquired images are recorded in an EXIF (exchangeable image file format) file, which incorporates information such as shutter speed, apertures, ISO (International Organization for Standardization) and GPS (global positioning system) coordinates (x, y and z).
In addition to the UAV surveys, four bionomic transects (i.e., perpendicular to the coastline), each approximately 200 m in length, were examined across the four monitored areas (Figure 1c) in order to identify the composition of the benthic community and validate the remote sensing data.

Flight Planning
The surveys were designed using DJI GS PRO [54], a flight planner distributed by DJI for use on an iPad device. This application enables detailed planning of all the aspects of the UAV mission: generating flight paths, setting camera parameters and directly checking data acquisition on the iPad display. For all the surveys, we set fixed flight and image acquisition parameters (Table 1). Additionally, in order to collect relevant data for temporal comparison, we selected similar environmental conditions for all surveys. During the 2017 and 2018 surveys (Figure 2), we always flew during low tide to reduce water column effects and with a low sun position in the horizon (half an hour after the sunrise or one hour before the sunset) to avoid the glint effect in images or inconvenient shadows for data processing [24]. Moreover, we chose days with no or little wind (less than 5 m/s) and, consequently, a calm sea state: waves strongly influence the ability to recognize and map benthic communities [49]. Such conditions are common during the dry season (around December to April) and occur irregularly during the wet season (around May to November). The described expedients allowed us to minimize image adjustments during the post-processing and comparison of the output data.
Remote Sens. 2020, 12, x 6 of 26 Additionally, in order to collect relevant data for temporal comparison, we selected similar environmental conditions for all surveys. During the 2017 and 2018 surveys (Figure 2), we always flew during low tide to reduce water column effects and with a low sun position in the horizon (half an hour after the sunrise or one hour before the sunset) to avoid the glint effect in images or inconvenient shadows for data processing [24]. Moreover, we chose days with no or little wind (less than 5 m/s) and, consequently, a calm sea state: waves strongly influence the ability to recognize and map benthic communities [49]. Such conditions are common during the dry season (around December to April) and occur irregularly during the wet season (around May to November). The described expedients allowed us to minimize image adjustments during the post-processing and comparison of the output data.

SfM Processing and Georeferencing
Agisoft Metashape [55], a commercial structure from motion (SfM) software, was used to process UAV images of each area. The software was chosen for the affordable price of the license and its spontaneous workflow, user-friendly interface and excellent quality of the outputs [56]. These features have made it a widely used tool by the scientific community [38,57,58]. The images were processed following four main steps: alignment of the photos using a high-accuracy setting and the creation of a sparse point cloud; then, generation of a dense point cloud with an aggressive filter setting; and finally, creation of a digital terrain model (DTM) from the dense cloud. As final outputs, high-resolution orthomosaics were obtained from the DTM. Processing parameters and output details are provided in Table 1. For additional information on the SfM procedure, see Verhoeven (2011) and Brunier et al. (2016) [59,60].
In order to accurately georeference the orthomosaics resulting from SfM processing, we collected ground control points (GCPs) in all the monitored areas. At least eight GCPs per area were identified over structures that were easily visible and detectable in the UAV imagery ( Figure 3). Global

SfM Processing and Georeferencing
Agisoft Metashape [55], a commercial structure from motion (SfM) software, was used to process UAV images of each area. The software was chosen for the affordable price of the license and its spontaneous workflow, user-friendly interface and excellent quality of the outputs [56]. These features have made it a widely used tool by the scientific community [38,57,58]. The images were processed following four main steps: alignment of the photos using a high-accuracy setting and the creation of a sparse point cloud; then, generation of a dense point cloud with an aggressive filter setting; and finally, creation of a digital terrain model (DTM) from the dense cloud. As final outputs, high-resolution orthomosaics were obtained from the DTM. Processing parameters and output details are provided in Table 1. For additional information on the SfM procedure, see Verhoeven (2011) and Brunier et al. (2016) [59,60].
In order to accurately georeference the orthomosaics resulting from SfM processing, we collected ground control points (GCPs) in all the monitored areas. At least eight GCPs per area were identified over structures that were easily visible and detectable in the UAV imagery ( Figure 3). Global navigation satellite system (GNSS) coordinates were recorded with Emlid Reach RS © [61], low-cost single-frequency RTK (real-time kinematic) GNSS receivers (a base station and a rover). This GNSS is able to operate in long range radio (LoRa 868/915 MHz) mode, with one receiver working as a base station that sends real-time correction to the second receiver, which acts as a rover [62]. The base station was set in a fixed position for all the surveys, and the position accuracy was calculated using the PPP correction service offered by The Canadian Geodetic Survey of Natural Resources Canada [63]. Considering this single frequency GNSS in PPP base correction, the accuracy of horizontal and vertical components of the collected points is in the order of a few centimeters. After the photo-alignment in Photoscan, markers were placed over the GCPs, and corresponding RTK coordinates were imported. This process allowed us to estimate the XYZ accuracy of the orthomosaics as the root-mean-square error (RMSE).
Remote Sens. 2020, 12, x 7 of 26 navigation satellite system (GNSS) coordinates were recorded with Emlid Reach RS © [61], low-cost single-frequency RTK (real-time kinematic) GNSS receivers (a base station and a rover). This GNSS is able to operate in long range radio (LoRa 868/915 MHz) mode, with one receiver working as a base station that sends real-time correction to the second receiver, which acts as a rover [62]. The base station was set in a fixed position for all the surveys, and the position accuracy was calculated using the PPP correction service offered by The Canadian Geodetic Survey of Natural Resources Canada [63]. Considering this single frequency GNSS in PPP base correction, the accuracy of horizontal and vertical components of the collected points is in the order of a few centimeters. After the photoalignment in Photoscan, markers were placed over the GCPs, and corresponding RTK coordinates were imported. This process allowed us to estimate the XYZ accuracy of the orthomosaics as the rootmean-square error (RMSE).

OBIA and Classification of High-Resolution Orthomosaics
The analysis and classification of the orthomosaics obtained from UAV surveys were performed through Trimble eCognition Developer 9.4 [64]. The use of this software allowed us to carry out a robust image analysis in order to create comparable benthic coverage maps. Analysis (OBIA) algorithms followed the methods described in previous studies for satellite images [65] on satellite images and on UAV data [26]. A multiresolution segmentation algorithm based on homogeneity criteria was applied to the high-resolution orthomosaics (Figure 4b,b'). The image layer weights were set equal for all three bands, and the optimal scale parameter found was 150. For the homogeneity criteria, we established a shape value of 0.1 and compactness of 0.5. After the segmentation process, nearest neighbor classification (a supervised classification technique) was used to select well-defined samples as training areas. The accuracy of the samples was established by bionomic transects collected by snorkelers ( Figure 1c). Classification takes advantage of the features calculated for each object, such as the spectral value, but also size, shape, texture and proximity [66]. Using the thresholds between these parameters, we defined a ruleset for the classification algorithms that allowed us to classify and to label the three main substrate types that characterize the four study areas: hard coral (including live and dead corals), sand, and coral rubble (Figure 4c,c'). The sand class was classified on the basis of the spectral difference (mainly in brightness level and band ratio) between coral rubble and hard coral. Instead, hard coral and coral rubble were classified considering the difference in brightness, band ratio, standard deviation, saturation and the form and texture of the segments. Unfortunately, it was not possible to discriminate between dead and alive colonies from only the RGB channels because of the inaccuracy of an analysis based only on visible colors.

OBIA and Classification of High-Resolution Orthomosaics
The analysis and classification of the orthomosaics obtained from UAV surveys were performed through Trimble eCognition Developer 9.4 [64]. The use of this software allowed us to carry out a robust image analysis in order to create comparable benthic coverage maps. Analysis (OBIA) algorithms followed the methods described in previous studies for satellite images [65] on satellite images and on UAV data [26]. A multiresolution segmentation algorithm based on homogeneity criteria was applied to the high-resolution orthomosaics (Figure 4b,b'). The image layer weights were set equal for all three bands, and the optimal scale parameter found was 150. For the homogeneity criteria, we established a shape value of 0.1 and compactness of 0.5. After the segmentation process, nearest neighbor classification (a supervised classification technique) was used to select well-defined samples as training areas. The accuracy of the samples was established by bionomic transects collected by snorkelers ( Figure 1c). Classification takes advantage of the features calculated for each object, such as the spectral value, but also size, shape, texture and proximity [66]. Using the thresholds between these parameters, we defined a ruleset for the classification algorithms that allowed us to classify and to label the three main substrate types that characterize the four study areas: hard coral (including live and dead corals), sand, and coral rubble (Figure 4c,c'). The sand class was classified on the basis of the spectral difference (mainly in brightness level and band ratio) between coral rubble and hard coral. Instead, hard coral and coral rubble were classified considering the difference in brightness, band ratio, standard deviation, saturation and the form and texture of the segments. Unfortunately, it was not possible to discriminate between dead and alive colonies from only the RGB channels because of the inaccuracy of an analysis based only on visible colors.

DTM Analysis and Structure Complexity Evaluation
From the derived DTMs of Area 2 and Area 4, we selected two sub-areas far from the edges and close to GCPs. The reason for this is that during the SfM process, a slight broad-scale deformation can be produced that could affect the accuracy of the z value [67]. Moreover, from the visual comparison of the orthomosaics over time (2017 to 2018), we observed an evident change in the substrate composition in both areas. These two aspects make the two selected sub-areas suitable sites to be tested to evaluate our workflow designed to detect changes in structural complexity.
The models were imported in ArcGIS 10.6 and analyzed by using 3D Analyst Tools and applying the benthic terrain modeler (BTM) algorithm (i.e., vector ruggedness measure- [68]). The DTMs were firstly resampled to a uniform resolution of 5 cm/pix (using Resample Tool in ArcGIS 10.6) and then a low-pass filter (3 × 3) was applied. On the processed models, we calculated the linear rugosity and the average slope along five virtual transects and terrain ruggedness over the area around these

DTM Analysis and Structure Complexity Evaluation
From the derived DTMs of Area 2 and Area 4, we selected two sub-areas far from the edges and close to GCPs. The reason for this is that during the SfM process, a slight broad-scale deformation can be produced that could affect the accuracy of the z value [67]. Moreover, from the visual comparison of the orthomosaics over time (2017 to 2018), we observed an evident change in the substrate composition in both areas. These two aspects make the two selected sub-areas suitable sites to be tested to evaluate our workflow designed to detect changes in structural complexity. The models were imported in ArcGIS 10.6 and analyzed by using 3D Analyst Tools and applying the benthic terrain modeler (BTM) algorithm (i.e., vector ruggedness measure- [68]). The DTMs were firstly resampled to a uniform resolution of 5 cm/pix (using Resample Tool in ArcGIS 10.6) and then a low-pass filter (3 × 3) was applied. On the processed models, we calculated the linear rugosity and the average slope along five virtual transects and terrain ruggedness over the area around these transects in order to estimate changes in structure complexity between years. The transects were virtually laid parallel to the reef crest and started at randomly selected points within the target area.
The linear rugosity was calculated as the exact distance of a transect, considering changes in vertical surfaces (see Appendix A, Figures A1 and A2), divided by the linear distance [34,69]: The closer that this value is to 1, the flatter the surface is considered to be. The average slope along the transects was estimated using functional surface algorithms included in 3D analysis tools. Higher slope values indicate a steep area, while small slope values reflect smooth terrain [34,70]. The measure of area ruggedness was obtained using a vector ruggedness measure tool included in BTM. This system of measurement was used to quantify the surface complexity (rugosity) by assessing the 3D dispersion of vectors orthogonal to the surface of the DTM. Vectors orthogonal to each grid cell are decomposed into their x, y, and z components. A resulting vector is estimated within a 3 × 3 cell frame centered on each cell. A zero value indicates a flat terrain, while larger values indicate highly complex areas [9,68].

Map Accuracy Assessment
Map accuracy was assessed by comparing the outputs with 732 randomly distributed points that proportionally cover all the monitored areas. The points were manually classified, and the accuracy of the maps was evaluated visually using the orthomosaics; this was possible because of the high resolution of the UAV images (1.5 cm/pix), which allowed an on-screen check to be conducted by an expert in coral reefs with good knowledge of the study area [26,71]. Then, the accuracy was evaluated by comparing the outputs with the reference data in a confusion matrix to estimate the user and producer accuracy and the overall accuracy of the maps and to calculate the kappa index [72]. Finally, the benthic cover maps from 2017 and 2018 were compared in ArcGis 10.6 to highlight the differences between the two years and calculate the gain and loss in the three classes.

Map Accuracy
The benthic assemblage maps had an overall accuracy of 79% and a kappa index of 56% ( Table 2), values that represent a moderate agreement level, according to Congalton and Green (2008) [72]. Individual classes had a user accuracy of 85% for sand, 66% for coral rubble and 86% for hard coral. The lower accuracy of the coral rubble is mainly due to the absence of precise boundaries between the other two classes. In some areas, the mixture of sand with coral rubble formed a transition zone, which may have led to misclassification. This also applies to the margin between hard coral and rubble. In addition, during the map processing in eCognition, the jagged margins of the segments were smoothed. This process may have created some centimetric imprecision in the boundaries between classes that, in several cases, may have contributed to a reduction in classification accuracy.

High-Resolution Orthomosaics, Substrate Type Maps and Structural Complexity Variation
The processing of the UAV images led to the realization of high-resolution orthomosaics (1.5 cm/px) of 16 hectares of shallow-water coral reef environments around the island of Magoodhoo for the years 2017 and 2018. During this process, the average positional accuracy obtained for the models was 0.21 ± 0.09 m RMSE. The achieved detail of the models provides an overview of the dynamic of the benthic communities with low-cost and time-saving surveys. Moreover, the quality of the data highlights the difference in benthic assemblage extension between the two years with a visual on-screen comparison of the orthomosaics (Figure 4a,a').
The benthic assemblage maps realized from the OBIA analysis define the coverage of hard coral, coral rubble and sand in selected portions of the four monitored areas a few months after the 2016 bleaching events and more than two years later (Figure 4c,c'). The comparison of the maps revealed a general loss in the extent of hard coral framework coverage, precisely 11-15% of the total mapped area, with a total loss of 6119 m 2 ( Figure 5). Regarding the other two classes, we noticed an increase in coral rubble (a class closely linked to the degradation and fragmentation of the coral colonies) in Area 1, Area 2 and Area 4, the area most affected by the reduction in hard coral.

High-Resolution Orthomosaics, Substrate Type Maps and Structural Complexity Variation
The processing of the UAV images led to the realization of high-resolution orthomosaics (1.5 cm/px) of 16 hectares of shallow-water coral reef environments around the island of Magoodhoo for the years 2017 and 2018. During this process, the average positional accuracy obtained for the models was 0.21 ± 0.09 m RMSE. The achieved detail of the models provides an overview of the dynamic of the benthic communities with low-cost and time-saving surveys. Moreover, the quality of the data highlights the difference in benthic assemblage extension between the two years with a visual onscreen comparison of the orthomosaics (Figure 4a,a').
The benthic assemblage maps realized from the OBIA analysis define the coverage of hard coral, coral rubble and sand in selected portions of the four monitored areas a few months after the 2016 bleaching events and more than two years later (Figure 4c,c'). The comparison of the maps revealed a general loss in the extent of hard coral framework coverage, precisely 11-15% of the total mapped area, with a total loss of 6119 m 2 ( Figure 5). Regarding the other two classes, we noticed an increase in coral rubble (a class closely linked to the degradation and fragmentation of the coral colonies) in Area 1, Area 2 and Area 4, the area most affected by the reduction in hard coral.  In Figures 6-8, changes in substrate type composition can be directly linked to the integrity of the reef system. In particular, the comparison between 2017 and 2018 DTMs (Figures 6b-e and 8b-e) showed a general flattening of the monitored areas with a consequent loss of rugosity (Figures 6c-f and 8c-f). Therefore, both the linear rugosity and average slope values (Tables 3 and 4), measured along randomly chosen transects over the selected areas, suggest a general deterioration in structural complexity in the 2018 models. The linear rugosity index decreases close to a value of 1, indicating an increase in the extent of flat terrains in all the analyzed transects. Moreover, a reduction in the average slope values is noticeable and confirms the trend between the two years (Tables 3 and 4).
in structural complexity in the 2018 models. The linear rugosity index decreases close to a value of 1, indicating an increase in the extent of flat terrains in all the analyzed transects. Moreover, a reduction in the average slope values is noticeable and confirms the trend between the two years (Table 3 and  Table 4).
The terrain ruggedness maps, calculated for the areas around the transects (Figure 6c-f and Figure 8c-f), showed a pronounced drop in high ruggedness values. The comparison between 2017 and 2018 highlighted a significant increase in areas with values between 0 and 0.02 (Figure 7 and Figure 9), indicating a progressive reef flattening and an associated extensive decline in the structural complexity. The overall processing times for all the workflow steps, from image acquisition to DTM and structural complexity evaluation, are summarized in Appendix B (Table A1).  Table 3). (c,f) Images representing the rugosity maps of the area around transects.  Table 3).
(c,f) Images representing the rugosity maps of the area around transects.
Remote Sens. 2020, 12, x 12 of 26 Table 3. Variation in linear rugosity and average slope between years along the transects T1, T2 and T3 in Area 4 ( Figure 6).      Table 4). (c,f) Images representing the rugosity maps of the area around the transects.  Table 4). (c,f) Images representing the rugosity maps of the area around the transects. Table 3. Variation in linear rugosity and average slope between years along the transects T1, T2 and T3 in Area 4 ( Figure 6). The terrain ruggedness maps, calculated for the areas around the transects (Figures 6c-f and 8c-f), showed a pronounced drop in high ruggedness values. The comparison between 2017 and 2018 highlighted a significant increase in areas with values between 0 and 0.02 (Figures 7 and 9), indicating a progressive reef flattening and an associated extensive decline in the structural complexity. The overall processing times for all the workflow steps, from image acquisition to DTM and structural complexity evaluation, are summarized in Appendix B (Table A1).

UAV Surveys and Image Processing
The use of UAV data combined with SfM algorithms has increased in the last several years in different fields of environmental research, such as plastic and marine litter monitoring [73,74], marine habitat mapping [75,76], marine megafauna surveys [76][77][78], coastal geomorphology [79,80], structural geology [58,81] and forestry sciences [71,82]. However, from the first description of the great potential of these coupled techniques in monitoring ecology and geomorphology of the coral reefs by Casella et al. (2016) [24], the studies that have improved this methodology in this field are still few in number [48,83,84].
In our study, we used for the first time UAV images coupled with SfM and OBIA algorithms to monitor and map temporal changes in coral reef benthic assemblages following a mass bleaching event. The high-resolution orthomosaics confirmed the potential of these low-cost and versatile tools for the monitoring of shallow-water reef environments. Consumer-grade UAVs, such as DJI Phantom 4, can take off from the shore or a small boat. This flexibility allows for the easy monitoring of reef structures, such as patch reefs in lagoons, far from the coast. Of course, UAVs cannot be compared to a satellite in terms of spatial cover: satellite imagery analysis can be used to map an entire reef system [22,[85][86][87]. However, even if the spatial and temporal resolutions of the satellite sensors have improved in the last decade (Sentinel 2 offers 10 m/pix in visible bands and a revisit time of 5 days), they are not comparable to the centimetric resolution obtained from UAV surveys. This level of detail

UAV Surveys and Image Processing
The use of UAV data combined with SfM algorithms has increased in the last several years in different fields of environmental research, such as plastic and marine litter monitoring [73,74], marine habitat mapping [75,76], marine megafauna surveys [76][77][78], coastal geomorphology [79,80], structural geology [58,81] and forestry sciences [71,82]. However, from the first description of the great potential of these coupled techniques in monitoring ecology and geomorphology of the coral reefs by Casella et al. (2016) [24], the studies that have improved this methodology in this field are still few in number [48,83,84].
In our study, we used for the first time UAV images coupled with SfM and OBIA algorithms to monitor and map temporal changes in coral reef benthic assemblages following a mass bleaching event. The high-resolution orthomosaics confirmed the potential of these low-cost and versatile tools for the monitoring of shallow-water reef environments. Consumer-grade UAVs, such as DJI Phantom 4, can take off from the shore or a small boat. This flexibility allows for the easy monitoring of reef structures, such as patch reefs in lagoons, far from the coast. Of course, UAVs cannot be compared to a satellite in terms of spatial cover: satellite imagery analysis can be used to map an entire reef system [22,[85][86][87]. However, even if the spatial and temporal resolutions of the satellite sensors have improved in the last decade (Sentinel 2 offers 10 m/pix in visible bands and a revisit time of 5 days), they are not comparable to the centimetric resolution obtained from UAV surveys. This level of detail places the monitoring activities of marine environments with UAVs on a scale that lies between a satellite and a snorkeling/diving survey [24,26].
On the other hand, the main limitations of this type of UAV are the low accuracy of the built-in GPS and the restriction of data analysis to RGB images. Regarding the precision of the GPS, the unmodified geolocation is expected to have an error of up to 3 m [88]. This level of inaccuracy does not permit a comparison of the same area, monitored during two different flights, without precise GCPs. To overcome this limitation, it is possible to use a UAV with an integrated RTK (such as the recent DJI Phantom 4 RTK quadcopter) or georeference the models with RTK GCPs collected in all the areas. Instead, problems in the classification system for high-resolution UAV RGB imagery (e.g., the extensive time spent performing on-screen manual classification) can be overcome with the application of OBIA algorithms. OBIA processing considers not only the spectral information but also the geometrical and spatial relation between pixels or groups of them [89], and this allowed us to manage RGB spatial data with low spectral separation between the classes (e.g., hard coral and coral rubble). Moreover, the use of OBIA helps to speed up the classification process: once the ruleset has been developed, it can be applied to all the orthomosaics that represent the same environments with slight changes in feature parameters.
Another critical factor that may restrict the potential of the UAVs in benthic habitat surveys is the influence of environmental conditions on the water column. Thus, it is essential to monitor weather conditions and oceanographic parameters before the flight. Therefore, we planned the surveys to minimize class misclassification errors, and data acquisition was performed in the same optimal conditions for all the monitored areas in 2017 and 2018.

Classification and Map Comparison
The classification maps produced from the processing of orthomosaics with OBIA had good overall accuracy, 79%, which confirmed the validity of the classification process ( Table 2). The workflow adopted allows for the classification of benthic substrate types (sand, coral rubble and hard coral) with higher accuracy than maps realized from free satellite images (Sentinel 2) with the same methodology [65,86]. The maps represent the composition of the substrate that characterizes the surrounding reef rim a few months after the 2016 coral bleaching event and more than two years later. Before the bleaching, the reef flat, the back reef and other minor reefs comprising the reef system of the Magoodhoo island were composed of flourished frameworks (Figure 10a,b) of Acropora spp. colonies (mainly branching but with some tabular) with an increase in massive Porites spp. colonies moving from the shore towards the reef crest. During the bleaching, the most affected corals were tabular and branching Acropora spp., the primary habitat-forming species of shallow reef environments in the Maldives [44]. Mortality rates in this genus reached 80% [90], and signs of colony desegregation were already recorded a few months after the event [91]. The benthic maps (Figure 4b-d) describe an environment in which most of the extension of the hard coral class consists of dead Acropora still in the living position (Figure 10c,d). In addition, this category was composed of 35% dead coral, 51% sand, rubble and other organisms and only 14% live corals [92]. Unfortunately, it was not possible to discriminate the status of corals just from UAV images, so direct snorkeling observations were necessary. The comparison of the maps between 2017 and 2018 showed a substantial reduction in the hard coral class ( Figure 5). The dead colonies gradually lost their structure because of the impairment caused by the action of bioeroders and the mechanical action of currents and waves. Therefore, a shrinking of the coral frameworks was evident in all the monitored areas ( Figure 5). Most of the decline was driven by the degradation of branching and tabular colonies that, in some areas, entirely The comparison of the maps between 2017 and 2018 showed a substantial reduction in the hard coral class ( Figure 5). The dead colonies gradually lost their structure because of the impairment caused by the action of bioeroders and the mechanical action of currents and waves. Therefore, a shrinking of the coral frameworks was evident in all the monitored areas ( Figure 5). Most of the decline was driven by the degradation of branching and tabular colonies that, in some areas, entirely disappeared from the reef (Figures 6a-d and 11). Consequently, an expansion of coral rubble was recorded. The extension of these degraded areas can lead to an increase in algal cover and the formation of unstable substrate unsuitable for coral larvae settlement [93].
Remote Sens. 2020, 12, x 16 of 26 Figure 11. Close-range view of the three zones indicated with yellow arrows (A-C) in Figure 6.
Unfortunately, we were not able to compare the entirety of point clouds that generated all the DTMs and their terrain ruggedness because of the slight broad-scale geometrical distortion (i.e., doming effect [67,96]) that affected the final DTMs, especially along the edges and at major distances from collected GCPs. For this reason, we had to focus our analysis on a selected portion of the orthomosaics, where we observed that the accuracy of the variability of the topographic surface described by the surface roughness at the scale chosen by our analysis window (3 x 3 pixel) was not affected by the broad-scale error produced by the slight doming effect that could have affected the accuracy of the obtained DTMs. In future studies, this limitation can be overcome with the collection of additional water depth GCPs, at least 12 points regularly distributed over the monitored area, to scale and calibrate the Z value correctly in all the models. In Area 3, the trends of the benthic assemblages, from 2017 to 2108, were slightly different from the other three ( Figure 5). This is due to the sheltered position of the reef inside a small lagoon, and the mechanical action of waves and currents on colonies is likely less pronounced. Nevertheless, the general effect observed in all the mapped environments unveiled a significant transformation in substrate type, with a resulting deterioration in the integrity of the hard coral frameworks composing the surveyed reef areas. Rugosity is a crucial ecological factor since biodiversity is strongly correlated with habitat complexity [70,94,95]. The negative variation in this parameter due to the transition from healthy coral reefs to post-bleaching degraded environments can have significant implications on reef-associated organisms such as fish. Moreover, the described changes in the structural complexity of the fore reef can lead to an increase in the rates of coastline erosion as a result of the loss of their coastal protection functions [1,14]. The comparison of selected portions of the DTMs obtained from the SfM workflow let us detect this decline (Figures 6b-d and 8b-d). Linear rugosity and average slope, measured along random virtual transects, showed a decrease between 2017 and 2018 in both indexes (Tables 3 and 4). These values, which are associated with the analysis of the terrain ruggedness, provide a detailed overview of the variation (Figures 6c-f and 8c-f) and reveal shifts in the substrate type structure. Unfortunately, we were not able to compare the entirety of point clouds that generated all the DTMs and their terrain ruggedness because of the slight broad-scale geometrical distortion (i.e., doming effect [67,96]) that affected the final DTMs, especially along the edges and at major distances from collected GCPs. For this reason, we had to focus our analysis on a selected portion of the orthomosaics, where we observed that the accuracy of the variability of the topographic surface described by the surface roughness at the scale chosen by our analysis window (3 x 3 pixel) was not affected by the broad-scale error produced by the slight doming effect that could have affected the accuracy of the obtained DTMs. In future studies, this limitation can be overcome with the collection of additional water depth GCPs, at least 12 points regularly distributed over the monitored area, to scale and calibrate the Z value correctly in all the models.

Conclusions
Because of the increase in the frequency of extreme heating events driven by climate change, the future of the coral reefs in the Anthropocene is unknown. Estimation of coral recovery after a bleaching event is definitely challenging [1]. In the Maldives, a full recovery from the 2016 mortality could take decades [97][98][99] and might be adversely affected by the increase in coral reef exploitation from fishery activities, tourism and land reclamation [2,42,100]. In this scenario, affordable monitoring techniques and reproducible protocols are essential to improve high-resolution data collection.
Our study demonstrated the enormous potential of UAV data coupled with SfM and OBIA analysis as a tool to monitor shallow-water coral reef assemblages over time. Although with some limitation, our protocol showed that high-resolution imagery collected with UAVs can be efficiently analyzed to monitor changes in substrate type and structural complexity. The spatial resolution achieved from the orthomosaics allowed us to classify benthic community types with good accuracy, and the comparison of the maps over time clearly shows the deterioration of the reef area surveyed, confirmed by the loss in structural complexity that is estimated by the documented variation in terrain ruggedness.
Limitations, such as the difficulty in the calibration of DTMs in order to calculate accurate elevation data over the entire area and over time, could be easily overcome in the future with a better GCP collection strategy and the continuous increase in the efficiency of UAV platforms and cameras to avoid broad-scale errors in the provided DTMs (i.e., the doming effect). However, the maps and the DTMs generated in this work prove that quantitative data can be effectively produced for a long-term and cost-effective monitoring program to check the dynamics of shallow-water coral assemblages in this new human-dominated era. Appendix B Table A1. Workstation set-up and a summary of the entire workflow time for all the monitored areas. The processing time is closely linked to the configuration of the computer and the number of collected images.