Coral Reef Health Monitoring around Vanuatu’s Coastline

This notebook provides a step-by-step guide for identifying coral reefs and assessing their health over time using Sentinel-2 satellite imagery. Coral ecosystems are highly sensitive to environmental change, and satellite data offers a powerful way to monitor reefs at scale.

We will walk through a remote sensing workflow, including:

1. Obtain Top of Atmosphere Sentinel-2 imagery: We will use provinces to get buffered coastlines as geographic selection criteria for a STAC search.

2. Obtain derived Sea Surface Temperature (SST) data from Sentinel-3: Use the same seach criteria to obtain SST data and interpolate it to the grid of the Sentinel-2 data.

3. Cloud and land masking: Remove clouds, shadows, and non-water pixels to isolate the marine environment.

4. Atmospheric correction: Convert raw Level-1C Sentinel-2 imagery to surface reflectance (Level-2A) using physics-based correction techniques.

5. Index calculation: Derive the Normalized Blue-Green Index (NBGI), which is sensitive to benthic features like live coral, algae, and sand.

6. Training a classifier: Use labeled coral reef polygons from a vector dataset to train a Random Forest classifier on the spectral and index data.

7. Time series analysis: Apply the trained model across the time series to generate predictions, and track trends in reef extent over time.

Define AOI — Vanuatu Coastal Buffer¶

Buffer one province for demonstration. We’ll use a subset of a province for this lesson, as we’ll be using a time series which consumes more memory than a single timestamp.

Search STAC for Sentinel-2 L1C Imagery¶

Load Imagery¶

Search STAC for derived Sea Surface Temperature data¶

Sea surface temperature is one of the most important indicators for coral reef health because corals are highly sensitive to even small changes in temperature. Corals have a narrow thermal tolerance, and even small deviations above normal can trigger physiological stress. Prolonged heat stress—typically when SST exceeds the long-term maximum summer average by 1 °C or more for several weeks—can cause mass coral bleaching, a process in which corals expel their symbiotic algae, leading to reduced growth, increased susceptibility to disease, and potentially mortality (NOAA Coral Reef Watch, n.d.; Hoegh-Guldberg, 1999). Elevated SST is widely recognized as the primary driver of recent global bleaching events, fundamentally altering reef community structure (Hughes et al., 2018). Consequently, consistent SST monitoring enables both early warning and targeted management responses to protect vulnerable reef systems (Maynard et al., 2008). In short, SST acts like a “fever thermometer” for reefs — it tells you when the system is overheating, can predict stress events, and provides a long-term health indicator.

References

• NOAA Coral Reef Watch. (n.d.). Coral Bleaching. Retrieved from https://coralreefwatch.noaa.gov/satellite/education/tutorial/bleaching.php

• Hoegh-Guldberg, O. (1999). Climate change, coral bleaching and the future of the world’s coral reefs. Marine Pollution Bulletin, 38(2), 39–44. https://doi.org/10.1016/S0025-326X(99)00002-0

• Hughes, T. P., Kerry, J. T., Álvarez-Noriega, M., et al. (2018). Global warming transforms coral reef assemblages. Nature, 556(7702), 492–496. Hughes et al. (2017)

• Maynard, J. A., Anthony, K. R., Marshall, P. A., & Masiri, I. (2008). Major bleaching events can be predicted weeks in advance. Estuarine, Coastal and Shelf Science, 79(1), 173–180. Walker et al. (2008)

We’ll start by obtaining the dates we need SST data for. We’ll get those from the L1C time series.

We also need a center lat/lon for our area of interest to filter for intersecting scenes.

With that, we’ll search for an intersecting SST image for each timestamp.

We should have 5 SST items returned, which would equal how many unique timestamps we have in our L1C time series.

Checkout the metadata for the first SST item.

Mask land and clouds in the Sentinel-2 data¶

Now, we’ll come back to the Sentinel-2 data to mask away land and clouds. The buffer extent subtracting the original island boundary polygon (the coastal offshore remainder) is what we will use to designate the valid region.

That above WKT describes the EPSG:32759 coordinate reference system, which is:

WGS 84 / UTM zone 59S

• Datum: WGS84

• Projection: Transverse Mercator

