Entre télédétection et intelligence artificielle, le Deep Learning ouvre de nouvelles perspectives dans l’analyse d’images satellites. Dans cet article, je vous montre comment un simple script QGIS peut exploiter un modèle de segmentation U-Net pour identifier automatiquement des zones d’intérêt sur des images multibandes, sans quitter votre environnement SIG.
Principe du script et fonctionnement général
Le script présenté ici s’appuie sur le moteur de traitement de QGIS (Processing) et sur la bibliothèque PyTorch, largement utilisée pour les modèles de Deep Learning. Il a été conçu pour exécuter, directement depuis QGIS, un modèle de segmentation d’images — typiquement un U-Net, très répandu dans les tâches de classification pixel par pixel.
Objectif
L’objectif du script est de permettre à un utilisateur de :
- charger une image satellite multibande (par exemple une image Sentinel-2) ;
- appliquer un modèle PyTorch pré-entraîné (fichier
.pth) ; - obtenir en sortie un raster de probabilité ou de classes, représentant les zones détectées par le modèle (eau, coraux, végétation, etc.).
Le traitement se fait de manière totalement intégrée à QGIS, sans ligne de commande externe.
Structure générale
Le script est défini comme un algorithme QGIS, accessible dans la boîte à outils Processing.
Il comprend trois grandes étapes :
- Chargement des données et du modèle
L’image est ouverte avec GDAL, et le modèle PyTorch est chargé en mémoire.
Le script détecte automatiquement le nombre de canaux d’entrée attendus par le modèle pour choisir les bandes appropriées (par exemple les bandes 4-3-2 de Sentinel-2 pour le RGB). - Prétraitement et masquage
Les valeurs de chaque bande sont normalisées entre 0 et 1.
Une option permet de calculer le NDWI (Normalized Difference Water Index) afin d’exclure les zones terrestres du traitement, utile pour les modèles spécialisés dans les environnements marins. - Traitement par blocs (tiling)
Pour gérer de grandes images, l’algorithme découpe la scène en blocs qui sont traités individuellement.
Chaque bloc est passé dans le modèle, et les résultats sont fusionnés à l’aide d’une fenêtre de Hanning qui adoucit les transitions entre blocs.
En sortie, le script génère un GeoTIFF géoréférencé, prêt à être superposé aux autres couches du projet QGIS.
Code complet du script
# -*- coding: utf-8 -*-
"""
Algorithme QGIS : Segmentation Coraux (U-Net) avec option masquage terrestre
Compatible QGIS 3.44+
"""
from qgis.core import (
QgsProcessing,
QgsProcessingAlgorithm,
QgsProcessingParameterRasterLayer,
QgsProcessingParameterFile,
QgsProcessingParameterEnum,
QgsProcessingParameterRasterDestination,
QgsProcessingParameterBoolean,
)
from qgis.PyQt.QtCore import QCoreApplication
import torch
import torch.serialization
import numpy as np
from osgeo import gdal
from scipy.signal import windows
class SegmentationCoraux(QgsProcessingAlgorithm):
"""Segmentation des images Sentinel-2 à 1m avec un modèle U-Net"""
BAND_OPTIONS = ["RGB (4-3-2)", "Fausse couleur (8-4-3)"]
def initAlgorithm(self, config=None):
self.addParameter(
QgsProcessingParameterRasterLayer(
"raster_input",
self.tr("Image Sentinel-2 (GeoTIFF multibandes)")
)
)
self.addParameter(
QgsProcessingParameterEnum(
"band_set",
self.tr("Sélection manuelle des bandes"),
options=self.BAND_OPTIONS,
defaultValue=0
)
)
self.addParameter(
QgsProcessingParameterFile(
"model_path",
self.tr("Modèle PyTorch (.pth)"),
extension="pth"
)
)
self.addParameter(
QgsProcessingParameterRasterDestination(
"output_raster",
self.tr("Raster de sortie (masque)")
)
)
# --- Nouvelle option : masquage terrestre ---
self.addParameter(
QgsProcessingParameterBoolean(
"mask_land",
self.tr("Masquer les zones terrestres (NDWI)"),
defaultValue=True
)
)
def processAlgorithm(self, parameters, context, feedback):
model_path = self.parameterAsFile(parameters, "model_path", context)
mask_land = self.parameterAsBool(parameters, "mask_land", context)
band_set = self.parameterAsEnum(parameters, "band_set", context)
feedback.pushInfo("Chargement du modèle PyTorch...")
try:
from segmentation_models_pytorch.decoders.unet.model import Unet
torch.serialization.add_safe_globals([Unet])
except Exception as e:
feedback.pushInfo(f"Avertissement : impossible d’ajouter Unet aux classes sûres : {e}")
model = torch.load(model_path, map_location=torch.device('cpu'), weights_only=False)
model.eval()
raster_input = self.parameterAsRasterLayer(parameters, "raster_input", context)
ds = gdal.Open(raster_input.source())
feedback.pushInfo(f"Ouverture du raster : {raster_input.source()}")
nrows, ncols = ds.RasterYSize, ds.RasterXSize
# --- Détecter automatiquement les canaux d'entrée ---
conv1_weights = model.encoder.conv1.weight.data
n_channels = conv1_weights.shape[1]
feedback.pushInfo(f"Le modèle attend {n_channels} canaux d'entrée.")
if n_channels == 3:
# Analyser l’importance des canaux
channel_mean = conv1_weights.abs().mean(dim=(0, 2, 3)).numpy()
sorted_idx = list(np.argsort(-channel_mean))
feedback.pushInfo(f"Classement des canaux par importance : {sorted_idx}")
band_selection = [4, 3, 2] # RGB classique
feedback.pushInfo(f"Utilisation automatique des bandes : {band_selection} (RGB)")
else:
band_selection = [4, 3, 2] if band_set == 0 else [8, 4, 3]
feedback.pushInfo(f"Utilisation manuelle des bandes : {band_selection}")
# --- Lecture des bandes sélectionnées ---
img_list = []
for b in band_selection:
band = ds.GetRasterBand(b)
arr = band.ReadAsArray().astype(np.float32)
img_list.append(arr)
img = np.stack(img_list, axis=0)
# Normalisation
feedback.pushInfo("Normalisation des valeurs...")
img = (img - img.min()) / (img.max() - img.min() + 1e-6)
# --- Calcul NDWI si option activée ---
if mask_land and ds.RasterCount >= 8:
try:
B3 = ds.GetRasterBand(3).ReadAsArray().astype(np.float32)
B8 = ds.GetRasterBand(8).ReadAsArray().astype(np.float32)
ndwi = (B3 - B8) / (B3 + B8 + 1e-6)
water_mask = ndwi > 0
feedback.pushInfo("NDWI calculé : les zones terrestres seront exclues.")
except Exception as e:
water_mask = np.ones((nrows, ncols), dtype=bool)
feedback.pushInfo(f"NDWI non calculable ({e}), toutes les zones seront traitées.")
else:
water_mask = np.ones((nrows, ncols), dtype=bool)
feedback.pushInfo("Masquage terrestre désactivé.")
# --- Traitement par blocs ---
block_size = 2048
overlap = 512
output_mask = np.zeros((nrows, ncols), dtype=np.float32)
weight = np.zeros((nrows, ncols), dtype=np.float32)
total_blocks = ((nrows // (block_size - overlap) + 1) *
(ncols // (block_size - overlap) + 1))
done = 0
feedback.pushInfo("Début du traitement par blocs...")
for y in range(0, nrows, block_size - overlap):
for x in range(0, ncols, block_size - overlap):
if feedback.isCanceled():
break
block = img[:, y:y+block_size, x:x+block_size]
if block.size == 0:
continue
mask_block_water = water_mask[y:y+block_size, x:x+block_size]
if not mask_block_water.any():
continue
# ======= BLOC ORIGINAL PRÉSERVÉ =======
with torch.no_grad():
block_tensor = torch.from_numpy(block).unsqueeze(0)
pred = model(block_tensor)
mask_block = pred.squeeze().numpy()
h, w = mask_block.shape
mask_block *= mask_block_water[:h, :w]
# Fenêtre Hanning
win_y = windows.hann(h)[:, None]
win_x = windows.hann(w)[None, :]
weight_block = win_y * win_x
output_mask[y:y+h, x:x+w] += mask_block * weight_block
weight[y:y+h, x:x+w] += weight_block
done += 1
progress = int(100 * done / total_blocks)
feedback.setProgress(progress)
# ======= FIN DU BLOC ORIGINAL =======
# Fusion finale
output_mask /= np.maximum(weight, 1e-6)
output_mask = np.clip(output_mask, 0, 1)
# --- Sauvegarde GeoTIFF ---
feedback.pushInfo("Sauvegarde du masque final...")
driver = gdal.GetDriverByName("GTiff")
out_path = self.parameterAsOutputLayer(parameters, "output_raster", context)
out_ds = driver.Create(out_path, ncols, nrows, 1, gdal.GDT_Float32)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(output_mask)
out_ds.FlushCache()
feedback.pushInfo("✅ Segmentation terminée avec succès.")
return {"output_raster": out_path}
# --- Métadonnées ---
def name(self):
return "segmentation_coraux"
def displayName(self):
return self.tr("Segmentation des images Sentinel 2 à 1m (U-Net)")
def group(self):
return self.tr("Deep Learning")
def groupId(self):
return "deeplearning"
def tr(self, string):
return QCoreApplication.translate("SegmentationCoraux", string)
def createInstance(self):
return SegmentationCoraux()
Opérations préalables
Nous allons créer un script de traitement QGis. Comme préalable il faut installer
python -m pip install torchgeo
python -m pip install torch torchvision
python -m pip install segmentation-models-pytorch
Puis nous devons télécharger le modèle souhaité. Pour cela, dans la console python de QGis entrez le script suivant:
import segmentation_models_pytorch as smp
import torch
# U-Net RGB pré-entraîné
model = smp.Unet(
encoder_name="resnet34",
encoder_weights="imagenet",
in_channels=3,
classes=1 # corail/non-corail
)
# Sauvegarde pour QGIS
torch.save(model, "c:/models/unet_coraux.pth")
Dans cet exemple on sauvegarde les modèles téléchargés dans un répertoire c:/models
Mise en place du script de traitement
Dans QGIS → Boîte à outils de traitement → Scripts → Nouveau script, collez le code ci-dessus.
Structure et paramètres du script dans QGIS
Le script est entièrement intégré à QGIS via la boîte à outils Processing, ce qui permet de l’utiliser sans quitter l’interface graphique. Il est conçu pour être flexible et compatible avec différents modèles U-Net ou autres modèles de segmentation PyTorch.
Paramètres principaux
- Image d’entrée (
raster_input)- Type : raster multibande (GeoTIFF recommandé)
- Exemple : image Sentinel-2 avec 10 bandes.
- Le script lit automatiquement les bandes nécessaires selon le nombre de canaux attendu par le modèle, mais un choix manuel est également possible.
- Sélection des bandes (
band_set)- Type : énumération (RGB classique 4-3-2 ou fausse couleur 8-4-3)
- Le script peut détecter automatiquement les canaux d’entrée du modèle et sélectionner les bandes les plus pertinentes.
- Modèle PyTorch (
model_path)- Type : fichier
.pth - Contient le réseau pré-entraîné (U-Net ou autre) destiné à la segmentation.
- Compatible avec des modèles personnalisés si le nombre de canaux correspond à l’image.
- Type : fichier
- Masquage des zones terrestres (
mask_land)- Type : case à cocher (booléen)
- Si activé, le script calcule le NDWI pour exclure les zones terrestres du traitement, utile pour les modèles ciblant uniquement l’eau.
- Raster de sortie (
output_raster)- Type : raster (GeoTIFF)
- Résultat final de la segmentation, géoréférencé et prêt à être visualisé dans QGIS.
Points clés du fonctionnement
- Traitement par blocs (tiling) : chaque bloc est traité indépendamment pour réduire la charge mémoire et permettre des images de grande taille.
- Fenêtre de Hanning : permet de fusionner les blocs sans créer de bords nets ou d’artefacts visibles.
- Flexibilité du modèle : tant que le modèle PyTorch accepte le nombre de canaux de l’image, il peut être utilisé avec ce script.
Cette structure rend le script polyvalent, adapté à la segmentation de coraux, de zones végétalisées, ou de tout autre objet détectable sur une image multibande.
Fonctionnement interne : traitement par blocs et prédiction PyTorch
Le script traite l’image en blocs (tiles) pour gérer efficacement de grandes images Sentinel-2 sans saturer la mémoire. Chaque bloc suit le même processus :
- Extraction du bloc
- L’image multibande est découpée en blocs de taille fixe (2048×2048 par défaut) avec un recouvrement (overlap) de 512 pixels pour lisser les transitions.
- Chaque bloc contient uniquement les bandes sélectionnées selon le modèle.
- Conversion en tenseur PyTorch
with torch.no_grad():block_tensor = torch.from_numpy(block).unsqueeze(0) pred = model(block_tensor) mask_block = pred.squeeze().numpy()torch.from_numpy(block)convertit le bloc en tenseur PyTorch.unsqueeze(0)ajoute une dimension batch nécessaire au modèle.model(block_tensor)applique le modèle pré-entraîné pour produire un masque de prédiction.squeeze()supprime la dimension batch pour obtenir un tableau 2D.
- Application du masque eau (optionnel)
- Si le masquage des zones terrestres est activé, le script applique le NDWI pour ne garder que les pixels d’eau :
mask_block *= mask_block_water[:h, :w] - Cela permet au modèle de ne pas se tromper sur les zones terrestres.
- Si le masquage des zones terrestres est activé, le script applique le NDWI pour ne garder que les pixels d’eau :
- Fenêtre de Hanning pour la fusion des blocs
- Chaque bloc est pondéré avec une fenêtre Hanning pour lisser les bords :
win_y = windows.hann(h)[:, None] win_x = windows.hann(w)[None, :] weight_block = win_y * win_x - Cette pondération réduit les artefacts sur les limites des blocs lors de la reconstruction de l’image complète.
- Chaque bloc est pondéré avec une fenêtre Hanning pour lisser les bords :
- Fusion pondérée dans l’image finale
- Le masque pondéré est ajouté à l’image de sortie :
output_mask[y:y+h, x:x+w] += mask_block * weight_block weight[y:y+h, x:x+w] += weight_block - La division finale par le poids total normalise le masque sur toute l’image :
output_mask /= np.maximum(weight, 1e-6) output_mask = np.clip(output_mask, 0, 1)
- Le masque pondéré est ajouté à l’image de sortie :
- Résultat final
- L’image complète est sauvegardée en GeoTIFF, avec géoréférencement intact.
- Les blocs recouverts et pondérés garantissent un masque continu et homogène, prêt pour la visualisation ou l’analyse dans QGIS.
Cette méthode est générique : elle fonctionne avec n’importe quel modèle PyTorch compatible, tant que le nombre de canaux d’entrée correspond aux bandes de l’image. Elle combine efficacité mémoire, précision de la segmentation et lissage des artefacts de bloc, ce qui la rend particulièrement adaptée aux images satellites haute résolution.
Détection automatique des bandes et importance des canaux
L’un des points forts du script est sa capacité à adapter automatiquement les bandes Sentinel-2 utilisées selon le modèle PyTorch. Cela permet de travailler avec différents modèles U-Net sans modifier le code manuellement.
1. Identification du nombre de canaux d’entrée
Chaque modèle U-Net possède une première couche de convolution (conv1) qui définit le nombre de canaux attendus. Le script analyse cette couche pour détecter le nombre de canaux :
conv1_weights = model.encoder.conv1.weight.data
n_channels = conv1_weights.shape[1]
feedback.pushInfo(f"Le modèle attend {n_channels} canaux d'entrée.")
- Si le modèle attend 3 canaux, le script choisira par défaut les bandes RGB classiques (B4-Rouge, B3-Vert, B2-Bleu).
- Si le modèle utilise plus ou moins de 3 canaux, le script propose une option manuelle, ou peut être ajusté pour des modèles multispectraux.
2. Importance des canaux
Pour optimiser la qualité de la segmentation, le script peut évaluer l’importance de chaque canal sur la première couche de convolution :
channel_mean = conv1_weights.abs().mean(dim=(0,2,3)).numpy()
sorted_idx = list(np.argsort(-channel_mean))
feedback.pushInfo(f"Classement des canaux par importance : {sorted_idx}")
- Cette approche identifie les canaux les plus “activés” par le modèle.
- On peut ensuite associer chaque canal à la bande Sentinel-2 la plus pertinente, par exemple pour un modèle RGB ou une combinaison fausse couleur (B8-NIR, B4-Rouge, B3-Vert).
3. Choix automatique des bandes
- Pour les modèles 3 canaux classiques :
[B4, B3, B2](RGB) ou[B8, B4, B3](faux couleurs). - Pour les modèles multispectraux : le script peut être adapté pour mapper les canaux du modèle aux bandes Sentinel-2 disponibles.
Cette logique permet au script de s’adapter à différents modèles U-Net, rendant l’outil polyvalent et simple à utiliser pour différents projets de segmentation d’images satellite dans QGIS.
Masquage des zones terrestres avec NDWI
Pour des applications comme la segmentation des coraux ou tout autre objet aquatique, il est souvent utile d’ignorer les zones terrestres. Le script intègre cette fonctionnalité grâce à l’indice NDWI (Normalized Difference Water Index).
1. Calcul de l’indice NDWI
Le NDWI est calculé à partir de deux bandes Sentinel-2 :
- B3 (vert)
- B8 (proche infrarouge, NIR)
La formule est simple :
.
\(\text{NDWI} = \frac{B3 – B8}{B3 + B8 + \epsilon}\)
.
où (\(\epsilon\)) est un petit nombre pour éviter la division par zéro.
B3 = ds.GetRasterBand(3).ReadAsArray().astype(np.float32)
B8 = ds.GetRasterBand(8).ReadAsArray().astype(np.float32)
ndwi = (B3 - B8) / (B3 + B8 + 1e-6)
water_mask = ndwi > 0
water_mask: masque booléen oùTruecorrespond à l’eau etFalseà la terre.
2. Option dans le script
Le script offre une case à cocher pour activer ou désactiver ce masquage :
mask_land = self.parameterAsBool(parameters, "mask_land", context)
- Si activé, seules les zones aquatiques sont traitées par le modèle.
- Si désactivé, toutes les zones du raster sont considérées, utile si le modèle est générique ou si l’on souhaite segmenter aussi la terre.
3. Application au traitement par blocs
Lors du traitement en blocs du raster, chaque bloc est multiplié par le masque aquatique avant la prédiction du modèle :
mask_block *= mask_block_water[:h, :w]
- Cela évite que le modèle tente de segmenter des zones terrestres, améliorant la précision et réduisant le bruit.
- Couplé à la fenêtre de pondération Hanning, le masque aquatique assure une transition douce entre blocs, limitant les artefacts sur les bords.
Cette approche rend le script plus robuste pour l’imagerie maritime ou lacustre, tout en restant flexible pour d’autres types de segmentation.
Traitement par blocs et pondération Hanning
Lorsque l’on travaille avec des images satellites de grande taille, comme les images Sentinel‑2, il est souvent impossible de traiter l’intégralité de l’image en une seule passe à cause des contraintes de mémoire. Pour contourner ce problème, le script de segmentation utilise un traitement par blocs.
Le principe est simple :
- Découpage de l’image en blocs : l’image est divisée en petites fenêtres de taille fixe (par exemple 2048×2048 pixels) avec un recouvrement entre les blocs (par exemple 512 pixels). Ce recouvrement permet d’éviter les effets de bord lors de la prédiction.
- Application du modèle sur chaque bloc : chaque bloc est converti en tenseur PyTorch et passé dans le modèle U‑Net pour produire un masque de probabilité.
- Masquage des zones non pertinentes : si l’option est activée, un masque NDWI permet d’exclure les zones terrestres afin que le modèle ne prédit que sur l’eau.
- Pondération Hanning pour fusionner les blocs : les prédictions des blocs se chevauchant sont combinées à l’aide d’une fenêtre Hanning, qui attribue moins de poids aux pixels proches des bords du bloc et plus de poids aux pixels centraux. La formule appliquée est :
. \(\text{mask_final}[y:y+h, x:x+w] += \text{mask_bloc} \times \text{Hanning}(h, w) \\
\) \(\text{weight}[y:y+h, x:x+w] += \text{Hanning}(h, w)
\).
- Fusion finale : après traitement de tous les blocs, le masque final est obtenu en divisant la somme pondérée par la somme des poids, ce qui assure une transition douce entre les blocs :
.
. \(\text{mask_final} = \frac{\text{mask_final}}{\text{weight}}
\).
Ce mécanisme garantit que le modèle peut traiter des images très grandes tout en conservant la précision des prédictions et en minimisant les artefacts aux frontières des blocs.
Traitement par blocs et pondération Hanning
Pour traiter des images Sentinel‑2 complètes, qui peuvent mesurer plusieurs milliers de pixels de côté, le script découpe l’image en blocs. Chaque bloc est envoyé au modèle U‑Net pour générer un masque de segmentation local. Cette approche présente plusieurs avantages :
- Gestion de la mémoire : les blocs de taille raisonnable (par exemple 2048 × 2048 pixels) permettent de ne pas saturer la mémoire RAM ou GPU lors du traitement.
- Parallélisation possible si nécessaire pour accélérer le calcul.
- Fusion plus précise : les blocs peuvent se recouvrir pour éviter les artefacts aux bords.
Recouvrement des blocs et pondération
Les blocs se superposent sur une zone d’overlap (par exemple 512 pixels). Pour combiner les prédictions des blocs qui se recouvrent, le script utilise une fenêtre de pondération Hanning.
La fenêtre Hanning lisse les valeurs des pixels aux bords des blocs afin de réduire les discontinuités lors de la fusion. Mathématiquement, pour un bloc de hauteur (h) et largeur (w) :
.
\(\text{weight_block} = \text{hann}(h) \otimes \text{hann}(w)\)
.
où (\otimes) représente le produit extérieur des deux vecteurs de pondération 1D de Hanning pour les lignes et les colonnes.
Chaque bloc traité fournit un masque prédictif (mask_block) que l’on multiplie par la pondération :
.
\(\text{output_mask}[y:y+h, x:x+w] += mask_block \times weight_block\)
.
Les poids sont également accumulés pour normaliser correctement la fusion :
.
\(\text{weight}[y:y+h, x:x+w] += weight_block\)
.
Une fois tous les blocs traités, la fusion finale se fait par division des valeurs accumulées par la somme des poids :
.
\(\text{output_mask} = \frac{\text{output_mask}}{\max(\text{weight}, \varepsilon)}\)
.
Cette approche garantit un masque continu, lisse et homogène, sans artefacts liés aux bords des blocs.
Sauvegarde du masque GeoTIFF et utilisation dans QGIS
Une fois la segmentation complète et la fusion des blocs effectuée, le script crée un fichier raster GeoTIFF correspondant au masque prédictif final. Ce fichier conserve :
- la projection spatiale de l’image d’origine,
- la géoréférence (transformée géométrique) pour assurer un alignement parfait avec les autres couches dans QGIS,
- une seule bande contenant les valeurs de probabilité de présence de la classe d’intérêt (eau, coraux, végétation, etc.), normalisées entre 0 et 1.
Le script utilise GDAL pour cette étape :
driver = gdal.GetDriverByName("GTiff")
out_ds = driver.Create(out_path, ncols, nrows, 1, gdal.GDT_Float32)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(output_mask)
out_ds.FlushCache()
Visualisation et exploitation dans QGIS
Dans QGIS, il suffit de charger le fichier GeoTIFF généré :
- Les pixels proches de 1 correspondent à des zones fortement prédites par le modèle.
- Les pixels proches de 0 correspondent à des zones absentes de la classe ciblée.
Pour faciliter l’interprétation, on peut :
- appliquer une rampes de couleurs (par exemple du bleu clair au bleu foncé pour les coraux),
- créer un masque binaire en appliquant un seuil de probabilité,
- superposer le masque sur l’image d’origine ou d’autres couches SIG pour analyse spatiale.
Ainsi, le script transforme automatiquement les prédictions du modèle U‑Net en une couche raster prête à l’analyse SIG, tout en conservant la précision géographique.
