index.Rmd

---
title: "Caractérisation des formes du relief à l'échelle de bassins-versants"
subtitle: "Analyse quantitative des formes du relief *via* l'algorithme *Geomorphons* pour trois bassins-versants montagnards"
# date: "2022-09-30"
date: "2023-05-02"
author: 
 - name: Olivier Theureaux
   affiliation: UMR LADYSS, CNRS
 - name: Paul Passy
   affiliation: UMR PRODIG - Géotéca, Université Paris-Cité 
 - name: Thierry Feuillet
   affiliation: UMR IDEES, Université Caen-Normandie
 - name: Déborah Birre
   affiliation: EA PLEIADE, Université Sorbonne Paris Nord
image: "treeline.png"
logo: "figures/rzine.png" 
output: 
  rzine::readrzine:
    highlight: kate
    number_sections: true
csl: Rzine_citation.csl
bibliography: biblio.bib
nocite : |
  @*
link_citations: true
github: "rzine-reviews/fiche_rzine_geomorphon"
# gitlab: "gitlab.huma-num.fr/author/repository"
doi: "10.48645/s8je-kz92"
licence: "by-sa"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Photographie : Pyrénées ariégeoises (cliché D. Birre, 2021)

<ul style="background-color: #dcdcdc ; border-color:navajowhite3; padding: 0.8em; font-size:115% ; color: CCCCCC">
**Prérequis** : La reproduction du code présenté dans cette fiche nécessite l'installation des logiciels QGIS et GRASS sur votre machine. 
</ul>

<p style="background-color: #dcdcdc ; border-color:navajowhite3; padding: 1em; font-size:115% ; color: CCCCCC">
Cette fiche présente une analyse quantitative des formes du relief réalisée à l'aide de l'algorithme *r.geomorphon*, disponible via le logiciel GRASS. Cet algorithme est accessible dans R grâce au package `qgisprocess`, dont la mise à disposition sur le CRAN semble être à l'étude (cf [issue 46](https://github.com/r-spatial/qgisprocess/issues/46){target="_blank"}). </p>


# Importance du relief dans la compréhension des processus bio-physiques

La caractérisation des formes du relief pour une région donnée est une tâche importante pour la compréhension des processus bio-physiques agissant sur ce territoire. Tout d'abord, les formes actuelles du relief résultent de processus s'étant déroulés dans le passé sur des échelles spatio-temporelles plus ou moins larges. Des territoires marqués par des épisodes glaciaires, des transgressions marines ou de longues périodes d'aridité présenteront des formes bien caractéristiques. Les lire nous renseigne immédiatement sur leurs histoires. Par exemple, un territoire sur lequel le géomorphologue identifie la présence de *drumlins* est un territoire qui a été englacé dans un passé plus ou moins lointain. Ces formes peuvent également nous renseigner sur les dynamiques actuelles. Ainsi, un banc de sable ou de galets au milieu d'une rivière nous indique que ce cours d'eau est sujet à des épisodes de hautes eaux présentant une énergie suffisante pour entraîner des sédiments plus ou moins gros.

Les formes du relief ont aussi un impact fort sur les processus hydrologiques. L'eau, par définition, suit les plus fortes pentes au sein d'une région donnée. Repérer les formes en creux sur un versant, ainsi que leurs pentes, nous permet d'identifier les chemins empruntés par les flux d'eau. Ces chemins, en retour, seront plus sujets à l'érosion hydrique que les régions voisines et pourvoiront les flux d'eau en sédiments et nutriments. Au contraire, les zones de replat ou de cuvette seront des zones d'étalement et de pièges à sédiments et à nutriments. La présence d'eau sur des terrains pentus pourra également influer sur les mouvements de terrain, tels que les éboulements ou les coulées de boues.

La géomorphologie a également un impact sur le climat local. Les zones sommitales ou les lignes de crêtes sont plus sujettes aux vents forts que les bas de versants par exemple. De même, les creux de versants ou les dépressions pourront connaître des températures plus basses que les zones alentours du fait d'une moins forte exposition aux rayons solaires ou de la subsidence de l'air froid. Cet impact sur le micro-climat aura, à son tour, des répercussions sur la biogéographie et la répartition des espèces. Certaines d'entre elles se développeront mieux dans les zones ensoleillées que les zones en ombre. D'autres, au contraire, profiteront de la position d'abri qu'offrent les creux de versants. La caractérisation des formes du relief d'un territoire donné intervient également dans les logiques de localisation de certains aménagements, tels que les éoliennes ou les infrastructures de transport.

