2.1 Techniques basées sur des fenêtres glissantes

La plupart des techniques historiques de caractérisation du relief sont basées d’une part sur des calculs morphométriques (pente, orientation, exposition, profil) et d’autre part sur des analyses de voisinage liées au le concept de fenêtres glissantes. Ces méthodes de classification automatique du relief peuvent reposer sur la valeur locale de la pente et de la courbure (convexité) du terrain (Dikau, 1991 ; Wood, 1996) ou en comparant ces valeurs au voisinage.

L’analyse de voisinage utilise la technique des fenêtres glissantes (Figure 1). La valeur de chaque pixel est comparée aux valeurs des pixels adjacents. Selon la taille de la fenêtre glissante définie par l’utilisateur, les pixels adjacents peuvent être ceux touchant directement le pixel initial (fenêtre 3×3) ou tous les pixels situés à moins de deux pixels du pixel initial (fenêtre 5×5).

Ces méthodes ont l’avantage d’être relativement simples à mettre en place, mais ont l’inconvénient d’être gourmandes en ressources. De plus, elles sont très contraintes par la résolution des MNT d’entrée et par la taille choisie des fenêtres glissantes. Ainsi, pour un même terrain, le résultat pourra être très différent selon la résolution du MNT d’entrée et selon le choix de la taille de fenêtre glissante sélectionnée par l’utilisateur. Cette approche présente aussi le désavantage de mettre en évidence les petits reliefs, de l’ordre de taille de la fenêtre glissante, au détriment des grandes formes de relief du paysage.

2.2 Techniques par identification de geomorphons

L’algorithme Geomorphons, pour geomorphologic phonotypes, est une méthode qui identifie les formes du relief en utilisant un schéma de reconnaissance développé par Jasiewicz et Stepinski (2013). Ces formes de relief sont appelées geomorphons par les auteurs. L’avantage de cette méthode est qu’elle n’est pas figée sur une valeur donnée de fenêtre glissante mais adapte la distance à prendre en compte grâce au concept de ligne de mire.

Ici, nous appliquerons cette méthode à nos trois bassins-versants. Un des premiers usages de l’algorithme est la création de la carte géomorphologique de la Pologne (Jasiewicz et Stepinski, 2013). Il est aussi utilisé pour décrire le relief du Monténégro (Frankl et al., 2016). L’algorithme permet aussi d’identifier certaines formes spécifiques du paysage, comme des drumlins (Sarasan et al., 2018) ou des surfaces plates (Veselsky et al., 2015). Bandura (2016) l’a aussi utilisé pour identifier et délimiter des sommets dans les montagnes slovaques. Cette méthode a également été mobilisée pour détecter des objets à une échelle plus fine : des formes à la surface d’un glacier (Gawrysiak, 2020), des traces de barrages organisés par des castors (Swift, 2021) ou encore la détection de formes fluviales côtières (Gioia, 2021).

Les geomorphons sont déterminés par une méthode de classification semi-automatique du MNT. Cette méthode ne prend pas en compte les valeurs des pixels voisins mais s’inspire du principe de “ligne de mire”. Elle est flexible et s’adapte à différents terrains car le critère de la ligne de mire est plus souple que le critère de la fenêtre glissante fixe que l’on trouve dans les autres méthodes de classification. Elle adapte la distance à prendre en compte selon la notion de terrain openess. Il s’agit de calculer le rapport entre l’angle au zénith (Φ : zenith angle) et l’angle au nadir (Ψ : nadir angle) au sein d’un rayon de recherche défini par l’utilisateur (lookup distance). Le rayon est égal à x fois la largeur d’un pixel. Par exemple, choisir une lookup distance de 10 pour un MNT d’une résolution de 5 mètres revient à choisir une lookup distance de 50 mètres. De plus, le paramètre threshold angle permet de définir le niveau minimum de planéité souhaité permettant de gommer les aspérités microtopographiques.

Cette méthode compare chaque pixel du MNT aux huit valeurs des voisins visibles les plus lointains et au maximum jusqu’à la lookup distance définie par l’utilisateur. La figure 2 présente les 10 résultats types d’une telle analyse. Les altitudes équivalentes à la valeur de la cellule centrale sont représentées par un point vert, celles plus basses en bleu et celles plus hautes en rouge.

