Blog SIG & Territórios

Em dois artigos anteriores, explorámos como passar de imagens Sentinel-2 padrão (resolução de 10 m) para uma versão super-resolvida de 1 m usando o módulo S2DR3, desenvolvido para a pesquisa sobre corais nas Maurícias. Esta etapa permite obter um nível de detalhe sem precedentes a partir de dados abertos, revelando estruturas finas muitas vezes invisíveis a 10 m: leitos de ervas marinhas, platôs, zonas arenosas ou mesmo transições subtis entre a vegetação e o recife.

Tradicionalmente, o cálculo dos índices espectrais (NDVI, NDWI, MNDWI, BSI, NDBI, etc.) é feito no QGIS, uma vez que as imagens são descarregadas e processadas. No entanto, muitas vezes é mais eficiente gerar esses índices diretamente na fonte, ou seja, assim que as imagens Sentinel-2 de alta resolução são produzidas.

Esta abordagem apresenta várias vantagens:

• Automatização completa: os índices são calculados assim que o conjunto de dados é criado, sem processamento manual adicional no QGIS.
• Coerência espacial: os índices utilizam exatamente as mesmas correções radiométricas e a mesma super-resolução que as bandas originais.
• Ganho de tempo: sem exportação intermédia nem configuração manual de algoritmos raster.
• Preparação para a análise: os resultados são imediatamente utilizáveis no QGIS ou em qualquer outra ferramenta SIG.

O script apresentado abaixo automatiza todo o fluxo:

1. download e super-resolução da cena Sentinel-2,
2. geração de uma máscara de nuvens inteligente (com exclusão das áreas de rebentação),
3. cálculo direto dos principais índices espectrais,
4. criação de um arquivo multibanda pronto para análise.

Essa integração simplifica consideravelmente a transição da imagem bruta para os mapas de interpretação, garantindo a reprodutibilidade do processamento científico.

Índices calculados

NDVI — Índice de Vegetação por Diferença Normalizada

• Objetivo: medir a densidade e a saúde da vegetação.
• Fórmula:

.

.

• NIR = banda do infravermelho próximo (B08)
• Red = banda vermelha (B04)
• Valores: -1 a 1; >0,2 vegetação, <0,1 solo nu ou água.
• Utilização: monitorização da vegetação, estimativa da biomassa, agricultura.

NDWI — Índice de Diferença Normalizada da Água

• Objetivo: detetar a presença de água na paisagem.
• Fórmula:

.

.

• Verde = faixa verde (B03)
• NIR = faixa infravermelha próxima (B08)
• Valores: positivos = água, negativos = solo/vegetação.
• Utilização: cartografia de lagos, rios, zonas húmidas.

MNDWI — NDWI modificado

• Objetivo: separar melhor a água das áreas urbanas ou construídas.
• Fórmula:

.

.

• SWIR1 = banda de comprimento de onda curto (B11)
• Valores: positivos = água; melhora a detecção em relação ao NDWI.
• Utilização: monitorização de massas de água em zonas urbanas ou costeiras.

BSI — Índice de solo nu

• Objetivo: detetar solos nus e zonas desnudadas.
• Fórmula:

.

.

• Azul = banda azul (B02)
• Red = B04, NIR = B08, SWIR1 = B11
• Valores: quanto mais alto o BSI, mais descoberta é a superfície.
• Utilização: monitorização de solos nus, zonas urbanas, zonas minerais.

NDBI — Normalized Difference Built-up Index (Índice de diferença normalizada de áreas construídas)

• Objetivo: detetar áreas construídas ou impermeabilizadas.
• Fórmula:

.

.

• Valores: positivos = construído, negativos = vegetação ou água.
• Utilização: cartografia urbana, monitorização da expansão das cidades.

EVI — Índice de Vegetação Aprimorado

• Objetivo: melhorar a sensibilidade do NDVI em áreas com vegetação muito densa e corrigir o efeito atmosférico.
• Fórmula:

.

.

• Azul = B02, Vermelho = B04, NIR = B08
• Valores: -1 a 1, semelhante ao NDVI, mas mais robusto em áreas densas.
• Utilização: florestas, áreas tropicais, monitorização da saúde da vegetação.

SAVI — Índice de Vegetação Ajustado ao Solo

• Objetivo: medir a vegetação corrigindo o efeito do solo (útil em áreas esparsas).
• Fórmula:

.

.

onde:

• NIR = banda de infravermelho próximo (Sentinel-2 B08)
• Red = banda vermelha (B04)
• L = fator de correção do solo, geralmente 0,5
• Valores: -1 a 1, como o NDVI, mas menos sensível ao solo nu.
• Utilização: vegetação fraca ou esparsa (pradarias, savanas, zonas costeiras com baixa cobertura vegetal).

UI — Índice Urbano (ou Índice de Urbanização)

• Objetivo: detetar e quantificar as zonas urbanas.
• Fórmula simples:

.

.

• SWIR = banda de comprimento de onda curto (Sentinel-2 B11 ou B12)
• NIR = infravermelho próximo (B08)
• Valores: positivos para áreas impermeabilizadas ou construídas, baixos/negativos para vegetação ou água.
• Utilização: monitorização da urbanização, deteção de edifícios ou superfícies artificiais.

RDI — Índice de diferença de vermelhidão

• Objetivo: detetar solos nus ou áreas ricas em óxido de ferro («vermelho») na paisagem.
• Fórmula:

.

.

• Vermelho = banda vermelha (B04), Verde = banda verde (B03)
• Valores: quanto mais «vermelho» for o solo, mais elevado será o RDI.
• Utilização: monitorização da erosão, solos nus, zonas minerais ou áridas.

Tabela resumida completa

Paleta de índices Sentinel

Cada índice é apresentado com:

