Spatial Variation in Coral Diversity and Reef Complexity in the Galápagos: Insights from Underwater Photogrammetry and New Data Extraction Methods

1. Introduction

The Galapagos archipelago, located in the Eastern Tropical Pacific approximately 1000 km west of Ecuador, is globally renowned for its biodiversity and unique ecological systems [1].

Among its marine ecosystems, coral communities range from those growing on volcanic rocks to fully developed biogenic reefs, such as Wellington Reef in Darwin [2,3].

Despite their relatively small size, these reefs play a crucial role in increasing regional marine connectivity and supporting diverse marine life  [4]. However, Galápagos corals face numerous threats including climate change, oscillations in ocean currents, cold-water stress [5], warm-water bleaching events [6,7,8], urban development, tourism impacts [1], intense grazing pressure [9], coral diseases, pollution, habitat destruction [8], and invasive species [10].

Several studies have examined the coral communities in the Galápagos and their responses to environmental disturbances [11]. Over the past four decades, extreme climatic events such as the 1982–1983 and 1997–1998 El Niño events have caused severe coral bleaching and mortality, leading to long-term declines in reef cover and biodiversity [6,7,12]. Additionally, cold-water upwelling events, such as the 2007 La Niña, have further impacted these ecosystems [13,14,15].

More recently, marine heatwaves linked to consecutive El Niño years have raised concerns about the ability of Galápagos corals to withstand recurrent disturbances [13,16].

Despite several studies documenting these impacts, many knowledge gaps remain regarding coral recovery dynamics, resilience mechanisms, and spatial variation in reef complexity across the archipelago [13].

Moreover, much remains to be learned on how the recent marine heat waves caused by consequent el Nino years [17] will affect these heavily disturbed reefs.

To bridge these gaps we started a long-term monitoring effort at seven sites around the archipelago using underwater photogrammetry (image-based mapping). Photogrammetry is a form of 3D imaging that has been used extensively by marine ecologists in the past decade (e.g., [18,19,20]). Using a series of subsequent images as input, Structure-from-Motion photogrammetry builds a 3D model of the scene and estimates the relative camera locations. The 3D models capture the community structure together with the structural complexity of the site and are thus useful for ecological studies.

Traditionally, measuring coral reef structural complexity and biodiversity was done in situ by divers swimming along transect tapes [21,22,23]. These were later replaced and supplemented with photographic methods such as photo quadrats [24,25] and video surveys [23]. With the advent of photogrammetry, dozens of works have dealt with extracting community composition and structural metrics from 3D models of coral reefs [26].

Structural complexity was historically measured in situ using a chain draped over the reef surface and a measuring tape to compare surface contour versus planar distance within a quadrat [21,27]. This foundational approach led to the chain intercept transect method, which extended the same ratio-based principle across longer transects  [23,28,29,30]. Although these manual approaches provided critical baseline data on reef complexity, photogrammetry has largely supplanted them by enabling high-resolution 3D reconstructions of reef environments [31].

The first predominant works on measuring structural complexity from 3D models of the seabed replicated the rugosity metric in 2D and 3D [32,33,34], followed by more metrics such as shelter space, vector dispersion, and fractal dimension [35,36,37,38].

Fractal dimension (FD) is a unitless measure of how a structure fills space, based on how spatial details change with measurement scale [39]. Fine-scale measuring tools capture more features in complex surfaces, leading to a higher FD values. While FD has been measured in situ [30,40] and from photogrammetric models of coral colonies and reefs [37,41,42,43,44,45], earlier approaches often relied on top-down projections that neglect overhangs and crevices. Since many cryptic reef organisms inhabit these hidden or vertical surfaces [22,46], our method extends FD calculation to account for such complex habitats. An important work in this context measured fractal dimension from digital rugosity [45] similar to our work. the main difference is that they used top down projection to calculate the chain while we use directed geodesic walks making our method more sensitive to overhangs and crevices, accounting for full 3D (see Figure 4).

Similarly, approaches to measuring biodiversity have evolved from diver-based surveys [22,47], photo quadrats [24,25], and video surveys [23] to photogrammetric surveys, which are increasingly used to analyze benthic communities [18,44,48,49]. Many studies utilize photomosaics derived from 3D models to examine coral community structure [26]; however, because photomosaics are 2D images generated from a single point of view, they inevitably lose some of the 3D information. Despite offering large spatial coverage (commonly on the order of 100 m² at sub-centimeter resolution) and enabling advanced analyses via semantic segmentation [50,51], orthomosaics overlook the three-dimensional complexity of reefs— particularly overhangs and crevices— where numerous cryptic organisms reside [46].

