Classification paysagère basée sur CLC+ Backbone 2018
Author
Eric Delmelle (pour Paris 8)
Published
December 15, 2025
1 Introduction
Ce document présente une analyse complète de la couverture du sol pour les communes où il y a eu un signalement de piqûres de tiques, basée sur les données CLC+ Backbone 2018 à une résolution de 10 mètres. Elle inclut :
Extraction des pourcentages par classe de couverture du sol
Analyse en Composantes Principales (ACP)
Classification K-means
Cartographie des composantes et des clusters
2 Classes CLC+ Backbone
Le référentiel CLC+ Backbone comprend 11 classes de couverture du sol :
Classe 1 : Surfaces artificielles
Classe 2 : Forêts de feuillus
Classe 3 : Forêts de conifères
Classe 4 : Forêts mixtes
Classe 5 : Végétation arbustive/herbacée naturelle
Classe 6 : Prairies
Classe 7 : Cultures
Classe 8 : Zones humides intérieures
Classe 9 : Zones humides côtières
Classe 10 : Plans d’eau
Classe 11 : Espaces marins
3 1. Configuration et chargement des données
3.1 1.1 Installation et chargement des packages
Code
### MODIFIE LOCALEMENT #### ATTENTION aux chemins des fichiers !!! #### MODIFIE LOCALEMENT #### ---------------------------------------------------------# 1. Installation et chargement des packages# ---------------------------------------------------------pkgs <-c("sf", "terra", "exactextractr", "FactoMineR", "cluster","factoextra", "dplyr", "tidyr", "ggplot2", "tmap","knitr", "patchwork")to_install <-setdiff(pkgs, rownames(installed.packages()))if (length(to_install)) {install.packages(to_install, repos ="https://cloud.r-project.org")}library(sf)library(terra)library(exactextractr)library(FactoMineR)library(factoextra)library(cluster)library(dplyr)library(tidyr)library(ggplot2)library(tmap)library(knitr)library(patchwork)
3.2 1.2 Chargement du raster CLC+ (EPSG:3035)
Le raster CLC+ est en projection EPSG:3035 (ETRS89 / LAEA Europe), optimisée pour l’analyse à l’échelle européenne.
Code
# chemin vers le raster### MODIFIE LOCALEMENT ###raster_path <-"C:/Users/nizar/Downloads/Projet Tutoré/CLMS_CLCplus_RASTER_2018_010m_eu_03035_V1_1.tif"### MODIFIE LOCALEMENT ###raster_clc <-rast(raster_path)# Vérifier le CRS actuelprint(crs(raster_clc, describe =TRUE))
name authority code area extent
1 ETRS89-extended / LAEA Europe EPSG 3035 <NA> NA, NA, NA, NA
Code
# Projection standard et non personnaliséecrs(raster_clc) <-"EPSG:3035"# verifcat("📊 Résolution:", res(raster_clc), "m\n")
# Lecture du shapefile des communes### MODIFIE LOCALEMENT ###communes <-st_read("C:/Users/nizar/Downloads/Projet Tutoré/communes piqure.shp", quiet =TRUE)### MODIFIE LOCALEMENT #### Vérifier le CRS actuelprint(st_crs(communes))
Coordinate Reference System:
User input: ETRS89-extended / LAEA Europe
wkt:
PROJCRS["ETRS89-extended / LAEA Europe",
BASEGEOGCRS["ETRS89",
ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
MEMBER["European Terrestrial Reference Frame 1989"],
MEMBER["European Terrestrial Reference Frame 1990"],
MEMBER["European Terrestrial Reference Frame 1991"],
MEMBER["European Terrestrial Reference Frame 1992"],
MEMBER["European Terrestrial Reference Frame 1993"],
MEMBER["European Terrestrial Reference Frame 1994"],
MEMBER["European Terrestrial Reference Frame 1996"],
MEMBER["European Terrestrial Reference Frame 1997"],
MEMBER["European Terrestrial Reference Frame 2000"],
MEMBER["European Terrestrial Reference Frame 2005"],
MEMBER["European Terrestrial Reference Frame 2014"],
MEMBER["European Terrestrial Reference Frame 2020"],
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[0.1]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["Europe Equal Area 2001",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",52,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",10,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",4321000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",3210000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (Y)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (X)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Statistical analysis."],
AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
BBOX[24.6,-35.58,84.73,44.83]],
ID["EPSG",3035]]
Code
# Transformer vers ETRS89 / LAEA Europe (EPSG:3035)communes <-st_transform(communes, 3035)# Vérifier le CRS après transformationprint(st_crs(communes))
Coordinate Reference System:
User input: EPSG:3035
wkt:
PROJCRS["ETRS89-extended / LAEA Europe",
BASEGEOGCRS["ETRS89",
ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
MEMBER["European Terrestrial Reference Frame 1989"],
MEMBER["European Terrestrial Reference Frame 1990"],
MEMBER["European Terrestrial Reference Frame 1991"],
MEMBER["European Terrestrial Reference Frame 1992"],
MEMBER["European Terrestrial Reference Frame 1993"],
MEMBER["European Terrestrial Reference Frame 1994"],
MEMBER["European Terrestrial Reference Frame 1996"],
MEMBER["European Terrestrial Reference Frame 1997"],
MEMBER["European Terrestrial Reference Frame 2000"],
MEMBER["European Terrestrial Reference Frame 2005"],
MEMBER["European Terrestrial Reference Frame 2014"],
MEMBER["European Terrestrial Reference Frame 2020"],
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[0.1]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["Europe Equal Area 2001",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",52,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",10,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",4321000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",3210000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (Y)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (X)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Statistical analysis."],
AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
BBOX[24.6,-35.58,84.73,44.83]],
ID["EPSG",3035]]
Code
cat("🏘️ Nombre de communes:", nrow(communes), "\n")
🏘️ Nombre de communes: 2039
3.4 1.4 Découpage du raster CLC+ aux seules communes avec signalements
Code
# Conversion en vecteur terrav_communes <-vect(communes)# réduire le raster à l'emprise des communesraster_clc_crop <-crop(raster_clc, v_communes)
# Vérification visuelleplot(raster_clc_mask, main ="Raster CLC+ restreint aux communes avec signalements")
4 2. Extraction des statistiques CLC+ par commune
(exactextractr : pourcentage de chaque classe CLC+)
Code
# Vérifier les CRSprint(crs(raster_clc_mask))
[1] "PROJCRS[\"ETRS89-extended / LAEA Europe\",\n BASEGEOGCRS[\"ETRS89\",\n ENSEMBLE[\"European Terrestrial Reference System 1989 ensemble\",\n MEMBER[\"European Terrestrial Reference Frame 1989\"],\n MEMBER[\"European Terrestrial Reference Frame 1990\"],\n MEMBER[\"European Terrestrial Reference Frame 1991\"],\n MEMBER[\"European Terrestrial Reference Frame 1992\"],\n MEMBER[\"European Terrestrial Reference Frame 1993\"],\n MEMBER[\"European Terrestrial Reference Frame 1994\"],\n MEMBER[\"European Terrestrial Reference Frame 1996\"],\n MEMBER[\"European Terrestrial Reference Frame 1997\"],\n MEMBER[\"European Terrestrial Reference Frame 2000\"],\n MEMBER[\"European Terrestrial Reference Frame 2005\"],\n MEMBER[\"European Terrestrial Reference Frame 2014\"],\n MEMBER[\"European Terrestrial Reference Frame 2020\"],\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[0.1]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4258]],\n CONVERSION[\"Europe Equal Area 2001\",\n METHOD[\"Lambert Azimuthal Equal Area\",\n ID[\"EPSG\",9820]],\n PARAMETER[\"Latitude of natural origin\",52,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",10,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"False easting\",4321000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",3210000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"northing (Y)\",north,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"easting (X)\",east,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Statistical analysis.\"],\n AREA[\"Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State.\"],\n BBOX[24.6,-35.58,84.73,44.83]],\n ID[\"EPSG\",3035]]"
Code
print(st_crs(communes))
Coordinate Reference System:
User input: EPSG:3035
wkt:
PROJCRS["ETRS89-extended / LAEA Europe",
BASEGEOGCRS["ETRS89",
ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
MEMBER["European Terrestrial Reference Frame 1989"],
MEMBER["European Terrestrial Reference Frame 1990"],
MEMBER["European Terrestrial Reference Frame 1991"],
MEMBER["European Terrestrial Reference Frame 1992"],
MEMBER["European Terrestrial Reference Frame 1993"],
MEMBER["European Terrestrial Reference Frame 1994"],
MEMBER["European Terrestrial Reference Frame 1996"],
MEMBER["European Terrestrial Reference Frame 1997"],
MEMBER["European Terrestrial Reference Frame 2000"],
MEMBER["European Terrestrial Reference Frame 2005"],
MEMBER["European Terrestrial Reference Frame 2014"],
MEMBER["European Terrestrial Reference Frame 2020"],
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[0.1]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["Europe Equal Area 2001",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",52,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",10,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",4321000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",3210000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (Y)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (X)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Statistical analysis."],
AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
BBOX[24.6,-35.58,84.73,44.83]],
ID["EPSG",3035]]
Code
results_df <- exactextractr::exact_extract( raster_clc_mask, communes,fun =function(values, coverage_fractions) {# création d'un dataframe valeurs x poids df <-data.frame(val = values, w = coverage_fractions)# Filtrer uniquement les classes CLC+ valides (1-11) df <- df[!is.na(df$val) & df$val %in%1:11, ]# Si aucune donnée, retourner des zérosif (nrow(df) ==0) return(rep(0, 11))# Agréger par classe et calculer les pourcentages agg <-aggregate(w ~ val, df, sum) pct <-numeric(11) pct[agg$val] <-100* agg$w /sum(agg$w)return(pct) },progress =FALSE)
5 3. Assemblage des résultats
Code
# Filtrer uniquement les classes CLC+ valides (1-11)if (is.matrix(results_df) &&nrow(results_df) ==11&&ncol(results_df) ==nrow(communes)) { results_df <-t(results_df)}# conversion en dataframeresults_df <-as.data.frame(results_df)colnames(results_df) <-paste0("class_", 1:11)rownames(results_df) <-NULL# vérif des dimension (optionnel)stopifnot(nrow(results_df) ==nrow(communes))# fusion avec les données communalescommunes_clc <- dplyr::bind_cols(communes, results_df)# aplatissement des col si nécessaire (format liste)class_cols <-paste0("class_", 1:11)for (cc in class_cols) { x <- communes_clc[[cc]]if (is.list(x)) x <-unlist(x, recursive =TRUE, use.names =FALSE) x <-suppressWarnings(as.numeric(x))if (length(x) <nrow(communes_clc)) { x <-rep_len(x, nrow(communes_clc)) } communes_clc[[cc]] <- x}cat("✅ Extraction terminée pour", nrow(communes_clc), "communes\n")
✅ Extraction terminée pour 2039 communes
5.1 3.2 Statistiques descriptives
Code
# aperçu des premières communes (tab)kable(head( communes_clc %>%st_drop_geometry() %>%select(NOM_COM, dplyr::starts_with("class_")),10 ),digits =2,caption ="Pourcentage de couverture par classe pour les 10 premières communes")
Pourcentage de couverture par classe pour les 10 premières communes
NOM_COM
class_1
class_2
class_3
class_4
class_5
class_6
class_7
class_8
class_9
class_10
class_11
Beaumont Saint-Cyr
5.81
7.05
27.94
0
2.82
21.27
32.17
0
0.15
2.78
0
Amboise
11.43
2.28
59.78
0
0.96
15.59
6.22
0
0.18
3.55
0
Athis
2.61
0.21
17.08
0
0.34
8.94
67.09
0
0.85
2.87
0
Chissay-en-Touraine
4.16
1.76
27.77
0
3.05
19.88
41.84
0
0.10
1.44
0
Château-l’Évêque
2.46
16.90
47.27
0
1.71
25.13
5.90
0
0.06
0.57
0
Saint-Nabord
8.40
30.56
32.78
0
0.18
24.47
2.30
0
0.36
0.94
0
Crusnes
10.08
1.43
23.86
0
0.22
19.66
44.73
0
0.02
0.00
0
Saint-Étienne-lès-Remiremont
6.07
55.67
18.96
0
0.15
16.97
1.42
0
0.10
0.66
0
Tucquegnieux
11.74
0.76
25.09
0
0.07
26.12
35.65
0
0.01
0.57
0
Bono
15.93
5.26
25.92
0
0.12
40.44
6.85
0
0.11
5.38
0
Code
# stat globales (moyenne par classe)cat("\n📈 Statistiques par classe (% moyen):\n")
Cette section présente la répartition spatiale de chaque classe de couverture du sol.
Code
# dictionnaire des noms de classesclc_labels <-c("class_1"="1. Surfaces arti.","class_2"="2. Feuillus","class_3"="3. Conifères","class_4"="4. Mixtes","class_5"="5. Arbust./herbacée","class_6"="6. Prairies","class_7"="7. Cultures","class_8"="8. Humides int.","class_9"="9. Humides côt.","class_10"="10. Plans d'eau","class_11"="11. Marins")# format long pour ggplotcommunes_long <- communes_clc %>%select(all_of(names(clc_labels)), NOM_COM, geometry) %>%pivot_longer(cols =all_of(names(clc_labels)),names_to ="classe", values_to ="pourcentage" ) %>%mutate(classe =factor(classe,levels =names(clc_labels),labels = clc_labels))# carte facettée des 11 classesp_classes <-ggplot(communes_long) +geom_sf(aes(fill = pourcentage), color =NA) +scale_fill_viridis_c(option ="magma", direction =-1, limits =c(0, 100) ) +facet_wrap(~classe, ncol =4) +labs(title ="Répartition spatiale des 11 classes CLC+ par commune",subtitle ="CLC+ Backbone 2018 (résolution 10 m) - Région parisienne",caption ="Projection : ETRS89 / LAEA Europe (EPSG:3035)" ) +theme_minimal(base_size =11) +theme(legend.position ="right",axis.text =element_blank(),axis.title =element_blank(),panel.grid =element_blank(),strip.text =element_text(face ="bold", size =9),plot.title =element_text(face ="bold", size =13),plot.subtitle =element_text(size =10) )print(p_classes)
Faire une carto dynamique ???
7 8. Analyse en Composantes Principales
Code
# extraction des données pour l'ACP (sans géométrie)data_pca <- communes_clc %>%st_drop_geometry() %>%select(all_of(class_cols))# ACP avec centrage et réductionpca_res <-PCA(data_pca, scale.unit =TRUE, graph =FALSE)# résumé de la variance expliquéecat("🔍 Variance expliquée par les 5 premiers axes (%):\n")
# Éboulis des valeurs propresp_scree <-fviz_eig( pca_res,addlabels =TRUE,barfill ="#2E86AB",barcolor ="#1B4F72",linecolor ="grey40",ncp =7,main ="Éboulis des valeurs propres",xlab ="Composante principale",ylab ="% de variance expliquée")print(p_scree)
Code
# Cercle des corrélationsp_circle <-fviz_pca_var( pca_res,col.var ="contrib",gradient.cols =c("#00AFBB", "#E7B800", "#FC4E07"),repel =TRUE,title ="Cercle des corrélations (axes 1 et 2)")print(p_circle)
Code
# Contributions aux axes 1 et 2p_contrib1 <-fviz_contrib(pca_res, choice ="var", axes =1, top =11) +labs(title ="Contributions à l'axe 1")p_contrib2 <-fviz_contrib(pca_res, choice ="var", axes =2, top =11) +labs(title ="Contributions à l'axe 2")print(p_contrib1 | p_contrib2)
Code
# Tableau des loadings (coordonnées des variables)loadings <-as.data.frame(round(pca_res$var$coord[, 1:5], 3))loadings$Variable <-rownames(loadings)loadings <- loadings[, c(6, 1:5)]kable( loadings,caption ="Corrélations entre les variables et les 5 premiers axes principaux",col.names =c("Variable", "Axe 1", "Axe 2", "Axe 3", "Axe 4", "Axe 5"))
Corrélations entre les variables et les 5 premiers axes principaux
Variable
Axe 1
Axe 2
Axe 3
Axe 4
Axe 5
class_1
class_1
0.100
0.369
-0.201
0.777
-0.070
class_2
class_2
0.448
-0.466
-0.078
-0.206
0.120
class_3
class_3
-0.122
-0.668
-0.334
0.146
0.405
class_4
class_4
0.460
0.399
-0.363
-0.375
0.052
class_5
class_5
0.551
0.370
-0.336
-0.359
0.020
class_6
class_6
0.032
-0.329
0.280
-0.239
-0.821
class_7
class_7
-0.648
0.503
0.374
-0.313
0.219
class_8
class_8
0.000
0.000
0.000
0.000
0.000
class_9
class_9
0.594
0.015
0.587
0.125
0.110
class_10
class_10
0.197
0.176
-0.084
0.428
-0.312
class_11
class_11
0.407
-0.009
0.656
0.120
0.325
7.0.1 interprétation des corrélations entre les classes CLC+ et les 5 premiers axes principaux :
7.1Axe 1 - Gradient d’artificialisation vs agricole
Fortement positif :
Classe 9 (Zones humides côtières) : 0.602
Classe 5 (Végétation arbustive/herbacée naturelle) : 0.582
Classe 4 (Forêts mixtes) : 0.490 → Ces milieux naturels structurent fortement l’axe 1
Fortement négatif :
Classe 7 (Cultures) : -0.594 → Opposition nette entre zones naturelles et agricoles
7.2Axe 2 - Gradient forestier vs prairial
Positif :
Classe 7 (Cultures) : 0.502
Classe 1 (Surfaces artificielles) : 0.409 → Pôle agricole et urbain
Négatif :
Classe 3 (Forêts de conifères) : -0.644
Classe 2 (Forêts de feuillus) : -0.508 → Pôle forestier très marqué
7.3Axe 3 - Gradient zones humides/marines
Très positif :
Classe 11 (Espaces marins) : 0.644
Classe 9 (Zones humides côtières) : 0.572 → Axe dédié aux écosystèmes aquatiques
Négatif :
Classe 3 (Forêts de conifères) : -0.351
Classe 4 (Forêts mixtes) : -0.340 → Opposition avec les milieux forestiers
7.4Axe 4 - Gradient urbain/naturel
Très positif :
Classe 1 (Surfaces artificielles) : 0.744
Classe 10 (Plans d’eau) : 0.441 → Axe dominé par l’artificialisation
Négatif :
Classe 4 (Forêts mixtes) : -0.379
Classe 5 (Végétation arbustive) : -0.375 → Opposition avec les milieux naturels
7.5Axe 5 - Spécialisation prairiale
Extrêmement négatif :
Classe 6 (Prairies) : -0.821 → Axe presque exclusivement défini par les prairies
Positif :
Classe 3 (Forêts de conifères) : 0.398
Classe 11 (Espaces marins) : 0.337 → Opposition forêts de conifères/espaces marins vs prairies
7.6Observations importantes :
Classe 8 (Zones humides intérieures) : Corrélations nulles sur tous les axes → Cette classe est mal représentée dans l’analyse ou très rare dans les communes étudiées
Structure paysagère : L’analyse révèle 5 dimensions indépendantes qui structurent le paysage :
Nature vs Agriculture
Forestier vs Agricole-Urbain
Aquatique vs Terrestre
Artificialisé vs Naturel
Prairial vs Conifères/Marin
Classes les plus discriminantes :
Prairies (Axe 5)
Surfaces artificielles (Axe 4)
Forêts de conifères (Axe 2)
Cultures (Axe 1 et 2)
8 9. Cartes des scores PCA
Code
# Ajout des scores PCA aux communescommunes_pca <-bind_cols( communes_clc,as.data.frame(pca_res$ind$coord[, 1:5]) %>%setNames(paste0("PCA", 1:5)))# Échelle commune pour comparaison entre cartesscore_min <-min( communes_pca$PCA1, communes_pca$PCA2, communes_pca$PCA3, communes_pca$PCA4, communes_pca$PCA5)score_max <-max( communes_pca$PCA1, communes_pca$PCA2, communes_pca$PCA3, communes_pca$PCA4, communes_pca$PCA5)# Fonction générique pour créer une carte PCAplot_pca_map <-function(data, var, title) {ggplot(data) +geom_sf(aes(fill = .data[[var]]), color ="grey70", size =0.15) +scale_fill_viridis_c(option ="plasma",direction =-1,limits =c(score_min, score_max),name ="Score" ) +labs(title = title) +theme_minimal(base_size =10) +theme(legend.position ="right",axis.text =element_blank(),axis.title =element_blank(),panel.grid =element_blank(),plot.title =element_text(face ="bold", size =11) )}p1 <-plot_pca_map(communes_pca, "PCA1","PCA1 — Gradient naturel vs agricole")p2 <-plot_pca_map(communes_pca, "PCA2","PCA2 — Gradient forestier vs agricole-urbain")p3 <-plot_pca_map(communes_pca, "PCA3","PCA3 — Gradient aquatique vs forestier")p4 <-plot_pca_map(communes_pca, "PCA4","PCA4 — Gradient artificialisé vs naturel")p5 <-plot_pca_map(communes_pca, "PCA5","PCA5 — Gradient prairial vs conifères/marin")print(p1); print(p2); print(p3); print(p4); print(p5)
9 10. Classification K-means
Code
# Données pour K-meansdata_kmeans <-as.data.frame(pca_res$ind$coord[, 1:5])colnames(data_kmeans) <-paste0("PCA", 1:5)# Méthode du coudep_elbow <-fviz_nbclust(data_kmeans, kmeans, method ="wss", k.max =10) +geom_vline(xintercept =5, linetype =2, color ="red") +labs(title ="Méthode du coude",x ="Nombre de clusters (k)",y ="Somme des carrés intra-cluster" ) +theme_minimal()# Méthode de la silhouettep_silhouette <-fviz_nbclust(data_kmeans, kmeans, method ="silhouette", k.max =10) +labs(title ="Analyse de la silhouette",x ="Nombre de clusters (k)" ) +theme_minimal()print(p_elbow | p_silhouette)
Code
# Application du K-means (k = 5)set.seed(1234)k_opt <-5km_res <-kmeans(data_kmeans, centers = k_opt, nstart =25)# Ajout des clusters aux donnéescommunes_km <- communes_pcacommunes_km$cluster <-factor(km_res$cluster)cat("🎯 Classification réalisée en", k_opt, "groupes\n")
💯 Qualité du clustering (between_SS / total_SS): 52.87 %
Code
# Profils moyens des clusters dans l'espace PCAcluster_profiles <-aggregate(data_kmeans, by =list(Cluster = km_res$cluster), mean)kable( cluster_profiles,digits =2,caption ="Profil moyen des clusters dans l'espace des 5 PCA")
Profil moyen des clusters dans l’espace des 5 PCA
Cluster
PCA1
PCA2
PCA3
PCA4
PCA5
1
0.43
0.91
-0.49
2.11
-0.34
2
0.29
-0.53
0.59
-0.34
-1.12
3
3.38
2.55
-2.12
-2.17
0.19
4
0.15
-1.14
-0.52
0.03
0.67
5
-1.09
0.87
0.57
-0.46
0.40
Code
# Statistiques descriptives par cluster (classes CLC+ originales)communes_km_summary <- communes_km %>%st_drop_geometry() %>%select(cluster, starts_with("class_")) %>%group_by(cluster) %>%summarise(across(starts_with("class_"), mean, .names ="{.col}"))kable( communes_km_summary,digits =1,caption ="Pourcentage moyen de chaque classe CLC+ par cluster")
Pourcentage moyen de chaque classe CLC+ par cluster
cluster
class_1
class_2
class_3
class_4
class_5
class_6
class_7
class_8
class_9
class_10
class_11
1
44.8
1.7
20.5
0.2
1.3
18.8
7.4
0
1.0
4.5
0.0
2
7.7
7.8
22.0
0.1
1.7
44.6
13.6
0
1.4
1.2
0.1
3
9.9
8.7
11.6
11.8
31.8
17.8
5.3
0
1.3
1.8
0.0
4
7.3
14.3
44.4
0.1
1.6
21.4
9.6
0
0.4
0.8
0.0
5
7.5
1.1
19.0
0.0
0.7
20.3
50.4
0
0.3
0.7
0.0
Code
# Carte des typologies paysagèresp_km_map <-ggplot(communes_km) +geom_sf(aes(fill = cluster), color ="grey70", size =0.15) +scale_fill_manual(values =c("#dbffe2", "#FF0014", "#F9FF00", "#18a843", "#095FF1"),name ="Cluster",labels =paste("Cluster", 1:5) ) +labs(title ="Typologie paysagère des communes",subtitle ="Classification K-means sur les 5 axes PCA (CLC+ Backbone 2018)",caption ="Projection : ETRS89 / LAEA Europe (EPSG:3035)" ) +theme_minimal(base_size =12) +theme(axis.text =element_blank(),axis.title =element_blank(),panel.grid =element_blank(),plot.title =element_text(face ="bold", size =14),legend.position ="right" )print(p_km_map)
Code
# Visualisation des clusters dans l'espace PCA (axes 1–2)p_km_biplot <-fviz_cluster( km_res,data = data_kmeans,geom ="point",ellipse.type ="norm",palette =c("#dbffe2", "#FF0014", "#F9FF00", "#18a843", "#095FF1"),ggtheme =theme_minimal(),main ="Clusters dans l'espace PCA (axes 1-2)")print(p_km_biplot)
---title: "Analyse de la Couverture du Sol des Communes"subtitle: "Classification paysagère basée sur CLC+ Backbone 2018"author: "Eric Delmelle (pour Paris 8)"date: todayformat: html: toc: true toc-location: left toc-depth: 3 number-sections: true code-fold: show code-tools: true theme: cosmo self-contained: trueeditor: visualexecute: warning: false message: false---# IntroductionCe document présente une analyse complète de la couverture du sol pour les communes où il y a eu un signalement de piqûres de tiques, basée sur les données **CLC+ Backbone 2018** à une résolution de 10 mètres. Elle inclut :1. Extraction des pourcentages par classe de couverture du sol\2. Analyse en Composantes Principales (ACP)\3. Classification K-means\4. Cartographie des composantes et des clusters------------------------------------------------------------------------# Classes CLC+ BackboneLe référentiel CLC+ Backbone comprend 11 classes de couverture du sol :- **Classe 1** : Surfaces artificielles- **Classe 2** : Forêts de feuillus- **Classe 3** : Forêts de conifères- **Classe 4** : Forêts mixtes- **Classe 5** : Végétation arbustive/herbacée naturelle- **Classe 6** : Prairies- **Classe 7** : Cultures- **Classe 8** : Zones humides intérieures- **Classe 9** : Zones humides côtières- **Classe 10** : Plans d’eau- **Classe 11** : Espaces marins------------------------------------------------------------------------# 1. Configuration et chargement des données## 1.1 Installation et chargement des packages```{r}### MODIFIE LOCALEMENT #### ATTENTION aux chemins des fichiers !!! #### MODIFIE LOCALEMENT #### ---------------------------------------------------------# 1. Installation et chargement des packages# ---------------------------------------------------------pkgs <-c("sf", "terra", "exactextractr", "FactoMineR", "cluster","factoextra", "dplyr", "tidyr", "ggplot2", "tmap","knitr", "patchwork")to_install <-setdiff(pkgs, rownames(installed.packages()))if (length(to_install)) {install.packages(to_install, repos ="https://cloud.r-project.org")}library(sf)library(terra)library(exactextractr)library(FactoMineR)library(factoextra)library(cluster)library(dplyr)library(tidyr)library(ggplot2)library(tmap)library(knitr)library(patchwork)```## 1.2 Chargement du raster CLC+ (EPSG:3035)Le raster CLC+ est en projection **EPSG:3035** (ETRS89 / LAEA Europe), optimisée pour l’analyse à l’échelle européenne.```{r}# chemin vers le raster### MODIFIE LOCALEMENT ###raster_path <-"C:/Users/nizar/Downloads/Projet Tutoré/CLMS_CLCplus_RASTER_2018_010m_eu_03035_V1_1.tif"### MODIFIE LOCALEMENT ###raster_clc <-rast(raster_path)# Vérifier le CRS actuelprint(crs(raster_clc, describe =TRUE))# Projection standard et non personnaliséecrs(raster_clc) <-"EPSG:3035"# verifcat("📊 Résolution:", res(raster_clc), "m\n")cat("🗺️ Projection:", crs(raster_clc, proj =TRUE), "\n")```## 1.3 Chargement des communes```{r}# Lecture du shapefile des communes### MODIFIE LOCALEMENT ###communes <-st_read("C:/Users/nizar/Downloads/Projet Tutoré/communes piqure.shp", quiet =TRUE)### MODIFIE LOCALEMENT #### Vérifier le CRS actuelprint(st_crs(communes))# Transformer vers ETRS89 / LAEA Europe (EPSG:3035)communes <-st_transform(communes, 3035)# Vérifier le CRS après transformationprint(st_crs(communes))cat("🏘️ Nombre de communes:", nrow(communes), "\n")```## 1.4 Découpage du raster CLC+ aux seules communes avec signalements```{r}# Conversion en vecteur terrav_communes <-vect(communes)# réduire le raster à l'emprise des communesraster_clc_crop <-crop(raster_clc, v_communes)# masquer pour obtenir seulement les pixels dans les polygonesraster_clc_mask <-mask(raster_clc_crop, v_communes)# Vérification visuelleplot(raster_clc_mask, main ="Raster CLC+ restreint aux communes avec signalements")```# 2. Extraction des statistiques CLC+ par commune``` (exactextractr : pourcentage de chaque classe CLC+)``````{r}# Vérifier les CRSprint(crs(raster_clc_mask))print(st_crs(communes))results_df <- exactextractr::exact_extract( raster_clc_mask, communes,fun =function(values, coverage_fractions) {# création d'un dataframe valeurs x poids df <-data.frame(val = values, w = coverage_fractions)# Filtrer uniquement les classes CLC+ valides (1-11) df <- df[!is.na(df$val) & df$val %in%1:11, ]# Si aucune donnée, retourner des zérosif (nrow(df) ==0) return(rep(0, 11))# Agréger par classe et calculer les pourcentages agg <-aggregate(w ~ val, df, sum) pct <-numeric(11) pct[agg$val] <-100* agg$w /sum(agg$w)return(pct) },progress =FALSE)```# 3. Assemblage des résultats```{r}# Filtrer uniquement les classes CLC+ valides (1-11)if (is.matrix(results_df) &&nrow(results_df) ==11&&ncol(results_df) ==nrow(communes)) { results_df <-t(results_df)}# conversion en dataframeresults_df <-as.data.frame(results_df)colnames(results_df) <-paste0("class_", 1:11)rownames(results_df) <-NULL# vérif des dimension (optionnel)stopifnot(nrow(results_df) ==nrow(communes))# fusion avec les données communalescommunes_clc <- dplyr::bind_cols(communes, results_df)# aplatissement des col si nécessaire (format liste)class_cols <-paste0("class_", 1:11)for (cc in class_cols) { x <- communes_clc[[cc]]if (is.list(x)) x <-unlist(x, recursive =TRUE, use.names =FALSE) x <-suppressWarnings(as.numeric(x))if (length(x) <nrow(communes_clc)) { x <-rep_len(x, nrow(communes_clc)) } communes_clc[[cc]] <- x}cat("✅ Extraction terminée pour", nrow(communes_clc), "communes\n")```## 3.2 Statistiques descriptives```{r}# aperçu des premières communes (tab)kable(head( communes_clc %>%st_drop_geometry() %>%select(NOM_COM, dplyr::starts_with("class_")),10 ),digits =2,caption ="Pourcentage de couverture par classe pour les 10 premières communes")# stat globales (moyenne par classe)cat("\n📈 Statistiques par classe (% moyen):\n")moyennes <- communes_clc %>%st_drop_geometry() %>%select(dplyr::starts_with("class_")) %>%summarise(across(everything(), mean, na.rm =TRUE))kable(t(moyennes), col.names ="Moyenne (%)", digits =2)```# 4. Cartographie des 11 classes CLC+Cette section présente la répartition spatiale de chaque classe de couverture du sol.```{r}# dictionnaire des noms de classesclc_labels <-c("class_1"="1. Surfaces arti.","class_2"="2. Feuillus","class_3"="3. Conifères","class_4"="4. Mixtes","class_5"="5. Arbust./herbacée","class_6"="6. Prairies","class_7"="7. Cultures","class_8"="8. Humides int.","class_9"="9. Humides côt.","class_10"="10. Plans d'eau","class_11"="11. Marins")# format long pour ggplotcommunes_long <- communes_clc %>%select(all_of(names(clc_labels)), NOM_COM, geometry) %>%pivot_longer(cols =all_of(names(clc_labels)),names_to ="classe", values_to ="pourcentage" ) %>%mutate(classe =factor(classe,levels =names(clc_labels),labels = clc_labels))# carte facettée des 11 classesp_classes <-ggplot(communes_long) +geom_sf(aes(fill = pourcentage), color =NA) +scale_fill_viridis_c(option ="magma", direction =-1, limits =c(0, 100) ) +facet_wrap(~classe, ncol =4) +labs(title ="Répartition spatiale des 11 classes CLC+ par commune",subtitle ="CLC+ Backbone 2018 (résolution 10 m) - Région parisienne",caption ="Projection : ETRS89 / LAEA Europe (EPSG:3035)" ) +theme_minimal(base_size =11) +theme(legend.position ="right",axis.text =element_blank(),axis.title =element_blank(),panel.grid =element_blank(),strip.text =element_text(face ="bold", size =9),plot.title =element_text(face ="bold", size =13),plot.subtitle =element_text(size =10) )print(p_classes)```Faire une carto dynamique ???# 8. Analyse en Composantes Principales```{r}# extraction des données pour l'ACP (sans géométrie)data_pca <- communes_clc %>%st_drop_geometry() %>%select(all_of(class_cols))# ACP avec centrage et réductionpca_res <-PCA(data_pca, scale.unit =TRUE, graph =FALSE)# résumé de la variance expliquéecat("🔍 Variance expliquée par les 5 premiers axes (%):\n")print(round(pca_res$eig[1:5, 2], 2))cat("\n📊 Variance cumulée (5 axes):",round(pca_res$eig[5, 3], 2), "%\n")# Éboulis des valeurs propresp_scree <-fviz_eig( pca_res,addlabels =TRUE,barfill ="#2E86AB",barcolor ="#1B4F72",linecolor ="grey40",ncp =7,main ="Éboulis des valeurs propres",xlab ="Composante principale",ylab ="% de variance expliquée")print(p_scree)# Cercle des corrélationsp_circle <-fviz_pca_var( pca_res,col.var ="contrib",gradient.cols =c("#00AFBB", "#E7B800", "#FC4E07"),repel =TRUE,title ="Cercle des corrélations (axes 1 et 2)")print(p_circle)# Contributions aux axes 1 et 2p_contrib1 <-fviz_contrib(pca_res, choice ="var", axes =1, top =11) +labs(title ="Contributions à l'axe 1")p_contrib2 <-fviz_contrib(pca_res, choice ="var", axes =2, top =11) +labs(title ="Contributions à l'axe 2")print(p_contrib1 | p_contrib2)# Tableau des loadings (coordonnées des variables)loadings <-as.data.frame(round(pca_res$var$coord[, 1:5], 3))loadings$Variable <-rownames(loadings)loadings <- loadings[, c(6, 1:5)]kable( loadings,caption ="Corrélations entre les variables et les 5 premiers axes principaux",col.names =c("Variable", "Axe 1", "Axe 2", "Axe 3", "Axe 4", "Axe 5"))```### interprétation des corrélations entre les classes CLC+ et les 5 premiers axes principaux :## **Axe 1 - Gradient d'artificialisation vs agricole**- **Fortement positif** : - Classe 9 (Zones humides côtières) : 0.602 - Classe 5 (Végétation arbustive/herbacée naturelle) : 0.582 - Classe 4 (Forêts mixtes) : 0.490\ *→ Ces milieux naturels structurent fortement l'axe 1*- **Fortement négatif** : - Classe 7 (Cultures) : -0.594\ *→ Opposition nette entre zones naturelles et agricoles*## **Axe 2 - Gradient forestier vs prairial**- **Positif** : - Classe 7 (Cultures) : 0.502 - Classe 1 (Surfaces artificielles) : 0.409\ *→ Pôle agricole et urbain*- **Négatif** : - Classe 3 (Forêts de conifères) : -0.644 - Classe 2 (Forêts de feuillus) : -0.508\ *→ Pôle forestier très marqué*## **Axe 3 - Gradient zones humides/marines**- **Très positif** : - Classe 11 (Espaces marins) : 0.644 - Classe 9 (Zones humides côtières) : 0.572\ *→ Axe dédié aux écosystèmes aquatiques*- **Négatif** : - Classe 3 (Forêts de conifères) : -0.351 - Classe 4 (Forêts mixtes) : -0.340\ *→ Opposition avec les milieux forestiers*## **Axe 4 - Gradient urbain/naturel**- **Très positif** : - Classe 1 (Surfaces artificielles) : 0.744 - Classe 10 (Plans d'eau) : 0.441\ *→ Axe dominé par l'artificialisation*- **Négatif** : - Classe 4 (Forêts mixtes) : -0.379 - Classe 5 (Végétation arbustive) : -0.375\ *→ Opposition avec les milieux naturels*## **Axe 5 - Spécialisation prairiale**- **Extrêmement négatif** : - Classe 6 (Prairies) : -0.821\ *→ Axe presque exclusivement défini par les prairies*- **Positif** : - Classe 3 (Forêts de conifères) : 0.398 - Classe 11 (Espaces marins) : 0.337\ *→ Opposition forêts de conifères/espaces marins vs prairies*## **Observations importantes :**1. **Classe 8 (Zones humides intérieures)** : Corrélations nulles sur tous les axes → Cette classe est mal représentée dans l'analyse ou très rare dans les communes étudiées2. **Structure paysagère** : L'analyse révèle 5 dimensions indépendantes qui structurent le paysage : - Nature vs Agriculture - Forestier vs Agricole-Urbain - Aquatique vs Terrestre - Artificialisé vs Naturel - Prairial vs Conifères/Marin3. **Classes les plus discriminantes** : - Prairies (Axe 5) - Surfaces artificielles (Axe 4) - Forêts de conifères (Axe 2) - Cultures (Axe 1 et 2)# 9. Cartes des scores PCA```{r}# Ajout des scores PCA aux communescommunes_pca <-bind_cols( communes_clc,as.data.frame(pca_res$ind$coord[, 1:5]) %>%setNames(paste0("PCA", 1:5)))# Échelle commune pour comparaison entre cartesscore_min <-min( communes_pca$PCA1, communes_pca$PCA2, communes_pca$PCA3, communes_pca$PCA4, communes_pca$PCA5)score_max <-max( communes_pca$PCA1, communes_pca$PCA2, communes_pca$PCA3, communes_pca$PCA4, communes_pca$PCA5)# Fonction générique pour créer une carte PCAplot_pca_map <-function(data, var, title) {ggplot(data) +geom_sf(aes(fill = .data[[var]]), color ="grey70", size =0.15) +scale_fill_viridis_c(option ="plasma",direction =-1,limits =c(score_min, score_max),name ="Score" ) +labs(title = title) +theme_minimal(base_size =10) +theme(legend.position ="right",axis.text =element_blank(),axis.title =element_blank(),panel.grid =element_blank(),plot.title =element_text(face ="bold", size =11) )}p1 <-plot_pca_map(communes_pca, "PCA1","PCA1 — Gradient naturel vs agricole")p2 <-plot_pca_map(communes_pca, "PCA2","PCA2 — Gradient forestier vs agricole-urbain")p3 <-plot_pca_map(communes_pca, "PCA3","PCA3 — Gradient aquatique vs forestier")p4 <-plot_pca_map(communes_pca, "PCA4","PCA4 — Gradient artificialisé vs naturel")p5 <-plot_pca_map(communes_pca, "PCA5","PCA5 — Gradient prairial vs conifères/marin")print(p1); print(p2); print(p3); print(p4); print(p5)```# 10. Classification K-means```{r}# Données pour K-meansdata_kmeans <-as.data.frame(pca_res$ind$coord[, 1:5])colnames(data_kmeans) <-paste0("PCA", 1:5)# Méthode du coudep_elbow <-fviz_nbclust(data_kmeans, kmeans, method ="wss", k.max =10) +geom_vline(xintercept =5, linetype =2, color ="red") +labs(title ="Méthode du coude",x ="Nombre de clusters (k)",y ="Somme des carrés intra-cluster" ) +theme_minimal()# Méthode de la silhouettep_silhouette <-fviz_nbclust(data_kmeans, kmeans, method ="silhouette", k.max =10) +labs(title ="Analyse de la silhouette",x ="Nombre de clusters (k)" ) +theme_minimal()print(p_elbow | p_silhouette)# Application du K-means (k = 5)set.seed(1234)k_opt <-5km_res <-kmeans(data_kmeans, centers = k_opt, nstart =25)# Ajout des clusters aux donnéescommunes_km <- communes_pcacommunes_km$cluster <-factor(km_res$cluster)cat("🎯 Classification réalisée en", k_opt, "groupes\n")cat("📊 Répartition des communes par cluster:\n")print(table(communes_km$cluster))cat("\n💯 Qualité du clustering (between_SS / total_SS):",round(km_res$betweenss / km_res$totss *100, 2), "%\n")# Profils moyens des clusters dans l'espace PCAcluster_profiles <-aggregate(data_kmeans, by =list(Cluster = km_res$cluster), mean)kable( cluster_profiles,digits =2,caption ="Profil moyen des clusters dans l'espace des 5 PCA")# Statistiques descriptives par cluster (classes CLC+ originales)communes_km_summary <- communes_km %>%st_drop_geometry() %>%select(cluster, starts_with("class_")) %>%group_by(cluster) %>%summarise(across(starts_with("class_"), mean, .names ="{.col}"))kable( communes_km_summary,digits =1,caption ="Pourcentage moyen de chaque classe CLC+ par cluster")# Carte des typologies paysagèresp_km_map <-ggplot(communes_km) +geom_sf(aes(fill = cluster), color ="grey70", size =0.15) +scale_fill_manual(values =c("#dbffe2", "#FF0014", "#F9FF00", "#18a843", "#095FF1"),name ="Cluster",labels =paste("Cluster", 1:5) ) +labs(title ="Typologie paysagère des communes",subtitle ="Classification K-means sur les 5 axes PCA (CLC+ Backbone 2018)",caption ="Projection : ETRS89 / LAEA Europe (EPSG:3035)" ) +theme_minimal(base_size =12) +theme(axis.text =element_blank(),axis.title =element_blank(),panel.grid =element_blank(),plot.title =element_text(face ="bold", size =14),legend.position ="right" )print(p_km_map)# Visualisation des clusters dans l'espace PCA (axes 1–2)p_km_biplot <-fviz_cluster( km_res,data = data_kmeans,geom ="point",ellipse.type ="norm",palette =c("#dbffe2", "#FF0014", "#F9FF00", "#18a843", "#095FF1"),ggtheme =theme_minimal(),main ="Clusters dans l'espace PCA (axes 1-2)")print(p_km_biplot)```# 11. Synthèse des typologies```{r}summary_clusters <- communes_km %>%st_drop_geometry() %>%group_by(cluster) %>%summarise(N =n(),Urban =mean(class_1, na.rm =TRUE),Foret =mean(class_2 + class_3 + class_4, na.rm =TRUE),Agriculture =mean(class_7, na.rm =TRUE),Prairie =mean(class_6, na.rm =TRUE),.groups ="drop" )kable( summary_clusters,digits =1,col.names =c("Cluster", "N communes", "% Urbain", "% Forêt", "% Agriculture", "% Prairie"),caption ="Synthèse des typologies paysagères")```# 12. Export des données (optionnel)```{r}# Export du shapefile avec clusters# (décommenter si souhaité)#st_write(communes_km, "communes_typologie_paysagere.shp", delete_dsn = TRUE)# Export des scores PCA + clusters en CSV#write.csv(# communes_km %>%# st_drop_geometry() %>%# select(INSEE_COM, NOM_COM, cluster,# starts_with("PCA"), starts_with("class_")),# "communes_scores_pca_clusters.csv",# row.names = FALSE# )```# 13. Informations de session```{r}sessionInfo()```#
