# ---------------------------------------------------
# TP Analyse d'une image Sentinel-2
# ---------------------------------------------------

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import sys

# -------------------
# Ouverture de l’image
# -------------------
fichier = "data/2025-07-25-00_Sentinel-2_L2A_B4B3B2B8.tiff"
ds = gdal.Open(fichier, gdal.GA_ReadOnly)

num_band = 1
band1 = ds.GetRasterBand(num_band)
array = band1.ReadAsArray()

# -------------------
# Étape 1 - Métadonnées et histogrammes
# -------------------
print("Taille (X,Y) :", ds.RasterXSize, ds.RasterYSize)
print("Résolution pixel (m) :", ds.GetGeoTransform()[1])
print("Système de coordonnées :", ds.GetProjection())
print("Type de données :", gdal.GetDataTypeName(band1.DataType))

plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.imshow(array, cmap="gray", vmin=10, vmax=100)
plt.title("Bande {} brute".format(num_band))
plt.colorbar(label="DN") # DN: Digital Number entre 0 et 255 en 8 bits

plt.subplot(1,2,2)
plt.hist(array.flatten(), bins=256, color='black', alpha=0.7)
plt.title("Histogramme - Bande {} brute".format(num_band))
plt.xlabel("Valeurs DN")
plt.ylabel("Fréquence")

# -------------------
# Étape 2 - Correction de contraste (étirement linéaire)
# -------------------
# On ignore les valeurs extrêmes (e.g. 2ème - 98ème percentiles)

def contrast_stretch(band):
    p2, p98 = np.percentile(band, (2, 98))
    band_stretch = np.clip((band - p2) / (p98 - p2), 0, 1)
    return band_stretch

array_stretch = contrast_stretch(array)

plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.imshow(array_stretch, cmap="gray")
plt.title("Bande {} etiree".format(num_band))
plt.colorbar(label="Valeurs normalisées [0-1]")

plt.subplot(1,2,2)
plt.hist(array_stretch.flatten(), bins = 100, color='black', alpha=0.7)
plt.title("Histogramme après correction de contraste")
plt.xlabel("Valeurs normalisées")
plt.ylabel("Fréquence")

# -------------------
# Étape 3 - Seuillage et Binarisation
# -------------------
# Hypothèse : eau = faibles valeurs radiométriques
# On choisit un seuil simple = moyenne - 0.4*écart-type

threshold = np.mean(array_stretch) - 0.4*np.std(array_stretch)
print(threshold)
binaire = (array_stretch > threshold).astype(np.uint8)

plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.imshow(array_stretch, cmap="gray")
plt.title("Bande 1 corrigée")

plt.subplot(1,2,2)
plt.imshow(binaire, cmap="gray")
plt.title(f"Seuillage (seuil = {threshold:.2f})\nEau=0, Terre=1")

# -------------------
# Étape 4 - Image miroir
# -------------------

def get_mirror(band, Ysize, direction="vertical"):

    if direction == "vertical":
        flip = np.flipud(band)
    elif direction == "horizontal":
        flip = np.fliplr(band)
    else:
        raise ValueError("Direction inconnue. Utilisez 'vertical' ou 'horizontal'.")

    mirror = np.concatenate((band[:int(Ysize / 2), :],
                             flip[int(Ysize / 2):, :]), axis=0)
    return mirror

mirror = get_mirror(array_stretch, ds.RasterYSize)

# --- affichage ---
fig, ax1 = plt.subplots(1, 1, figsize=(8,7))
ax1.imshow(mirror, cmap="gray")
ax1.set_title("Image miroir")
ax1.axis('off')

# -------------------
# Étape 5 - Composition Colorée
# -------------------

# --- lecture des bandes ---
red = ds.GetRasterBand(1).ReadAsArray().astype(float)    # B4
green = ds.GetRasterBand(2).ReadAsArray().astype(float)  # B3
blue = ds.GetRasterBand(3).ReadAsArray().astype(float)   # B2
nir = ds.GetRasterBand(4).ReadAsArray().astype(float)    # B8

# --- correction contraste pour chaque bande ---
red_stretch   = contrast_stretch(red)
green_stretch = contrast_stretch(green)
blue_stretch  = contrast_stretch(blue)
nir_stretch   = contrast_stretch(nir)

# --- composition colorée vraie couleur (RGB) ---
rgb = np.dstack((red_stretch, green_stretch, blue_stretch))
print(np.shape(rgb))

# --- composition colorée infrarouge (NIR, R, G) ---
false_color = np.dstack((nir_stretch, red_stretch, green_stretch))

# --- affichage final ---
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14,7))

ax1.imshow(rgb)
ax1.set_title("Composition colorée - Vraie couleur (corrigée gamma)")
ax1.axis('off')

ax2.imshow(false_color)
ax2.set_title("Composition colorée - Infrarouge (corrigée gamma)")
ax2.axis('off')

plt.show()