Alternative approaches for analyzing coral communities from 3D models include 3D point annotations [52,53], and several works that employed deep learning for full segmentation of 3D models [54,55,56,57]. In this work we perform 3D annotations from multiple views using a synthetic view on each annotation point providing the annotator 3D context (See Methods Section 2.5), which to the best of our knowledge was not done before.

In this study we applied underwater photogrammetry to assess coral community structure and reef complexity at seven sites across the Galápagos archipelago. Specifically, we introduce two novel analytical approaches: a 3D point-annotation technique for benthic classification and a fractal dimension assessment using directed geodesic walks, which improves upon traditional rugosity and fractal dimension measurements by accounting for the Z axis and incorporating overhangs and cryptic spaces. By integrating these techniques, we seek to provide a baseline for long-term coral monitoring, contribute to the understanding of reef resilience mechanisms, and inform conservation strategies to prioritize structurally diverse reefs that may serve as climate refugia.

2. Materials and Methods

2.1. Study Sites

The locations are distributed throughout the archipelago (Figure 1,Figure 2) and were selected based on the Charles Darwin Foundation team’s long-term monitoring sites [13] and earlier research in the area [2,4,8,9,12,58]. At each location, we began by examining the area to find a $10\times 10m2$ section featuring the highest coral density and species diversity at a depth of 8–12 m. Balancing these criteria is challenging, as the area with the most coral coverage is often not the most diverse. Additionally, the reef depth varies at each site, as outlined in Table 2. Once the plot was selected and its boundaries marked using a transect tape and a compass we installed permanent metal stakes at each corner of the plot.

2.2. Image Acquisition and Scene Set-Up

Before image acquisition we placed a set of photogrammetric markers and color cards in the plot. These are used for scaling the map and assessing the quality of the images in terms of sharpness and illumination during image acquisition and also in post-processing. Moreover, the targets are placed in a way that helps the diver swim in a systematic manner while maintaining overlap between reciprocal legs.

Image acquisition was done using a Nikon D850 DSLR camera with a 35 mm Nikkor lens and two BigBlue video lights. Images were acquired from a 2-3 m distance to the substrate at one frame per second while the diver holding the camera maintained a slow and steady swimming pace and a downward-looking angle. The diver swam in a boustrophedonic pattern (lawn-mower pattern) across the plot to complete image acquisition from the first distance. Then the diver changed depth to a ∼ 5m distance from the substrate and imaged the plot again. Imaging from two distances helps the image alignment process of the 3D reconstruction [20].

2.3. 3D Processing:

Agisoft Metashape software (version 2.12) was used to reconstruct the 3D models from the series of subsequent images. We used the parameters described in (Table 1).

We scaled the models using the photogrammetric targets as scale bars. 10-15 scale bars were defined per model in order to minimize the scaling error. For each scale bar we marked the points of known distance in several input images and entered the known distance value. We marked these points in several images until the scale errors were reduced to sub-centimeter values.

2.4. 3D Model Alignment and Export

After the models were constructed and scaled we aligned them to the XY axis. This step is important for the structural complexity measurements as we don’t want the effect of aspect (slope of the reef— differences in depth between corners of the plot) to affect the measurement. We first fit a plane to each model using the open3D python [59] package. The normal of this plane is used for finding the 3D rotation that aligns the model on the XY plane. We find the rotation matrix of the normal and rotate the model by the inverse. We then use the model’s bounding box to find a translation to the center of the axis (this is detailed in the supplementary code "orient model and set region limits").

After rotating and translating the model, we build a point cloud from the 3D mesh using a distance-based sampling at 1 cm point density in Metashape. Finally, a $10{\mathrm{m}}^3$ region around the center of the axes was defined and the point clouds were exported for further analysis (cube counting and directed geodesic walks).

The point clouds are released as part of the dataset, please see section Supp. 1.2 and Figure Supp. 2 for more details.

2.5. 3D Annotation

We developed a novel approach for annotating 3D models using the best camera view and a synthetic image. First, we use the sparse point cloud generated in image alignment and read it using the open3D library [59]. We build a mesh from these points using the open3D alpha shapes surface reconstruction method. We sample points on this mesh using the open3D sample_points_poisson_disk method where each point has approximately the same distance to the neighbor points. Then we use a KD-tree to find the closest point in the sparse cloud to each of the points. We then read these points back as markers on the 3D model. The number of points per model is determined by measuring the surface area of the model in the region and multiplying it by ten in order to achieve approximately ten points per m².

