Blog GIS & Territories

In two previous articles, we explored how to convert standard Sentinel-2 images (10 m resolution) to a super-resolution 1 m version using the S2DR3 module, developed for coral research in Mauritius. This step provides an unprecedented level of detail from open data, revealing fine structures that are often invisible at 10 m: seagrass beds, reef flats, sandy areas, and subtle transitions between vegetation and reef.

Traditionally, spectral indices (NDVI, NDWI, MNDWI, BSI, NDBI, etc.) is done in QGIS, once the images have been downloaded and processed. However, it is often more efficient to generate these indices directly at the source, i.e., as soon as the super-resolution Sentinel-2 images are produced.

This approach has several advantages:

• Complete automation: the indices are calculated as soon as the dataset is created, without additional manual processing in QGIS.
• Spatial consistency: the indices use exactly the same radiometric corrections and super-resolution as the original bands.
• Time savings: no intermediate exports or manual configuration of raster algorithms.
• Preparation for analysis: the results can be used immediately in QGIS or any other GIS tool.

The script presented below automates the entire workflow:

1. downloading and super-resolution of the Sentinel-2 scene,
2. generation of an intelligent cloud mask (excluding wave-overlap areas),
3. direct calculation of the main spectral indices,
4. creation of a multiband file ready for analysis.

This integration greatly simplifies the transition from raw image to interpretation maps, while ensuring the reproducibility of scientific processing.

The calculated indices

NDVI — Normalized Difference Vegetation Index

• Purpose: to measure vegetation density and health.
• Formula:

.

.

• NIR = near-infrared band (B08)
• Red = red band (B04)
• Values: -1 to 1; >0.2 vegetation, <0.1 bare soil or water.
• Use: vegetation monitoring, biomass estimation, agriculture.

NDWI — Normalized Difference Water Index

• Purpose: detect the presence of water in the landscape.
• Formula:

.

.

• Green = green band (B03)
• NIR = near-infrared band (B08)
• Values: positive = water, negative = soil/vegetation.
• Use: mapping lakes, rivers, wetlands.

MNDWI — Modified NDWI

• Purpose: to better separate water from urban or built-up areas.
• Formula:

.

.

• SWIR1 = short wavelength band (B11)
• Values: positive = water; improves detection compared to NDWI.
• Use: monitoring water bodies in urban or coastal areas.

BSI — Bare Soil Index

• Purpose: to detect bare soil and denuded areas.
• Formula:

.

.

• Blue = blue band (B02)
• Red = B04, NIR = B08, SWIR1 = B11
• Values: the higher the BSI, the more open the surface.
• Use: monitoring bare soil, urban areas, mineral areas.

NDBI — Normalized Difference Built-up Index

• Purpose: to detect built-up or impervious areas.
• Formula:

.

.

• Values: positive = built-up, negative = vegetation or water.
• Use: urban mapping, monitoring of urban sprawl.

EVI — Enhanced Vegetation Index

Purpose: to improve the sensitivity of NDVI in areas with very dense vegetation and correct for atmospheric effects.

Formula:

.

.

• Blue = B02, Red = B04, NIR = B08
• Values: -1 to 1, similar to NDVI but more robust in dense areas.
• Usage: forests, tropical areas, monitoring vegetation health.

SAVI — Soil Adjusted Vegetation Index

• Purpose: to measure vegetation by correcting for the effect of soil (useful in sparse areas).
• Formula:

.

.

where:

• NIR = near-infrared band (Sentinel-2 B08)
• Red = red band (B04)
• L = soil correction factor, often 0.5
• Values: -1 to 1, like NDVI, but less sensitive to bare soil.
• Use: sparse or low vegetation (grasslands, savannas, coastal areas with low vegetation cover).

UI — Urban Index (or Urbanization Index)

• Purpose: to detect and quantify urban areas.
• Simple formula:

.

.

• SWIR = short wavelength band (Sentinel-2 B11 or B12)
• NIR = near infrared (B08)
• Values: positive for impervious or built-up areas, low/negative for vegetation or water.
• Use: monitoring urbanization, detecting buildings or artificial surfaces.