• seu nome e acrônimo,
• tipo de superfície visada (vegetação, água, solo, urbano),
• valores típicos/escala,
• cor representativa para visualização no mapa (por exemplo, NDVI → verde, NDWI → azul, BSI → marrom claro, etc.).

Por que usar uma máscara de nuvens?

As nuvens e suas sombras são um problema na imagem de satélite porque:

1. Elas obscurecem as superfícies: solo, vegetação, água, corais.
2. Elas perturbam os índices (NDVI, NDWI…) porque os valores radiométricos não correspondem à superfície real.
3. As áreas cobertas por nuvens podem gerar valores atípicos: muito altos ou muito baixos, tornando os cálculos de tendências distorcidos.

Portanto, para a análise de índices e a detecção de características no solo ou na água, é importante mascarar esses pixels.

Máscara de nuvens para Sentinel-2 a 10 m

O Sentinel-2 oferece várias opções para detectar nuvens em sua resolução nativa (10 m, 20 m ou 60 m, dependendo da banda):

a) Máscara oficial fornecida pela ESA

• Arquivo QA60 nos produtos Sentinel-2 L2A: máscara binária de 10 m indicando nuvens e cirros.
• Vantagens: confiável e fácil de aplicar.
• Desvantagens: máscara grosseira → nem sempre diferencia bem nuvens finas ou sombras.

b) Métodos baseados em índices espectrais

• Faixa SWIR e NIR: as nuvens são muito refletivas no visível, mas menos no SWIR.
• Exemplo de índice simples:

.

.

• Pixels claros no visível e fracos no SWIR → provavelmente nuvens.
• Permite criar uma máscara de nuvens adaptada ao contexto, útil se se quiser detectar nuvens e cirros locais.

c) Algoritmos automáticos

• Fmask ou s2cloudless: detecção sofisticada utilizando machine learning e espectros multiespectrais.
• Prático para grandes séries temporais.

Como procedemos no script Colab

No script, não utilizamos o QA60 oficial, mas sim uma detecção por índices espectrais e textura local, adaptada às zonas costeiras e à super-resolução de 1 m:

Etapas principais:

1. Normalização das banda :
   blue_n = blue / max
green_n = green / max
red_n = red / max
nir_n = nir / max
swir1_n = swir1 / max → odos os valores entre 0 e 1 para homogeneidade..
2. Condição “nuvem inicial” baseada em 4 critérios :
   • Muito claro no visível (albedo_vis > 0,35) SWIR1 elevado (swir1_n > 0,15) Relação azul/vermelho elevada (blue_n / red_n > 1,2) NIR elevado (nir_n > 0,25)

   → os pixels que respeitam pelo menos 3/4 dos critérios são considerados nuvens..

3. Máscara “espuma” para as zonas marinhas :

   • Detecta-se a espuma/rebentação das ondas através do desvio padrão local na banda azul e no albedo visível..
   • Esses pixels são removidos da máscara de nuvens :
     cloud_final = cloud_init & (~foam_mask)
   • Isso evita mascarar a área de corais ou o mar agitado.

4. Filtragem de pequenos componentes :

   • Remove nuvens muito pequenas que poderiam ser ruído → mantemos apenas objetos > 500 pixels.ls.

5. Resultado final :

   • cloud_mask = máscara de nuvens limpa, pronta para ser usada para excluir os pixels nublados durante o cálculo dos índicess.

Por que esse método?

• A cena tem alta resolução (1 m) → o QA60 oficial não existe nessa resolução.
• A área é costeira/coralina → a espuma e as ondas devem ser diferenciadas das nuvens.
• A abordagem espectral + textura local permite:

  • detectar as nuvens
  • não ocultar a água ou a espuma
  • trabalhar com imagens de alta resolução geradas pelo S2DR3

Modo de operação

Abra o Google Colabb: https://colab.research.google.com

No Collab, vincule seu notebook ao seu Google Drive

Agora, para que o processamento não demore muito, você deve ativar uma GPU no ambiente de execução. Olhe no canto inferior direito do seu notebook, você verá algo como “Python 3” com um ícone de chip. Clique nele e selecione Alterar o tipo de ambiente de execução. Uma janela será aberta, na qual você deverá procurar por Acelerador de hardware e selecionar GPU T4. Salve as alterações e o notebook será reiniciado automaticamente com a nova configuração.

Copie este código en una nueva celda de su cuaderno y ejecútelo:

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é.")

Depois de executado, você terá os arquivos .tif no seu Google Drive.

Conclusão: rumo a um processamento integrado das imagens Sentinel-2

A integração do cálculo dos índices espectrais diretamente no processo de super-resolução representa uma evolução natural do processamento dos dados Sentinel-2.
Ao automatizar tanto a geração da imagem em 1 m quanto a produção dos índices, obtém-se um conjunto de dados homogêneo, completo e imediatamente utilizável para análise espacial.

Aplicado aos ambientes costeiros de Maurício e Rodrigues, esse método abre novas perspectivas:

• monitoramento da vegetação litorânea e dos manguezais,
• observação de recifes e ervas marinhas,
• detecção de zonas de turbidez ou assoreamento,
• cartografia das mudanças no uso do solo em zonas costeiras sensíveis.

Ao reunir todas essas etapas em um único fluxo reproduzível, esse script se torna uma ferramenta de pesquisa e monitoramento ambiental particularmente adequada para ilhas tropicais.
Ele também demonstra que as ferramentas livres (Sentinel-2, Google Colab, QGIS) permitem hoje atingir um nível de análise espacial antes reservado a plataformas profissionais.

Próxima etapa: integrar esses índices em uma análise cronológica com várias datas, a fim de acompanhar a evolução sazonal ou pós-ciclônica dos habitats costeiros.

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