For each marker, we extract the image with the smallest reprojection error. Reprojection error is measured by projecting the point from the 3D model (origin point) to the image and back to the 3D model (target point) and calculating the distance between them. After selecting the image with minimal error, it is cropped around the point in 1600 X 1600 pixels. Eventually we resized the images by half to enable their upload resulting in 800X800 pixels per image. We also extract a novel view for each point- a synthetic image of the 3D model using a top-down projection. This image is called synthetic because it was never taken by a camera. It is a rendering of the 3D model from a given viewpoint. We used Metashape’s functions to produce these images with each marker as the center of each image and a fixed size of 1600 X 1600 pixels (Figure 3). The advantage of using both images is that the close up camera image helps to identify the class of the object while the synthetic view helps to put that in context. This is helpful for distinguishing corals such as Porites and Montipora which look similar from a close up view but have different 3D structure in the colony.

The pair of images (best camera view and synthetic image) per point was uploaded to labelbox [60], an online platform for image annotation. The camera view is a close up on the point while the synthetic image shows the 3D context of the point. Therefore, using the image pairs helps to annotate the points rapidly and accurately.

The annotations are released as part of the dataset, please see section Supp. 1.2 and Figure Supp. 2 for more details.

2.6. Classification

A total of 17,807 image pairs were annotated in 12 classes (results are detailed in Table 2): Pavona, Pocillopora, Porites, Psammocora, Rock, Sand, Tubastraea, Caulerpa, Dead Coral-Framework, Dead Coral-Rubble, Other, and Unknown. The class Dead Coral-Framework was used to describe biogenic substrate found mostly in Darwin. The class Unknown was used for points where the images (best camera view and synthetic image) did not show the same point. The class Other was used for points that were not in any category.

The images were annotated by three members of the marine monitoring team of the Charles Darwin Foundation (W.I, F.T, and W.B.S) including reviewing all image pairs to validate the classifications.

We also used a classification: "Projection quality", to describe the fit between the synthetic image and the best camera view. When both images showed the same point, we classify it as good, and when the synthetic view does not exactly align with the best camera view we classify it as bad. This can happen if the 3D point is obscured by an overhang, since the synthetic views are generated only from a top down view. We used this classification becuase in the future we aim to use these synthetic images together with the best camera view to train a neural network for coral identification, and using the synthetic views can increases the amount of training data almost twofold. Please see section Supp. 1.1 and Figure Supp. 1 for more information on using this data for training a classifier.

Figure 3. Method for 3D annotation using the best camera and synthetic views: each row is from a different site. The left images show a 2X1 m close-up of the 3D model with the markers and their identities. The right images show the pair of images that were used to annotate a marker. For each point, two images are extracted: a 1600-pixel crop around the annotation point on the best image that depicts the point, and a novel synthetic image from a top-down view of the point.

Figure 3. Method for 3D annotation using the best camera and synthetic views: each row is from a different site. The left images show a 2X1 m close-up of the 3D model with the markers and their identities. The right images show the pair of images that were used to annotate a marker. For each point, two images are extracted: a 1600-pixel crop around the annotation point on the best image that depicts the point, and a novel synthetic image from a top-down view of the point.

2.7. Community Metrics

We calculated the relative abundance of each class per site by dividing the number of observations from each genus by the total number of observations in the model. We calculated the Shannon diversity only for coral classes: Pavona, Pocillopora, Porites, and Psammocora. We used R Statistical Software (v4.1.2; R Core Team 2021) and the R package vegan [61,62]. Shannon diversity is calculated as:

$$H=-\sum _{i=1}^S{p}_{i}log{p}_i$$

(1)

where $p_i$ is the proportion of species i, and S is the number of species.

We calculated percent coral cover calculating the proportions of four coral classes— Pavona, Pocillopora, Porites, and Psammocora, out of all classes.

We used Principal Component Analysis (PCA) to visualize the difference in community structure on two axes. We ommited the classes Unknown, Other, and Caulerpa, and from this analysis, thus reducing nine dimensions to two (Figure 5). PCA is a linear method that identifies the axes (principal components) along which the data has the greatest variance, projecting the data onto these components.

2.8. Structural Complexity Metrics

2.8.1. Fractal Dimension from Cube Counting

We previously worked on cube-counting for assessing fractal dimension of coral reefs and detailed the method extensively [37]. In short, the cube-counting method [63] enables to calculate the fractal dimension of objects as:

$$FD=\underset{ϵ\to 0}{lim}{\displaystyle \frac{logN\left(ϵ\right)}{log(1/(ϵ\left)\right)}},$$

(2)

where $ϵ$ is the length of the cube and N($ϵ$) is the number of cubes required to cover the object at a given length. N($ϵ$) is calculated for several cube-lengths ($ϵ$) and the fractal dimension is calculated as the slope of the fitted line between log(N($ϵ$)) and log(1/$ϵ$).

We take each pointcloud and bound it with a minimal bounding cube ($\sim 10$ m³). The cube is then divided into 8 equal cubes by dividing each axis. In each iteration, cubes that contain a part of the mesh are counted (Eq. 2) and used in the next iteration.

