*2.3. Satellite-Derived Bathymetry*

#### 2.3.1. Ratio Transform Method

The SDB ratio transform algorithm provided by ENVI 5.3 (L3Harris Geospatial Solutions, Broomfield, CO, USA) was used to retrieve the relative DDM of the study area. The DDM is based on the relationship between reflectance and bathymetry, which is described by [23]. The semi-empirical model provides values of relative bathymetry by computing the logarithmic ratio of the reflectance of two spectral bands from a MS imagery. First, the 2 m pixel size imagery was converted into TOA reflectance values and was geometrically projected to the WGS84/UTM38S. At this point, the spatial resolution was enhanced using the Gram-Schmidt pan-sharpening method [42,43]. Second, the MS image was cropped with a spatial subset tool to isolate the geographical area of interest. Finally, the ratio transform was implemented with the algorithm developed by [23]. A map of the relative water depth (i.e., log ratio of the spectral bands) was derived from the log ratio between the green and blue spectral bands [44–49].

This method is one of the most commonly used methods in SDB studies, as it proved to provide accurate results and does not require many points for the calibration phase [18,23]. The advantages are that only two parameters need to be set and it works on all types of albedos. This method is also mainly suited for clear case 1 water, which is the case in this study [18,23].

Working with spectral band ratios is a way to compensate for the variability of ocean bottom type, since changes in the albedo values will affect approximately equally both spectral bands. On the contrary, a variation of depth has a higher impact on the spectral band, which is the most intensely absorbed in the water column. Therefore, depth is expected to be retrieved by this method independently of bottom albedo and can be obtained by inverting the radiative transfer equation as follows [23]:

$$z = m\_1 \frac{\ln \left( nR\_{\text{av}} \left( \lambda\_j \right) \right)}{\ln \left( nR\_{\text{av}} \left( \lambda\_i \right) \right)} - m\_{0\prime} \tag{4}$$

where *z* is the depth, *n* is a constant needed for the ratio to stay positive, *Rw* is the reflectance of the water, *m*<sup>0</sup> is the offset for a depth of 0 m, and *m*<sup>1</sup> is the gain coefficient.

Each pixel of the MS imagery was assigned a value between 0 and 1. The final DDM was obtained by calibrating the relative bathymetry product with field-based depth measurements.

The final DDM was produced by finding the equation that best fits (i.e., lowest RMSE and higher R2 with a simple equation formula) the relationship between the relative bathymetric values and the ground truth depth measurements. If bathymetric soundings measured by ICESat-2 are reliable and accurate enough, they could be used as a calibration/validation dataset to produce bathymetry solely from satellite observations.

#### 2.3.2. Calibration with ICESat-2 Soundings

During the calibration phase, pixels from the relative DDM were matched to bathymetric points measured by ICESat-2.