• Zone: 59 south (central meridian at 171°E)

• Units: meters (Easting/Northing)

• False easting: 500,000 m

• False northing: 10,000,000 m (because it’s in the southern hemisphere)

Our Sentinel-2 data is in EPSG:32759, which means the x and y coordinates are in meters east/north from the zone’s origin.

Now we will mask clouds. The Sentinel-2 L2A comes with an extra “band” called SCL which stands for Scene Classification Layer. The values in this band denote various high level land cover types, such as water and vegetation, but more importantly, there are also classes for cloud types. We will select those to use for masking away clouds.

Mask the L1C data, as that is what we will use for downstream processing.

Top of Atmosphere (TOA) vs. Surface Reflectance & Atmospheric Correction¶

TOA reflectance refers to the reflectance values measured by a satellite sensor before any correction for atmospheric effects. It represents the radiance reaching the sensor at the top of the atmosphere, and includes contributions from:

• Direct solar reflectance from the surface,

• Scattering by atmospheric molecules (Rayleigh scattering),

• Scattering by aerosols and haze,

• Reflection from clouds and surrounding areas.

Sentinel-2 Level-1C products are provided in TOA reflectance.

Surface Reflectance is the fraction of incoming solar radiation reflected by the Earth’s surface, as it would appear if the atmosphere were not present. It represents the true reflectance of land, water, or vegetation and is more reliable for quantitative analysis.

Sentinel-2 Level-2A products are atmospherically corrected to surface reflectance.

Sentinel-2 Level-2A data is produced using the Sen2Cor processor, which performs atmospheric correction to derive surface reflectance from top-of-atmosphere (TOA) reflectance. It was designed primarily for land surfaces — not water. We’ll come back to this in a moment.

Atmospheric correction is the process of removing atmospheric effects from TOA reflectance to estimate surface reflectance. It compensates for:

• Scattering (blue light by molecules, haze),

• Absorption (by gases like water vapor, ozone, CO₂),

• Aerosols (dust, smoke, sea salt).

There are two main types:

• Empirical methods: e.g., Dark Object Subtraction (DOS), based on scene content.

• Physics-based methods: e.g., Sen2Cor, 6S, or ACOLITE, which use radiative transfer models and ancillary data (e.g., aerosol optical thickness).

Water bodies have unique optical properties that make atmospheric correction far more complex:

1. Low Reflectance (Dark Target Problem)

• Water reflects very little sunlight, especially in the visible and near-infrared.

• This makes it hard to distinguish water signal from atmospheric scattering, leading to unstable or noisy surface reflectance values.

• Overcorrected water pixels may even appear as negative reflectance, which is physically meaningless.

2. High Sensitivity to Atmospheric Effects

• Even small amounts of aerosols, thin clouds, or Rayleigh scattering can dominate the signal received from water.

• Sen2Cor is tuned for land reflectance properties, and may overcorrect these subtle signals in aquatic environments.

3. Incorrect Assumptions

• Sen2Cor assumes a Lambertian (diffuse) surface, which is true for land but not for water, which reflects sunlight specularly (like a mirror).

• It also uses land-based visibility and elevation models that don’t apply well to open water.

4. No Water-Specific Tuning

• Sen2Cor doesn’t use water-leaving radiance models or bidirectional reflectance functions (BRDF) specific to aquatic systems.

• It lacks aerosol correction schemes optimized for open water, unlike ocean color processors (e.g., ACOLITE or C2RCC).

What To Use Instead Over Water?¶

For more reliable atmospheric correction over water, use:

• Dark Object Subtraction (DOS) – A simple, empirical method suitable for ocean scenes with few clouds and minimal aerosols.

• ACOLITE – Optimized for coastal and inland waters.

Today, we’ll use DOS as it doesn’t require us to download any external software. It isn’t the MOST accurate method we could use for this task, as we’ll explain next, but it does well enough to be suitable for this a lightweight and flexible approach.

Dark Object Subtraction¶

Dark Object Subtraction (DOS) is a simple and widely used atmospheric correction method that estimates and removes the effects of atmospheric scattering in satellite imagery. It is particularly useful for correcting Sentinel-2 Level-1C (Top-of-Atmosphere reflectance) imagery in oceanic or coastal environments.