Ces 10 geomorphons types sont identifiables sur les rasters produits par l’algorithme (Figure 3).

4.1 L’algorithme Geomorphons

L’algorithme Geomorphons est disponible dans QGIS via les logiciels GRASS GIS, SAGA GIS ou encore WhiteboxTools.

Nous avons décidé d’employer le module proposé par GRASS GIS. Pour ce faire, nous exploitons la librairie QGIS process qui permet d’utiliser facilement des algorithmes de QGIS depuis R. Ainsi, Geomorphons est appelé dans R via la librairie qgisprocess.

La première chose à faire est d’installer cette librairie. Elle n’est pas disponible sur les dépôts officiels, il faut donc le faire depuis son dépôt sur GitHub. Pour ce faire, il faut suivre la démarche proposée ici.

if (!require("remotes")) install.packages("remotes")
 remotes::install_github("r-spatial/qgisprocess")

Pour vérifier sa bonne installation, il est possible d’afficher la configuration de l’installation, comme présenté ci-dessous.

library(qgisprocess)
qgis_configure()

# Pour pouvoir utiliser l'algorithme r.geomorphon :
# Activer le plugin grassprovider qui nous permettra de l'appeler dans R. 
# Pour ce faire utiliser la fonction :
qgis_enable_plugins()

4.2 Autres packages nécessaires

Six packages sont requis pour reproduire les analyses présentées dans cette fiche :

• sf pour la manipulation de couches vecteurs;
• terra pour la manipulation de rasters;
• tmap pour la création de cartes;
• tidyverse pour la manipulation de dataframe;
• reshape2 pour modifier la structure de dataframe;
• ggplot2 pour les graphiques.

Le code ci-dessous vérifie que toutes les librairies requises sont bien installées. Si ce n’est pas le cas, les installer puis les charger dans la session de travail.

# Noms de packages nécessaires
my_packages <- c("sf", "terra", "tmap", "tidyverse", "reshape2", "ggplot2")

# Vérifier si ces packages sont installés
missing_packages <- my_packages[!(my_packages %in% installed.packages()[,"Package"])]

# Installation des packages manquants depuis le CRAN
if(length(missing_packages)) install.packages(missing_packages, 
                                              repos = "http://cran.us.r-project.org")

# Chargement des packages nécessaires
lapply(my_packages, library, character.only = TRUE)

4.3 Données mobilisées

Dans ce projet, nous avons besoin de deux types de données :

• les limites des bassins-versants
• les MNT des bassins-versants

Télécharger les données

Pour rejouer l’analyse, décompressez les données contenues dans l’archive data.zip

5.1 Préparation des données

# Chargement de la couche des bassins de toute la France que nous avons téléchargé au format gpkg
bv <- st_read("data/BassinVersantTopographique_FXX.gpkg")

# Extraction du bassin du Ribérot selon son nom
riberot <- bv %>% filter(bv$TopoOH == "Le Ribérot")
# Export du bassin du Ribérot en une couche autonome (facultatif)
st_write(riberot, "data_processed/bv_riberot.gpkg", delete_layer = TRUE)

# Extraction du bassin de la Séveraisse
severaisse <- bv %>% filter(bv$TopoOH == "La Séveraisse de sa source au torrent de Navette inclus")

# Extraction du bassin de la Tarentaine 
tarentaine <- bv %>% filter(bv$TopoOH == "La Tarentaine (Trentaine) de sa source au confluent de l'Eau Verte")

# Enregistrement de deux couches géographiques en format GeoPackage
st_write(severaisse, "data_processed/bv_severaisse.gpkg", delete_layer = TRUE)
st_write(tarentaine, "data_processed/bv_tarentaine.gpkg", delete_layer = TRUE)

library(qgisprocess)

# Pour le département de l'Ariège (09)

# Indication du répertoire contenant les dalles de MNT de l'Ariège
path_dem_09 <- "/data/RGE_5m_Dep09/"

