🔍 Explication des packages :
raster : permet de manipuler des données raster (images satellites, modèles numériques, etc.).
sf : utilisé pour gérer les données spatiales vectorielles selon les standards modernes (format simple features).
ggplot2 : l’un des outils les plus puissants pour créer des graphiques élégants et personnalisés.
terra : successeur du package raster, il offre des fonctions plus rapides et modernes pour les traitements raster. Ici, on l’utilise notamment pour trim(), qui supprime les bordures nulles d’un raster.
viridis : propose des palettes de couleurs adaptées aux graphiques scientifiques, avec une bonne lisibilité, même en impression ou pour les personnes daltoniennes.
# Chargement des bibliothèques nécessaires
library(raster)Loading required package: splibrary(sf)Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(ggplot2)
library(terra) # Pour trim()terra 1.8.42library(viridis)Loading required package: viridisLiteDans cette première étape, nous préparons les données géographiques nécessaires pour effectuer une analyse de la végétation à partir d’images satellites Sentinel-2. Tout d’abord, nous chargeons une couche vectorielle des communes (au format GeoPackage) en utilisant le package sf. Cette couche est initialement en système de coordonnées géographiques (EPSG:4326), c’est-à-dire avec des coordonnées latitude/longitude. Pour garantir une bonne compatibilité avec les rasters Sentinel-2, cette couche est reprojetée dans le système de coordonnées projetées EPSG:32631 (UTM zone 31N, utilisé couramment pour la France métropolitaine). Comme certaines fonctions du package raster nécessitent encore un format géospatial plus ancien, nous convertissons la couche sf en objet de type Spatial.
Ensuite, une fonction nommée charger_bandes est définie pour faciliter le chargement des deux bandes spectrales nécessaires à l’analyse : la bande B4 (rouge) et la bande B8 (proche infrarouge). Ces deux bandes sont indispensables pour le calcul du NDVI (Normalized Difference Vegetation Index), un indicateur essentiel pour évaluer la santé et la densité de la végétation. Grâce à cette fonction, on peut facilement charger les paires de bandes pour différentes dates ou années.
# Charger la couche des communes (EPSG:4326)
communes_sf <- st_read("/Users/admin/Documents/communes.gpkg")Reading layer `communes' from data source `/Users/admin/Documents/communes.gpkg' using driver `GPKG'
Simple feature collection with 7 features and 26 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 1.933051 ymin: 48.71395 xmax: 2.106251 ymax: 48.80044
Geodetic CRS: WGS 84# Reprojeter en EPSG:32631 (SCR des rasters Sentinel-2)
communes_sf_proj <- st_transform(communes_sf, crs = 32631)
# Convertir la couche sf reprojetée en Spatial (nécessaire pour 'mask')
communes_sp <- as(communes_sf_proj, "Spatial")
# Fonction pour charger les bandes d'une image Sentinel-2
charger_bandes <- function(b4_path, b8_path) {
list(
b4 = raster(b4_path),
b8 = raster(b8_path)
)
}# Charger les images par année
bandes_2015 <- charger_bandes(
b4 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20150716T105026_N0500_R051_T31UDQ_20231009T223128.SAFE/GRANULE/L2A_T31UDQ_A000334_20150716T105024/IMG_DATA/R10m/T31UDQ_20150716T105026_B04_10m.jp2",
b8 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20150716T105026_N0500_R051_T31UDQ_20231009T223128.SAFE/GRANULE/L2A_T31UDQ_A000334_20150716T105024/IMG_DATA/R10m/T31UDQ_20150716T105026_B08_10m.jp2"
)
bandes_2017 <- charger_bandes(
b4 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20170618T105621_N0500_R094_T31UDQ_20231012T193147.SAFE/GRANULE/L2A_T31UDQ_A010387_20170618T110415/IMG_DATA/R10m/T31UDQ_20170618T105621_B04_10m.jp2",
b8 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20170618T105621_N0500_R094_T31UDQ_20231012T193147.SAFE/GRANULE/L2A_T31UDQ_A010387_20170618T110415/IMG_DATA/R10m/T31UDQ_20170618T105621_B08_10m.jp2"
)
bandes_2019 <- charger_bandes(
b4 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2B_MSIL2A_20190514T105629_N0500_R094_T31UDQ_20221204T195224.SAFE/GRANULE/L2A_T31UDQ_A011417_20190514T105747/IMG_DATA/R10m/T31UDQ_20190514T105629_B04_10m.jp2",
b8 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2B_MSIL2A_20190514T105629_N0500_R094_T31UDQ_20221204T195224.SAFE/GRANULE/L2A_T31UDQ_A011417_20190514T105747/IMG_DATA/R10m/T31UDQ_20190514T105629_B08_10m.jp2"
)
bandes_2021 <- charger_bandes(
b4 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20210614T105031_N0500_R051_T31UDQ_20230317T045627.SAFE/GRANULE/L2A_T31UDQ_A031222_20210614T105443/IMG_DATA/R10m/T31UDQ_20210614T105031_B04_10m.jp2",
b8 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20210614T105031_N0500_R051_T31UDQ_20230317T045627.SAFE/GRANULE/L2A_T31UDQ_A031222_20210614T105443/IMG_DATA/R10m/T31UDQ_20210614T105031_B08_10m.jp2"
)
bandes_2023 <- charger_bandes(
b4 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20230528T105621_N0510_R094_T31UDQ_20240908T090134.SAFE/GRANULE/L2A_T31UDQ_A041418_20230528T105909/IMG_DATA/R10m/T31UDQ_20230528T105621_B04_10m.jp2",
b8 = "/Users/admin/Documents/ETUDES/Projet_IDGEO/vegetation/S2A_MSIL2A_20230528T105621_N0510_R094_T31UDQ_20240908T090134.SAFE/GRANULE/L2A_T31UDQ_A041418_20230528T105909/IMG_DATA/R10m/T31UDQ_20230528T105621_B08_10m.jp2"
)Cette section du script est consacrée au calcul du NDVI (Normalized Difference Vegetation Index), un indicateur couramment utilisé pour mesurer la vigueur de la végétation à partir des images satellites. Le NDVI est calculé à l’aide de deux bandes spectrales : la bande rouge (B4) et la bande proche infrarouge (B8). La fonction calc_ndvi applique la formule standard du NDVI :
NDVI=(B8−B4) / (B8+B4)
Pour garantir des résultats fiables, la fonction intègre deux filtres :
Elle remplace les valeurs infinies (résultant d’une division par zéro) par des valeurs manquantes (NA).
Elle élimine également les valeurs aberrantes qui ne se situeraient pas dans l’intervalle théorique du NDVI (entre -1 et 1).
Une fois la fonction définie, elle est utilisée pour calculer le NDVI sur les cinq années disponibles (2015, 2017, 2019, 2021 et 2023), en appliquant les bandes chargées précédemment.
Ensuite, un seuil est appliqué au NDVI pour identifier les zones fortement végétalisées. Le critère choisi ici est NDVI ≥ 0.4, ce qui correspond généralement à une couverture végétale dense ou saine. Une fonction seuillage est définie à cet effet et appliquée à chaque raster NDVI, produisant ainsi des rasters booléens (valeurs TRUE pour NDVI ≥ 0.4, FALSE ou NA ailleurs), utiles pour l’analyse spatiale et la quantification des surfaces végétalisées.
# Fonction de calcul NDVI avec correction de zones sans données
calc_ndvi <- function(b4, b8) {
ndvi <- (b8 - b4) / (b8 + b4)
ndvi[is.infinite(ndvi)] <- NA # éviter les divisions par zéro
ndvi[ndvi < -1 | ndvi > 1] <- NA # filtrer les valeurs aberrantes
return(ndvi)
}
# Calcul NDVI
ndvi_2015 <- calc_ndvi(bandes_2015$b4, bandes_2015$b8)
ndvi_2017 <- calc_ndvi(bandes_2017$b4, bandes_2017$b8)
ndvi_2019 <- calc_ndvi(bandes_2019$b4, bandes_2019$b8)
ndvi_2021 <- calc_ndvi(bandes_2021$b4, bandes_2021$b8)
ndvi_2023 <- calc_ndvi(bandes_2023$b4, bandes_2023$b8)
# Seuillage NDVI ≥ 0.4
seuillage <- function(ndvi) ndvi >= 0.4
seuil_2015 <- seuillage(ndvi_2015)
seuil_2017 <- seuillage(ndvi_2017)
seuil_2019 <- seuillage(ndvi_2019)
seuil_2021 <- seuillage(ndvi_2021)
seuil_2023 <- seuillage(ndvi_2023)Une fois les rasters NDVI seuillés obtenus (c’est-à-dire identifiant les zones avec NDVI ≥ 0.4), cette partie du code applique une découpe spatiale à l’aide de la couche des communes pour se focaliser uniquement sur l’emprise géographique d’intérêt. Cela se fait grâce à la fonction mask() qui conserve uniquement les pixels situés à l’intérieur des limites communales.
Chaque raster seuillé pour les années 2015 à 2023 est ainsi masqué avec la géométrie des communes (communes_sp), produisant cinq nouveaux objets (coupe_2015, coupe_2017, etc.).
Ensuite, on utilise la fonction trim() sur chacun de ces rasters. Cette opération permet de supprimer les marges nullesautour des données (zones sans information), ce qui facilite un affichage plus propre et centré.
Pour la visualisation, la fonction par(mfrow = c(2, 3), mar = c(3, 3, 3, 1)) divise la fenêtre graphique en 2 lignes et 3 colonnes (format grille) avec des marges ajustées. Cela permet d’afficher en parallèle les différentes cartes des zones végétalisées pour les années 2015, 2017, 2019, 2021 et 2023. Chaque plot() affiche l’évolution de la végétation avec en titre l’année correspondante, facilitant ainsi la lecture comparative visuelle de la dynamique de végétalisation sur le territoire.
# Découpe avec la couche des communes
coupe_2015 <- mask(seuil_2015, communes_sp)
coupe_2017 <- mask(seuil_2017, communes_sp)
coupe_2019 <- mask(seuil_2019, communes_sp)
coupe_2021 <- mask(seuil_2021, communes_sp)
coupe_2023 <- mask(seuil_2023, communes_sp)
# Trim pour enlever les marges nulles
trim2015 <- trim(coupe_2015)
trim2017 <- trim(coupe_2017)
trim2019 <- trim(coupe_2019)
trim2021 <- trim(coupe_2021)
trim2023 <- trim(coupe_2023)
# Ajustement marges graphiques
par(mfrow = c(2, 3), mar = c(3, 3, 3, 1))
# Affichage des images NDVI seuillées
plot(trim2015, main = "NDVI ≥ 0.4 - 2015")
plot(trim2017, main = "NDVI ≥ 0.4 - 2017")
plot(trim2019, main = "NDVI ≥ 0.4 - 2019")
plot(trim2021, main = "NDVI ≥ 0.4 - 2021")
plot(trim2023, main = "NDVI ≥ 0.4 - 2023")
🔲 Affichage final et réinitialisation graphique :
Après avoir affiché les 5 cartes NDVI seuillées, une case vide est ajoutée avec plot.new() et text() pour équilibrer la disposition en grille 2x3. Ce petit détail esthétique permet de compléter visuellement la dernière colonne restée vide. Ensuite, la configuration graphique est réinitialisée avec par(mfrow = c(1, 1)) pour revenir à un affichage par défaut.
📐 Calcul de surface végétalisée :
Une fonction personnalisée calc_surface() est ensuite définie. Elle permet de calculer la surface totale végétalisée à partir des rasters seuillés. Cela se fait en :
additionnant les valeurs des cellules valides (valeurs TRUE),
puis en les multipliant par la résolution des pixels (aire d’un pixel), grâce à prod(res(r)).
📊 Création du tableau statistique (data.frame) :
Les superficies végétalisées ainsi calculées pour chaque année sont stockées dans une table structurée (superficie_df)avec deux colonnes :
Annee : contenant les années 2015, 2017, 2019, 2021, 2023,
Superficie_Vegetation : correspondant à la surface végétalisée en mètres carrés pour chaque année.
# Dernière case vide
plot.new()
text(0.5, 0.5, "Fin des images", cex = 1.5)
# Réinitialisation graphique
par(mfrow = c(1, 1))
# Fonction de calcul de surface végétalisée
calc_surface <- function(r) {
cellStats(r, stat = "sum", na.rm = TRUE) * prod(res(r))
}
# DataFrame des superficies en m²
superficie_df <- data.frame(
Annee = c(2015, 2017, 2019, 2021, 2023),
Superficie_Vegetation = c(
calc_surface(coupe_2015),
calc_surface(coupe_2017),
calc_surface(coupe_2019),
calc_surface(coupe_2021),
calc_surface(coupe_2023)
)
)Le graphique ci-dessus illustre l’évolution de la superficie végétalisée au sein des communes étudiées, sur la période 2015-2023. La végétation est ici définie par un indice NDVI (Normalized Difference Vegetation Index) supérieur ou égal à 0,4, seuil couramment utilisé pour identifier une couverture végétale dense.
Chaque barre verticale représente la superficie totale en kilomètres carrés (km²) où le NDVI dépasse ce seuil, calculée à partir d’images satellite Sentinel-2 pour chaque année d’étude. Les valeurs sont affichées directement au-dessus des barres pour faciliter la lecture.
La coloration dégradée des barres, obtenue via la palette Viridis, met en évidence les différences de surface végétalisée d’une année sur l’autre, permettant d’identifier visuellement les tendances de croissance ou de déclin de la couverture végétale.
Ce graphique offre ainsi un outil visuel simple et efficace pour suivre les changements spatio-temporels de la végétation dans la zone d’étude, contribuant à mieux comprendre les dynamiques environnementales locales.
# Graphique évolution superficie végétalisée
# S'assurer que la colonne est bien numérique
superficie_df$Superficie_Vegetation <- as.numeric(superficie_df$Superficie_Vegetation)
# Graphique de l'évolution de la superficie végétalisée
ggplot(superficie_df, aes(
x = as.factor(Annee),
y = Superficie_Vegetation / 1e6,
fill = Superficie_Vegetation / 1e6 # Remplissage avec une variable numérique
)) +
geom_col(width = 0.6, show.legend = FALSE) +
geom_text(aes(label = paste0(round(Superficie_Vegetation / 1e6, 2), " km²")),
vjust = -0.5, size = 4) +
labs(
title = "Évolution de la superficie végétalisée (NDVI ≥ 0.4)",
x = "Année",
y = "Superficie (km²)"
) +
theme_minimal() +
scale_fill_viridis(option = "C")