The core idea behind DOS is that certain “dark objects” in the scene—such as deep, clear water or dense vegetation—should theoretically have near-zero reflectance in some bands (especially the blue and shortwave bands). Any non-zero signal observed in these dark areas is attributed to path radiance (light scattered by the atmosphere).

How It Works¶

1. Identify dark pixels: Find the minimum (or low-percentile, e.g., 1%) reflectance values over dark surfaces, usually in the blue band (Band 2) or coastal aerosol band (Band 1).

2. Estimate path radiance: Assume this value represents atmospheric scattering (haze, Rayleigh scattering).

3. Subtract offset: Subtract this value from all pixels in the scene for each band to estimate surface reflectance:

   ρ_surface = ρ_TOA - ρ_min

4. Clip negatives: Any resulting negative reflectance values are set to zero.

Application Over the Ocean¶

• Clear deep water serves as an ideal dark object.

• Works best in clear-sky ocean conditions where water reflectance is minimal.

• Helps reduce haze and Rayleigh scattering impacts in shorter wavelengths (e.g., blue, green).

• Less effective in turbid or super shallow waters where bottom reflectance or suspended particles are present.

• Simple, fast, and does not require external atmospheric data.

• Useful for preprocessing imagery in remote areas with limited ancillary data.

Limitations¶

• Assumes perfect dark targets exist in the scene.

• May overcorrect in bright or turbid water conditions.

• Ignores adjacency effects and variable atmospheric thickness.

• Does not correct for absorption effects (e.g., water vapor, aerosols).

xarray tries to infer the time coordinate values from the datasets being concatenated.

But since we created each slice by selecting individual timesteps (isel(time=t)), the individual datasets might lose the original time coordinate values or have them replaced by default integer indexes (0,1,2...).

So the concatenated dataset ds_corrected may end up with a generic integer index for time instead of the original timestamp values.

Visualize the L1C true color imagery before atmospheric correction¶

Visualize the L1C true color imagery after atmospheric correction¶

Add the Sea Surface Temperature data to the corrected L1C imagery¶

Now we match the SST data to the L1C data by latitude/longitude coordinates.

• ds_l2a['x'] is the Sentinel-2 longitude coordinate array.

• ds_l2a['y'] is the Sentinel-2 latitude coordinate array.

• meshgrid expands them into 2-D arrays (lon2d, lat2d) so there’s a lat/lon pair for every pixel in the Sentinel-2 scene.

• points=(data['longitude'].values, data['latitude'].values) defines the known SST measurement locations in lon/lat space.

• xi=(lon2d, lat2d) says: interpolate onto these exact lat/lon positions.

• It’s doing the interpolation in geographic space, not in raster grid space.

• Each Sentinel-2 pixel gets the SST value of the nearest SST observation point (because method='nearest'), based purely on the geographic coordinates.

Check the value range for SST in each timestamp.

Compute Reef Health Indicators¶

The Blue-Green Index (BGI) is a spectral index used primarily in aquatic and coastal remote sensing to help discriminate benthic habitats (like coral reefs, seagrass, sand) and assess water properties (Bannari et al., 2022). It leverages the difference between the blue and green spectral bands.

Formula¶

• Calculated using surface reflectance values.

Interpretation of BGI values¶

$D$ = difference

• Positive $D>0$: Blue reflectance > Green → More blue reflectance → indicates clear water, possibly healthy coral or bright sandy substrate.

• Negative $D<0$: Green > Blue → often shallow water, vegetation, or sediment that boosts green. May suggest turbid water, algal blooms, or seagrass.

• Zero $D=0$: equal blue/green reflectance.

How BGI helps with coral reefs¶

• Differentiates benthic substrates: Coral, seagrass, algae, sand, and bare substrate often have distinct blue and green reflectance signatures, so BGI can help separate these classes.

• Highlights shallow water features: Coral reefs usually occur in shallow, clear water where blue and green light penetrates well, making BGI sensitive to benthic composition.

• Tracks changes over time: By analyzing BGI time series, you can detect changes in reef cover, algal blooms, sedimentation, or coral bleaching events affecting reflectance.

Limitations and considerations¶