RDI — Redness Difference Index

• Purpose: to detect bare soil or areas rich in iron oxide (“red”) in the landscape.
• Formula:

.

.

• Red = red band (B04), Green = green band (B03)
• Values: the more “red” the soil, the higher the RDI.
• Use: monitoring erosion, bare soil, mineral or arid areas.

Complete summary table

Sentinel Index Palette

Each index is presented with:

• its name and acronym,
• target surface type (vegetation, water, soil, urban),
• typical values/scale,
• representative color for map visualization (e.g., NDVI → green, NDWI → blue, BSI → light brown, etc.).

Why use a cloud mask?

Clouds and their shadows pose problems in satellite imagery because:

1. They obscure surfaces: soil, vegetation, water, coral.
2. They distort indices (NDVI, NDWI, etc.) because the radiometric values do not correspond to the actual surface.
3. Areas covered by clouds can generate outliers: very high or very low values, skewing trend calculations.

Therefore, for index analysis and the detection of features on land or in water, it is important to mask these pixels.

Cloud mask for Sentinel-2 at 10 m

Sentinel-2 offers several options for detecting clouds at its native resolution (10 m, 20 m, or 60 m depending on the band):

a) Official mask provided by ESA

• QA60 file in Sentinel-2 L2A products: 10 m binary mask indicating clouds and cirrus.
• Advantages: reliable and easy to apply.
• Disadvantages: coarse mask → does not always differentiate well between thin clouds and shadows.

b) Methods based on spectral indices

• SWIR and NIR bands: clouds are highly reflective in the visible spectrum but less so in SWIR.
• Example of a simple index:

.

.

• Bright pixels in the visible spectrum and faint pixels in SWIR → probably clouds.
• Allows you to create a cloud mask adapted to the context, useful if you want to detect local clouds and cirrus clouds.

c) Automatic algorithms

• Fmask or s2cloudless: sophisticated detection using machine learning and multispectral spectra.
• Useful for large time series.

How we proceeded in the Colab script

In the script, we did not use the official QA60, but rather detection by spectral indices and local texture, adapted to coastal areas and 1 m super-resolution:

Main steps:

1. Band normalization :
   blue_n = blue / max
green_n = green / max
red_n = red / max
nir_n = nir / max
swir1_n = swir1 / max → all values between 0 and 1 for homogeneity.
2. “Initial cloud” condition based on 4 criteria :
   • TVery clear in the visible (albedo_vis > 0.35) High SWIR1 (swir1_n > 0.15) High blue/red ratio (blue_n / red_n > 1.2)High NIR (nir_n > 0.25)

   → pixels meeting at least 3/4 criteria are considered clouds..

3. “Foam” mask for marine areas :

   • We detect foam/breaking waves via the local standard deviation on the blue band and visible albedo.
   • These pixels are removed from the cloud mask :
     cloud_final = cloud_init & (~foam_mask)
   • This avoids masking coral reefs or rough seas..

4. Filtering small components: :

   • Removes very small clouds that could be noise → only objects > 500 pixels are kept.

5. Final result :

   • cloud_mask = cleaned cloud mask, ready to be used to exclude cloudy pixels when calculating indicess.

Why this method?

• The scene is super-resolved (1 m) → the official QA60 does not exist at this resolution.
• The area is coastal/coral → foam and breaking waves must be distinguished from clouds.
• The spectral + local texture approach allows us to:

  • detect clouds
  • avoid masking water or foam
  • work on super-resolved images generated by S2DR3

Proceduree

Open Google Colab: https://colab.research.google.com

In Collab, link your notebook to your Google Drive

Now, to ensure that processing does not take too long, you need to enable a GPU in the runtime environment. Look in the lower right corner of your notebook, where you will see something like “Python 3” with a chip icon. Click on it and select Change Runtime Environment Type. A window will open where you will need to search for Hardware Accelerator and select GPU T4. Save the changes and the notebook will automatically restart with the new configuration.

Copy this code into a new cell in your notebook and run it:

script Colab

