TD : Faire une classification supervisée.

Vincent GODARD - V1 - 29/04/2025

Cours SIG et télédétection

Département de géographie - L3 - Université de Paris 8

Inspiré de :

• J. Ronald Eastman, 2015, Supervised Classification, in TerrSet Tutorial, 256-258 (https://clarklabs.org/download/)

• Machine Learning > Carte d’occupation des sols avec SENTINELLE 2 https://geoafrica.fr/data-science-carte-doccupation-des-sols-2018avec-sentinelle-2/

• Mastering Machine Learning based Land Use Classification with Python: A Comprehensive Guide! https://waleedgeo.medium.com/lulc-py-78cb954673d

Sources :

• données Landsat Thematic Mapper (TM) du 10 septembre 1987, ouest de Worcester, Massachusetts (issues du https://clarklabs.org/download/)

• scikit-image https://scikit-image.org/

• matplotlib https://matplotlib.org/

Voir où est-ce : https://www.peakbagger.com/peak.aspx?pid=43307

    Téléchargement des documents nécessaires :

    Dossier compressé à télécharger => ici.

1. Chargement des données et des bibliothèques

1.1 Chargement des données

Créez dans le répertoire TD16 un sous répertoire "data" où vous déposerez (bouton Upload Files), depuis votre zone de stockage, les fichiers :

○ "how87tm1.rst et rdc" ; "how87tm2.rst et rdc" ; etc.

○ "sample_ods_6cl.tif" qui contient les catégories d'occupation du sol en mode raster

Ce noteboock et le répertoire "data" seront au même niveau dans l'arborescence, tous les deux sous TD16.

1.2 Chargement des bibliothèques (library)

Si dans le JupyterLab la bibliothèque pandas est déjà installée, la bibliothèque skimage ne l'est pas, il faut l'installer, idem pour le Jupyter Notebook.

Il faut ensuite les charger (fonction import).

## Si vous n'avez jamais installé skimage sur votre ordinateur, il est probablement absent des bibliothèques par défaut. 
## Pour l'installer, effacez le # (avant le "!").

#!pip install scikit-image
#!pip install rasterio

## Import des bibliothèques
import numpy as np
from pathlib import Path
from skimage.io import imread
from sklearn.ensemble import RandomForestClassifier
import rasterio
from rasterio.transform import Affine
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

Il se peut que l'import de la version de skimage génère un message d'erreur, comme une incompatibilité avec NumPy, par exemple. Une des bibliothèques est trop ancienne par rapport à l'autre. Il faut donc faire une mise à jour.

## SciPy est un projet qui fédère des bibliothèques Python à usage scientifique.
## Scipy utilise les tableaux et matrices du module NumPy (https://fr.wikipedia.org/wiki/SciPy).

## Mise à jour de scipy, effacez le # (avant le "!") si nécessaire.
#!pip install --upgrade scipy

## Ensuite, il faut l'importer
#import scipy

## Puis vérifier que cette nouvelle version est compatible
#scipy.__version__

On pourra avantageusement regarder à quoi servent les bibliothèques https://scikit-learn.org/stable/ ou https://scikit-image.org/ dont il existe parfois des Wiki en français https://fr.wikipedia.org/wiki/Scikit-learn ou https://fr.wikipedia.org/wiki/Scikit-image.

2. Lecture des bandes de l'image TM et des vérités terrain

Lecture d'un fichier image centré sur le secteur des Howe Hill, à l'ouest de Worcester dans le Massachusetts (USA, https://www.peakbagger.com/peak.aspx?pid=43307). Ce fichier image est composé de 7 canaux issus des bandes TM1 (bande du bleu), TM2 (bande du vert), TM3 (bande du rouge), TM4 (bande du PIR), TM5 (bande du MIR-1), TM6 (bande de l'IR thermique) et TM7 (bande du MIR-2) du Landsat Thematic Mapper (TM) du 10 septembre 1987.

2.1. Information sur les métadonnées, dont la taille des images

Les images rasters sont issues d'Idrisi, maintenant TerrSet/liberaGIS, elles sont en binôme : un fichier rst (qui contient l'image) et un fichier rdc qui contient les métadonnées.

# Affichage de la taille de l'image .rst en lisant les métadonnées du fichier .rdc (cf. columns et rows)
# Vérifier le contenu réel du fichier .rdc
print(Path("data/h87tm1.rdc").read_text())

# Affichage de la taille de l'image TIFF (celle des échantillons de vérités terrain)
tif_path = Path("data/sample_ods_6cl.tif")
if tif_path.exists():
    img = imread(tif_path)
    print(f"{tif_path.name} : dimensions {img.shape}")
else:
    print(f"{tif_path.name} : fichier absent")

Si les fichiers ont la même taille [nombre de lignes (rows) et de colonnes (columns)] identiques, alors on poursuit.

2.2. Lecture des canaux de l'image TM

Lecture des seuls fichiers .rst car, dans cet exemple, le ref. system est "plane". Pas de système géodésique ici.

## Chemin des fichiers
data_dir = Path("data")
output_path = "results/h87_fr_6cl.tif"

## Paramètres fixes (d'après h87tm1.rdc)
## Taille des images lue dans les métadonnées du fichier .rdc (cf. columns et rows)
ROWS = 86
COLS = 72
RESOLUTION = 30.0  # 30m/pixel
MIN_X, MAX_Y = 3930.0, 3600.0  # Coin supérieur gauche

## === Lecture des bandes RST ===
bands = []
for i in range(1, 8):
    rst_file = data_dir / f"h87tm{i}.rst"
    arr = np.fromfile(rst_file, dtype=np.uint8).reshape(ROWS, COLS)
    bands.append(arr)

X = np.dstack(bands).reshape(-1, 7)  # Shape (86*72, 7)

2.3. Chargement de la vérité terrain

Par convention standard en *machine learning* : X majuscule pour les features, y minuscule pour les labels.

# === Chargement vérité terrain ===
def load_ground_truth():
    gt = imread(data_dir / "sample_ods_6cl.tif")
    return gt.flatten()

y = load_ground_truth()

2.4. Entraînement et prédiction

## Réalisation de la classification
mask = y > 0
clf = RandomForestClassifier(n_estimators=100).fit(X[mask], y[mask])
classification = clf.predict(X).reshape(ROWS, COLS)

2.5. Création du profil géographique

# === Création du profil géographique ===
transform = Affine(RESOLUTION, 0, MIN_X, 0, -RESOLUTION, MAX_Y)
profile = {
    'driver': 'GTiff',
    'height': ROWS,
    'width': COLS,
    'count': 1,
    'dtype': 'uint8',
    'crs': '',  # À vérifier dans vos métadonnées. Ne rien mettre entre '' si plane !
    # EPSG:32619 - WGS 84 / UTM zone 19N (utilisé pour les données globales et GPS)
    # EPSG:26919 - NAD83 / UTM zone 19N (utilisé pour les données nord-américaines)
    'transform': transform,
    'nodata': 0
}

2.6. Sauvegarde du résultat

# === Sauvegarde avec rasterio (si besoin de géoréférencement) ===
with rasterio.open(output_path, 'w', **profile) as dst:
    dst.write(classification.astype(np.uint8), 1)

print(f"✅ Carte sauvegardée : {output_path}")

2.7. Affichage des résultats bruts

# === Visualisation brute ===
plt.figure(figsize=(10, 8))
img = plt.imshow(classification, cmap='tab10', extent=[MIN_X, MIN_X + COLS*RESOLUTION, MAX_Y - ROWS*RESOLUTION, MAX_Y])
plt.colorbar(label='Classes')
plt.title("Classification de l'occupation du sol par Random Forest Classifier")
plt.xlabel('Coordonnée X (m)')
plt.ylabel('Coordonnée Y (m)')
plt.show()

2.8. Affichage des résultats légendés et recolorés

## Amélioration des couleurs et des étiquettes de légende
#import matplotlib.pyplot as plt
#import matplotlib.colors as mcolors

# Palette personnalisée pour 6 classes (voir ci-dessous pour d'autres couleurs)
cmap = mcolors.ListedColormap([
    '#01579b',  # Classe 1 - Eaux peu profondes
    '#08306b',  # Classe 2 - Eaux profondes
    '#f7dc6f',  # Classe 3 - Agriculture
    '#e31a1c',  # Classe 4 - Urbain
    '#28b463',  # Classe 5 - Feuillus
    '#196f3d',  # Classe 6 - Conifères
])

# Norme pour mapper 1-6 directement
norm = mcolors.BoundaryNorm(boundaries=[0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5], ncolors=6)

# Pour l'affichage géoréférencé
# (MIN_X, MAX_Y, COLS, ROWS, RESOLUTION doivent être définis dans ton script principal)
extent = [
    MIN_X,                   # gauche (X min)
    MIN_X + COLS * RESOLUTION,  # droite (X max)
    MAX_Y - ROWS * RESOLUTION,  # bas (Y min)
    MAX_Y                    # haut (Y max)
]

plt.figure(figsize=(10, 8))
img = plt.imshow(classification, cmap=cmap, norm=norm, extent=extent, origin='upper')
plt.title("Howe Hill 1987, MA, USA\n Classification de l'occupation du sol par Random Forest Classifier (6 classes)")
cbar = plt.colorbar(img, ticks=[1, 2, 3, 4, 5, 6])
cbar.ax.set_yticklabels([
    'Eaux peu profondes',
    'Eaux profondes',
    'Agriculture',
    'Urbain',
    'Feuillus',
    'Conifères'
])
plt.xlabel('Coordonnée X (m)')
plt.ylabel('Coordonnée Y (m)')
plt.tight_layout()
plt.show()

Liste non exhaustive de nuanciers complets avec les codes hexadécimaux à utiliser ensuite dans une *ListedColormap*, quelques ressources en ligne :

• HTML Color Codes (https://htmlcolorcodes.com/fr/) : propose un sélecteur de couleur interactif, des tableaux de couleurs populaires, et des palettes prêtes à l’emploi avec les codes hexadécimaux, RGB et les noms de couleurs. 

• Color Hex (https://www.color-hex.com/) : permet de rechercher une couleur, d’obtenir son code hexadécimal, et de générer des palettes harmonieuses (triadiques, complémentaires, etc.). Il affiche aussi des palettes des exemples CSS.

• ColorSpace et d’autres générateurs de palettes (https://graphiste.com/blog/choisir-palette-couleurs/) : entrez une couleur de base et l’outil propose des combinaisons variées, avec tous les codes hex nécessaires.

• Liste des couleurs Pantone et équivalents hex (https://www.agence-casanova.fr/professionnel/identitaire/pantone.html) : palettes inspirées du nuancier Pantone avec équivalents hexadécimaux.

• Webdi Couleurs web (https://webdi.fr/couleur-hexa.php) : tableaux de couleurs avec noms et codes hexadécimaux.

• Adobe Color (https://color.adobe.com/fr/create/color-wheel) : color.adobe.com pour créer, explorer et exporter des palettes hex.

Mais quelle est la qualité du résultat ? Comment l'évaluer ?

3. Évaluation de la classification

Classiquement, l'échantillon des "vérités terrain" est scindé en deux, 70p.100 des pixels sont utilisés pour initier les traîtements et 30p.100 pour les contrôler. Une visualisation est réalisée dans une matrice de confusion. Les pixels sur la diagonale sont déclarés bien classés, ceux qui s'en écartent, ne le sont pas et indiquent d'où ils viennent ou vers quel thème ils ont migré. L'indice de Kappa est une indice de performance globale.

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix, classification_report, cohen_kappa_score
import matplotlib.pyplot as plt
import seaborn as sns

# === 1. Filtrer les pixels non étiquetés (classe 0) AVANT la séparation ===
mask = y > 0
X_filtered = X[mask]
y_filtered = y[mask]

# === 2. Séparation des données en conservant toutes les classes ===
X_train, X_test, y_train, y_test = train_test_split(
    X_filtered, 
    y_filtered,
    test_size=0.3, # taille à faire varier si résultat trop bon (cf. infra)
    random_state=42,
    stratify=y_filtered  # Maintient la distribution des classes
)

# === 3. Entraînement et prédiction ===
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)

# === 4. Définir explicitement les classes attendues (1-6) ===
classes = [1, 2, 3, 4, 5, 6]
class_names = [
    'Eaux peu profondes', 
    'Eaux profondes', 
    'Agriculture', 
    'Urbain', 
    'Feuillus', 
    'Conifères'
]

# === 5. Matrice de confusion complète ===
cm = confusion_matrix(y_test, y_pred, labels=classes)

plt.figure(figsize=(10,8))
sns.heatmap(
    cm, 
    annot=True, 
    fmt='d', 
    cmap='Blues',
    xticklabels=class_names,
    yticklabels=class_names
)
plt.title('Matrice de confusion (6 classes)')
plt.xlabel('Prédictions')
plt.ylabel('Vérité terrain')
plt.show()

# === 6. Rapport de classification avec gestion des classes manquantes ===
print(classification_report(
    y_test, 
    y_pred, 
    labels=classes,
    target_names=class_names,
    zero_division=0  # Supprime les avertissements
))

# === 7. Coefficient Kappa ===
kappa = cohen_kappa_score(y_test, y_pred)
print(f"Coefficient Kappa de Cohen : {kappa:.3f}")

Cas particulier détecté :

Votre sortie montre seulement 50 échantillons de test (support total si test = 0.3). Avec si peu de données :

• Un Random Forest peut atteindre 100% de précision sur de petits jeux de test

• Le coefficient Kappa devient instable

3.3. Vérification des dimensions

print(f"Train shape: {X_train.shape}, Test shape: {X_test.shape}")
# Doit retourner (ex) : Train (70, 7), Test (30, 7)

3.4. Analyse de la distribution des classes

Où l'on voit que certaines classes ont très peu de pixels tests !

Il faut accroître le nombre de pixels test => test_size=0.5, # 50% en test

from collections import Counter
print("Distribution train:", Counter(y_train))
print("Distribution test:", Counter(y_test))

3.5. Test de permutation (pour éliminer le sur-apprentissage) :

from sklearn.model_selection import cross_val_score
scores = cross_val_score(clf, X, y, cv=5)
print(f"Accuracy cross-val: {scores.mean():.2f} ± {scores.std():.2f}")

Interprétation des composantes

a) Précision moyenne (0.98)

Indique une performance globale excellente

Dans la plupart des contextes, >0.9 est considéré comme très bon

b) Écart-type (0.04)

Mesure la variabilité entre les folds

Une valeur faible (<0.05) suggère une bonne stabilité du modèle

Signifie que le modèle ne dépend pas fortement de découpages aléatoires spécifiques

3.6. Études complémentaires recommandées

a) Vérifier l'équilibre des classes

from collections import Counter
print(Counter(y))

from collections import Counter
print("Distribution train:", Counter(y_train))
print("Distribution test:", Counter(y_test))

b) Comparer avec un hold-out set (jeu test à part)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
clf.fit(X_train, y_train)
print("Accuracy hold-out:", clf.score(X_test, y_test))

Accuracy = nombre total de prédictions / nombre de prédictions correctes. Si => 1, c'est parfait (ou louche, cf. supra).

c) Analyser d'autres métriques

from sklearn.metrics import f1_score, roc_auc_score
print("F1-score:", cross_val_score(clf, X, y, cv=5, scoring='f1_macro').mean())

Le F1-score est une métrique d’évaluation des modèles de classification qui combine la précision (proportion de bonnes prédictions parmi les positifs prédits) et le rappel (proportion de vrais positifs détectés parmi tous les positifs réels) en une seule valeur, via leur moyenne harmonique.

F1 = 2 × (précision x rappel) / (précision + rappel)

Le F1-score varie de 0 à 1 :

• 1 indique un modèle parfait (précision et rappel à 1)

• 0 indique de très mauvaises performances

d) Tester la significativité statistique

from scipy.stats import ttest_1samp
_, p_value = ttest_1samp(scores, popmean=0.95)
print(f"p-value vs 95%: {p_value:.4f}")

Interprétation directe

1. Comparaison au seuil :

    • p-value (0.1550) > α (0.05)

    • Conclusion : Vous ne pouvez pas rejeter l'hypothèse nulle (H₀).

    • La différence/relation observée dans vos données n'est pas statistiquement significative au seuil de 5%.

2. Signification concrète :

    • Il y a 15.5% de probabilité d'observer un résultat aussi extrême (ou plus) que celui de votre étude si l'hypothèse nulle était vraie.

    • Ce niveau de risque est jugé trop élevé pour conclure à un effet significatif.

Recommandations

1. Ne pas conclure à "pas d'effet" :

    • Une p-value élevée ne prouve pas que H₀ est vraie, seulement que les données ne permettent pas de la rejeter.

    • Exemple : Si vous testiez l'efficacité d'un médicament, cela ne signifie pas qu'il est inefficace, mais que votre étude n'a pas détecté d'effet significatif.

2. Vérifier la puissance statistique :

    • Un manque de puissance (échantillon trop petit) peut expliquer un résultat non significatif.

    • Calculez la taille d'effet observée et estimez la puissance a posteriori.

3. Analyser les intervalles de confiance :

    • Si l'intervalle de confiance à 95% inclut la valeur nulle (ex: 0 pour une différence de moyennes), cela confirme la non-significativité.

4. Test avec SVM (Support Vector Machine)

À des fins de comparaison avec Random Forest.