For flat shapes, the size of the cube should not affect the measurement while in complex shapes more surface detail is revealed with smaller cube sizes, and the number of cubes increase exponentially when decreasing the cube-length. flat reefs are expected to have lower values close to two, and reefs that are rich in structural features— i.e., structurally complex reefs, are expected to have values closer to three [37].

2.8.2. Fractal Dimension from Directed Geodesic Walks

This is a novel method for calculating the fractal dimension that elaborates the classic rugosity metric and applies a similar measurement over multiple scales of measurement across the full model.

$$\mathtt{Fractaldimensionfromgeodesicwalks}=\underset{ϵ\to 0}{lim}{\displaystyle \frac{log\left(walkRatio\right)}{log(1/(stepSize\left)\right)}}$$

(3)

$$\mathtt{Walkratio}\left(\mathtt{Stepsize}\left(\mathtt{r}\right)\right)={\displaystyle \frac{\mathtt{Distance}\mathtt{travelled}}{\mathtt{Geodesic}\mathtt{Distance}}}$$

(4)

We first construct a mesh from the point cloud using the Open3D Python package. The mesh is then cropped into 5 cm-wide slices with 5 cm intervals along both axes, resulting in approximately 200 slices per site. To measure the structural complexity of these slices we simulate a walk along each slice over a range of step sizes. The starting and ending points of the walk are defined as the minimal and maximal coordinates along the slice. Using PyVista [64] python package we calculate the geodesic path: the shortest path along the mesh surface—between these points. This path determines the order of points for subsequent walks therefore we call this method directed geodesic walks. We use the recorded vertices of the geodesic path as a baseline for building a new point cloud by densifying these points. This enables us to maintain the order of points rather than their coordinates for the walk algorithm.

It is crucial to balance the slice width: too wide, and the walk avoids the actual terrain, underestimating complexity; too thin, and the mesh fragments, making the walk ineffective. Moreover, we had several cases (57 out of 29,343) where the walk did not work because the algorithm was unable to find points in the specified step size. We excluded these from the FD assessment.

The algorithm follows two key rules: avoid revisiting previous points (using point indices rather than coordinates), and choose the minimal point on the axis of the walk, for example if the walk is along the x axis, choose the point with the minimal x coordinate (see Figure 4). For each slice, the walk is conducted twice per step-size : once from start to finish and then in reverse. The average number of steps is recorded, and the walk ratio is calculated as the distance traveled divided by the geodesic path length. A linear model is then fitted to the log relationship between walk ratio and step size for each slice, with the overall score per site being the mean of these measurements.

This method captures structural complexity by observing how the walk ratio decreases with increasing step size (Figure 4). For flat surfaces, step size has little effect, as both small and large steps cover similar distances. For complex surfaces, however, smaller steps traverse more intricate paths, highlighting structural complexity. This difference provides an estimate of the dimension and complexity of the reef structure. The idea is that for complex structures the walk ratio will decrease when increasing the step size. For flat surfaces, the step size r has no influence on the distance of the walk. To portray this, imagine an ant, a horse, and an elephant crossing a surface from start to end. If the surface is flat they will all travel the same distance. If the surface is wavy, rugged, and complex, the ant will travel a larger distance than the horse and elephant. Measuring how the distance traveled changes by step size gives an estimation of the dimension of the mesh and the structural complexity of the reef.

Figure 4. Method for calculating fractal dimension from geodesic walks: A) A slice is extracted from the 3D model (left) and used to calculate the geodesic path between start and end points. This process is repeated at 5 cm intervals to extract approximately 100 slices per side (top, right). The model is then rotated 90° around the z-axis, and the process is repeated to capture slices on both axes. B) A profile of a geodesic path extracted from the 3D model. c) a close up on the slice of the mesh and the geodesic path. The geodesic path links vertices and enables us to assign the walk order by their index (ordinal) and not only their position. The algorithm works by finding points in radius (step-size) and choosing the point with the minimal coordinate on the walk axis. The bottom image in (C) shows the steps of the walk and a red circle with radius step size to illustrate how steps are chosen. the right image is a close up of the secyion outlined in dashed blue that shows step n and possible points within the radius(marked with asterisk).The chosen point n+1 is the point with minimal x coord. D) a walk ratio is determined for each step-size and used to calculate FD for the slice of the mesh. On complex paths, increase in step size causes a decrease in walk ratio. The overall score per site is the average of 200 slices (A, right).