# =========================================================
#   S2DR3 + Calcul indices Sentinel-2 + Masque nuages
#   Étude des coraux à Blue Bay / Île aux Aigrettes
# =========================================================

# --- 1. Monter Google Drive ---
from google.colab import drive
drive.mount('/content/drive')

# Crée un dossier pour les résultats
!mkdir -p /content/drive/MyDrive/Sentinel2_Coraux_S2DR3/output
!ln -s /content/drive/MyDrive/Sentinel2_Coraux_S2DR3/output /content/output

# --- 2. Installer le paquet S2DR3 ---
!pip -q install https://storage.googleapis.com/0x7ff601307fa5/s2dr3-20250905.1-cp312-cp312-linux_x86_64.whl

# --- 3. Importer le module principal ---
import s2dr3.inferutils
import os, glob, numpy as np, rasterio, scipy.ndimage as ndi
from skimage.transform import resize

# --- 4. Définir la zone d'intérêt et la date ---
lonlat = (57.73, -20.44)
date = '2025-10-09'

# --- 5. Lancer le traitement S2DR3 ---
s2dr3.inferutils.test(lonlat, date)

# =========================================================
#   6. Lecture du fichier _MS.tif
# =========================================================
root_dir = '/content/output'
files = sorted(glob.glob(os.path.join(root_dir, '**/*_MS.tif'), recursive=True),
               key=os.path.getmtime, reverse=True)
if not files:
    raise FileNotFoundError("❌ Aucun fichier *_MS.tif trouvé dans /content/output")
ms_path = files[0]
outdir = os.path.dirname(ms_path)
print(f" Fichier trouvé : {os.path.basename(ms_path)}")
print(f" Dossier : {outdir}")

# Lecture des bandes
with rasterio.open(ms_path) as src:
    profile = src.profile
    n_bands = src.count
    print(f"Nombre de bandes détectées : {n_bands}")
    bands = [src.read(i).astype('float32') for i in range(1, n_bands+1)]

def get_band(index, name):
    try:
        return bands[index-1]
    except IndexError:
        print(f" Bande {name} manquante, remplacée par zéro.")
        return np.zeros_like(bands[0])

# Attribution bandes principales Sentinel-2
blue  = get_band(1, "B02")
green = get_band(2, "B03")
red   = get_band(3, "B04")
nir   = get_band(4, "B08")
swir1 = get_band(9, "B11")
swir2 = get_band(10, "B12")

# =========================================================
#   7. Calcul des indices
# =========================================================
def safe_div(a, b, mask=None):
    res = np.where((b!=0) & (~np.isnan(a)) & (~np.isnan(b)), a/b, 0)
    if mask is not None:
        res[mask] = np.nan
    return res

L = 0.5  # facteur SAVI

indices = {
    'NDVI' : safe_div(nir - red, nir + red),
    'NDWI' : safe_div(green - nir, green + nir),
    'MNDWI': safe_div(green - swir1, green + swir1),
    'BSI'  : safe_div((swir1 + red) - (nir + blue), (swir1 + red) + (nir + blue)),
    'NDBI' : safe_div(swir1 - nir, swir1 + nir),
    'EVI'  : 2.5 * safe_div((nir - red), (nir + 6*red - 7.5*blue + 1)),
    'SAVI' : safe_div((nir - red) * (1 + L), (nir + red + L)),
    'UI'   : safe_div(swir1 - nir, swir1 + nir),
    'RDI'  : safe_div(red - green, red + green)
}

# =========================================================
#   8. Masque nuages léger adapté 1 m
# =========================================================
eps = 1e-8
blue_n, green_n, red_n, nir_n, swir1_n = [x/(np.nanmax(x)+eps) for x in [blue, green, red, nir, swir1]]
albedo_vis = (blue_n + green_n + red_n)/3

cond_bright = albedo_vis > 0.35
cond_swir = swir1_n > 0.15
cond_ratio = (blue_n/(red_n+1e-6)) > 1.2
cond_nir = nir_n > 0.25
cloud_init = (cond_bright.astype(int) + cond_swir.astype(int) + cond_ratio.astype(int) + cond_nir.astype(int)) >= 3