Dans cette fiche nous présentons une analyse morphométrique de trois bassins-versants de montagne réalisée à partir d'outils géomatiques et de *modèles numériques de terrain* (MNT) à haute résolution. Le bassin-versant, défini comme l'espace géographique alimentant un cours d'eau et drainé par lui, est un niveau d'analyse pertinent car les actions bio-physiques qui s'y déroulent dépendent fortement de la géomorphologie. L'analyse morphométrique, à savoir l'analyse quantitative des formes du relief et des processus qui les ont produits, est réalisée ici à travers le langage R et différentes librairies d'analyse spatiale. Ce développement permet (i) de représenter la répartition géographique des formes du relief au sein de chacun des trois bassins, (ii) de procéder à des comparaisons morphométriques, (iii) de présenter des distinctions entre des bassins-versants de moyenne montagne et de haute montagne.

# Techniques de caractérisation du relief basées sur les MNT

La caractérisation automatique des formes du relief se fait dans la plupart des cas sur des rasters de MNT. Sur un MNT, à chaque pixel est associée une valeur représentant son altitude. Ainsi, en comparant les valeurs d'un pixel aux valeurs des pixels voisins, il est possible d'en dériver des caractéristiques géomorphologiques.

## 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*).

<figure class="center">
  <img src="figures/fenetres_glissantes.png" width="500">
  <figcaption style="font-size:13px;">*Figure 1 : Analyse de voisinage avec une fenêtre glissante de 3×3 pixels (à gauche) et avec une fenêtre glissante de 5×5 pixels (à droite)*.
  </figcaption>
</figure>

</br>

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.

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

<figure class="center">
  <img src="figures/geomorphon_10valeurs.png" width="600">
  <figcaption style="font-size:13px;">*Figure 2 : Les 10 geomorphons identifiés par la méthode. Par exemple pour le Pic (Peak), tous les pixels voisins sont plus bas que le pixel central, leurs altitudes relatives sont donc représentées en bleu. Pour le fond de vallée (Valley), deux pixels sont d'altitudes équivalentes donc représentés en vert et les autres sont plus hauts, donc représentés en rouge (figure reprise de Jasiewicz et Stepinski, 2012).*  
  </figcaption>
</figure>  

</br>

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

<figure class="center">
  <img src="figures/geomorphon_10formes_v2.png" width="500">
  <figcaption style="font-size:13px;">*Figure 3 : Exemple d’un raster présentant les 10 geomorphons types tels que calculés par l’algorithme (figure reprise de Jasiewicz et Stepinski, 2012). A noter que "summit" correspond au terme "peak" et que "depression" correspond au terme "pit" de la figure 2.*  
  </figcaption>
</figure>

# Objectifs de cette fiche

L'objectif de cette analyse est de caractériser et de quantifier les formes de relief grâce à la technique des *geomorphons* pour trois bassins-versants montagnards. Les trois bassins choisis sont le [Ribérot](https://fr.wikipedia.org/wiki/Riberot){target="_blank"} dans les Pyrénées, la [Séveraisse](https://fr.wikipedia.org/wiki/S%C3%A9veraisse){target="_blank"} dans les Alpes et la [Tarentaine](https://fr.wikipedia.org/wiki/Tarentaine){target="_blank"} dans le Massif Central. Ces trois bassins sont faiblement anthropisés et présentent des formes de reliefs et des altitudes variées. Pour chacun de ces bassins, la proportion de chaque classe de relief sera calculée. Nous tenterons ainsi de voir si les trois bassins partagent des traits communs et les possibles spécificités topographiques de chacun.

# Outil utilisé et données

## L'algorithme *Geomorphons*

