Estimating benthic reflectance of deep coral atoll lagoons from satellite imagery and bathymetry - Analysis code and case studies (NESP MaC 2.3, AIMS)

This code repository and dataset details a method for determining benthic reflectance from a combination of satellite imagery and bathymetry. Its key benefit is that it can map benthic reflectance up to 50 - 60 m in depth in clear waters, when using Sentinel 2 B2 channel combined with matching bathymetry data. Benthic reflectance is a measure of how much light the seafloor reflects and is useful for distinguishing areas that are sand (high reflectance) or vegetation such as seagrass, algae and coral (low reflectance). 

The proposed method allows the reflectance to be estimated much deeper than existing multi-spectral approaches that rely solely on the satellite imagery. These typically only work to depths of 10 - 20 m as they require reflected light across a wide spectral range to disentangle the depth and reflectance. In deeper areas only the blue end of the spectrum remains, making it ambiguous whether an area is dark because it is deep or dark because the substrate is dark. We resolve this ambiguity by using an independent bathymetry digital elevation model. 

This method only works in regions where all the following conditions are true:
1. The water is clear enough that the benthic features of interest are visible in the imagery. 
2. There is a matching bathymetry dataset (with similar resolution and coverage to the satellite imagery and it is not derived from satellite imagery itself).
3. There are enough relatively flat light and dark areas that they can be used to fit the relationship between bathymetry, depth and reflectance.

This dataset contains case studies for three reef complexes in the Coral Sea (Flinders Reefs, Holmes Reefs and Lihou Reefs) and a small part of the GBR north east of Townsville including Davies Reefs, Grub Reef, Chicken Reef, and Bowl Reef. It also contains additional tests to determine how sensitive the results are to degraded input data, including using a bathymetry dataset that is 1/10th the resolution (100 m) of the satellite imagery (10 m) and the effect of not performing sun glint correction prior to the benthic reflection estimation. 

This method was developed to assist in the mapping of oceanic vegetation on seafloor of coral atolls in the Coral Sea. For most of these atolls the available bathymetry is too low resolution and so we need to rely on manual visual mapping. This dataset serves as a visual reference where the vegetation can be mapped with greater confidence due to the estimated benthic reflectance.

In this dataset we perform the same processing for each of the first four colour bands of Sentinel 2 imagery (B1 UV - 443 nm, B2 Blue - 490 nm, B3 Green - 560 nm, B4 Red - 665 nm). This is done to assess the relative effectiveness of each band and to what depth they can be used for mapping benthic reflectance.

Methods:
The benthic reflectance is estimated from satellite imagery and bathymetry. The satellite image brightness is scaled pixel by pixel between two thresholds corresponding to the expected brightness for low and high benthic reflectance. These thresholds are estimated from the bathymetry and a model of the depth verses brightness and benthic reflectance. This model is parameterised in each area using sampling points, chosen across the depth gradient, that have been classified as high or low benthic reflectance. The reflected light is determined by the attenuation coefficient, which is the fraction of light lost in each metre of water, the amount of scattered reflected light from deep water, and the relative brightness of the wet substrate with no water cover. This is based on the model developed by Jupp 1988.

