Coral Sea features satellite imagery and raw depth contours (Sentinel 2 and Landsat 8) 2015 – 2021 (AIMS)

This dataset contains Sentinel 2 and Landsat 8 cloud free composite satellite images of the Coral Sea reef areas and some parts of the Great Barrier Reef. It also contains raw depth contours derived from the satellite imagery. This dataset was developed as the base information for mapping the boundaries of reefs and coral cays in the Coral Sea. It is likely that the satellite imagery is useful for numerous other applications. The full source code is available and can be used to apply these techniques to other locations.

This dataset contains two sets of raw satellite derived bathymetry polygons for 5 m, 10 m and 20 m depths based on both the Landsat 8 and Sentinel 2 imagery. These are intended to be post-processed using clipping and manual clean up to provide an estimate of the top structure of reefs. This dataset also contains select scenes on the Great Barrier Reef and Shark bay in Western Australia that were used to calibrate the depth contours. Areas in the GBR were compared with the GA GBR30 2020 (Beaman, 2017) bathymetry dataset and the imagery in Shark bay was used to tune and verify the Satellite Derived Bathymetry algorithm in the handling of dark substrates such as by seagrass meadows. This dataset also contains a couple of small Sentinel 3 images that were used to check the presence of reefs in the Coral Sea outside the bounds of the Sentinel 2 and Landsat 8 imagery.