# Nous transformons ce chemin en chemin absolu
path_dem_09 <- paste(getwd(),"/data/RGE_5m_Dep09/", sep="")

# Nous définissons le nom et le chemin du raster virtuel qui sera produit
dem_09 <- "data_processed/RGE_5m_Dep09.vrt"

# Nous construisons la liste de toutes les dalles de MNT du département
ras_lst <- list.files(path_dem_09, full.names =TRUE, pattern = ".asc$", recursive=TRUE)

# Conversion vers le type qgis_list_input
ras_lst <- qgis_list_input(!!! ras_lst)

# Nous construisons le raster virtuel de l'Ariège via QGIS process
qgis_run_algorithm("gdal:buildvirtualraster",
                INPUT=ras_lst,
                SEPARATE=FALSE,
                ASSIGN_CRS='EPSG:2154',
                OUTPUT=dem_09)

# Pour le département des Hautes-Alpes (05)

# Indication du répertoire contenant les dalles de MNT des Hautes-Alpes
path_dem_05 <- "/data/RGE_5m_Dep05/"

# Transformation de ce chemin en chemin absolu
path_dem_05 <- paste(getwd(),"/data/RGE_5m_Dep05/", sep="")

# Définition du nom et du chemin du raster virtuel qui sera produit
dem_05 <- "data_processed/RGE_5m_Dep05.vrt"

# Construction de la liste de toutes les dalles de MNT du département
ras_lst <- list.files(path_dem_05, full.names = TRUE, pattern = ".asc$", recursive=TRUE)

# Conversion vers le type qgis_list_input
ras_lst <- qgis_list_input(!!! ras_lst)

# Construction du raster virtuel des Hautes-Alpes via QGIS process
qgis_run_algorithm("gdal:buildvirtualraster",
                INPUT=ras_lst,
                SEPARATE=FALSE,
                ASSIGN_CRS='EPSG:2154',
                OUTPUT=dem_05)

# Pour le département du Puy-de-Dôme (63)

# Indication du répertoire contenant les dalles de MNT du Puy-de-Dôme
path_dem_63 <- "/data/RGE_5m_Dep63/"

# Transformation de ce chemin en chemin absolu
path_dem_63 <- paste(getwd(),"/data/RGE_5m_Dep63/", sep="")

# Définition du nom et du chemin du raster virtuel qui sera produit
dem_63 <- "data_processed/RGE_5m_Dep63.vrt"

# Construction de la liste de toutes les dalles de MNT du département
ras_lst <- list.files(path_dem_63, full.names = TRUE, pattern = ".asc$", recursive=TRUE)

# Conversion vers le type qgis_list_input
ras_lst <- qgis_list_input(!!! ras_lst)

# Construction du raster virtuel du Puy-de-Dôme via QGIS process
qgis_run_algorithm("gdal:buildvirtualraster",
                INPUT=ras_lst,
                SEPARATE=FALSE,
                ASSIGN_CRS='EPSG:2154',
                OUTPUT=dem_63)

5.2 Découpage des MNT selon l’emprise des bassins-versants

Nous avons maintenant les limites de nos trois bassins, récupérées grâce à la BD TOPAGE, ainsi que les MNT de leurs départements d’appartenance. Le but est maintenant d’extraire les MNT de nos bassins et seulement de nos bassins.

# Si les limites des bassins ne sont plus chargées dans le projet, nous pouvons les recharger
# en réactivant les trois lignes suivantes (c'est-à-dire en supprimant le #) :
# Riberot <- st_read("data_processed/bv_riberot.gpkg")
# Severaisse <- st_read("data_processed/bv_severaisse.gpkg")
# Tarentaine <- st_read("data_processed/bv_tarentaine.gpkg")

# Chargement du MNT virtuel de l'Ariège
dem_09 <- rast("data_processed/RGE_5m_Dep09.vrt")

# Conversion des limites du Ribérot en un objet SpatVector utilisable par terra
riberot <- vect(riberot)

# Découpage de ce MNT selon l'étendue spatiale des limites du Ribérot
dem_riberot <- crop(dem_09, riberot)

# Masque du MNT selon les limites du Ribérot
dem_riberot <- mask(dem_riberot, riberot)