To estimate the benthic reflectance we:
1. Start with an cloud free, clear water, low noise Sentinel 2 image composite of a region and matching high resolution Digital Elevation Model (DEM) for the same region. In our case we test the approach in the Coral Sea and the GBR using Sentinel 2 image composites prepared and described in the Marine satellite image test collections (AIMS) (Lawrey and Hammerton, 2024) and the GBR 30 m and 100 m bathymetry (Beaman 2017, Beaman 2020) data sets. The GBR 30 m bathymetry covers the GBR and some of the atolls close to the GBR, including North Flinders and Holmes reefs that we use in this study. The GBR 100 m, provides coverage of the whole Coral Sea, but at a significantly lower resolution, and quality (greater levels of interpolation).
2. We manually select and record locations in the image that visually correspond to high and low benthic reflectance. In deep areas the contrast between areas that are high benthic reflectance (sand) and low benthic reflectance (vegetation or coral) is low. In these areas it is often ambiguous whether the area is dark due to depth or substrate reflectance. We therefore select sampling points where the ambiguity is low. This occurs where there is a transition between dark vegetation and sand in an area that is likely to be flat or gently sloping.
3. Each Sentinel 2 swath is imaged by 12 detectors that are staggered and overlap (MSI Overview, n.d.). In clear water areas there is a noticeable brightness difference between successive detectors in the image. Some combination of the slight differences in the sensors and parallax angle between odd and even detectors result in light and dark banding in the Sentinel 2 imagery. With the large amount of contrast stretching that is needed to estimate the benthic reflectance in deep waters, these small brightness differences can result in very large errors in the final benthic reflectance. To compensate for this we divide the Sentinel 2 imagery into sections corresponding to each of the staggered detectors in the Sentinel 2 MSI instrument, rather than whole Sentinel 2 image scenes. The depth verses image brightness modelling is performed independently on each detector.
4. We extract from the satellite imagery and the DEM triplets of bathymetry, brightness and reflectance classification. 
5. We then visualise and review the depth verses brightness curves to look for outliers, then review the underlying imagery and bathymetry data for the potential reason for large deviation. To correct the outliers we would typically move the point a locations less affected the cause of the error. In most cases the errors were due to the limited resolution of the bathymetry and its errors in very shallow areas (&lt; 1-2 m deep). We also reviewed the number of sample points in each band of depths, adding new points to improve the coverage of all depths.
6. We then fit a model for each reflectance level, and detector segment, mapping the relationship between depth and image brightness. We use a simple model, based on Jupp 1988, that assumes that the brightness of the reflected signal is the addition of scattered light plus the reflected light that exponentially decays with depth. We parameterise this model for each detector swath area based on the data established in step 4 and least squares fitting, using scipy.optimize.curve_fit. This model estimates the depth averaged attenuation coefficient for the downwelling and upwelling light for each Sentinel 2 band and the background scattered light, which largely matches the brightness of open ocean water. This model assumes that the attenuation coefficient and background scattered light is constant over the detector segment. While this model doesn't fully parameterise all the inherent optical properties of the water it does provide a very good fit under most conditions. 
7. We then synthesise an estimate for the upper (corresponding to the high reflectance model) and lower (corresponding to the low reflectance model) brightness expected for each location based on the bathymetry. 
8. The brightness of the original satellite imagery is then scaled between these limits, estimating the reflectance for each pixel. In deep areas the contrast is greatly enhanced scaling the high and low reflectance areas to match the contrast of shallow areas. Small errors in the model (offsets in estimates of the high and low reflectance) get magnified by the amount of contrast enhancement. To limit these errors we constrain the maximum contrast enhancement to that needed to normalise to a depth of approximately 50 m. This threshold was determined experimentally. 
9. To reduce the amount of noise in the image we apply a small gaussian filter, with a filter sigma radius of 15 m (1.5 pixels in the final image).

Limitations of the data:

This approach only works in areas where the water is clear enough to see the benthic features and there is an independent source of high quality bathymetry. The bathymetry needs to be close in resolution to the imagery, otherwise it introduces significant errors in the conversion of satellite imagery to reflectance. In the tests performed in this study we achieve very good results using bathymetry that was three times lower resolution (30 m) than the satellite imagery (10 m). Tests using bathymetry at 1/10th the satellite imagery showed significant problems. Each area being mapped requires manually selecting points of high and low reflectance to perform the conversion. In some regions there may be insufficient high and low reflectance areas at each depth level to create the curves needed to fit the models. 

If the water constituents are consistent across the imagery then the scattering from suspended sediment, raising the brightness, or increases in the light absorption from coloured dissolved organic matter (CDOM), lowering the brightness, will be compensated for by the empirical sampling and modelling of the depth verses brightness for the image. If however, there is a variation in the water across the scene, such as a turbidity gradient or plumes of high CDOM coming off reefs and marine vegetation, then these brighter or darker regions will be misinterpreted by the algorithm as changes in the benthic reflectance rather than changes in water conditions. This means that any areas with a higher concentration of CDOM will be darker and falsely interpreted as having a lower benthic reflectance. Any areas that are affected by suspended sediment will be brighter and falsely interpreted by the algorithm as a higher benthic reflectance.  

The conversion from the satellite imagery to benthic reflectance is calibrated by manual sampling of locations across the image. These locations are classified as high or low benthic reflectance based on visual inspection of the imagery. In practice no parts of the images correspond to pure white (benthic reflectance of 1) or pure black (benthic reflectance of 0). The points labelled as high benthic reflectance typically correspond to sandy areas and low benthic reflectance areas correspond to algae or reef areas. These areas correspond to intermediate benthic reflectance values. In this study we assume that the high reflectance sampling locations correspond to a benthic reflectance of 0.8 and the low reflectance locations corresponding to a benthic reflectance of 0.4. These values are only approximate and were chosen to ensure the resulting benthic reflectance data had sufficient contrast to assist in vegetation mapping. The resulting maps are thus closer to a relative estimate of benthic reflectance than a calibrated estimate of benthic reflectance from 0 to 1.