The Sentinel 2 and Landsat 8 imagery was prepared using the Google Earth Engine, followed by post processing in Python and GDAL. The processing code is available on GitHub (https://github.com/eatlas/CS_AIMS_Coral-Sea-Features_Img).

This collection contains composite imagery for Sentinel 2 tiles (59 in Coral Sea, 8 in GBR) and Landsat 8 tiles (12 in Coral Sea, 4 in GBR and 1 in WA). For each Sentinel tile there are 3 different colour and contrast enhancement styles intended to highlight different features. These include:
- `TrueColour` - Bands: B2 (blue), B3 (green), B4 (red): True colour imagery. This is useful to identifying shallow features are and in mapping the vegetation on cays.
- `DeepFalse` - Bands: B1 (ultraviolet), B2 (blue), B3 (green): False colour image that shows deep marine features to 50 - 60 m depth. This imagery exploits the clear waters of the Coral Sea to allow the ultraviolet band to provide a much deeper view of coral reefs than is typically achievable with true colour imagery. This imagery has a high level of contrast enhancement applied to the imagery and so it appears more noisy (in particular showing artefact from clouds) than the TrueColour styling.
- `Shallow` - Bands: B5 (red edge), B8 (Near Infrared) , B11 (Short Wave infrared): This false colour imagery focuses on identifying very shallow and dry regions in the imagery. It exploits the property that the longer wavelength bands progressively penetrate the water less. B5 penetrates the water approximately 3 - 5 m, B8 approximately 0.5 m and B11 < 0.1 m. Features less than a couple of metres appear dark blue, dry areas are white. This imagery is intended to help identify coral cay boundaries.

For Landsat 8 imagery only the `TrueColour` and `DeepFalse` stylings were rendered.

All Sentinel 2 and Landsat 8 imagery has Satellite Derived Bathymetry (SDB) depth contours.
- `Depth5m` - This corresponds to an estimate of the area above 5 m depth (Mean Sea Level).
- `Depth10m` - This corresponds to an estimate of the area above 10 m depth (Mean Sea Level).
- `Depth20m` - This corresponds to an estimate of the area above 20 m depth (Mean Sea Level).

For most Sentinel and some Landsat tiles there are two versions of the DeepFalse imagery based on different collections (dates). The R1 imagery are composites made up from the best available imagery while the R2 imagery uses the next best set of imagery. This splitting of the imagery is to allow two composites to be created from the pool of available imagery. This allows any mapped features to be checked against two images. Typically the R2 imagery will have more artefacts from clouds. In one Sentinel 2 tile a third image was created to help with mapping the reef platform boundary.

The satellite imagery was processed in tiles (approximately 100 x 100 km for Sentinel 2 and 200 x 200 km for Landsat 8) to keep each final image small enough to manage. These tiles were not merged into a single mosaic as it allowed better individual image contrast enhancement when mapping deep features. The dataset only covers the portion of the Coral Sea where there are shallow coral reefs and where their might have been potential new reef platforms indicated by existing bathymetry datasets and the AHO Marine Charts. The extent of the imagery was limited by those available through the Google Earth Engine.

# Methods:

The Sentinel 2 imagery was created using the Google Earth Engine. The core algorithm was:
1. For each Sentinel 2 tile, images from 2015 – 2021 were reviewed manually after first filtering to remove cloudy scenes. The allowable cloud cover was adjusted so that at least the 50 least cloud free images were reviewed. The typical cloud cover threshold was 1%. Where very few images were available the cloud cover filter threshold was raised to 100% and all images were reviewed. The Google Earth Engine image IDs of the best images were recorded, along with notes to help sort the images based on those with the clearest water, lowest waves, lowest cloud, and lowest sun glint. Images where there were no or few clouds over the known coral reefs were preferred. No consideration of tides was used in the image selection process. The collection of usable images were grouped into two sets that would be combined together into composite images. The best were added to the R1 composite, and the next best images into the R2 composite. Consideration was made as to whether each image would improve the resultant composite or make it worse. Adding clear images to the collection reduces the visual noise in the image allowing deeper features to be observed. Adding images with clouds introduces small artefacts to the images, which are magnified due to the high contrast stretching applied to the imagery. Where there were few images all available imagery was typically used.
2. Sunglint was removed from the imagery using estimates of the sunglint using two of the infrared bands (described in detail in the section on Sun glint removal and atmospheric correction).
3. A composite image was created from the best images by taking the statistical median of the stack of images selected in the previous stage, after masking out clouds and their shadows (described in detail later).
4. The brightness of the composite image was normalised so that all tiles would have a similar average brightness for deep water areas. This correction was applied to allow more consistent contrast enhancement. Note: this brightness adjustment was applied as a single offset across all pixels in the tile and so this does not correct for finer spatial brightness variations.
5. The contrast of the images was enhanced to create a series of products for different uses. The `TrueColour` colour image retained the full range of tones visible, so that bright sand cays still retain detail. The `DeepFalse` style was optimised to see features at depth and the `Shallow` style provides access to far red and infrared bands for assessing shallow features, such as cays and island.
6. The various contrast enhanced composite images were exported from Google Earth Engine and optimised using Python and GDAL. This optimisation added internal tiling and overviews to the imagery. The depth polygons from each tile were merged into shapefiles covering the whole for each depth.

## Cloud Masking

Prior to combining the best images each image was processed to mask out clouds and their shadows.

The cloud masking uses the COPERNICUS/S2_CLOUD_PROBABILITY dataset developed by SentinelHub (Google, n.d.; Zupanc, 2017). The mask includes the cloud areas, plus a mask to remove cloud shadows. The cloud shadows were estimated by projecting the cloud mask in the direction opposite the angle to the sun. The shadow distance was estimated in two parts.

A low cloud mask was created based on the assumption that small clouds have a small shadow distance. These were detected using a 40% cloud probability threshold. These were projected over 400 m, followed by a 150 m buffer to expand the final mask.

A high cloud mask was created to cover longer shadows created by taller, larger clouds. These clouds were detected based on an 80% cloud probability threshold, followed by an erosion and dilation of 300 m to remove small clouds. These were then projected over a 1.5 km distance followed by a 300 m buffer.

The buffering was applied as the cloud masking would often miss significant portions of the edges of clouds and their shadows. The buffering allowed a higher percentage of the cloud to be excluded, whilst retaining as much of the original imagery as possible.

The parameters for the cloud masking (probability threshold, projection distance and buffer radius) were determined through trial and error on a small number of scenes. The algorithm used is significantly better than the default Sentinel 2 cloud masking and slightly better than the COPERNICUS/S2_CLOUD_PROBABILITY cloud mask because it masks out shadows, however there is potentially significant improvements that could be made to the method in the future.

Erosion, dilation and buffer operations were performed at a lower image resolution than the native satellite image resolution to improve the computational speed. The resolution of these operations were adjusted so that they were performed with approximately a 4 pixel resolution during these operations. This made the cloud mask significantly more spatially coarse than the 10 m Sentinel imagery. This resolution was chosen as a trade-off between the coarseness of the mask verse the processing time for these operations. With 4-pixel filter resolutions these operations were still using over 90% of the total processing resulting in each image taking approximately 10 min to compute on the Google Earth Engine.

One problem with the cloud masking was that very bright sand cays were persistently interpreted as clouds resulting in transparent holes in the final composite images. To counter this the cloud masked composite was overlaid on a composite made with no cloud masking. This would infill the holes created by false positive masking. Since the imagery was selected to have minimal clouds over reef and cay areas the under lay introduced little or no clouds into the final composites.

## Sun glint removal and atmospheric correction.

Sun glint was removed from the images using the infrared B8 band to estimate the reflection off the water from the sun glint. B8 penetrates water less than 0.5 m and so in water areas it only detects reflections off the surface of the water. The sun glint detected by B8 correlates very highly with the sun glint experienced by the ultra violet and visible channels (B1, B2, B3 and B4) and so the sun glint in these channels can be removed by subtracting B8 from these channels.

This simple sun glint correction fails in very shallow and land areas. On land areas B8 is very bright and thus subtracting it from the other channels results in black land. In shallow areas (< 0.5 m) the B8 channel detects the substrate, resulting in too much sun glint correction. To resolve these issues the sun glint correction was adjusted by transitioning to B11 for very shallow areas. B11 penetrates the water much less than B8. We don't use B11 everywhere because it is half the resolution of B8.

Land areas need their tonal levels to be adjusted to match the water areas after sun glint correction. Ideally this would be achieved using an atmospheric correction that compensates for the contrast loss due to haze in the atmosphere. Complex models for atmospheric correction involve considering the elevation of the surface (higher areas have less atmosphere to pass through) and the weather conditions. Since this dataset is focused on coral reef areas, elevation compensation is unnecessary due to the very low and flat land features being imaged. Additionally the focus of the dataset it on marine features and so only a basic atmospheric correction is needed. Land areas (as determined by very bright B8 areas) were assigned a fixed smaller correction factor to approximate atmospheric correction. This fixed atmospheric correction was determined iteratively so that land areas matched the tonal value of shallow and water areas.

## Depth contours

The depth contours in this dataset were derived from the composite imagery. Satellite Derive Bathymetry (SDB) requires reference bathymetry to tune transformation from brightness in the imagery to its associated depth. No high resolution shallow bathymetry was available in the Coral Sea and so the SDB algorithms were calibrated against reefs in the outer Great Barrier Reef where reference bathymetry, GA GBR30 2020 (Beaman, 2017) was available.

In this dataset we focused on generating specific depth contours rather than a fully linear bathymetry dataset.

The 5 m and 10 m depth contours were estimated using the following algorithm: (ln(Green)/ln(Blue-BlackLevel))*UnitScalar-UnitOffset, where Green corresponds to the green channel (B3) and Blue corresponds to the blue channel (B2). The BlackLevel adjusts the relative contrast and colour balance between the blue and green channel and was tuned to provide the best compensation for dark substrates such as seagrass at 8 - 10 m. This tuning was performed using the imagery from Shark Bay, as this region has shallow clear waters with a mix of sand and seagrass. The UnitScalar and UnitOffset values transform the log ratio into units of metres. These were tuned by starting with a rough estimate then successively tuning the values so that the 5 m and 10 m contours would best match the GBR30 dataset on the Great Barrier Reef, and the results obtained with Landsat 8 and Sentinel 2 would align in the Coral Sea. To further calibrate and validate these contours they were compared to the deepest features visible in the red (B4) and far red (B5) channels of Sentinel 2. These channels are limited in their depth penetration, ~5 m for B5 and ~6-8 m for B4. This depth limit is due to the absorption of light in the red and far red part of the spectrum and thus they represent an excellent independent reference depth.

The SDB algorithm performs poorly for depths below 12-15 m and so could not be used for creating the 20 m depth contour product. Instead the 20 m contour was generated by applying a brightness threshold to the green channel (B3). The threshold was calibrated against multiple Sentinel 2 scenes on the GBR and the GBR30 dataset. This technique relies on the fact that the green only penetrates to a depth of 28 - 30 m in clear water and so features deeper than this are guaranteed to not create false positives. At depths close to the limit (such as 20 m) the contrast is reduced leading to less susceptibility to variations in substrate brightness. See change log 2023-02-06.

Spatial filtering was applied to reduce the noise in all depth contours products. This means that many of the small reef features (< ~50 m) are not represented in the contours.

The resulting contours were: -5 m MSL +- 1.5 m, -10 m +- 2 m MSL and -20 m MSL +- 6 m.

Additional details about the depth calibration can be found in the Google Earth Engine source code on GitHub.

# Format:

- Sentinel 2 - GeoTiff - LZW compressed, 8 bit channels, 0 as NoData, Imagery as values 1 - 255. Internal tiling and overviews. Average size: 11500 x 11500 pixels and 300 MB per image.
- Landsat 8 - GeoTiff - LZW compressed, 8 bit channels, 0 as NoData, Imagery as values 1 - 255. Internal tiling and overviews. Average size: 7750 x 7780 pixels and 35 - 60 MB per image.
- Raw 5 m, 10 m, 20 m Depth Contour polygons - Shapefiles - Coral Sea 94 MB, GBR 60 MB
- Preview imagery - JPEG - 50% resolution versions of the Sentinel 2 and Landsat 8 images (1 - 5 MB each)

## Filename conventions

The images in this dataset are all named using a naming convention. An example file name is `CS_AIMS_Coral-Sea-Features_Imagery_S2_R1_DeepFalse_54LZP.tif`. The name is made up from:
- Dataset name (`CS_AIMS_Coral-Sea-Features_Imagery`). This is the imagery used in the digitisation of the `CS_AIMS_Coral-Sea-Features` dataset.
- Satellite source: `L8` for Landsat 8 or `S2` for Sentinel 2.
- Best imagery (`R1`) or second reference imagery (`R2`), `R3` for third reference image,
- Colour and contrast enhancement applied (`DeepFalse`, `TrueColour`,`Shallow`,`Depth5m`,`Depth10m`,`Depth20m`),
- Image tile (example: Sentinel 2 `54LZP`, Landsat 8 `091086`)

# Limitations:

In some scenes (56KQB and 56KPC) very few Sentinel 2 images were available leading to poor image composites. In these cases Landsat 8 images were developed and used for reef boundary mapping.

Due to the high contrast enhancement applied in the `DeepFalse` composites masked out clouds create significant visual artefacts that need to be considered when mapping from this imagery. Additionally slight differences in the sensitivity and angle of the Sentinel 2 MSI imager sensors results in uncorrected diagonal tonal bands in the imagery. The `Slope` composites show the cloud masks a rings in the images, and care needs to be made not to interpret these as reef features.

Only simple atmospheric correction was applied to land areas and as a result the imagery only approximates the bottom of atmosphere reflectance.

For the sentinel 2 imagery the sun glint correction algorithm transitions between different correction levels from deep water (B8) to shallow water (B11) and a fixed atmospheric correction for land (bright B8 areas). Slight errors in the tuning of these transitions can result in unnatural tonal steps in the transitions between these areas.

For the Landsat 8 image the sun glint correction did not consider land areas and so these appear as black in the imagery.

The depth contours are estimated using satellite derived bathymetry that is subject to errors caused by cloud artefacts, substrate darkness, calibration issues and uncorrected tides. The depth contours in this dataset are RAW and contain many false positives due to clouds. They should not be used without additional dataset cleanup. The intention is that areas outside the reef boundaries should be clipped (to remove cloud artefacts) and where a reef is disrupted by clouds, such as Kenn Reef in the Sentinel 2 imagery, then the Landsat 8 depth contours can be used instead.

The raw depth contours were generated in both the GBR and Coral Sea region, however depth contours in the GBR are only valid in areas where the water clarity is high. The depth contours of all inshore areas will be WRONG due to lower water clarity.

The depth contours products should NOT be used for navigation. The algorithm smooths out small features to lower noise in the contours, resulting in small reefs (< 50 m) being removed. As a result there will be many shipping hazards not detected in these contours. Additionally the depth contours are affected by substrate reflectivity and so darker substrates need to be shallower to be included in the contour.

The depth contour products are very slow to render as the entire regions (GBR and Coral Sea) are represented as a single multi-polygon. Speed issues could be resolved by subsequent post processing of the dataset split it into multiple polygons.

# Change log

2022-11-12 Depth contour products for 5 m, 10 m and 20 m, were refined and better calibrated. Imagery and depth contours from the GBR used for the depth calibration were added to the dataset, along with the imagery for Shark Bay. These were added for completeness.

2023-02-06 Comparison between the 20 m contour and the GBR bathymetry 30 m 2020 dataset on North and South Flinders island show that the contour corresponds to deeper waters over clear sandy areas and shallower over reef substrate. These are due to differences in the reflectivity of the substrate. In sandy areas on North and South Flinders reefs the contour closely matches 23 m depth whilst over reef substrate the contour corresponds to 17 m depth contour. On Saumarez Reefs the sandy areas match 23 m depth, whilst reef substrate corresponds to 14 m depth.

# References:

Beaman, R.J. 2017. High-resolution depth model for the Great Barrier Reef - 30 m. Geoscience Australia, Canberra. http://dx.doi.org/10.4225/25/5a207b36022d2

Google (n.d.) Sentinel-2: Cloud Probability. Earth Engine Data Catalog. Accessed 10 April 2021 from https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S...

Zupanc, A., (2017) Improving Cloud Detection with Machine Learning. Medium. Accessed 10 April 2021 from https://medium.com/sentinel-hub/improving-cloud-detection-with-machine-l...

International Hydrographic Organization, Intergovernmental Oceanographic Commission, The IHO-IOC GEBCO Cook Book, IHO Publication B-11, Monaco, Oct. 2019, 493pp - IOC Manuals and Guides 63, France, Oct. 2019, 493pp

# Data Location:

This dataset is filed in the eAtlas enduring data repository at: data\2020-2029-AIMS\CS_AIMS_Coral-Sea-Features_Img
The source code is available on [GitHub](https://github.com/eatlas/CS_AIMS_Coral-Sea-Features_Imagery).