Very High-Resolution Satellite-Derived Bathymetry and Habitat Mapping Using Pleiades-1 and ICESat-2

: Accurate and reliable bathymetric data are needed for a wide diversity of marine research and management applications. Satellite-derived bathymetry represents a time saving method to map large shallow waters of remote regions compared to the current costly in situ measurement techniques. This study aims to create very high-resolution (VHR) bathymetry and habitat mapping in Mayotte island waters (Indian Ocean) by fusing 0.5 m Pleiades-1 passive multispectral imagery and active ICESat-2 LiDAR bathymetry. ICESat-2 georeferenced photons were ﬁltered to remove noise and corrected for water column refraction. The bathymetric point clouds were validated using the French naval hydrographic and oceanographic service Litto3D ® dataset and then used to calibrate the multispectral image to produce a digital depth model (DDM). The latter enabled the creation of a digital albedo model used to classify benthic habitats. ICESat-2 provided bathymetry down to 15 m depth with a vertical accuracy of bathymetry estimates reaching 0.89 m. The benthic habitats map produced using the maximum likelihood supervised classiﬁcation provided an overall accuracy of 96.62%. This study successfully produced a VHR DDM solely from satellite data. Digital models of higher accuracy were further discussed in the light of the recent and near-future launch of higher spectral and spatial resolution satellites.