# Export du MNT (facultatif)
terra::writeRaster(dem_riberot,
      "data_processed/dem_Riberot.tif",
      filetype="GTiff", 
      overwrite=TRUE)

# Chargement du MNT virtuel des Hautes-Alpes
dem_05 <- rast("data_processed/RGE_5m_Dep05.vrt")

# Conversion des limites de la Séveraisse en un objet SpatVector utilisable par terra
severaisse <- vect(severaisse)

# Découpage de ce MNT selon l'étendue spatiale des limites de la Severaisse
dem_severaisse <- crop(dem_05, severaisse)

# Masque du MNT selon les limites de la Severaisse
dem_severaisse <- mask(dem_severaisse, severaisse)

# Export du MNT (facultatif)
terra::writeRaster(dem_severaisse,
      "data_processed/dem_Severaisse.tif",
      filetype="GTiff",
      overwrite=TRUE)

# Chargement du MNT virtuel du Puy-de-Dôme
dem_63 <- rast("data_processed/RGE_5m_Dep63.vrt")

# Conversion des limites de la Tarentaine en un objet SpatVector utilisable par terra
tarentaine <- vect(tarentaine)

# Découpage de ce MNT selon l'étendue spatiale des limites de la Tarentaine
dem_tarentaine <- crop(dem_63, tarentaine)

# Masque du MNT selon les limites de la Tarentaine
dem_tarentaine <- mask(dem_tarentaine, tarentaine)

# Export du MNT (facultatif)
terra::writeRaster(dem_tarentaine,
      "data_processed/dem_Tarentaine.tif",
      filetype="GTiff",
      overwrite=TRUE)

Nous disposons maintenant des MNT découpés selon chacun de nos bassins-versants. Nous pouvons passer au traitement des geomorphons.

5.3 Le traitement Geomorphons : détection de 10 formes de relief

À partir du MNT de chaque bassin-versant, nous pouvons lancer l’analyse des formes de relief grâce à l’algorithme Geomorphons. Nous utilisons son implémentation dans GRASS GIS à l’aide de la librairie qgisprocess. Nous obtenons le raster des 10 formes de relief telles que répertoriées par Geomorphons pour chacun de nos bassins.

L’utilisateur doit spécifier une distance de recherche, c’est-à-dire la distance maximale sur laquelle une forme de relief doit être recherchée. Ici, nous pouvons la spécifier à 20 pixels, ce qui équivaut à une distance de recherche de 100 mètres, car un pixel mesure 5 mètres de côté.

# Analyse des geomorphons du Ribérot
# Définition du chemin et du nom du raster de geomorphons en sortie
geom_riberot <- "data_processed/geomorphons_Riberot.tif"

qgis_run_algorithm("grass7:r.geomorphon",
                  elevation = dem_riberot,
                  search = 20, # lookup distance/rayon de recherche
                  flat = 1,
                  forms = geom_riberot)


# Analyse des geomorphons de la Séveraisse
# Définition du chemin et du nom du raster de geomorphons en sortie
geom_severaisse <- "data_processed/geomorphons_Severaisse.tif"

qgis_run_algorithm("grass7:r.geomorphon",
                  elevation = dem_severaisse,
                  search = 20, # lookup distance/rayon de recherche
                  flat = 1,
                  forms = geom_severaisse)

# Analyse des geomorphons de la Tarentaine
# Définition du chemin et du nom du raster de geomorphons en sortie
geom_tarentaine <- "data_processed/geomorphons_Tarentaine.tif"

qgis_run_algorithm("grass7:r.geomorphon",
                  elevation = dem_tarentaine,
                  search = 20, # lookup distance/rayon de recherche
                  flat = 1,
                  forms = geom_tarentaine)

Nous disposons à ce niveau de trois rasters présentant les 10 formes de relief pour chacun de nos bassins-versants. Selon l’algorithme Geomorphons, ces 10 classes correspondent au tableau suivant.

6.1 Représentation cartographique

Nous pouvons représenter la répartition géographique des formes de relief de nos trois bassins en les cartographiant. Pour ce faire nous commençons par placer les noms et les couleurs des formes du relief disponibles dans geomorphons dans un vecteur noms pour les noms et couleurs pour les couleurs.