Figure 4. Method for calculating fractal dimension from geodesic walks: A) A slice is extracted from the 3D model (left) and used to calculate the geodesic path between start and end points. This process is repeated at 5 cm intervals to extract approximately 100 slices per side (top, right). The model is then rotated 90° around the z-axis, and the process is repeated to capture slices on both axes. B) A profile of a geodesic path extracted from the 3D model. c) a close up on the slice of the mesh and the geodesic path. The geodesic path links vertices and enables us to assign the walk order by their index (ordinal) and not only their position. The algorithm works by finding points in radius (step-size) and choosing the point with the minimal coordinate on the walk axis. The bottom image in (C) shows the steps of the walk and a red circle with radius step size to illustrate how steps are chosen. the right image is a close up of the secyion outlined in dashed blue that shows step n and possible points within the radius(marked with asterisk).The chosen point n+1 is the point with minimal x coord. D) a walk ratio is determined for each step-size and used to calculate FD for the slice of the mesh. On complex paths, increase in step size causes a decrease in walk ratio. The overall score per site is the average of 200 slices (A, right).

3. Results

We studied seven reefs using 3D point annotations and structural complexity metrics. All results are summarized in Table 2:

Table 2. Summary table of all counts, structural complexity measurements, GPS coordinates, and depth per site

site name	Caulerpa	DeadCoral-Framwork	DeadCoral-Rubble	Other	Pavona	Pocillopora	Porites	Rock	Sand	Tubastraea
Darwin- Wellington Reef	105	66	1	46	566	29	1108	855	177	9
Española- Xarifa	1	31	65	7	0	0	584	353	800	0
Floreana- Tres Cuevitas	0	11	12	12	731	1	0	1409	463	1
Floreana-Punta Cormorant	0	16	0	13	351	0	16	1317	485	1
Marchena -Roca Espejo	0	90	38	18	35	613	361	1063	265	1
Pinta - Cabo Ibetson	0	332	39	18	3	674	10	1423	64	0
Wolf - Corales	0	40	0	12	812	7	961	580	93	2
Site name	Unknown	Psammocora	Cube Counting FD	mean_slope	lacunarity	SurfaceArea	Lat	Lon	Depth
Darwin - Wellington Reef	129	0	2.267	0.108951	1.495	309	-91.995824	1.67824	15
Española - Xarifa	8	1	2.138	0.035663	1.626	185	-89.644622	-1.357863	5
Floreana- Tres Cuevas	130	8	2.229	0.101956	1.998	225	-90.408089	-1.236136	10
Floreana- Punta Cormorant	51	1	2.202	0.074877	1.615	276	-90.419325	-1.223949	12
Marchena - Roca Espejo	74	19	2.212	0.062882	1.458	258	-90.401241	0.312162	8
Pinta - Cabo Ibetson	71	18	2.208	0.062978	1.536	265	-90.720945	0.544216	8
Wolf - Corales	101	0	2.225	0.096196	1.684	261	-91.815861	1.386624	12

Table 2. Summary table of all counts, structural complexity measurements, GPS coordinates, and depth per site

site name	Caulerpa	DeadCoral-Framwork	DeadCoral-Rubble	Other	Pavona	Pocillopora	Porites	Rock	Sand	Tubastraea
Darwin- Wellington Reef	105	66	1	46	566	29	1108	855	177	9
Española- Xarifa	1	31	65	7	0	0	584	353	800	0
Floreana- Tres Cuevitas	0	11	12	12	731	1	0	1409	463	1
Floreana-Punta Cormorant	0	16	0	13	351	0	16	1317	485	1
Marchena -Roca Espejo	0	90	38	18	35	613	361	1063	265	1
Pinta - Cabo Ibetson	0	332	39	18	3	674	10	1423	64	0
Wolf - Corales	0	40	0	12	812	7	961	580	93	2
Site name	Unknown	Psammocora	Cube Counting FD	mean_slope	lacunarity	SurfaceArea	Lat	Lon	Depth
Darwin - Wellington Reef	129	0	2.267	0.108951	1.495	309	-91.995824	1.67824	15
Española - Xarifa	8	1	2.138	0.035663	1.626	185	-89.644622	-1.357863	5
Floreana- Tres Cuevas	130	8	2.229	0.101956	1.998	225	-90.408089	-1.236136	10
Floreana- Punta Cormorant	51	1	2.202	0.074877	1.615	276	-90.419325	-1.223949	12
Marchena - Roca Espejo	74	19	2.212	0.062882	1.458	258	-90.401241	0.312162	8
Pinta - Cabo Ibetson	71	18	2.208	0.062978	1.536	265	-90.720945	0.544216	8
Wolf - Corales	101	0	2.225	0.096196	1.684	261	-91.815861	1.386624	12

3.1. Community Metrics from 3D Annotations