Introduction
Mapping coastal areas is essential to tackle a broad range of environmental and social issues [1][2][3]. Therefore, a wide variety of scientific research disciplines could benefit from a better knowledge of this interface, especially regarding the monitoring and protection of coral reefs in archipelagos or the production of navigational charts [4].
Numerous reliable and accurate techniques exist to acquire bathymetric soundings. Data are often obtained through marine surveys equipped with multibeam or single-beam echosounders [5]. However, these approaches are usually impracticable in remote and shallow areas as well as time consuming and limited in terms of spatial coverage, and therefore remain costly. As an alternative, remote sensing is increasingly used to retrieve coastal bathymetry. Airborne data acquired with bathymetric LiDAR are useful to map larger areas but remain costly and limited spatially [6,7].
Over the past few years, satellite-derived bathymetry (SDB) has been increasingly used as it offers a more affordable and time saving alternative. Scientific studies have demonstrated the possibility of obtaining reliable bathymetric data through hyperspectral and multispectral (MS) imagery at various spatial resolutions, due to a correlation between water depth and reflectance data [8][9][10][11]. Nevertheless, SDB mostly relies on passive imagery, coefficient of 490 nm measured at 4 km resolution by the moderate-resolution imaging spectroradiometer (MODIS-Aqua, publicly accessible from https://oceancolor.gsfc.nasa.gov/l3/, last accessed: 28 December 2021) [32]. A diffuse attenuation coefficient value of 0.0615 m −1 was obtained for the month of May 2020 for the study site, indicating a very clear water type. Previous studies using ICESat-2 for bathymetric estimations had a diffuse attenuation coefficient ranging from 0.032 m −1 for the Virgin Island to 0.123 m −1 for the Bahamas, both known to be areas with very clear water [21,33].

Litto3D ® Reference Dataset
The French Oceanographic and Hydrographic Marine Service (SHOM) and the French National Institute for Geographical and Forest Information (IGN) conducted a joint altimetric and hydrographic survey of Mayotte from 2003 to 2010. Most of the island was mapped using the airborne topographic and bathymetric LiDAR and multibeam echosounder. The resulting Litto3D ® product provides soundings located in a three-dimensional geometric reference system with high spatial resolution and a land-sea continuum (data are available for free from https://diffusion.shom.fr/, last accessed: 28 December 2021). Data extracted from this dataset, corresponding to the ICESat-2 ground track, and used for comparison, include bathymetric points acquired using the bathymetric Li-DAR and multibeam echosounder.
Litto3D ® soundings are provided in a cartesian coordinate system in the horizontal plane and with orthometric heights. The point cloud density is constrained by the acquisition method over a specific area and the gridded model is provided with a spatial spacing of either 1 m or 5 m. Specification regarding the positioning and the geodesy of the Water clarity is a key parameter in SDB estimation. Clarity is related to light ray penetration in the water column, thus impacting the quality and quantity of the available bathymetric soundings [27][28][29][30][31]. This variable can be estimated using a diffuse attenuation coefficient of 490 nm measured at 4 km resolution by the moderate-resolution imaging spectroradiometer (MODIS-Aqua, publicly accessible from https://oceancolor.gsfc.nasa. gov/l3/, last accessed: 28 December 2021) [32]. A diffuse attenuation coefficient value of 0.0615 m −1 was obtained for the month of May 2020 for the study site, indicating a very clear water type. Previous studies using ICESat-2 for bathymetric estimations had a diffuse attenuation coefficient ranging from 0.032 m −1 for the Virgin Island to 0.123 m −1 for the Bahamas, both known to be areas with very clear water [21,33].

Litto3D ® Reference Dataset
The French Oceanographic and Hydrographic Marine Service (SHOM) and the French National Institute for Geographical and Forest Information (IGN) conducted a joint altimetric and hydrographic survey of Mayotte from 2003 to 2010. Most of the island was mapped using the airborne topographic and bathymetric LiDAR and multibeam echosounder. The resulting Litto3D ® product provides soundings located in a three-dimensional geometric reference system with high spatial resolution and a land-sea continuum (data are available for free from https://diffusion.shom.fr/, last accessed: 28 December 2021). Data extracted from this dataset, corresponding to the ICESat-2 ground track, and used for comparison, include bathymetric points acquired using the bathymetric LiDAR and multibeam echosounder.
Litto3D ® soundings are provided in a cartesian coordinate system in the horizontal plane and with orthometric heights. The point cloud density is constrained by the acquisition method over a specific area and the gridded model is provided with a spatial spacing of either 1 m or 5 m. Specification regarding the positioning and the geodesy of the Litto3D ® dataset are presented in Table 1 [34]. A local geoid, RGM04, was used as a reference for this dataset. This dataset was used during the validation phase to measure the accuracy of the SDB, but Litto3D ® data were not used as calibration points for the models.

ICESat-2 LiDAR Satellite Soundings
ICESat-2 is in a near-polar orbit at an altitude of 496 km and operates with a revisit period of 91 days over oceans [19,36]. ICESat-2 was mainly designed to measure icesheet topography, sea ice, and various inherent properties of the atmosphere and terrestrial vegetation, although ocean and inland surface waters are also observed. The Advanced Topographic Laser Altimeter (ATLAS), a photon-counting LiDAR, is the only sensor onboard the satellite, emitting a green laser beam at a wavelength of 532 nm. ATLAS enhances spatial sampling by splitting the laser beam into three pairs of beams separated by 3.3 km. Each pair, separated by 90 m, consists of a "weak" energy beam and a "strong" beam with a four-fold higher pulse energy [19,36]. ICESat-2 data can be downloaded with different degrees of processing, depending on the users' needs. This study uses the 3rd version of the L2 ATL03 georeferenced photons (data publicly available at https://search.earthdata.nasa.gov/search, last accessed: 28 December 2021) [37]. Data about each photon are provided with the latitude, the longitude, and the height relative to the WGS84 ellipsoid as well as other ancillary information. Considering that ICESat-2 was not designed to study the sub-surface water or the bottom topography, it is necessary to include in the analyses a correction for refraction bias induced by the water column.
We selected the ICESat-2 track acquired on the date closest to the acquisition date of the MS imagery. The two datasets were acquired 10 days, 10 h and 33 min apart. Then, the specific study area in Mayotte was selected based on the range of depths for which calibration data were available. ICESat-2 passed over Mayotte on 14 May 2020, at 20 h 51 min UTC, and collected bathymetric data down to a depth of 15 m.

Data Processing
Most of the ICESat-2 photons that reach the oceans penetrate into the water. However, compared to the water surface returns, only a small fraction is returned from the water column backscatter and bottom reflectance. Therefore, ICESat-2 signal photons correspond primarily to the water surface reflectance, water column backscatter, seabed reflectance, and noise.
ICESat-2 ATL03 data are provided with a preliminary classification of every photon regarding how likely it is to be signal or noise (confidence levels are: Noise, Low, Medium, High, and Buffer). Photons classified as "Buffer" are identified after all the signal photons are clustered. These are the photons for which doubt remains, which are at the limit to be identified as part of the signal. Therefore, this category has been created to ensure that all of the photons identified as signal are present in the corrected product [37]. Figure 2 presents the transect used in this study, indicating the original classification of the georeferenced photons. In this figure, photon positions (latitude and longitude) were projected onto a local geographic plane. Therefore, the horizontal axis corresponds to the along track distance. The origin point corresponds to the northernmost location of the trajectory. However, this clustering is not suited to underwater environments as it considers a considerable amount of the seafloor as noise. Therefore, all of the photons were considered and a modified density-based spatial clustering of application with noise (DBSCAN) algorithm was used to separate the photons characterizing the noise and the sea surface from those related to the seabed [17,[38][39][40].
ICESat-2 ATL03 data are provided with a preliminary classification of every photon regarding how likely it is to be signal or noise (confidence levels are: Noise, Low, Medium, High, and Buffer). Photons classified as "Buffer" are identified after all the signal photons are clustered. These are the photons for which doubt remains, which are at the limit to be identified as part of the signal. Therefore, this category has been created to ensure that all of the photons identified as signal are present in the corrected product [37]. Figure 2 presents the transect used in this study, indicating the original classification of the georeferenced photons. In this figure, photon positions (latitude and longitude) were projected onto a local geographic plane. Therefore, the horizontal axis corresponds to the along track distance. The origin point corresponds to the northernmost location of the trajectory. However, this clustering is not suited to underwater environments as it considers a considerable amount of the seafloor as noise. Therefore, all of the photons were considered and a modified density-based spatial clustering of application with noise (DBSCAN) algorithm was used to separate the photons characterizing the noise and the sea surface from those related to the seabed [17,[38][39][40].

Noise Removal and Detection of the Sea Surface
In the dataset, noise corresponds to sparse points with a low spatial density compared to the sea surface and seabed clusters. Georeferenced photons likely to be noise were removed, and photons associated with the sea surface were identified.
Here, a density-based spatial clustering method was used, which is an unsupervised learning method used to identify clusters in a dataset. The method is based on the premise that each cluster is defined as a region of points with a given density and spatially isolated from other groups by areas of lower density. The DBSCAN algorithm used scanned the entire dataset and established a search radius on each point successively. The point considered during a given step is a "core point". DBSCAN allows the users to specify a search radius size according to two criteria: The search circle radius ϵ and the minimum number of points MinPts. Once a criterion is no longer satisfied, the algorithm begins a new classification group [41].
Previous studies successfully implemented DBSCAN on the ICESat-2 dataset of islands located in the south of China and in the Bahamas. One particular study provides formulas to configure the MinPts and the ϵ radius parameters of the DBSCAN algorithm [17]. In the present research, the search radius was manually chosen by the user to guide

Noise Removal and Detection of the Sea Surface
In the dataset, noise corresponds to sparse points with a low spatial density compared to the sea surface and seabed clusters. Georeferenced photons likely to be noise were removed, and photons associated with the sea surface were identified.
Here, a density-based spatial clustering method was used, which is an unsupervised learning method used to identify clusters in a dataset. The method is based on the premise that each cluster is defined as a region of points with a given density and spatially isolated from other groups by areas of lower density. The DBSCAN algorithm used scanned the entire dataset and established a search radius on each point successively. The point considered during a given step is a "core point". DBSCAN allows the users to specify a search radius size according to two criteria: The search circle radius and the minimum number of points MinPts. Once a criterion is no longer satisfied, the algorithm begins a new classification group [41].
Previous studies successfully implemented DBSCAN on the ICESat-2 dataset of islands located in the south of China and in the Bahamas. One particular study provides formulas to configure the MinPts and the radius parameters of the DBSCAN algorithm [17]. In the present research, the search radius was manually chosen by the user to guide the clustering process and optimize the results. In addition, the MinPts parameter is defined by Equation (1) [17] (this formula is suited for a study of a water column whose depth is not expected to exceed 60 m). where SN 1 is the number of expected photons corresponding to signal and noise and defined by Equation (2): where N 1 is the total number of photons (both signal and noise), h is the vertical range and l is the along track range. SN 2 is the expected noise photons number and is defined by Equation (3): where N 2 corresponds to the number of photons in the layer with the fewer bathymetric photons, while h 2 is the height of the corresponding layer [17]. The variable MinPts is constrained to a value no lower than 3. If the previous formula provides a value lower than this threshold, then MinPts was set to 3 [17]. This algorithm might not be optimal in the present situation, since the dataset contains isolated photons from the seabed which could be identified as noise. Considering the small number of photons from the seabed, it was decided not to optimize the noise cleaning process, even if it meant that some manual cleaning had to be done. Therefore, the remaining noise points were removed manually using GlobalMapper software 22.1.0 (Blue Marble Geographics, Hallowell, ME, USA).

Detection of the Seabed
The sea surface is the cluster with the highest number of photons. It is a high-density group of photons spread over a continuous line, depending on the state of the sea. The sea surface cluster is clearly visible in blue in Figure 2.
According to [17], after removing the noise, the seabed is defined as every signal photon below a threshold value underneath the water surface. Therefore, every photon whose elevation is lower than LMS-3SV (where LMS is the Local Mean Sea level and SV the Surface Variance), was identified as a return signal from the seabed [17].

Correction for the Refraction Bias
Geolocated signal photons located below the water surface are not corrected for the refraction effect that redirects the light, and thus the LiDAR beams. It induces a positioning bias for photons in the water layer that would impact the bathymetry estimate. A relevant coordinate system is important to compute simple correction formulas. In this paper, the results presented were obtained using the coordinate system defined by [21] and the correction formulas were recomputed from this point. Corrections were applied in a satellite-centered coordinate system, where Z is the vertical direction (opposite to the direction of local gravity), and Y is orthogonal to Z (in the horizontal plane) and oriented along the azimuth of the pointing vector [21]. The resulting geometry is presented in Figure 3.

Validation of ICESat-2 Seabed Ellipsoidal Heights
The Litto3D ® dataset was used to validate the seabed photon ellipsoidal heights corrected from the refraction bias. While the ICESat-2 photons' geographic coordinates are projected onto a local tangent plane (ENU) during the refraction bias correction, the vertical references are different. Soundings measured by the SHOM bathymetric LiDAR are provided with orthometric heights. The latter were converted into ellipsoidal heights (relative to the IAG GRS 80 ellipsoid) using Circe 5.2.1 (IGN free software). Data visualization and extraction were conducted in GlobalMapper software, and the points were processed using Spyder 4.15 python interpreter (an open-source MIT environment) to qualify the accuracy. The point density from Litto3D ® is considerably higher compared to the ICESat-2 dataset, it was interpolated to allow a comparison with the ICESat-2 dataset. Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 24

Validation of ICESat-2 Seabed Ellipsoidal Heights
The Litto3D ® dataset was used to validate the seabed photon ellipsoidal heights corrected from the refraction bias. While the ICESat-2 photons' geographic coordinates are projected onto a local tangent plane (ENU) during the refraction bias correction, the vertical references are different. Soundings measured by the SHOM bathymetric LiDAR are provided with orthometric heights. The latter were converted into ellipsoidal heights (relative to the IAG GRS 80 ellipsoid) using Circe 5.2.1 (IGN free software). Data visualization and extraction were conducted in GlobalMapper software, and the points were processed using Spyder 4.15 python interpreter (an open-source MIT environment) to qualify the accuracy. The point density from Litto3D ® is considerably higher compared to the ICESat-2 dataset, it was interpolated to allow a comparison with the ICESat-2 dataset.

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][45][46][47][48][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

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][45][46][47][48][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]: where z is the depth, n is a constant needed for the ratio to stay positive, R w is the reflectance of the water, m 0 is the offset for a depth of 0 m, and m 1 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 R 2 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.

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. pected to be retrieved by this method independently of bottom albedo and can be obtained by inverting the radiative transfer equation as follows [23]: where z is the depth, n is a constant needed for the ratio to stay positive, Rw is the reflectance of the water, m0 is the offset for a depth of 0 m, and m1 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 R 2 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.

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. 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 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.

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, R 2 . 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.

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]: where A b is the bottom albedo, R w is the water column reflectance, R ∞ is the reflectance in deep water, and K d 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 K d 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 K d was computed for the three spectral bands in the visible range (Table 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). 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). 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). 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). 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).

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.