noms = c("", "flat (plan)", "peak/summit (sommet)", "ridge (crête)", "shoulder (haut de versant", 
         "spur (eperon)", "slope (pente)", "hollow (creux de versant)", "footslope (bas de versant)", 
         "valley (fond de vallée)", "pit/depression (dépression)")

couleurs = c("#ffffffff",
             "#dcdcdc", 
             "#380000", 
             "#c80000", 
             "#ff5014", 
             "#fad23c", 
             "#ffff3c", 
             "#b4e614", 
             "#3cfa96", 
             "#0000ff", 
             "#000038")

Puis nous utilisons les lignes de commande suivantes :

library(terra)
library(tmap)

# Cartographie du Ribérot
# Chargement du raster des geomorphons du Ribérot
geom_riberot <- rast("data_processed/geomorphons_Riberot.tif")

tm_shape(geom_riberot) + 
  tm_raster(style="cat",
            labels=noms,
            palette=couleurs,
            title="Geomorphons du Ribérot") + 
  tm_compass(type="4star", size=1, position=c("right", "top")) +
  tm_scale_bar(breaks=1) +
  tm_layout(legend.outside=TRUE)

# Cartographie de la Séveraisse
# Chargement du raster des geomorphons de la Séveraisse
geom_severaisse <- rast("data_processed/geomorphons_Severaisse.tif")

tm_shape(geom_severaisse) + 
  tm_raster(style="cat",
            labels=noms,
            palette=couleurs,
            title="Geomorphons de la Séveraisse") + 
  tm_compass(type="4star", size=1, position=c("right", "top")) +
  tm_scale_bar(breaks=1) +
  tm_layout(legend.outside=TRUE)

# Cartographie de la Tarentaine
# Chargement du raster des geomorphons de la Tarentaine
geom_tarentaine <- rast("data_processed/geomorphons_Tarentaine.tif")

tm_shape(geom_tarentaine) + 
  tm_raster(style="cat",
            labels=noms,
            palette=couleurs,
            title="Geomorphons de la Tarentaine") + 
  tm_compass(type="4star", size=1, position=c("right", "top")) +
  tm_scale_bar(breaks=1) +
  tm_layout(legend.outside=TRUE)

Sur les trois bassins, la forme de relief dominante est clairement la “Pente”, ce qui est tout à fait logique au vu du contexte montagnard de ces bassins. Ensuite, les grandes “Vallées” apparaissent clairement et laissent deviner le réseau hydrographique des bassins. Les grandes vallées sont séparées les unes des autres par des “Lignes de crête” et des “Éperons”. Entre les lignes de crête et les fonds de vallées, outre les pentes, nous retrouvons fréquemment des “Creux de versant”. Les autres classes sont, quant à elles, peu représentées.

6.2 Comparaison quantitative des formes de relief

L’idée est maintenant de regarder si ces trois bassins montagnards partagent des traits de reliefs en commun ou s’ils sont très différents. Nous pouvons également nous demander si les deux bassins de haute montagne (Ribérot et Séveraisse) se démarquent par rapport à celui de moyenne montagne (Tarentaine).

Pour cette comparaison, nous calculons la superficie relative occupée par chaque forme de relief dans nos trois bassins. Nous utilisons la fonction freq() de la librairie terra, qui permet de calculer le nombre de pixels de chaque classe d’un raster et de récupérer ce résultat dans un dataframe.

Nous commençons par calculer le nombre de pixels de chaque classe pour le raster de geomorphons du Ribérot.

relief_riberot <- freq(geom_riberot)

# Affichage de la table
knitr::kable(relief_riberot, format = "html")

Nous obtenons bien un dataframe à trois colonnes. La première colonne (layer) est utile seulement si nous travaillons sur un raster multibande, ce qui n’est pas le cas ici. La deuxième colonne (value) contient le nom des classes et la troisième (count) le nombre de pixels de chaque classe. Nous allons créer une nouvelle colonne nommée simplement Ribérot dans laquelle nous allons calculer la proportion de chaque classe, ce qui nous permettra de faire des comparaisons entre les bassins.