• Scale-sensitive: Sensitive to illumination, sensor gain, and shadows — not standardized across conditions. A difference of 0.1 at low reflectance (e.g. 0.15–0.05) is not the same “contrast” as 0.1 at high reflectance (0.80–0.70), but raw $D$ treats them identically.

• Usually masked to water pixels (e.g., using land/water masks or cloud masks).

• Water column effects: Water depth, turbidity, and dissolved materials affect blue and green light differently, potentially confounding BGI values.

• Atmospheric correction needed: Accurate surface reflectance is essential for reliable BGI values, especially over water.

• Not a direct health measure: BGI reflects substrate and water color but not coral health metrics like bleaching directly. It is an indirect indicator.

• Supplement with other indices and data: Combine BGI with other indices (e.g., NDVI for algae, bathymetry data, or hyperspectral data) for better reef monitoring.

Because the Blue Green index has such limitations, we can normalize it to make it robust to lighting changes, atmospheric variation, and sensor inconsistencies — making it more reliable for time-series analysis or cross-scene comparisons. Importantly, coral bleaching often leads to whitening, reducing pigment (chlorophyll) and increasing reflectance in the visible spectrum (especially blue and green). Chlorophyll strongly absorbs blue (~450 nm) and red (~670 nm) light for photosynthesis, so healthy coral with chlorophyll content shows low reflectance in blue and red bands. When chlorophyll content decreases (e.g., due to stress or bleaching), the absorption in the blue and red wavelengths lessens, so more light is reflected in those bands. Green reflectance is usually high already because chlorophyll reflects green (~550 nm) more than blue or red, but with reduced chlorophyll, the overall structure and pigments may change, sometimes leading to subtle changes in green reflectance as well.

The normalized blue green index uses a normalized difference of the blue and green bands to detect coral bleaching, which manifests as a spectral signature change in shallow waters. This is a form of a Normalized Difference Index (like NDVI), designed to highlight spectral changes in coral reefs:

• Healthy coral reflects more in the green spectrum and absorbs more blue light.

• Bleached coral loses pigmentation and structure, becoming brighter in blue, reducing the green–blue contrast.

Interpretation of NBGI values¶

• Positive ($0< \text{NBGI}\le1$): Blue dominates relative to total signal → likely clearer/deeper water or substrates with strong blue reflectance (e.g., bleached coral).

• Negative ($-1\le \text{NBGI}<0$): Green dominates → shallow water, benthic vegetation (seagrass, algae), or sediment.

• Zero ($\text{NBGI}=0$): equal contribution of blue and green.

• Higher index values → likely healthy coral (stronger green reflectance).

• Lower index values → potential bleached coral (higher blue reflectance).

Advantages¶

• Normalized: automatically adjusts for overall brightness—sun angle, water depth, sensor gain—so values are directly comparable across images and times.

• Bounded: easy to interpret thresholds (e.g., $\text{NBGI}<−0.2$ for dense vegetation, $\text{NBGI}>0.2$ for clear water).

• Contrast-enhancing: emphasizes relative difference rather than absolute magnitude.

Considerations:¶

• Only works in shallow, clear water — turbid water or shadows can affect accuracy.

• Requires masking clouds, deep water, and land first.

• You can couple it with red edge or NIR bands, or other indices like:

  • NDVI

  • Red-Edge Chlorophyll Index

  • Turbidity Index (Blue/SWIR)

Ground truth data integration and segmentation¶

We will use a reef map provided by Vanuatu Bureau of Statistics.

Let’s also download a subset from the Allen Coral Atlas for the same area.

The reef map is a vector polygon dataset. We are rasterizing it to use as extent labels for the MSI and SST pixels.

Live Coral vs. Algae/Sand¶

Live coral has unique reflectance in the red-edge and NIR bands. As well, a sudden rise in visible reflectance and drop in red-edge/NIR may indicate bleaching.

Useful Bands:

• Red Edge

• Near Infrared

• Short Wave Infrared (for benthic type separation)

Live coral is not just a bare rock or dead skeleton — it’s covered by zooxanthellae algae living in the coral tissue. These algae give coral its color and contribute to unique pigment absorption features:

• Strong chlorophyll absorption in the blue and red bands.

• Higher reflectance in the green (green gap between chlorophyll absorption peaks).