Correction of ICESat-2 Dataset
The ICESat-2 gt1l transect from the 2020 dataset corresponds to the strong beam and presents some variability in the depth range. The results obtained in this study, after removing the noise photons and correcting the signal for refraction bias, are presented in Figure 6. For comparison, the official signal detection and classification provided with the downloaded L2 ATL03 dataset were presented in Figure 2.
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).

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.

Correction of ICESat-2 Dataset
The ICESat-2 gt1l transect from the 2020 dataset corresponds to the strong beam and presents some variability in the depth range. The results obtained in this study, after removing the noise photons and correcting the signal for refraction bias, are presented in Figure 6. For comparison, the official signal detection and classification provided with the downloaded L2 ATL03 dataset were presented in Figure 2. In Figure 6, photon positions (latitude and longitude) were projected onto a local geographic plane. Therefore, the horizontal axis corresponds to the along track distance. The origin point corresponds to the northernmost location of the trajectory. For the current study, the analysis indicates that (1) points from the seabed detected in this study are following the bottom topography, distinguishable on the satellite imagery; (2) while the ATL03 data classified with a high and medium confidence level (corresponding to the blue and green points in Figure 2) are located at or close to the surface, our bathymetry algorithm correctly identifies the bottom topography from low and buffer confidence points, a much smaller portion of the returned signal; and finally (3) the ATL03 dataset removed only a small fraction of the noise photons in comparison to the results generated with the DBSCAN.
The DBSCAN algorithm was configured based on an empirical approach with the value = 0.65 m. This value, valid for our study area, allowed us to retrieve a majority of the seabed signal while eliminating most of the noise photons. The remaining noise points can be manually removed during the validation phase.