To produce a DDM, i.e., a map of the water height at the acquisition date of the satellite MS imagery, ICESat-2 vertical heights were converted into the appropriate datum. First, the ellipsoidal heights measured by the ICESat-2 satellite were referenced to the IAG GRS 80 ellipsoid and had to be referenced to the chart datum. The SHOM (https://data.shom.fr/, last accessed: 28 December 2021) provides accurate altimetric information over Mayotte island, including the distance between the ellipsoid and the chart datum in Dzaoudzi locality (distance of −21.74 m). Second, the water height above the chart datum, at the acquisition time of the MS satellite imagery, was added. The closest tide gauge from the study site was also located at Dzaoudzi and the measurements were available from the SHOM website. The tide gauge recorded a water height of 1.02 m above the chart datum on 25 May 2020, at 07 h 24 min UTC. Figure 4 illustrates the different variables involved to compute the bathymetry from ICESat-2 ellipsoidal heights.

**Figure 4.** Variables involved in the process of retrieving ICESat-2 bathymetry at the acquisition time of Pleiades-1A.

Relative bathymetry points from the MS imagery were collected at the exact same location as the measurement points of ICESat-2 and an equation linking the two datasets was determined. Finally, the equation was applied to the relative DDM using the ENVI band math tool to generate the final DDM.

#### 2.3.3. Digital Depth Model Validation

We compared different regression models. The aim was to find the model that best matches the bathymetry measurements of ICESat-2 with the remotely sensed spectral values of Pleiades-1A. The best regression was chosen based on the RMSE and based on the coefficient of determination, R2. The final DDM was validated in comparison to the Litto3D® reference dataset by computing the root mean square error (RMSE) and the maximum absolute error (MAE) statistical indicators.

#### *2.4. Benthic Habitats Mapping*

#### 2.4.1. Processing of the Multispectral Imagery

This paper further intends to classify the seabed into five general types: Sand, sand with coral rubble, rock with coral rubble, corals and algae, and deep water. The classification is based upon a DAM obtained from the MS imagery and the DDM.

In theory, it would be feasible to train the classification algorithm directly on the MS imagery, without any preliminary corrections. It would also be conceivable to add a fifth spectral band, corresponding to the bathymetry, to add extra information for the algorithm to get better results. However, this method was not optimal as classifying MS imagery without correcting the image for the decay of light rays in the water column might induce confusion between the spectral signatures of benthic habitats [50]. A better solution was to generate a DAM. Bathymetric data are necessary to quantify the loss of light in the water column, to compensate for this loss, and finally to obtain a DAM [25,50]. Benthic albedo values were obtained with Equation (5) [25]:

$$A\_{\flat} = (R\_{\mathcal{w}} - R\_{\infty})e^{2K\_d z} + R\_{\infty} \tag{5}$$

where *Ab* is the bottom albedo, *Rw* is the water column reflectance, *R*<sup>∞</sup> is the reflectance in deep water, and *Kd* is the diffuse attenuation coefficient.

The TOA reflectance value for each spectral band was obtained after processing. First, the 2 m pixel size MS imagery was cropped to the area of interest and then orthorectified. Second, the image was converted into bottom of atmosphere (BOA) reflectance values using the FLAASH algorithm (see [50] for further details). At this point, it was possible to enhance the spatial resolution using pan-sharpening. These reflectance values were applied in Equation (5) to remove the water column contribution and obtain a bottom of hydrosphere (BOH) reflectance imagery from the BOA reflectance imagery.

The diffuse attenuation coefficient *Kd* was estimated for every spectral band of the visible range using values from a previous study on case 1 waters [51]. Values were given for a wide range of wavelength and were weighed with the appropriate factor found according to the wavelength sensitivities of Pleiades-1A sensor. Finally, an average value of *Kd* was computed for the three spectral bands in the visible range (Table 2).

**Table 2.** Diffuse attenuation coefficient (*Kd*) for every spectral band of Pleiades-1A sensor in the visible range.


#### 2.4.2. Supervised Classification Process

Similar to most coastal areas worldwide, the study area is lacking high-resolution benthic habitat data. Therefore, the supervised classification process was based on a visual identification of marine habitats assisted by two local marine scientists familiar with the area. The visual mapping by these experts was performed using the 0.5 m MS BOH reflectance imagery and the 0.5 m MS TOA reflectance imagery in parallel. Five benthic habitats were identified on several areas of the MS imagery (see Figure 5).

**Figure 5.** Maps of the five benthic habitats identified on the 0.5 m Pleiades-1A TOA reflectance imagery.

Corals are often visible on the edges of patch reefs. They appear in the form of brown and dark slender spots on the MS imagery and they are often colonized by algae. Sand areas are very bright areas often found on the border of the patch reefs, alongside corals. Two other main categories of benthic habitats can be distinguished. The first one corresponds to a mix of sand and coral rubbles. It appears as bright areas with dark brown spots. The other group contains mainly rocks and coral rubbles and appears as dark areas. While coral colonies are found on the edges, these two areas are often located towards the inner parts of the patch reefs. Coral rubbles are often transported towards the inner part of patch reefs by waves. Finally, both deep water areas, where the bottom was not visible, and land areas were masked.

In an attempt to enhance the classification accuracy, the albedo imagery was separated according to depth ranges. A first mask was created to suppress the depth values higher than the highest depth value measured by ICESat-2 over the area (i.e., 15 m). Then, the extinction depth of every spectral band was computed based on [50,51] (see Table 3 for results).

**Table 3.** Extinction depths for every spectral band of Pleiades-1A in the visible range.


A first DAM was created for depths in the range of 0–3.8 m using three spectral bands (Red, Green, and Blue). A second DAM was created for depths in the range of 3.8–15 m using two spectral bands (Green and Blue), as the extinction depth of the red band was exceeded. Due to the fact that some habitats are not present in the studied depth range, two sets of regions of interest (ROI) were created for each DAM. The DAM with the lower depth range was classified using four ROIs, namely: Sand, coral and algae, sand and coral

rubble, and rocks and coral rubble. The DAM with the higher depth range was classified using three ROIs, namely: Sand, coral and algae, and deep water.

Moreover, a choice of three morphological predictors was made to complement the MS bands of the DAM in order to enhance the classification results. The first predictor added was the slope before adding the aspect and the profile convexity all together. Classifications were computed using a 3 × 3 pixel kernel size.

Three classification algorithms were compared in this study: Neural network (NN), maximum likelihood (ML), and support vector machine (SVM).

#### 2.4.3. Validation of the Supervised Classification

A validation dataset based on the MS imagery was produced with the same knowledge as for the calibration phase. A post-classification accuracy assessment using a confusion matrix provided information on overall accuracy (OA) and the kappa coefficient (κ) [52,53].

The ML and the SVM classifiers were set with the default parameters. A neural network classification was implemented with one and three hidden neurons in one hidden layer, in order to test the depth of the neuronal architecture.
