Automated Global Shallow Water Bathymetry Mapping Using Google Earth Engine

: Global shallow water bathymetry maps offer critical information to inform activities such as scientiﬁc research, environment protection, and marine transportation. Methods that employ satellite-based bathymetric modeling provide an alternative to conventional shipborne measurements, offering high spatial resolution combined with extensive coverage. We developed an automated bathymetry mapping approach based on the Sentinel-2 surface reﬂectance dataset in Google Earth Engine. We created a new method for generating a clean-water mosaic and a tailored automatic bathymetric estimation algorithm. We then evaluated the performance of the models at six globally diverse sites (Heron Island, Australia; West Coast of Hawai’i Island, Hawai’i; Saona Island, Dominican Republic; Punta Cana, Dominican Republic; St. Croix, United States Virgin Islands; and The Grenadines) using 113,520 ﬁeld bathymetry sampling points. Our approach derived accurate bathymetry maps in shallow waters, with Root Mean Square Error (RMSE) values ranging from 1.2 to 1.9 m. This automatic, efﬁcient, and robust method was applied to map shallow water bathymetry at the global scale, especially in areas which have high biodiversity (i.e., coral reefs).

Historically, methods for estimating shallow water bathymetry have suffered from a variety of trade-offs and limitations. Conventional methods such as shipborne sounding or airborne LiDAR have limited spatial coverage [21]. Such methods cannot be deployed at global scales since they are both time-consuming and labor-intensive [19,22]. As a result, alternative methods using satellite data to derive bathymetry have been developed [23][24][25][26]. Past studies have applied low-spatial-resolution (>30 m) satellite sensors including NASA MODIS, SeaWiFS, and VIIRS [27][28][29][30], as well as land-viewing satellite imagery including Landsat-8, Sentinel-2, and Worldview 2-4 [22,23,[31][32][33][34] to estimate bathymetry. However, these methods involve trade-offs between spatial resolution and satellite-based observation frequency. Low spatial-resolution sensors can lead to large uncertainties in land-water mixed pixels along coastlines, while land-viewing sensors can have a low temporal frequency that results in a more limited selection of cloud-free satellite images, especially in tropical coastal regions with frequent cloud cover [4]. Furthermore, satellite-based bathymetric algorithms often rely on field data calibration [19,20,22,35,36] or intensive calculations to physically simulate broad water column conditions [21,23,26,[31][32][33]37,38], and are often limited to a single scene of multispectral or hyperspectral satellite images. Overall, these many drawbacks limit the mapping of shallow water bathymetry at a global scale.
To overcome many of these challenges in estimating shallow water bathymetry, we investigate the use of Google Earth Engine (GEE), a powerful cloud-based computational platform that provides easy access to high-coverage global-scale analysis-ready satellite reflectance datasets and high-resolution innovation to bathymetry estimation. GEE has been used for a variety of global-scale products including land cover, forest change, water surface extent, and urban land use [39][40][41][42][43]. In particular, the Sentinel-2 surface reflectance dataset in GEE has been widely used in terrestrial environment studies [22,43], but it is not fully understood if the Sentinel-2 reflectance dataset in GEE can be directly used to generate shallow water bathymetry at the global scale. Overall, an automated (i.e., no field calibration) and computationally efficient shallow water bathymetry estimation algorithm is needed for the Sentinel-2 reflectance dataset in GEE at the global scale.
Here, we developed a new method to automatically map global shallow water depth using Google Earth Engine in both coastal and offshore marine environments. We used the Sentinel-2 surface reflectance dataset in the GEE to build a "clean-water" mosaic that minimizes cloud, cloud shadow, sun glint, and other disturbances. This clean-water mosaic enabled us to derive a 10 m spatial resolution bathymetry in shallow marine environments (< 20 m) without field data calibration. We evaluated the performance of our bathymetric estimation method using 113,520 field sampling points across six globally distributed sites that represent different depths and diverse benthic habitats. Our study strives to resolve the challenge of developing an automatic, efficient, and robust method to map bathymetry at the global scale.

Study Sites and Data
We tested our automatic bathymetric mapping method across broad geographic ranges, benthic types, and water column conditions at six study sites worldwide. Figure 1 presents Sentinel-2 true color mosaics and field sample locations that were used for our study sites in the Pacific Ocean and the Caribbean Sea: Heron Island, Australia; West Coast of Hawai'i Island, Hawai'i (hereafter "West Hawai'i"); Saona Island, Dominican Republic; Punta Cana, Dominican Republic; St. Croix, U.S. Virgin Islands; The Grenadines (Table 1, Figure 1). These sites represent a variety of geomorphic zones (reef crest, patch reef, lagoon, and open ocean) and benthic types (coral reef, seagrass, sand). These sites provide a wide range of bottom reflectance signatures and geographic conditions for examining our automatic shallow water bathymetry mapping method. . Sentinel-2 true color mosaics (i.e., 12-month mosaic) and field sample locations for the six study sites. We used lidar-derived bathymetry in the Big Island, Hawaiʻi (plotted in pink polygon). Other bathymetry samplings were measured at consecutive points. Sentinel-2 true color mosaics (i.e., 12-month mosaic) and field sample locations for the six study sites. We used lidar-derived bathymetry in the Big Island, Hawai'i (plotted in pink polygon). Other bathymetry samplings were measured at consecutive points.
We utilized field-measured bathymetry samples to verify satellite-derived bathymetry results (Table 1). In total, we collected 113,520 field sampling points to globally validate our bathymetry method. To supplement depth measurements, we also recorded benthic composition types in the Dominican Republic and the US Virgin Islands [25].

Clean Water Mosaic Generation Using Google Earth Engine
We developed a method to build a Sentinel-2 surface reflectance mosaic for shallow waters using the Reducer function in GEE. We selected satellite images with minimal cloud Remote Sens. 2021, 13, 1469 4 of 17 coverage, sun glint, and water turbidity over a selected period (e.g., 12 months). We then aggregated the input dataset into a single, clean-water mosaic output. This step is fundamental for building an automatic bathymetry mapping method at a global scale. A clean water mosaic with minimal water column attenuation enables our automatic bathymetry estimation algorithm to have reduced uncertainties caused by water column attenuations.
We applied different filters to mask low-quality (e.g., cloud, cloud shadow, wave breaks, sun glint, etc.) Sentinel-2 surface reflectance (ρ(λ)) in the input raster dataset. First, we used the QA60 band to exclude pixels with clouds. Next, we applied the Scene Classification map (SCL) band to mask non-water targets including cloud shadows, vegetation, bare soil, clouds of medium probability, clouds of high probability, and cirrus. We also used several band threshold values to mask high turbidity waters, sun glint and wave breaks after intensive testing (i.e., green band > 0.01, red edge 1 band < 0.1, NIR band < 0.03, 0.005 < water vapor band < 0.03). We further masked non-water targets by using a normalized difference water index (NDWI) with surface reflectance inputs [44], using only pixels with positive NDWI values in the input datasets: Finally, we aggregated the selected clean satellite reflectance by taking the median of the pixel stack (i.e., median statistic of the Reducer function). We built mosaics over different time ranges (i.e., three months (Jan 2019 to March 2019), six months (Jan 2019 to June 2019), 12 months (Jan 2019 to December 2019)) to compare our method by mosaic time period.

An Automatic Bathymetry Estimation Algorithm
We developed a new automatic bathymetry mapping method based on a previous single-scene adaptive bathymetry algorithm [25]. Our algorithm was tailored to the clean water mosaic built by GEE. Figure 2 illustrates the detailed steps of our approach. We first calculated remote sensing reflectance R rs from the mosaic surface reflectance ρ(λ) as [45]: Next, we derived below-surface remote sensing reflectance (r rs (λ)) from the R rs (λ) to remove the air-water surface effect [38]: We estimated shallow water bathymetry by quantifying different attenuation levels between the blue and green bands as [19]: The bathymetry estimation parameters (m 0 and m 1 ) were calculated using Chlorophylla (Chl-a) concentration values as [25] representative for clean offshore waters:

Clean Shallow Water Mosaic Created from Google Earth Engine
The mosaics over three different time ranges (i.e., 3, 6, and 12 months) were created using GEE and provided clean shallow water composites without clouds, breaking waves, and sun glint. Figure 3 illustrates our clean-water mosaic results across varying time ranges within three of the six sites. Overall, benthic habitats, such as coral reef, sandy bottom, hard bottom, and seagrass beds, are clearly identifiable in the mosaics (i.e., shallow waters of St. Croix, USVI).
Selecting scenes over a longer mosaic time range resulted in an overall higher quality mosaic (i.e., lower numbers of missing points) ( Figure 3). We found that the 12-month mosaic had many fewer missing point quantities than the other mosaics in all three sites, Figure 2. Google Earth Engine global bathymetry estimation diagram. Sentinel-2 satellite reflectance images were selected to build the clean water mosaic to derive bathymetry.
As noted previously, we only selected satellite images with low water turbidity to build the mosaic. Additionally, given that water mosaic values represent the median value over time (e.g., 12 months), we used a fixed Chl-a value (Chl-a = 0.5 mg m-3) to calculate m 0 and m 1 in our clean shallow water mosaic. This Chl-a value (Chl-a = 0.5 mg m-3) is a mean value calculated from GEE outputted 12-month clean water mosaic in 26 sites globally (Chl-a ranges in these 26 sites: 0.4 mg m-3 < Chl-a < 0.6 mg m-3

Clean Shallow Water Mosaic Created from Google Earth Engine
The mosaics over three different time ranges (i.e., 3, 6, and 12 months) were created using GEE and provided clean shallow water composites without clouds, breaking waves, and sun glint. Figure 3 illustrates our clean-water mosaic results across varying time ranges within three of the six sites. Overall, benthic habitats, such as coral reef, sandy bottom, hard bottom, and seagrass beds, are clearly identifiable in the mosaics (i.e., shallow waters of St. Croix, USVI).
where missing points are areas where no reflectance pixels were selected in the time range to build the mosaic. While image quality was similar when comparing between the sixmonth and three-month mosaics, performance was much improved going from six to twelve months. Missing points were primarily concentrated in near-coastal regions in the Dominican Republic and St. Croix, USVI, and in the latter, almost all missing points were located along the coastline. In West Hawaiʻi, missing points were located primarily in the open ocean. To examine mosaic reflectance, we plotted surface reflectance values of the 12-month Dominican Republic mosaic using different benthic targets (Figure 4), since the 12-month mosaic has the best quality. Using the field sample locations, we extracted mean spectral reflectance values of different benthic targets at similar depths (Gorgonian and Reef, Sand and Seagrass at 6 m depth, and Ocean deeper than 10 m depth, Figure 4a), and compared Selecting scenes over a longer mosaic time range resulted in an overall higher quality mosaic (i.e., lower numbers of missing points) ( Figure 3). We found that the 12-month mosaic had many fewer missing point quantities than the other mosaics in all three sites, where missing points are areas where no reflectance pixels were selected in the time range to build the mosaic. While image quality was similar when comparing between the sixmonth and three-month mosaics, performance was much improved going from six to twelve months. Missing points were primarily concentrated in near-coastal regions in the Dominican Republic and St. Croix, USVI, and in the latter, almost all missing points were located along the coastline. In West Hawai'i, missing points were located primarily in the open ocean.
To examine mosaic reflectance, we plotted surface reflectance values of the 12-month Dominican Republic mosaic using different benthic targets (Figure 4), since the 12-month mosaic has the best quality. Using the field sample locations, we extracted mean spectral reflectance values of different benthic targets at similar depths (Gorgonian and Reef, Sand and Seagrass at 6 m depth, and Ocean deeper than 10 m depth, Figure 4a), and compared them to the same benthic targets across different depths (Figure 4b). The spectral shape and corresponding values represented the signal variations across different benthic targets and varying depths. For instance, a sandy bottom was characterized as having the highest average reflectance values, while seagrass beds and ocean (dark targets) had the lowest average values. A decrease in surface reflectance values from the red to NIR bands was observed across all targets. As expected, we also observed a decreasing pattern in spectral values from shallow to deep water across the same targets (Gorgonian and Reef, Figure 4b).
Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 18 them to the same benthic targets across different depths (Figure 4b). The spectral shape and corresponding values represented the signal variations across different benthic targets and varying depths. For instance, a sandy bottom was characterized as having the highest average reflectance values, while seagrass beds and ocean (dark targets) had the lowest average values. A decrease in surface reflectance values from the red to NIR bands was observed across all targets. As expected, we also observed a decreasing pattern in spectral values from shallow to deep water across the same targets (Gorgonian and Reef, Figure 4b).  Figure 5 shows the results of bathymetry maps generated using 12-month mosaics that were created in GEE. These bathymetric models exhibited expected spatial patterns, such as bathymetry values increasing from shallow to deep waters in the direction from the coastline to the ocean. This pattern was consistent across all five land-linked study  Figure 5 shows the results of bathymetry maps generated using 12-month mosaics that were created in GEE. These bathymetric models exhibited expected spatial patterns, such as bathymetry values increasing from shallow to deep waters in the direction from the coastline to the ocean. This pattern was consistent across all five land-linked study areas. For example, in Punta Cana, Dominican Republic, shallow depths (<5 m) were observed nearshore, transitioning to medium depths (5-10 m) at a distance of~1 km from land, followed by deeper water (>12 m) in the open ocean. In Heron Island, Australia (nearly 80 km north east of Queensland, Australia), we found the expected depth pattern in a reef area. It showed the bathymetry increased from 3 to 15 m between the reef rim and open ocean. In the reef rim regions, depth was observed less than 3 m. By contrast, beyond the reef rim regions, depth sharply increased to over 15 m. Moreover, our bathymetry maps showed constant high values (e.g., depths > 15 m) in the optically deep ocean waters.

Bathymetry Spatial Variation Analysis
Each of the bathymetry maps that were created using the different temporal mosaic composites (i.e., 3-, 6-, and 12-month mosaic) were investigated to determine the effect the mosaic time period had on the estimated bathymetry values ( Figure 6). Consistent with the surface reflectance mosaic results, the 12-month mosaic provided the lowest number of missing data pixels and the cleanest bathymetric spatial patterns ( Figure 6).

Evaluation of Bathymetry Estimation
We compared depth point values from the satellite bathymetry map with field depth sample locations. Table 2 shows four error metrics capturing differences between the satellite-derived bathymetry and field-derived depth values. In the 12-month mosaic, our approach effectively generated depth values with a Root Mean Square Error (RMSE) ranging from 1.  (Table 2).    . Sentinel-2 bathymetric models based on different temporal mosaic composites. As the time frame selection is increased, the quality of the mosaic also increases, resulting in fewer data gaps and more accurate bathymetric models. The white color represents the missing points in the bathymetry results.

Evaluation of Bathymetry Estimation
We compared depth point values from the satellite bathymetry map with field depth sample locations. Table 2 shows four error metrics capturing differences between the satellite-derived bathymetry and field-derived depth values. In the 12-month mosaic, our approach effectively generated depth values with a Root Mean Square Error (RMSE) ranging from 1.  . Sentinel-2 bathymetric models based on different temporal mosaic composites. As the time frame selection is increased, the quality of the mosaic also increases, resulting in fewer data gaps and more accurate bathymetric models. The white color represents the missing points in the bathymetry results.
We compared 12-month mosaic-derived depth with field depth samples at six different sites (Figure 7). We found the satellite-derived depth was highly correlated with field depth samples across all sites. As shown in the plots, the majority of points follow the 1:1 line. Overall, R 2 values ranged from 0.78 to 0.98 (Table 2)

Bathymetry Estimation Performance Impacted by Mosaic, Depth and Bottom Conditions
We evaluated the performance of our approach across different parameters including mosaic time range (three months, six months and 12 months), depth range (0 to 5 m, 5 to 10 m, 10 to 15 m and 15 to 20 m), and bottom benthic types (coral reef, sand, seagrass, and ocean). We found that longer-time frame mosaics offered better results in all error metrics ( Table 2). 12-month mosaics had 20% lower RMSE values and better Bias, MNB and R 2 values compared with six-month and three-month mosaics. Performance (e.g., RMSE and R 2 values) did not dramatically improve between three-and six-month mosaics. Overall, a majority of the mosaic depth results had the same Bias and MNB trend (e.g., over-estimation: Saona Island; under-estimation: Heron, The Grenadines).
We compared satellite-derived depths with field-measured values along two transects in St. Croix and Heron Island (Figure 8). These transects traversed different geomorphologic regions and benthic types, and therefore provided variable environments for testing performance. Depth results from all period mosaics and study sites followed the spatial trend of field measurements. As expected, 12-month mosaics showed the best correlation with field measurements, while three-and six-month mosaics resulted in a lower correlation. For example, three-month mosaics slightly over-estimated (1 to 2 m higher) values at edge points of the transect in St. Croix, and underestimated a discrete section of the transect at Heron Island.

Bathymetry Estimation Performance Impacted by Mosaic, Depth and Bottom Conditions
We evaluated the performance of our approach across different parameters including mosaic time range (three months, six months and 12 months), depth range (0 to 5 m, 5 to 10 m, 10 to 15 m and 15 to 20 m), and bottom benthic types (coral reef, sand, seagrass, and ocean). We found that longer-time frame mosaics offered better results in all error metrics (Table 2). 12-month mosaics had 20% lower RMSE values and better Bias, MNB and R 2 values compared with six-month and three-month mosaics. Performance (e.g., RMSE and R 2 values) did not dramatically improve between three-and six-month mosaics. Overall, a majority of the mosaic depth results had the same Bias and MNB trend (e.g., over-estimation: Saona Island; under-estimation: Heron, The Grenadines).
We compared satellite-derived depths with field-measured values along two transects in St. Croix and Heron Island (Figure 8). These transects traversed different geomorphologic regions and benthic types, and therefore provided variable environments for testing performance. Depth results from all period mosaics and study sites followed the spatial trend of field measurements. As expected, 12-month mosaics showed the best correlation with field measurements, while three-and six-month mosaics resulted in a lower correlation. For example, three-month mosaics slightly over-estimated (1 to 2 m higher) values at edge points of the transect in St. Croix, and underestimated a discrete section of the transect at Heron Island.
To further analyze how our method performed at different depths, we calculated RMSE values of bathymetry results (12-month mosaics) in different depth groups (0 to 5 m, 5 to 10 m, 10 to 15 m and 15 to 20 m) according to the field depth records (Figure 9). We found that RMSE increases as depth increased. Our method performed better in shallow To further analyze how our method performed at different depths, we calculated RMSE values of bathymetry results (12-month mosaics) in different depth groups (0 to 5 m, 5 to 10 m, 10 to 15 m and 15 to 20 m) according to the field depth records (Figure 9). We found that RMSE increases as depth increased. Our method performed better in shallow (0 to 5 m, 5 to 10 m) and medium depth waters (10 to 15 m), and worse in deeper Multiple benthic habitats at the study sites allowed us to further analyze algorithm performance by combining mosaic time range and bottom reflectance types ( Figure 10). In St. Croix, USVI and the Dominican Republic, we compared different mosaic derived depths with field measured data according to different benthic types (coral reef, sand, seagrass, and ocean). Sand and seagrass showed stronger correlations between satellite and field data compared to other benthic habitats. Ocean points had the lowest correlated pattern of all benthic types. Additionally, increasing the mosaic time range from three to twelve months clearly resulted in more accurate depth points for coral reef, sand, and seagrass benthic types, based on a comparison of results to the 1:1 line. In the 12-month mosaics, we found the trend line of coral reef, sand, and seagrass followed the 1:1 line more accurately than the three-month and six-month mosaic. Increasing mosaic time ranges showed the lowest improvement in ocean points than in coral reef, sand and seagrass.  Multiple benthic habitats at the study sites allowed us to further analyze algorithm performance by combining mosaic time range and bottom reflectance types ( Figure 10). In St. Croix, USVI and the Dominican Republic, we compared different mosaic derived depths with field measured data according to different benthic types (coral reef, sand, seagrass, and ocean). Sand and seagrass showed stronger correlations between satellite and field data compared to other benthic habitats. Ocean points had the lowest correlated pattern of all benthic types. Additionally, increasing the mosaic time range from three to twelve months clearly resulted in more accurate depth points for coral reef, sand, and seagrass benthic types, based on a comparison of results to the 1:1 line. In the 12-month mosaics, we found the trend line of coral reef, sand, and seagrass followed the 1:1 line more accurately than the three-month and six-month mosaic. Increasing mosaic time ranges showed the lowest improvement in ocean points than in coral reef, sand and seagrass.

Discussion
Today, open available global bathymetry datasets (e.g., GEBCO's gridded bathymetric data set) are often designed for the deep ocean and have coarse spatial resolutions. In Figure 10. Results performance in coral reef, sand, seagrass and ocean in three different mosaic time frames. The trend line of 12-month mosaic followed the 1:1 line more accurately than the three-month and six-month mosaics.

Discussion
Today, open available global bathymetry datasets (e.g., GEBCO's gridded bathymetric data set) are often designed for the deep ocean and have coarse spatial resolutions. In this study, we developed an automatic bathymetric mapping approach for shallow water in coastal and offshore marine environments at the global scale using Sentinel-2 satellite images (10 m resolution). In particular, our method can efficiently derive water bathymetry in high biodiversity benthic habitats (i.e., coral reefs). Our approach is based on a new clean-water mosaic built in GEE, which allowed us to directly use the process-ready Sentinel-2 surface reflectance datasets in GEE [44]. Although GEE's Sentinel-2 reflectance datasets are mostly designed for land application, we demonstrated the novel application of Sentinel-2 to water-color studies by filtering out clouds, cloud shadow, break waves, and high-turbidity waters. We largely reduced the uncertainties from atmospheric correction for turbid waters, cloud shadow, and waves [45,46]. Once these errors were removed, we aggregated selected pixels into a single clean-water mosaic. Therefore, our method is suitable for the sites where the bottom is visible in the satellite images during a year. As demonstrated by the resulting clean surface reflectance mosaics ( Figure 3) and accurate reflectance spectral data for benthic targets (Figure 4), our GEE clean-water mosaic provided reliable Sentinel-2 imagery for bathymetric modeling. All processing was conducted in the GEE platform which took advantage of powerful cloud computation capabilities [47].
Our clean-water mosaic using Sentinel-2 provides a reliable method for creating highly detailed information in multiple types of coastal environment studies, including bathymetry to support mapping and monitoring of coastal benthic habitats (coral reef, seagrass, and algae cover), marine transportation navigation, and physical modeling of waves and flooding [4,19,20,23,24,32,[48][49][50][51][52][53]. The amount of time incorporated into the mosaics determined the image quality since a longer time period allows more clean-water pixels to enter the aggregation process. Therefore, 12-month mosaics showed the highest quality, especially near the coastline (wave-and bubble-enriched regions) and ocean (high cloud/cloud shadow regions) which offered the best bathymetry results [54][55][56]. However, three-month mosaics provided reasonable reflectance values and can be usefully applied to seasonal coastal environmental monitoring, such as tracking benthic habitat seasonal changes [57][58][59].
We generated bathymetry through a tailored clean water bathymetric estimation algorithm. In previous studies, the dynamics of water column attenuation have posed a challenge in bathymetry estimation. Thus, previous bathymetric algorithms have relied on field calibration or complex physical modeling to simulate a broad range of water conditions [19,51,60]. However, these methods are not sufficient to overcome the dynamics of water column attenuation at a global scale. Our new algorithm is applied to a cleanwater mosaic and was therefore developed for low-attenuation waters. We used a large number of field depth measurements (113,520) to validate our approach across a wide range of geographic locations and benthic habitat types. The field sampling points are in low latitude regions in which the benthic geomorphology showed distinct variations in shallow depth (i.e., from lagoon to reef crest). Therefore, the bathymetry variations can be validated in a diverse range of geomorphology. The high accuracy (RMSE ranges from 1.26 to 1.92 m) of our automatic method demonstrated that the fixed water attenuation condition index (Chl-a values related m 0 and m 1 ) is sufficient to generate bathymetry.
We evaluated the performance of our approach under different depth and bottom reflectance conditions to uncover strengths and weaknesses. Generally, our approach works best in shallow to moderate water depths (0 to 15 m). Our algorithm was designed based on the quantification of water attenuation differences between blue and green bands [25]. In the same water condition, high water attenuation in deep water (depth > 15 m) had weaker differences than shallow water (depth < 10 m) [23]. Users could easily tailor our approach for specific applications to further reduce uncertainty or constrain the estimation bounds. For example, this could include adjacency correction of depth values, or additional pixel masking during mosaic generation to exclude pixels with a low noise equivalent difference in reflectance [23,31]. Different benthic habitats have different strengths of bottom reflectance [24,59].

Conclusions
We developed an automated bathymetric mapping method using the Sentinel-2 surface reflectance dataset in Google Earth Engine. A new clean water mosaic creation method and a tailored bathymetry estimation algorithm were designed. We tested our method in six diverse globally distributed sites with abundant field bathymetry sampling points. Our automatic and efficient method could be applied to map shallow water bathymetry at a global scale.