TD : Analyse en composantes principales d'une image satellitale.

Vincent GODARD - V3 - 05/04/2024

Cours SIG et télédétection

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

Inspiré de :

• J. Ronald Eastman, 2015, Principal component analysis for multi-spectral imagery, in TerrSet Tutorial, 256-258 (https://clarklabs.org/download/)

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/

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 TD15 un sous répertoire "data" où vous déposerez (bouton Upload Files), depuis votre zone de stockage, les fichiers :

○ "h87tm1.xxx" à "h87tm7.xxx"

○ "viet1.xxx" à "viet7.xxx"

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

1.2 Chargement des bibliothèques (library)

Dans Jupyter, la bibliothèque Numpy est probablement déjà installée. Si elle ne l'est pas, comme la bibliothèque Skimage, lors de la première utilisation, il faut l'installer.

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 rasterio
#!pip install numpy
#!pip install scikit-image
#!pip install scikit-learn


## ipympl permet d'utiliser les fonctionnalités interactives de matplotlib dans les blocs-notes Jupyter, Jupyter Lab, Google Colab et VSCode.

# !pip install ipympl # si version antérieure à Python 3.01.2, sinon message d'erreur avec le widget (window gadget)

## Import des bibliothèques

import rasterio
import numpy as np
#import sklearn
from sklearn.decomposition import PCA
import pandas as pd
import matplotlib.pyplot as plt

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__

2. ACP d'un premier exemple

Exemple dans la continuité des précédents basés sur le tutoriel du ClarkLabs et ses données de la région de Worcester, MA.

2.1. Lecture des données

Lecture d'un fichier image centré sur le secteur des Howe Hill, à l'ouest de Worcester dans le Massachusetts (USA). Les données sont issues les bandes TM1 (bande du bleu), TM2 (bande du vert), TM3 (bande du rouge), TM4 (PIR, bande du proche infra-rouge), TM5 (MIR1, bande du moyen infra-rouge le plus court), TM6 (IRTh, bande de l'infra-rouge thermique) et TM7 (MIR2, bande du moyen infra-rouge le plus long) du Landsat Thematic Mapper (TM) du 10 septembre 1987.

## Chemins d'accès aux images raster

raster_paths = ['data/h87tm1.rst', 'data/h87tm2.rst', 'data/h87tm3.rst', 'data/h87tm4.rst', 'data/h87tm5.rst', 'data/h87tm6.rst', 'data/h87tm7.rst']

## Lire et empiler des bandes raster avec une double boucle

bands_list = [] # Crée une liste vide pour stocker chaque bande de l'image raster
for raster_path in raster_paths: # Parcourt le chemin pour chaque fichier d'image raster
    src = rasterio.open(raster_path) # Ouvre le fichier image raster en utilisant rasterio
    bands = [] # Crée une liste vide pour stocker les bandes de l'image raster
    for i in range(1, src.count + 1): # Parcourt chaque bande de l'image raster
        band = src.read(i) # Lit une bande de l'image raster
        bands.append(band) # Ajoute cette bande à la liste des bandes de l'image raster
    bands_list.append(np.stack(bands, axis=-1)) # Empile les bandes de l'image raster dans un tableau 3D à partir du dernier (axis=-1)

Une fois les boucles terminées, bands_list contiendra une liste de tableaux 3D, où chaque tableau représente toutes les bandes d'une image matricielle unique empilées ensemble.

## Affichage des valeurs du 1er pixel ("en haut à gauche") au dernier pixel ("en bas à droite") de la première bande (tm1)

print(bands_list[0:1]) # [0:2] afficherait également les valeurs de la 2ème bande

2.2. Transformation du tableau en matrice des données

## Préparer la matrice de données

data_matrix = np.concatenate(bands_list, axis=-1) # Concatène la liste des tableaux 3D (bands_list) le long du dernier axe (axis=-1).
rows, cols, bands_total = data_matrix.shape # Récupère la forme du tableau data_matrix, qui représente toutes les bandes de toutes les images raster empilées ensemble.

# Cette ligne remodèle la matrice de données 3D en une matrice de données 2D.
data_matrix = np.reshape(data_matrix, (rows * cols, bands_total)) # La fonction reshape modifie la forme du tableau tout en préservant le nombre total d'éléments.

Après avoir exécuté ces lignes de code, data_matrix sera un tableau NumPy 2D où chaque ligne représente un pixel et chaque colonne représente une bande des images raster d'origine.

## Affiche la valeur (en lignes) de chaque pixel pour chaque bande (en colonnes) pour les 3 premiers pixels

print("Valeurs successives des pixels de tm1 à tm7 (en colonnes) pour le 1er, le 2ème et le 3ème pixel (en lignes) :")
print(data_matrix[0:3]) # de tm1 à tm7 pour, successivement, le 1er, le 2ème et le 3ème pixel.

## Indique combien de il y a de lignes et de colonnes dans la matrice (tableau 2D)

print('[nombre de lignes (pixels), nombre de colonnes (bands)]:', {data_matrix.shape})

2.3. Exécution de l'ACP sur une matrice

## Lancement de l'ACP

n_components = 7  # Choix du nombre de composantes principales
pca = PCA(n_components=n_components) # Création d'un objet PCA à l'aide de la bibliothèque scikit-learn, avec spécification du nombre de composants à conserver comme défini précédemment.
pca.fit(data_matrix) # Ajustement du modèle PCA à la matrice de données remodelée. Étape qui calcule les composantes principales des données.

Il faut maintenant choisir les résultats à affficher.

2.4. Calcul et analyse des coefficients de corrélation

Pour savoir quels sont les canaux d'origines porteurs d'une information "originale" (non partagée avec les autres).

Le coefficient de corrélation est une mesure de la force et de la direction de la relation linéaire entre deux variables.

Il varie entre -1 et +1.

-1 : corrélation négative parfaite.
 0 : absence de corrélation.
+1 : corrélation positive parfaite.

## Calcul de la matrice de corrélation sur les données d'origine

# Calcul du coefficient de corrélation de Pearson entre chaque paire de variables
# sur la matrice de données transposée "data_matrix"
original_correlation_matrix = np.corrcoef(data_matrix.T) 

## Choix de la précision d'affichage (nb de décimales)

np.set_printoptions(precision=2) #  Configure l'affichage des nombres décimaux avec une précision de 2 chiffres après la virgule

## Affichage de le matrice des corrélations des canaux d'origine

print("Matrice des corrélations des canaux d'origine:")
print(original_correlation_matrix)

De gauche à droite et de haut en bas : tm1 à tm7. La matrice est symétrique par rapport à la diagonale (valeur 1., corrélation d'un canal avec lui-même). • 0.9 forte corrélation positive (corr. entre tm1 et tm2) • 0.22 faible corrélation positive (corr. entre tm4 et tm3)

Quelle est le canal le plus corrélé avec tm1 ? Avec tm4 ?

2.5. Part de variance (d'information) contenue dans chaque composante

Calcul des valeurs propres (eigenvalues) pour chaque composante (néocanal) pour déterminer leur importance.

En ACP (PCA), les valeurs propres représentent la variance capturée par chaque composante principale.

Des valeurs propres élevées indiquent une composante principale capturant une grande partie de la variance totale des données.

## Création d'un tableau avec les valeurs propres (eigenvalues) et la proportion de celles-ci
## (Proportion of Variance) pour chaque composante

## Calcul des valeurs propres, du rapport de variance expliqué et de la proportion de variance

eigenvalues = pca.explained_variance_ # Récupération des valeurs propres du modèle PCA ajusté (pca).

# Le ratio de variance expliquée est la proportion de la variance totale expliquée par chaque composante principale.
# Ratio calculé en divisant chaque valeur propre par la somme de toutes les valeurs propres.
explained_variance_ratio = pca.explained_variance_ratio_ 

# Calcul de la proportion de variance
# Calcul le même résultat que la ligne précédente (explained_variance_ratio), mais "manuellement".
proportion_variance = eigenvalues / np.sum(eigenvalues)


## Création un DataFrame pour stocker les valeurs propres et les proportions

# Création du DataFrame Pandas "eigen_df" pour stocker les résultats dans un format structuré.
# Contient 3 colonnes :
eigen_df = pd.DataFrame({
    'Principal Component': range(1, n_components + 1),
    'Eigenvalue': eigenvalues,
    'Proportion of Variance': proportion_variance
})

## Impression du DataFrame

print("Valeurs propres et proportions:")
print(eigen_df)

Si la part de variance de la première composante est de plus de 86%, quelle est la part de variance des composantes 1, 2 et 3 prises séparément ?

Quelle est la part de variance des deux premières composantes réunies (faire la somme) ?

L'examen du DataFrame "eigen_df" aide à décider du nombre de composantes principales à étudier en ne conservant que les composantes qui expliquent la part la plus importante de la variance avec le moins de dimensions (bandes).

2.7. Équations de transformation

Il est également possible de calculer les coefficients de transformation des vecteurs propres (eigenvectors) pour chaque composante (néocanal).

En ACP (PCA), les vecteurs propres représentent les directions de plus grande variance dans les données.

Ils définissent les nouveaux axes (composantes principales) le long desquels les données sont projetées après réduction de dimensionnalité.

## Calcul des vecteurs propres (eigenvectors)

eigenvectors = pca.components_ # Récupération des vecteurs propres du modèle PCA ajusté (pca).

## Création d'un DataFrame pour stocker les vecteurs propres

# Avec transposition de la matrice des vecteurs propres "eigenvectors.T"
# La transposition permute les lignes et les colonnes de la matrice (chaque colonne représente un seul vecteur propre).
eigenvectors_df = pd.DataFrame(eigenvectors.T, columns=[f'Eigenvector {i+1}' for i in range(len(eigenvectors))])

# Affichage du DataFrame
print("Vecteurs propres pour chaque composante principale:")

# Ajout d'une colonne avec les numéros des composantes principales (CP)
eigenvectors_df.insert(0, "CP", range(1, n_components + 1))

print(eigenvectors_df)

Les vecteurs propres ne sont pas directement interprétables, mais ils jouent un rôle crucial dans l’ACP. Ils impriment les nouvelles directions (composantes) dans l’espace de dimensionnalité réduite qui capturent le plus de variance dans les données d’origine.

Il est plus intéressant de calculer la contribution des canaux d'origine à chaque nouvelle composante.

2.8. Contributions des canaux d'origine aux nouvelles composantes de l'ACP

Pour calculer la contribution globale de chaque bande d'origine à chaque composante principal de l'ACP, il faut élever les poids (loadings) au carré ce qui souligne la magnitude de la contribution de chaque canal d'origine à chaque composante principale.

La somme en colonnes des contributions est égale à 1. La contibution moyenne est cette somme divisée par le nombre de bandes, soit ici environ 0,14. Toutes les valeurs supérieures à la contribution moyenne indiquent que le canal concerné est contributif à cette composante.

Si pour la première composante sont considérées comme contributives les bandes 6 (IR thermique, 0.55) et 4 (PIR, 0.25), dont la part de variance est de 86 p.100, qu'en est-il de la deuxième composante, dont la part de variance n'est plus que de 11 p.100 ?

Nota : Si on veut enregistrer la matrice des contributions pour faire le calcul de la somme des colonnes (contributions) et calculer la contribution moyenne dans un tableur :

## Enregistre la matrice des contributions dans un fichier CSV

#np.savetxt("data/global_contribution_df.csv", global_contribution_df, delimiter=",")

2.9. Affichage en "images" des nouvelles composantes de l'ACP

Pour être sûr d'avoir la "bonne" matrice, on refait tourner l'ACP comme à la section 2.3.

## Lancement de l'ACP

n_components = 7  # Choix du nombre de composantes principales
pca = PCA(n_components=n_components) # Création d'un objet PCA à l'aide de la bibliothèque scikit-learn, avec spécification du nombre de composants à conserver comme défini précédemment.
pca.fit(data_matrix) # Ajustement du modèle PCA à la matrice de données remodelée. Étape qui calcule les composantes principales des données.

## On réextrait les composantes principales

principal_components = pca.transform(data_matrix)

## Remodelage des composantes principales vers la forme d'origine de l'image

principal_components_images = np.reshape(principal_components, (rows, cols, n_components))

2.9.1. Affichage d'une composante

Première composante.

## Visualisation de la première composante principale

plt.imshow(principal_components_images[:, :, 0], cmap='gray') # /!\ 0 = CP 1
plt.colorbar() # Coordonnées des pixels sur le premier axe
plt.title('Première Composante Principale')
plt.show()

Dernière composante.

## Visualisation de la dernière composante principale

plt.imshow(principal_components_images[:, :, 6], cmap='gray') # /!\ 6 = CP 7
plt.colorbar() # Coordonnées des pixels sur le premier axe
plt.title('Dernière Composante Principale')
plt.show()

Comparez la dynamique des deux images avec leur légende respective (colorbar).

De quelle information est porteuse la dernière composante ?

2.9.2. Affichage de plusieurs composantes

Les trois premières composantes.

## Préparation de l'affichage

# Créez une figure avec trois composantes disposées horizontalement
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Afficher chaque image de composante principale dans une sous-image
for i in range(3):
    axes[i].imshow(principal_components_images[:, :, i], cmap='gray')
    axes[i].set_title(f'Composante Principale {i+1}')
    axes[i].axis('off')

# Afficher
plt.tight_layout()
plt.show()

Et maintenant avec les légendes pour comparer la dynamique (variance) des composantes.

## Préparation de l'affichage

# Créez une figure avec trois composantes disposées horizontalement
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Afficher chaque image de composante principale dans une sous-image
for i in range(3):
    im = axes[i].imshow(principal_components_images[:, :, i], cmap='gray')
    axes[i].set_title(f'Composante Principale {i+1}')
    axes[i].axis('off')
    plt.colorbar(im, ax=axes[i-3])

# Afficher
plt.tight_layout()
plt.show()

Les trois suivantes.

## Préparation de l'affichage

# Créez une figure avec trois composantes disposées horizontalement
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Afficher chaque image de composante principale dans une sous-image
for i in range(3, 6):
    im = axes[i-3].imshow(principal_components_images[:, :, i], cmap='gray')
    axes[i-3].set_title(f'Composante Principale {i+1}')
    axes[i-3].axis('off')
    plt.colorbar(im, ax=axes[i-3])

# Afficher
plt.tight_layout()
plt.show()

Que voit-on de plus ?

2.9.3. Affichage de toutes les composantes

## Création d'une figure avec sept composantes disposées horizontalement

fig, axes = plt.subplots(1, 7, figsize=(20, 5))

#  Afficher chaque image de composante principale dans une sous-image
for i in range(7):
    axes[i].imshow(principal_components_images[:, :, i], cmap='gray')
    axes[i].set_title(f'Composante Principale {i+1}')
    axes[i].axis('off')

# Afficher
plt.tight_layout()
plt.show()

3. À votre tour

Sauriez-vous refaire une ACP sur les images du Viet-Nam ?