First, we examined the relative abundance of each class. The results (Figure 5, A) show that Darwin (n=3091) is dominated by Porites (36%) and Pavona (18%). Española (n = 1850) is characterized by Sand (43%) and Porites (32%). Fl_cormorant_pt (n = 2251) is dominated by rock (58%) with contributions from Sand (21%) and Pavona (15%).Fl_tres_Cuevas (n=2778) shows a high abundance of rock (51%) and Pavona (26%). Marchena (n=2577) is characterized by rock (41%) and Pocillopora (24%), with higher evenness among coral classes. Pinta (n = 2652) has substantial rock cover (54%) and Pocillopora (25%). Wolf (n=2608) is dominated by Porites (37%) and Pavona (31%).

Coral diversity (Figure 5, B) was calculated based on the coral classes Pavona, Pocillopora, Porites, and Psammocora. The results show that Marchena has the highest diversity (0.87), followed by Darwin and Wolf (0.71). Española and Fl_tres_cuevas have the lowest diversity. When examining the relative abundance of coral classes, it shows that Española is dominated by Porites, while Fl_cormorant_pt and Fl_tres_cuevas are dominated by Pavona. In contrast, Marchena, Darwin, and Wolf exhibit higher evenness among coral classes.

The results of coral cover per site (Figure 5, C) varied from 14% in Fl_cormorant_pt to 42% in Wolf. Darwin had 38% coral cover, while Fl_tres_cuevas, Pinta, and Española had 21–24%. Marchena showed 29% coral cover.

3.2. Multidimensional Scaling

We performed a Principal Component Analysis (PCA) to reduce the dimensionality of the data (Figure 5, D) and visualize site similarity and class contributions. The first and second PCs explain 49% and 25% of the variation in the data. Caulerpa and Tubastraea classes were only apparent in Darwin and Unknown and Other classes are points that were not identified and therefore we dropped them from the PCA.

The PCA biplot (Figure 5, D) shows that Darwin and Wolf are close due to their higher proportions of Porites and Pavona. Floreana sites (Fl_tres_cuevas and Fl_cormorant_pt) cluster together and are affected by Pavona, Sand, and Rock. Pinta and Marchena group together due to higher proportions of Pocillopora, and Psammocora, Española is separate from all sites due to the high proportion of Sand.

Figure 5. A) Percent relative abundance of benthic classes across surveyed sites in the Galápagos Archipelago, highlighting site-specific trends in coral and non-coral categories. Darwin (n = 3091) is dominated by Porites and Pavona. Española (n = 1850) is characterized by Sand and Porites. Fl_cormorant_pt (n = 2251) is dominated by rock with contributions from Sand and Pavona. Fl_tres_cuevas (n = 2778) shows a high abundance of rock and Pavona. Marchena (n = 2577) is characterized by rock and Pocillopora, with higher evenness among coral classes. Pinta (n = 2652) has substantial rock cover and Pocillopora. Wolf (n = 2608) is dominated by Porites and Pavona. B) Shannon diversity index (H’) of coral genera, Marchena is the most diverse site followed by Darwin and Wolf, Española has the lowest diversity. C) Percentage of coral cover calculated as the proportion of annotations from Pavona, Pocillopora, Porites, Psammocora and Tubastraea out of all annotations, showing the highest cover in Wolf (42%) and Darwin (38%), and the lowest in Fl_cormorant_pt (14%). D) Principal Component Analysis (PCA) biplot illustrating the separation of reef sites based on their class compositions. Vectors represent the contributions of individual classes (e.g., Porites, Pavona, Pocillopora, etc.) to the variance in the data. The first and second PCs explain 49% and 25% of the variation in the data. Darwin and Wolf cluster close to one another, primarily influenced by the Porites and Pavona vectors. Floreana sites also group together, with the rock and Pavona vectors driving their composition. Pinta and Marchena are positioned further apart but are influenced by Psammocora and Pocillopora classes. Española is situated further from the other sites, with Sand contributing the most to its separation.

Figure 5. A) Percent relative abundance of benthic classes across surveyed sites in the Galápagos Archipelago, highlighting site-specific trends in coral and non-coral categories. Darwin (n = 3091) is dominated by Porites and Pavona. Española (n = 1850) is characterized by Sand and Porites. Fl_cormorant_pt (n = 2251) is dominated by rock with contributions from Sand and Pavona. Fl_tres_cuevas (n = 2778) shows a high abundance of rock and Pavona. Marchena (n = 2577) is characterized by rock and Pocillopora, with higher evenness among coral classes. Pinta (n = 2652) has substantial rock cover and Pocillopora. Wolf (n = 2608) is dominated by Porites and Pavona. B) Shannon diversity index (H’) of coral genera, Marchena is the most diverse site followed by Darwin and Wolf, Española has the lowest diversity. C) Percentage of coral cover calculated as the proportion of annotations from Pavona, Pocillopora, Porites, Psammocora and Tubastraea out of all annotations, showing the highest cover in Wolf (42%) and Darwin (38%), and the lowest in Fl_cormorant_pt (14%). D) Principal Component Analysis (PCA) biplot illustrating the separation of reef sites based on their class compositions. Vectors represent the contributions of individual classes (e.g., Porites, Pavona, Pocillopora, etc.) to the variance in the data. The first and second PCs explain 49% and 25% of the variation in the data. Darwin and Wolf cluster close to one another, primarily influenced by the Porites and Pavona vectors. Floreana sites also group together, with the rock and Pavona vectors driving their composition. Pinta and Marchena are positioned further apart but are influenced by Psammocora and Pocillopora classes. Española is situated further from the other sites, with Sand contributing the most to its separation.