library(tidyverse)
relief_riberot$Ribérot <- relief_riberot$count / sum(relief_riberot$count)

# Nous en profitons pour enlever les colonnes inutiles
relief_riberot <- relief_riberot %>% select(-one_of('layer', 'count'))

# Affichage de la table
knitr::kable(relief_riberot, format = "html")

Nous procédons de même pour les deux autres bassins-versants.

# Pour la Séveraisse
relief_severaisse <- freq(geom_severaisse)
relief_severaisse$Séveraisse <- relief_severaisse$count / sum(relief_severaisse$count)
relief_severaisse <- relief_severaisse %>% select(-one_of('layer', 'count')) 

# Affichage de la table
knitr::kable(relief_severaisse, format = "html")

# Pour la Tarentaine
relief_tarentaine <- freq(geom_tarentaine)
relief_tarentaine$Tarentaine <- relief_tarentaine$count / sum(relief_tarentaine$count)
relief_tarentaine <- relief_tarentaine %>% select(-one_of('layer', 'count'))

# Affichage de la table
knitr::kable(relief_tarentaine, format = "html")

Nous assemblons maintenant ces trois dataframes en un seul via une jointure attributaire sur le champ value. Nous pouvons ainsi comparer les proportions de chaque classe de relief dans nos bassins.

# Création d'une liste contenant les trois dataframes
relief_bassins <- list(relief_riberot, relief_severaisse, relief_tarentaine)

# Grâce à la fonction "reduce" de la librairie tidyverse
# Jointure attributaire sur les dataframes en une seule fois
relief_bassins <- relief_bassins %>% reduce(full_join, by='value')

# Nous renommons la colonne "value" en "relief"
relief_bassins <- rename(relief_bassins, 'relief' = 'value')

# Affichage de la table
knitr::kable(relief_bassins, format = "html")

À partir de ce dataframe groupant les proportions de formes de reliefs des trois bassins, nous traçons des graphiques en bâtons de ces proportions. Nous utilisons la librairie ggplot2. Cette librairie prend en entrée les dataframes au format long, c’est-à-dire des dataframes où les données sont organisées en lignes plutôt qu’en colonnes. Nous commençons par transformer le dataframe des formes de relief dans ce format grâce à la librairie reshape2 et à sa fonction melt.

relief_bassins_long <- melt(relief_bassins,
                            id.vars = c('relief'))

# Affichage de la table
knitr::kable(relief_bassins_long, format = "html")

Nous pouvons maintenant tracer le graphique.

ggplot(relief_bassins_long, aes(fill=variable, y=value, x=relief)) + 
  geom_bar(position='dodge', stat='identity') +
  theme_minimal() + 
  labs(x='Relief', y='Proportion', title='Proportion de relief par bassin') +
  theme(plot.title = element_text(hjust=0.5, size=20, face='bold')) +
  scale_fill_manual('Bassin-versant', values=c('coral2', 'steelblue', 'orange')) +
  theme(legend.position = c(.11, .85))

Au vu de ce graphique, les trois bassins sont largement similaires en termes de répartition de formes de reliefs. Cependant, en y regardant de plus près, le bassin de moyenne montagne, à savoir celui de la Tarentaine dans le Massif Central, présente une proportion de pentes (slope) moins importante que les bassins pyrénéen et alpin, respectivement représentés par le Ribérot et la Séveraisse. Le bassin de la Tarentaine présente 55 % de pentes contre un peu plus de 65 % pour les deux autres bassins. En contrepartie, le bassin de la Tarentaine présente un peu plus d’éperons (spur) que les deux autres bassins. La proportion d’éperons est de 17,4 % pour la Tarentaine contre à peu près 14 % pour le Ribérot et la Séveraisse. Il en est de même pour les fonds de vallées (valley), où la proportion atteint 6 % pour la Tarentaine et seulement 1,5 à 2 % pour les deux autres bassins. Ces résultats sont tout à fait logiques au regard du contexte des trois bassins. En moyenne montagne, les pentes sont moins abruptes et les fonds de vallées plus larges, ce qui se reflète bien dans cette classification.