MNDWI_calc = np.where((green + swir1)!=0, (green - swir1)/(green + swir1), -1)
water_mask = MNDWI_calc > 0

size = 7
mean = ndi.uniform_filter(blue_n, size=size)
mean_sq = ndi.uniform_filter(blue_n**2, size=size)
local_std = np.sqrt(np.maximum(0, mean_sq - mean**2))
foam_mask = water_mask & (albedo_vis>0.25) & (local_std>0.03)

cloud_mask = cloud_init & (~foam_mask)
label_im, nb = ndi.label(cloud_mask)
sizes = ndi.sum(cloud_mask, label_im, range(1, nb+1))
mask_keep = np.zeros(nb+1, dtype=bool)
mask_keep[1:] = sizes >= 500
cloud_mask = mask_keep[label_im].astype(bool)

print(f" Nuages finaux : {cloud_mask.sum()} pixels")

# =========================================================
#   9. Sauvegarde indices individuels
# =========================================================
for name, data in indices.items():
    out_path = os.path.join(outdir, f"{name}.tif")
    profile_idx = profile.copy()
    profile_idx.update(dtype='float32', count=1, compress='lzw', nodata=None)
    with rasterio.open(out_path, 'w', **profile_idx) as dst:
        dst.write(data.astype('float32'), 1)
    print(f" {name} enregistré → {out_path}")

# =========================================================
#   10. Création du stack multibande robuste
# =========================================================
all_layers = list(indices.values()) + [cloud_mask.astype('float32')]
band_names = list(indices.keys()) + ['CLOUD_MASK']

ref_shape = indices['NDVI'].shape
for i, layer in enumerate(all_layers):
    if layer.shape != ref_shape:
        print(f" Redimensionnement couche {band_names[i]} de {layer.shape} → {ref_shape}")
        all_layers[i] = resize(layer, ref_shape, order=0, preserve_range=True, anti_aliasing=False).astype('float32')

profile_stack = profile.copy()
profile_stack.update(
    count=len(all_layers),
    dtype='float32',
    compress='lzw',
    nodata=None,
    height=ref_shape[0],
    width=ref_shape[1]
)

stack_path = os.path.join(outdir, 'indices_stack.tif')
with rasterio.open(stack_path, 'w', **profile_stack) as dst:
    for i, layer in enumerate(all_layers, start=1):
        dst.write(layer.astype('float32'), i)
        dst.set_band_description(i, band_names[i-1])

if os.path.exists(stack_path):
    print(f"\n Fichier multibande créé avec succès : {stack_path}")
    print(" Bandes : " + ', '.join(band_names))
else:
    print(f"\n Échec de création de {stack_path}")

print("\n Traitement complet terminé.")

Once executed, you will have the .tif files in your Google Drive.

Conclusion: towards integrated processing of Sentinel-2 images

The integration of spectral index calculation directly into the super-resolution process represents a natural evolution in Sentinel-2 data processing.
By automating both the generation of the 1 m image and the production of indices, we obtain a homogeneous, complete, and immediately usable dataset for spatial analysis.

Applied to the coastal environments of Mauritius and Rodrigues, this method opens up new perspectives:

• monitoring of coastal vegetation and mangroves,
• observation of reef flats and seagrass beds,
• detection of areas of turbidity or siltation,
• mapping of land use changes in sensitive coastal areas.

By combining all these steps into a single reproducible workflow, this script becomes a research and environmental monitoring tool particularly suited to tropical islands.
It also demonstrates that open-source tools (Sentinel-2, Google Colab, QGIS) now make it possible to achieve a level of spatial analysis that was previously reserved for professional platforms.

Next step: integrate these indices into a multi-date chronological analysis to track seasonal or post-cyclonic changes in coastal habitats.

Si cet article vous a intéressé et que vous pensez qu'il pourrait bénéficier à d'autres personnes, n'hésitez pas à le partager sur vos réseaux sociaux en utilisant les boutons ci-dessous. Votre partage est apprécié !

Atilio Francois