Validation of ICESat-2 Data
The comparison between ICESat-2 ellipsoidal heights and the Litto3D ® ellipsoidal heights is presented in Figure 7.
ATL03 data classified with a high and medium confidence level (corresponding to the blue and green points in Figure 2) are located at or close to the surface, our bathymetry algorithm correctly identifies the bottom topography from low and buffer confidence points, a much smaller portion of the returned signal; and finally (3) the ATL03 dataset removed only a small fraction of the noise photons in comparison to the results generated with the DBSCAN.
The DBSCAN algorithm was configured based on an empirical approach with the value ϵ = 0.65 m. This value, valid for our study area, allowed us to retrieve a majority of the seabed signal while eliminating most of the noise photons. The remaining noise points can be manually removed during the validation phase.

Validation of ICESat-2 Data
The comparison between ICESat-2 ellipsoidal heights and the Litto3D ® ellipsoidal heights is presented in Figure 7.   Figure 7 between the reference dataset and the corrected ICESat-2 data close to the surface.

Digital Depth Model
ICESat-2 corrected and validated data were used to calibrate the relative DDM. The calibration required the identification of a good model to bound the ICESat-2 dataset to the relative water depth values from the 0.5 m MS imagery.
Scaled pixel values were extracted from the relative water depth map derived from the 0.5 m MS imagery at the same location than ICESat-2 bathymetric points using QGIS 3.18.3 (open-source geographic information system). Several regression models were   Figure 7 between the reference dataset and the corrected ICESat-2 data close to the surface.