This dataset only contains a limited study area in the Coral Sea (Flinders Reefs, Holmes Reefs, Lihou Reef) and a small part of the GBR north east of Townsville including Davies Reefs, Grub Reef, Chicken Reef, and Bowl Reef.

This method is appropriate for use in select locations where there is both good bathymetry and clear water. It is most useful for studying coral atolls.

In this method we cut up the Sentinel 2 imagery into narrow swaths corresponding to each of the staggered detectors of the MSI instrument. We independently model the depth verses image brightness and reflectance for each of these detector swaths to compensate for the slight brightness differences in each detector swath. This however requires that there are enough useful benthic sampling sites (high and low benthic reflection) in each modelled area. 

Output file descriptions:

This data repository contains the files used as part of the analysis including:
- new-data/Depth-Reflectance-Sampling-Points.shp: This file contains 1678 manually positioned point locations across the study areas, classified as high or low reflectance. These points are used to characterize the relationship between depth, image brightness, and benthic reflectance. The points were initially placed in areas with high confidence in benthic reflectance, such as edges between sandy and vegetated areas or around grazing halos. Outliers in the brightness versus depth curves were identified and adjusted based on classification mistakes or bathymetry issues, with points moved to flatter locations when necessary.
- new-data/Swath-analysis-areas-Poly.shp: This file contains areas where the brightness versus depth relationship is modeled independently. The Sentinel 2 imagery is divided into areas matching the width of each detector to account for their slightly different brightness offset characteristics, which are magnified by the benthic estimation processing. The depth versus brightness and benthic reflectance modeling is performed separately for each detector.
- new-data/metadata/Benth-Reflect_dataset-bounds.shp: This contains the boundary of the dataset and is used for the creation of the dataset metadata record. It represents the extent of the study area.

The following correspond to analysis results for 7 test case studies.
- output/55KFA-8: North Flinders reef composite of 8 images and GBR 2020 30 m bathymetry - Reference case with clearest imagery and good bathymetry.
- output/55KFA-1: North Flinders reef with single best image - To show much much difference using an image composite makes to the result.
- output/55KFA-8-gbr100: North Flinders reef with lower resolution 100 m bathymetry - To show the effect of lower resolution bathymetry.
- output/55KFA-8-NoSGC: North Flinders reef with an 8 image composite, but without sun glint correction - To show the benefit of sun glint correction.
- output/55KEB: Holmes Reefs GBR 2020 30 m bathymetry - For comparison with JCU drop camera benthic surveys.
- output/56KLF-7-gbr100: Lihou Reefs - To see the effectiveness when the bathymetry is limited.
- output/55KEV: GBR north east of Townsville including Davies Reefs, Grub Reef, Chicken Reef, and Bowl Reef - To assess the performance in waters with lower water clarity and high suspended sediment than the Coral Sea.

Data dictionary:

The following files are available for each of the test case studies.
output/*/02B-Depth_Reflect-class_S2-Bright.csv
- Latitude: Location of point sample
- Longitude: Location of point sample
- ID: Sequential counter of point sample
- Reflect: Substrate brightness, either 'High' or 'Low'
- SWATH_SEG: Integer 1 or 2, corresponding to two separate areas to repeat the modelling over.
- Depth_m: Bathymetry of the point sample in metres
- S2_R1_B1: Sentinel 2 image brightness band 1 (UV)
- S2_R1_B2: Sentinel 2 image brightness band 2 (Blue)
- S2_R1_B3: Sentinel 2 image brightness band 3 (Green)
- S2_R1_B4: Sentinel 2 image brightness band 4 (Red)

output/03C-benthic-reflect_SEG_{swath}.tif
- Estimated benthic reflectance scaled from 1 - 255. 0 is reserved for no data.

References:
Jupp, D. L. B., 1988. Background and extensions to Depth of Penetration (DOP)  mapping in shallow coastal waters. Symposium on Remote Sensing of the Coastal Zone. Gold Coast, Queensland, Session 4, Paper 2

MultiSpectral Instrument (MSI) Overview. (n.d.). Sentinel Online. Retrieved March 8, 2024, from https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/msi-instrument

eAtlas Processing:
No modifications were made to the data as part of publication.

Location of the data:
This dataset is filed in the eAtlas enduring data repository at: data\custodian\2022-2024-NESP-MaC-2\2.3_Improved-Aus-Marine-Park-knowledge\CS_NESP-MaC-2-3_AIMS_Benth-Reflect

Tags: BENTHIC HABITAT, Coral Sea, marine
Organization: Australian Ocean Data Network