3.3. Structural Complexity from Cube Counting and Walk Ratios

The results of these two metrics show similar trends. To compare them we normalized them using min-max scaling and plotted them against each other (Figure 6, C). Darwin, Fl_tres_cuevas, Wolf, and Española Maintain their ranks of structural complexity in both methods. However, the sites Marchena and Pinta have the same cube counting results and slightly different walk results. Fl_cormorant_pt has higher Geodesic walk scores than these two sites, but lower FD from cube counting.

These findings emphasize the importance of using complementary metrics to understand reef structural complexity.

4. Discussion

Galápagos coral reefs have experienced severe and recurrent environmental disturbances over the past four decades, yet the extent of their resilience remains uncertain. El Niño events of 1982 and 1997 [8,11] and a La Niña event in 2007 caused cold-water upwelling and bleaching [13]. More recently, the third global bleaching event (2014–2017) led to widespread coral mortality worldwide [65], although in the Galápagos reefs displayed some resilience, potentially due to localized upwelling [66] or possible shifts in symbiont associations [11]. The fourth global bleaching event was preceded by El Niño conditions [17] and declared in March 2024 [67]. Predictions indicated that the 2024 El Niño could transition to La Niña [68] possibly adding stress to Galápagos reefs through cold-water fluxes [5]. We initiated this program in 2024 to document the corals before the El Niño and during the transition to colder water from La Niña using new technologies. Our findings provide critical new insights into the spatial variation of coral cover, diversity, and structural complexity across seven reef sites, offering a valuable baseline for future ecological assessments.

This study demonstrates the advantages of photogrammetry over traditional chain transects and planar rugosity indices, which often underestimate structural complexity by failing to capture overhangs and crevices. The high-resolution 3D models produced here not only provide greater spatial accuracy but also enable automated analyses that can be replicated across time and space, facilitating long-term coral monitoring efforts.

High structural complexity is often correlated with increased biodiversity and reef resilience, as complex habitats provide refuges for juvenile corals and associated reef organisms. Our findings suggest that reefs with high fractal dimensions, such as those at Darwin and Wolf, may serve as ecological strongholds capable of supporting diverse marine life despite environmental stressors.

Protocols for annotating photogrammetric models from underwater environments usually use the photomosaic [50,51] which is a 2D representation of the 3D model. We aimed to segment the 3D model in order to obtain samples that are obscured in the photomosaic such as from overhangs and crevices. A similar method for pointcloud annotation was previously presented [52,69] however it requires additional software. Our workflow operates in Metashape and the annotations can be distributed online, e.g., on Labelbox [60]. Moreover, using the synthetic image (the 3D model view) improves the annotation accuracy and speed because it enables to view the point in context while the close-up images facilitate precise identification. Previous research on 3D semantic segmentation of coral reef models was done via deep learning and multiview images (e.g., [54,55,56,57]). In this work, we have taken the first step toward semantic segmentation of 3D models from the region by compiling a large dataset of annotations (see section Supp. 1.1 and Supp. 1.2). An important follow up study will focus on training a deep neural net based on the image pairs presented here and testing if the novel views increase the accuracy in predicting benthic classes.

Our results reveal significant variation in coral cover and diversity across sites. Northern sites, such as Darwin and Wolf, exhibited the highest coral cover, primarily dominated by Porites and Pavona species. In contrast, Marchena demonstrated the highest species diversity, with more even representation among Pavona, Pocillopora, Psammocora, and Porites. Floreana sites had relatively lower coral cover but exhibited high structural complexity, suggesting that these reefs may still provide essential habitat for reef-associated organisms despite their lower coral abundance. These patterns suggest that both geographical location and depth significantly influence coral community composition, although this was not directly tested. Our findings align with previous research, which documented higher coral cover in northern sites and patchier distributions in central and southern locations. However, our study expands upon prior work by incorporating 3D structural complexity assessments, which provide a more comprehensive understanding of reef habitat variation.