Digital Depth Model
ICESat-2 corrected and validated data were used to calibrate the relative DDM. The calibration required the identification of a good model to bound the ICESat-2 dataset to the relative water depth values from the 0.5 m MS imagery.
Scaled pixel values were extracted from the relative water depth map derived from the 0.5 m MS imagery at the same location than ICESat-2 bathymetric points using QGIS 3.18.3 (open-source geographic information system). Several regression models were tested, and the corresponding equations and their performance (in terms of the RMSE) are presented in Figure A1 (Appendix A). The regression chosen is a 2nd degree polynomial (Equation (6)). The expression of the latter is simple and performs well, with a RMSE of 0.895 m (see Figure 8). y = −45.87x 2 + 74.567x − 31.581 (6) Figure 9 shows the final DDM obtained using ENVI through the application of the model to the relative DDM. The hatched areas correspond to depths greater than the maximum calibration depth.
tested, and the corresponding equations and their performance (in terms of the RMSE) are presented in Figure A1 (Appendix A). The regression chosen is a 2nd degree polynomial (Equation (6)). The expression of the latter is simple and performs well, with a RMSE of 0.895 m (see Figure 8).  Figure 9 shows the final DDM obtained using ENVI through the application of the model to the relative DDM. The hatched areas correspond to depths greater than the maximum calibration depth.  Figure 9 shows the final DDM obtained using ENVI through the application of the model to the relative DDM. The hatched areas correspond to depths greater than the maximum calibration depth.

Digital Depth Model Validation
The reliability of the DDM was quantified by comparing the estimated bathymetry to the Litto3D ® reference dataset. Points from the bathymetric LiDAR point clouds were extracted, along the path of ICESat-2, from the DDM. The predicted RMSE was 0.895 m and the observed RMSE was 0.874 m along the ICESat-2 path. The predicted R 2 coefficient was 0.931 and the observed R 2 coefficient was 0.97. In addition, we report a MAE of 0.701 m.

Benthic Habitat Classification
Results from the four classifiers tested using different combinations of predictors are summarized in Figure 10, presenting both the overall accuracy and the kappa coefficient values.
The reliability of the DDM was quantified by comparing the estimated bathymetry to the Litto3D ® reference dataset. Points from the bathymetric LiDAR point clouds were extracted, along the path of ICESat-2, from the DDM. The predicted RMSE was 0.895 m and the observed RMSE was 0.874 m along the ICESat-2 path. The predicted R 2 coefficient was 0.931 and the observed R 2 coefficient was 0.97. In addition, we report a MAE of 0.701 m.

