## 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
import rasterio
from sklearn.cluster import KMeans
import glob
import pathlib
import matplotlib.pyplot as plt

## 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__

## Chargement des bandes TM et initialisation du profile
band_files = sorted(glob.glob('data/h87tm*.rst'))  # À adapter si besoin
bands = [rasterio.open(band) for band in band_files]
profile = bands[0].profile  # <-- INITIALISATION CRUCIALE ICI

# Création du stack multibande
stack = np.dstack([band.read(1) for band in bands])
height, width, n_bands = stack.shape
data = stack.reshape(-1, n_bands)

## Configuration K-means
n_clusters = 10  # Nb de classes de départ à ajuster selon les besoins
kmeans = KMeans(
    n_clusters=n_clusters, # Définit le nombre de clusters (classes) à générer.
    random_state=0, # Contrôle l’aléatoire de l’initialisation des centroïdes pour garantir la reproductibilité des résultats (0 au début, none ensuite)
    n_init=10 # Nombre de relances de l'algorithme avec des centroïdes initiaux différents. Le meilleur résultat (en termes d’inertie) est retenu.
)

## Mise en forme du traitement
labels = kmeans.fit_predict(data) # entraîne le modèle et attribue à chaque pixel une classe (cluster).
result = labels.reshape(height, width) #  reconstitue une image 2D des classes, pour visualisation ou sauvegarde.

## Création explicite du dossier de sortie avant l'écriture :
output_dir = pathlib.Path('results')
output_dir.mkdir(parents=True, exist_ok=True)

profile.update(
    driver='GTiff',  # Obligatoire si l'extension est .tif
    dtype=rasterio.uint8,
    count=1,
    compress='lzw'   # Compression recommandée
)

with rasterio.open('results/h87_k-means_10cl.tif', 'w', **profile) as dst:
    dst.write(result.astype(rasterio.uint8), 1)
    print("Fichier écrit avec succès dans :", dst.name)

#import matplotlib.pyplot as plt

# Ouvrir le raster de classification
with rasterio.open('results/h87_k-means_10cl.tif') as src:
    classification = src.read(1)
    # Récupérer les informations spatiales
    transform = src.transform
    width = src.width
    height = src.height

# Calculer les coordonnées d'emprise
MIN_X = transform.c
MAX_Y = transform.f
RESOLUTION = transform.a
COLS = width
ROWS = height

# Affichage
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 avec k-means")
plt.xlabel('Coordonnée X (m)')
plt.ylabel('Coordonnée Y (m)')
plt.show()

#import rasterio
#import matplotlib.pyplot as plt

# Ouvrir les deux fichiers raster
file1 = 'results/h87_k-means_10cl.tif'
file2 = '../tdpy16/results/h87_fr_6cl.tif'

with rasterio.open(file1) as src1, rasterio.open(file2) as src2:
    img1 = src1.read(1)
    img2 = src2.read(1)
    transform1 = src1.transform
    transform2 = src2.transform
    width1, height1 = src1.width, src1.height
    width2, height2 = src2.width, src2.height

# Calcul des coordonnées d'emprise pour chaque image
MIN_X1 = transform1.c
MAX_Y1 = transform1.f
RESOLUTION1 = transform1.a
COLS1 = width1
ROWS1 = height1

MIN_X2 = transform2.c
MAX_Y2 = transform2.f
RESOLUTION2 = transform2.a
COLS2 = width2
ROWS2 = height2

# --- Configuration de la figure ---
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
plt.subplots_adjust(wspace=0.4)  # Espace supplémentaire pour la colorbar

# Affichage image 1
img1_extent = [MIN_X1, MIN_X1 + COLS1*RESOLUTION1, MAX_Y1 - ROWS1*RESOLUTION1, MAX_Y1]
im1 = axs[0].imshow(img1, cmap='tab10', extent=img1_extent)
axs[0].set_title('Classif. k-means en 10 classes')
axs[0].set_xlabel('Coordonnée X (m)')
axs[0].set_ylabel('Coordonnée Y (m)')

# Affichage image 2
img2_extent = [MIN_X2, MIN_X2 + COLS2*RESOLUTION2, MAX_Y2 - ROWS2*RESOLUTION2, MAX_Y2]
im2 = axs[1].imshow(img2, cmap='tab10', extent=img2_extent)
axs[1].set_title('Classif. random forest en 6 classes')
axs[1].set_xlabel('Coordonnée X (m)')
axs[1].set_ylabel('Coordonnée Y (m)')

# --- Positionnement précis de la colorbar entre les deux axes ---
pos_ax1 = axs[0].get_position()
pos_ax2 = axs[1].get_position()

cax = fig.add_axes([
    pos_ax1.x1 + 0.02,  # Début après le premier axe
    pos_ax1.y0,         # Même hauteur verticale
    0.02,               # Largeur de la colorbar
    pos_ax1.height      # Hauteur identique aux images
])

cbar = plt.colorbar(im1, cax=cax)
cbar.set_label('Classes', labelpad=2)
cbar.set_ticks(range(10))  # Toutes les classes (0-9)

plt.suptitle('Classifications non supervisée (k-means) et supervisée (random forest) dans les Howe Hills, MA, USA, 1987', y=0.98)
plt.show()

TD : Faire une classification non supervisée.¶

Cours SIG et télédétection¶

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

Téléchargement des documents nécessaires :

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

1.1 Chargement des données¶

1.2 Chargement des bibliothèques (library)¶

2. Lecture des bandes de l'image TM¶

3. Création du stack multibande¶

4. Configuration K-means¶

5. Classification¶

6. Sauvegarde du résultat¶

7. Affichage des résultats bruts¶

8. Affichage des classifications supervisée et non supervisée côte à côte¶