Previous studies in the region showed that coral cover is highest at the northern sites and that diversity however was highest at Marchena to the south. Coral cover at Darwin was quantified very close to our site using 0.25 m² photo quadrats sampling 37.5 m² per site (15 transects × 10 photos × 0.25 m²) [70]. The authors reported a total coral cover of 21.1%, with Porites contributing 19.5%, and Pavona and Pocillopora spp. each contributing less than 0.4%. By comparison, we sampled a 100 m² area at Darwin and found 41.7% coral cover, with Porites contributing 35.8%, Pavona 18.4%, and Pocillopora 1%.

Point-intercept methods along several 100 m transects per site at Wolf, Darwin, and Marchena were previously conducted [13]. At Wolf (Bahía Tiburón B), the authors reported 48% live coral, 25% rock, and 5% sand, with Porites and Pavona as the most abundant corals. For comparison, we found 42% coral cover, also dominated by Porites and Pavona. At Darwin (Anchorage Point 4), they observed 53% live coral, 37% dead coral, and 11% sand, with Porites and Pavona as the dominant genera. Our results were similar, with 38% coral cover dominated by Porites and Pavona. At Marchena, they noted high diversity and recorded 7% live coral, 34% rock, and 25% sand, with Porites, Pavona, and Pocillopora as the most abundant genera. In comparison, we found 29% live cover, dominated by Pocillopora and Porites. Moreover, we found that Fl_tres_cuevas had low coral cover and diversity but exhibited high structural metrics suggesting that it provides a rich habitat for reef organisms via diverse structural features.

Calculating fractal dimension from geodesic walks is a novel method developed for analyzing wide-scale reef models. Cube counting for fractal dimension was implemented by us previously on smaller models [37], however, it’s less effective for large area models because of their extensive lateral extent and relatively shallow vertical features, which result in larger cube sizes being mostly empty. In Figure 6, A and B appear different primarily due to the range on the y-axis. Since cube counting operates in 3D space its slope is steeper compared to walk ratios which operate in 1D. Darwin, Fl_tres_cuevas, and Wolf exhibit the highest fractal dimension from box counting and geodesic walks. Española has the lowest slope, suggesting a more flat reef structure. Geodesic walks enable to depict structural complexity in large area models effectively. This is evident in the results (Figure 6c) where geodesic walks maintain similar rank orderings among sites as cube counting. Nevertheless, there are some discrepancies, most notably between Pinta/Marchena and Fl_cormorant_pt. Here, cube counting overestimates structural complexity. The first two are relatively flat models abundant in Pocillopora. On the other hand, Fl_cormorant_pt has an abundance of Pavonas colonies scattered across the substrate in different sizes (Figure 2). Cube counting rated these sites high probably because of differences in the small box sizes.

Our application of novel photogrammetric techniques represents a significant methodological advancement in Galápagos coral research. Traditional reef complexity measurements, such as chain transects and planar rugosity indices, often fail to capture the full three-dimensionality of reef structures, especially overhangs and crevices where cryptic species reside. By implementing fractal dimension analysis through directed geodesic walks, we introduce a refined approach that better quantifies fine-scale reef complexity across multiple spatial scales. This method enhances our ability to compare reefs with different structural compositions and provides a robust metric for monitoring long-term changes in habitat complexity. We release a new dataset: Galápagos_3D, including the 3D models used in this study together with 3D annotations and classified image pairs. This data is useful for studies in ecology and computer graphics such as 3D segmentation from sparse labels via region growing [71].

While our study provides a valuable baseline, several limitations should be considered. First, although photogrammetry allows for detailed 3D modeling, manual annotation of benthic classes remains labor-intensive. Future research should explore machine-learning approaches to automate annotation, potentially improving efficiency and scalability. Additionally, while we documented spatial variation in reef complexity, temporal monitoring is needed to assess long-term resilience. Repeated surveys, particularly following extreme climatic events, will be crucial in understanding whether certain coral communities are developing adaptive responses to recurring stressors.

Given the increasing frequency of marine heatwaves and El Niño events, understanding which coral communities are most resilient and why will be critical for conservation planning. Our findings suggest that structural complexity may play a key role in reef persistence by providing refuges for coral recruits and associated marine life. Conservation efforts should prioritize sites with high coral cover and structural heterogeneity, as these areas may serve as ecological strongholds in the face of climate change.

Future studies should focus on longitudinal assessments, combining photogrammetry with genetic and physiological analyses to identify potential resilience mechanisms in Galápagos corals. Additionally, integrating satellite remote sensing data with underwater photogrammetry could enhance large-scale monitoring efforts, providing a more comprehensive picture of coral ecosystem health. In conclusion, this study highlights the importance of high-resolution 3D mapping for understanding coral reef dynamics and sets a new standard for reef monitoring in the Galápagos. By continuing this long-term research effort, we can contribute to evidence-based conservation policies, ensuring the preservation of these ecologically vital ecosystems for future generations.

• The authors used chatGPT and the Metashape online forum to generate the code and figures in the mansucript.