Benthic Habitat Classification
Results from the four classifiers tested using different combinations of predictors are summarized in Figure 10, presenting both the overall accuracy and the kappa coefficient values.
(a) (b) The NN algorithm configured with three hidden neurons for one layer systematically provided the lowest accuracy and the lowest kappa score with an overall accuracy which does not exceed 33.33% and kappa coefficient values that are all null. The neural network using only one hidden neuron produced better, yet inconclusive results, except for the combination of predictors "2".
The ML and the SVM algorithms generally produced the best classification results. The ML algorithm allows for a global accuracy of 96.62% and a kappa coefficient of 0.94 when using two spectral bands (Green and Blue) and with the addition of the three geomorphic predictors.
The SVM results were not as affected by the absence of the red spectral band. The results are constant and reached an overall accuracy of 96.50% and a kappa coefficient of 0.95 when using a DAM with three spectral bands and the slope as the only geomorphic predictor.
Maps of the benthic habitats with the best classification results for each depth range, are presented in Figure 11. The NN algorithm configured with three hidden neurons for one layer systematically provided the lowest accuracy and the lowest kappa score with an overall accuracy which does not exceed 33.33% and kappa coefficient values that are all null. The neural network using only one hidden neuron produced better, yet inconclusive results, except for the combination of predictors "2".
The ML and the SVM algorithms generally produced the best classification results. The ML algorithm allows for a global accuracy of 96.62% and a kappa coefficient of 0.94 when using two spectral bands (Green and Blue) and with the addition of the three geomorphic predictors.
The SVM results were not as affected by the absence of the red spectral band. The results are constant and reached an overall accuracy of 96.50% and a kappa coefficient of 0.95 when using a DAM with three spectral bands and the slope as the only geomorphic predictor.
Maps of the benthic habitats with the best classification results for each depth range, are presented in Figure 11.
These results, compared to the 0.5 m MS imagery and the DAM, are consistent with the classification presented in Figure 5, except for a bias in the classification of corals and algae, appearing in Figure 11a. The green band on the right of the image is an error in the classification process. These results, compared to the 0.5 m MS imagery and the DAM, are consistent with the classification presented in Figure 5, except for a bias in the classification of corals and algae, appearing in Figure 11a. The green band on the right of the image is an error in the classification process.

Bathymetric Errors
RMSE values comparing the DDM to the Litto3D ® dataset were obtained along the ICESat-2 ground track, thus at the same place used for the calibration. Here, we discuss the impact of both the depth values and the location of the validation transect on the results. Figure 12 shows the DDM from the Litto3D ® dataset. While the DDM produced in our study reaches 33 m depth, the Litto3D ® survey of this area indicated depths reaching at least 82 m.