• Sharp reflectance changes in the red-edge region.

• Distinct absorption features in the SWIR due to water content in living tissue.

Dead coral, macroalgae, and sand each lack some of these features or have them in different magnitudes.

Red Edge (700–750 nm)

• The red-edge region captures the abrupt transition between strong chlorophyll absorption in the red and high reflectance in the NIR.

• Live coral’s pigment-rich tissue produces a less steep slope compared to terrestrial vegetation but still much steeper than bare sand.

• Algae also show a red-edge shift, but spectral shape and slope differ — useful for discrimination.

• Live coral → moderate red-edge slope, distinct curvature.

• Macroalgae → steeper slope, more vegetation-like.

• Sand → almost flat, no steep rise.

Near Infrared (NIR, ~800–900 nm underwater)

• In shallow water, NIR reflectance is very low because water absorbs NIR strongly.

• However, differences in how little NIR is reflected can still be meaningful.

  • Live coral: reflect slightly higher than sand (due to scattering from tissue and skeleton microstructure).

  • Sand: minimal NIR return, especially in deeper water.

  • Algae: slightly higher NIR than coral because of cell structure.

NIR underwater is weak, but the relative differences in shallow clear water still carry classification power.

Short Wave Infrared (SWIR, >1000 nm)

• SWIR is almost completely absorbed by water at even small depths, so it’s not for direct benthic reflectance in situ.

• However:

  • In above-water remote sensing (like airborne or satellite with sunglint correction), SWIR is used for atmospheric and sunglint correction.

  • Good for very shallow / above-water cases to separate moisture-rich from dry substrates.

  • It can also help separate benthic types in very shallow or exposed conditions (intertidal coral flats) because:

    • Water content in coral tissue absorbs differently than in macroalgae.

    • Sand has distinct SWIR slopes.

Summary:

• Red Edge + NIR: Distinguishes living photosynthetic tissue (live coral, algae) from non-living substrates (sand, dead coral). Helps separate algae and coral based on pigment density.

• SWIR (if applicable): Helps reject false positives from wet sand or turbidity artifacts.

Flatten the features and newly rasterized labels for use with a random forest classifier. We will use the first timestamp as features for model training.

Data Splitting¶

Now that we have the arrays flattened, we can split the datasets into training and testing partitions. We will reserve 80 percent of the data for training, and 20 percent for testing.

Ensure all labels are in each partition.

Now we will set up a small random forest classifider with 10 trees. We use a seed (random_state) to ensure reproducibility. Calling the .fit() method on the classifier will initiate training.

Once the classifier is finished training, we can use it to make predictions on our test dataset.

It’s important to know how well our classifier performs relative to the true labels (y_test). For this, we can calculate the accuracy metric to measure agreement between the true and predicted labels.

We can also produce a classification report to check the precision, recall and F1 scores for each class.

We can also plot a confusion matrix to explore per-class performance.

Save the model to file so that it can be loaded and reused again without needing to repeat training.

If we want to generate predictions for the entire dataset in order to plot a map of predicted reefs for the entire area of interest, we can do this using the test province dataset.

Vectorize and save the predicted reef map to a geojson file.

Time series analysis and change detection¶

Trend analysis¶

Now let’s compute the per-pixel linear trend over time from the time series prediction maps.

• time_num is a numeric time axis (e.g., [0, 1, 2, ...]) corresponding to each time step in the data.

• The inner function fit_linear_trend(y) takes a 1D time series (one pixel’s values across time) and:

  • Masks out NaNs with np.isfinite.

  • Fits a linear model (np.polyfit) only if there are at least 2 valid observations.

  • Returns the slope of the best-fit line, representing the rate of change over time.

  • Returns NaN if the fit fails or insufficient data exists.

• xr.apply_ufunc applies fit_linear_trend to every pixel in the dataset using xarray’s vectorized and parallelized infrastructure, preserving spatial dimensions and returning a new DataArray of trend values.

This function is useful for detecting temporal changes (e.g., coral bleaching trends, coral reef recession or expansion) on a per-pixel basis.

Threshold the change to detect significant reef change¶

Track loss in reflectance in blue/green bands or increase in brightness. Bleached coral appears whiter and brighter due to algae loss.