L'algorithme *Geomorphons* est disponible dans QGIS *via* les logiciels [GRASS GIS](https://grass.osgeo.org/grass82/manuals/r.geomorphon.html){target="_blank"}, [SAGA GIS](https://saga-gis.sourceforge.io/saga_tool_doc/7.3.0/ta_lighting_8.html){target="_blank"} ou encore [WhiteboxTools](https://www.whiteboxgeo.com/manual/wbt_book/available_tools/geomorphometric_analysis.html){target="_blank"}. 

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](https://github.com/r-spatial/qgisprocess){target="_blank"}.

```{r message=FALSE, eval=TRUE, warning=FALSE, echo=TRUE, results='hide'}
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.

```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
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()

```

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

```{r message=FALSE, warning=FALSE, echo=TRUE, results=FALSE}
# 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)
```

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

<br/>

<p class="center">[<span style="font-size: 230%;" class="glyphicon glyphicon-download-alt"></span> <br/> Télécharger les données](https://rzine.fr/docs/20230425_geomorphon/data.zip)</p>

<br/>

Pour rejouer l'analyse, décompressez les données contenues dans l'archive [data.zip](https://rzine.fr/docs/20230425_geomorphon/data.zip)

### Les bassins-versants issus de la BD Topage® de l'IGN-F

Les limites des bassins-versants proviennent de la [base de données (BD) Topage®](https://www.sandre.eaufrance.fr/actualite/diffusion-de-la-1%C3%A8re-version-de-la-bd-topage%C2%AE-m%C3%A9tropole){target="_blank"}, référentiel hydrographique français, en remplacement du référentiel BD Carthage®. Au sein de cette base de données, nous téléchargeons la couche *Bassins-versant topographiques - Métropole 2019 - BD Topage®* au format *GeoPackage*. Cette couche se nomme *BassinVersantTopographique_FXX.gpkg*. Nous y trouvons l'attribut *TopoOH* qui contient le toponyme de chaque bassin. Nous nous basons sur cet attribut pour extraire les bassins de cette étude.

Pour rappel, nous travaillons ici sur trois bassins-versants : 

* un bassin-versant pyrénéen : le Ribérot (49 km<sup>2</sup>), situé en Ariège. Il prend sa source au mont Valier (2838 m).
* un bassin-versant alpin : la Séveraisse, du torrent de Navette au Drac (107 km<sup>2</sup>). Elle prend sa source dans les glaciers des Écrins.
* un bassin-versant du Massif Central : la Tarentaine, de sa source au confluent de l'Eau Verte (44 km<sup>2</sup>). Elle prend sa source à 1785 m dans les monts Dore.

### Le MNT RGE de l'IGN à 5 m

Le calcul des *geomorphons* se fait à partir d'un modèle numérique de terrain. Aujourd'hui, des MNT de très haute résolution sont couramment produits et sont disponibles relativement facilement. L'IGN met à disposition gratuitement des MNT de très haute résolution (1 m, 5 m) sur [l'intégralité du territoire français](https://geoservices.ign.fr/rgealti){target="_blank"}. Nous travaillons ici sur le MNT à la résolution de 5 m nommé *Référentiel à grande échelle (RGE) ALTI® 5 m*. À l'échelle des bassins étudiés, cette résolution est tout à fait suffisante. Notons que le RGE est concrètement un *modèle numérique d'élévation*, prenant donc en compte notamment les bâtiments et les infrastructures. Néanmoins, les bassins étudiés étant peu anthropisés, nous pouvons nous contenter de cette donnée.

# Mise en application

Nous commençons par préparer nos données afin de disposer d'un MNT à 5 m par bassin-versant. Ensuite, nous calculons les *geomorphons* sur chacun de nos trois MNT. Enfin, nous calculons quelques statistiques descriptives du relief des trois bassins.

## Préparation des données

### Extraction des bassins de la BD TOPAGE®

Nous commençons par charger la couche des bassins-versants à l'échelle de la France métropolitaine issue de la [BD TOPAGE®](https://www.sandre.eaufrance.fr/atlas/srv/fre/catalog.search#/metadata/8885e6fc-9c87-49e7-bd2f-d91f6632e44c){target="_blank"}. Nous la chargeons sous la forme d'un objet *`sf`*.

```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
# 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")

```

Cette couche contient un attribut *TopoOH* qui identifie le nom des bassins-versants. Dans cette couche, les trois bassins sont identifiés avec les noms suivants :

* Ribérot : *Le Ribérot*
* Séveraisse : *La Séveraisse de sa source au torrent de Navette inclus*
* Tarentaine : *La Tarentaine (Trentaine) de sa source au confluent de l'Eau Verte*

Par une sélection attributaire, nous individualisons nos trois bassins. Pour ce faire, nous utilisons la fonction *`filter`* de la librairie *`dplyr`*.

```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
# 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)
```

### Création du raster virtuel de MNT

Une fois les limites physiques des bassins-versants individualisées, nous construisons les MNT correspondants. Pour chaque bassin-versant, il est nécessaire de télécharger le MNT de chaque département depuis le site des [Géoservices de l'IGN-F](https://geoservices.ign.fr/rgealti){target="_blank"}. Pour nos MNT nous partirons sur un pas de 5 mètres. Attention le téléchargement peut être long.

Nos trois bassins se trouvent respectivement dans les départements suivants :

* Ribérot : Ariège (09)
* Séveraisse : Hautes-Alpes (05)
* Tarentaine : Puy-de-Dôme (63)

Pour chaque département, le MNT est découpé en dalles carrées. Une fois l'ensemble des dalles téléchargées, il faut créer un raster virtuel par département, comme présenté dans le code ci-après. Ce raster virtuel est découpé selon les limites des bassins.


```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
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)
```

Nous procédons de même pour le département des Hautes-Alpes.

```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
# 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)
```

Enfin, nous procédons pareillement pour le département du Puy-de-Dôme.

```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
# 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)
```

À ce stade, nous disposons des MNT de chacun des départements contenant nos bassins. La prochaine étape va consister à extraire seulement les MNT selon les limites de nos trois bassins.

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

```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
# 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*.