Bathymetric Errors
RMSE values comparing the DDM to the Litto3D ® dataset were obtained along the ICESat-2 ground track, thus at the same place used for the calibration. Here, we discuss the impact of both the depth values and the location of the validation transect on the results. Figure 12 shows the DDM from the Litto3D ® dataset. While the DDM produced in our study reaches 33 m depth, the Litto3D ® survey of this area indicated depths reaching at least 82 m.
The map of the differences between the Litto3D ® and the DDM calibrated with ICESat-2 dataset is shown in Figure 13. The extrapolation works very well over the range of depths used in the calibration. However, the error increases at deeper depths. Once again, in this figure, the hatched areas correspond to depths greater than the maximum calibration depth of the model.
Large error values seem to appear in deeper waters (over 15 m depth), probably due to the fact that the bathymetry was calibrated with a limited depth range (the ICESat-2 dataset does not exceed a depth of 15 m). On the other hand, the fact data sampling was restricted to the ICESat-2 track is not ideal for the calibration. Most of the soundings represent a depth lower than 5 m (this concerns nearly 85% of the total amount of points for the 2020 dataset). The dataset only has a few points representing higher depth values. This was confirmed by the study of three other transects, taken at different places over the area (a transect along the ICESat-2 swath, a transect perpendicular to the swath, and a transect extracted far from the swath). The results of those tests are visible in Figures A2-A4   The map of the differences between the Litto3D ® and the DDM calibrated with I Sat-2 dataset is shown in Figure 13. The extrapolation works very well over the rang depths used in the calibration. However, the error increases at deeper depths. Once ag in this figure, the hatched areas correspond to depths greater than the maximum calib tion depth of the model.   The map of the differences between the Litto3D ® and the DDM calibrated with ICE-Sat-2 dataset is shown in Figure 13. The extrapolation works very well over the range of depths used in the calibration. However, the error increases at deeper depths. Once again, in this figure, the hatched areas correspond to depths greater than the maximum calibration depth of the model.  The 10-year difference in the time of acquisition between Litto3D ® and ICESat-2 could also induce a systematic error, due to a change in the bottom topography caused, for instance, by erosion or a change in the mean sea level (MSL), although those changes are likely to be well beyond the vertical accuracy provided by the method.
Mayotte has been prone to a succession of earthquakes since May 2018. The origin of these earthquakes is located to the east of the island. There are four permanent GPS stations in Mayotte and their rigorous monitoring has allowed experts to observe a displacement of all the stations by several centimeters towards the East and a subsidence of several centimeters since the beginning of the events [54][55][56].
In addition, the conversion of datums, using CIRCE software, could be a source of error. The grid used for the calculation is the GGM04V1 and the resulting vertical accuracy is estimated by the software at 10/20 cm.
On the other hand, during the correction of the refraction effect, n 1 and n 2 refractive indices were assumed and could therefore contribute to a small bias.
The method used to retrieve the map of the relative water depth could be improved to obtain more accurate DDM by implementing more recent and innovative approaches, such as IMBR, OBRA, MODPA or SMART-SDB [11,24,[57][58][59].

Impact of the Spatial Resolution of the Multispectral Imagery
In this study, the DDM was generated at the VHR of 0.5 m. However, other studies used sensors providing a spatial resolution of 10 m (Sentinel-2) or 30 m (Landsat-8) [17,39,60]. The spatial resolution drawn from the Pleiades-1A sensor was degraded in order to compare the RMSE. This process was done in ENVI with the "Resize Data" tool. The pixel size of the output was set according to the desired spatial resolution. The results are presented in Figure 14 and Table 4.
values, the results are similar, regardless of the transect location.
The 10-year difference in the time of acquisition between Litto3D ® and ICESat-2 could also induce a systematic error, due to a change in the bottom topography caused, for instance, by erosion or a change in the mean sea level (MSL), although those changes are likely to be well beyond the vertical accuracy provided by the method.
Mayotte has been prone to a succession of earthquakes since May 2018. The origin of these earthquakes is located to the east of the island. There are four permanent GPS stations in Mayotte and their rigorous monitoring has allowed experts to observe a displacement of all the stations by several centimeters towards the East and a subsidence of several centimeters since the beginning of the events [54][55][56].
In addition, the conversion of datums, using CIRCE software, could be a source of error. The grid used for the calculation is the GGM04V1 and the resulting vertical accuracy is estimated by the software at 10/20 cm.
On the other hand, during the correction of the refraction effect, n1 and n2 refractive indices were assumed and could therefore contribute to a small bias.
The method used to retrieve the map of the relative water depth could be improved to obtain more accurate DDM by implementing more recent and innovative approaches, such as IMBR, OBRA, MODPA or SMART-SDB [11,24,[57][58][59].

Impact of the Spatial Resolution of the Multispectral Imagery
In this study, the DDM was generated at the VHR of 0.5 m. However, other studies used sensors providing a spatial resolution of 10 m (Sentinel-2) or 30 m (Landsat-8) [17,39,60]. The spatial resolution drawn from the Pleiades-1A sensor was degraded in order to compare the RMSE. This process was done in ENVI with the "Resize Data" tool. The pixel size of the output was set according to the desired spatial resolution. The results are presented in Figure 14 and Table 4.    Figure 14 highlights the importance of the imagery spatial resolution on the accuracy of the bathymetry. RMSE remained under 1 m when the spatial resolution remained below 5 m, but increased rapidly after.
In other studies, the RMSE reached between 1.5 and 2 m for the Yongle atoll, located in South China [17]. It was 1.2 m on average in the Acklins islands in the Bahamas based on the Sentinel-2 MS satellite with a 10 m spatial resolution [17]. The RMSE was 0.96 m with the MS satellite Sentinel-2B (10 m spatial resolution) and 1.54 m using the Landsat-8 satellite (30 m spatial resolution) in the Virgin Islands [33]. Moreover, these study areas had a diffuse attenuation coefficient (K d ) similar to the Mayotte study area.
One of these studies obtained different results using Sentinel-2 observations and ICESat-2 observations from multiple swaths [60]. The DDM was produced with an extrapolation process conducted over the entire area with a RMSE of 3.36 m. However, when the study area was constrained between the two ICESat-2 swaths, the RMSE decreased to 0.35 m. This study area had a higher turbidity of K d = 1.68 m −1 [60].
During the acquisition time of some of these studies, meteorological events such as hurricanes occurred and could have impacted the topography of the bottom and affected the results. However, the evaluation of possible episodic events was not reported for this study or investigated.

Benthic Classification
Classification algorithms performed very differently. It is difficult to assess the impact of the geomorphic predictors on the results. It seems that adding extra information did not impact the SVM and ML classifications, but could have degraded the NN (+1HL) classification. A large amount information could have undermined the results due to a redundancy in the information. The major change seems to be related to the use of the red spectral band. The results were sometimes better without the red spectral band, probably due to the fact that the corresponding maps were in the depth range of 3.7-15 m, for which benthic classes, such as sand and coral rubble and rocks and coral rubble, are not present. A reduced number of groups tends to enhance the algorithm performance.
The benthic classification is based on a visual recognition of general benthic classes based on experts' knowledge. Although, commonly done, identifying benthic classes on MS imagery is not as reliable as direct underwater observations. Living corals could have been confused for dead corals colonized by algae. As a matter of fact, the classification presented in Figure 11a presented very good results, while having a major bias in the classification of coral and algae in areas of deep water. Some regions selected both to train the algorithm and for further validation presented corals which were distinguishable but very dark, due to the depth. Those were confused with deep and dark water areas.
Moreover, Mayotte is a complex area. The tidal range reaches 4 m. Therefore, when the tide is the lowest, corals can be above the water surface and bleached. Dead corals are theoretically recognizable by their bright color, but they can get darker as they are often colonized by algae. The winds and the waves have the effect to break coral colonies and to create coral rubble areas which are difficult to identify, as they get mixed with sand and rocks and can be mixed up with areas of isolated corals.

Conclusions
This study aimed to evaluate the quality of VHR DDM and DAM generated from satellite data. A DDM calibrated with data from the satellite ICESat-2 presented a RMSE of 0.89 m along ICESat-2 ground track, i.e., around 6% of the maximum depth retrieved by ICESat-2. Bathymetric results were generally satisfying down to a depth of around 15 m, which is close to the maximum depth of the calibration data used. Marine habitat classification results were very heterogeneous, depending on the number of predictors used, the type of predictors, and the algorithm used. However, some combinations of parameters provided satisfactory results. The classification with the ML classification using Blue and Green spectral bands with the three geomorphic predictors provided an overall accuracy of 96.62% and a κ coefficient of 0.94. In addition, the SVM classification using Blue, Green, and Red spectral bands with the addition of the slope geomorphic predictor presented an overall accuracy of 96.50% and a κ coefficient of 0.95. This approach can be of strong interest to map coastal areas lacking bathymetry and marine habitat maps and for which field observations are difficult.
While the quality of the results obtained in this study can support coastal management and conservation, the accuracy of bathymetry predictions remains limited for applications, such as navigation, that require higher spatial accuracy. It would be interesting to pursue this research to get more accurate DDMs.
Further work, implementing this method on diverse study sites, would confirm the robustness of the method implemented. In the prospect of future studies, it would be relevant to consider several ICESat-2 ground tracks from the area of interest and even to add the points from other times that ICESat-2 surveyed the area. This would provide a better variability of depths and a better spatial distribution of the data for the calibration process. Moreover, this increase in the number of points opens prospects for the use of deep learning methods to generate DDMs.
Developing an algorithm dedicated to the processing of seafloor data generated from ICESat-2 datasets would be important. The correction for the refraction effect has proven necessary and reliable, but could be further enhanced. The water column properties are changing with depth and the refraction correction should also adapt according to the water column properties.
It would be relevant to also improve the seabed signal correction by considering the state of the sea (for instance, presence of waves on the water surface), in helping to develop a method that could be used in less sheltered areas [17].
In this study, the results presented were obtained using a MS imagery acquired by the Pleiades-1A sensor with four spectral bands and a VHR of 0.5 m using the panchromatic band. The correlation between spatial resolution and the quality of the resulting bathymetry has been demonstrated in this paper. Therefore, future studies could consider generating better quality DDMs using the WV3 sensor (eight spectral bands at 0.30 m using the panchromatic band) or even the sensor of the new Pleiades Neo constellation launched in early 2021 (six spectral bands at 0.3 m with the panchromatic band).
The ICESat-2 products produced by NASA are constantly enhanced and one can be very optimistic regarding the future quality of DDMs and by-products obtained using ICESat-2 measurements.