## 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é.

```{r message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
# 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.

| Classe de relief | Valeur de pixel | Code couleur par défaut |
|------------------|-----------------|-------------------------|
| Plan             | 1               | #dcdcdc                 |
| Sommet           | 2               | #380000                 |
| Crête            | 3               | #c80000                 |
| Haut de versant  | 4               | #ff5014                 |
| Éperon           | 5               | #fad23c                 |
| Pente            | 6               | #ffff3c                 |
| Creux de versant | 7               | #b4e614                 |
| Bas de versant   | 8               | #3cfa96                 |
| Fond de vallée           | 9               | #0000ff                 |
| Dépression       | 10              | #000038                 |

# Résultats

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

```{r, message=FALSE, warning=FALSE, echo=TRUE, results='hide'}
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 :

```{r, message=FALSE, warning=FALSE, echo=TRUE, results='hide'}

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.


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

```{r message=FALSE, warning=FALSE}
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.

```{r message=FALSE, warning=FALSE}
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.

```{r message=FALSE, warning=FALSE}
# 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.

```{r message=FALSE, warning=FALSE}
# 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`*.

```{r message=FALSE, warning=FALSE}
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.

```{r message=FALSE, warning=FALSE}
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.

# Récapitulatif de la méthode proposée

Nous avons mis au point une chaîne de traitement permettant de classifier et de comparer les formes de relief pour différents bassins-versants. La méthode présentée ici part de l'extraction de bassins choisis par l'utilisateur à partir de la couche hydrographique nationale de référence que représente la BD TOPAGE®. Cette chaîne de traitements l'IGN. Ensuite, pour chaque bassin, à partir de son MNT, est calculé le raster des formes de relief grâce à l'algorithme *geomorphons* implémenté dans GRASS GIS, interfacé par QGIS et appelé *via* la librairie *`qgisprocess`* de R. Pour finir, les résultats ont été représentés sous forme de cartes et de graphiques en bâtons groupés.






# Bibliographie {-}

<div id="refs"></div>



# Annexes {-}

## Infos session  {-}

```{r session_info, echo=FALSE}
kableExtra::kable_styling(knitr::kable(rzine::sessionRzine()[[1]], row.names = F))
kableExtra::kable_styling(knitr::kable(rzine::sessionRzine()[[2]], row.names = F))
```

## Citation {-}

```{r Citation, echo=FALSE}
rref <- bibentry(
   bibtype = "misc",
   title = "Caractérisation des formes du relief à l’échelle de bassins-versants",
   subtitle = "Analyse quantitative des formes du relief via l’algorithme Geomorphons pour trois bassins-versants montagnards",
   author = c("Olivier Theureaux", "Paul Passy", "Thierry Feuillet", "Déborah Birre"),
   doi = "10.48645/s8je-kz92",
   url = "https://rzine.fr/publication_rzine/20230425_geomorphon/",
   keywords ="FOS: Other social sciences",
   language = "fr",
   publisher = "FR2007 CIST",
   year = 2023,
   copyright = "Creative Commons Attribution Share Alike 4.0 International")
``` 

`r capture.output(print(rref))`


### BibTex : {-}

```{r generateBibTex, echo=FALSE}
writeLines(toBibtex(rref), "cite.bib")
toBibtex(rref)
``` 

<br/>