Tutoriel: Calculer l’indice d’humidité topographique à partir d’un MNT

Site: OpenCourseWare for GIS
Cours: Tutoriels avancés de QGIS
Livre: Tutoriel: Calculer l’indice d’humidité topographique à partir d’un MNT
Imprimé par: Guest user
Date: vendredi 29 mars 2024, 15:07

1. Introduction

L’indice topographique d’humidité (Topographic Wetness Index - TWI) est un modèle utile pour estimer où l’eau s’accumulera dans une zone avec des différences d’altitude. Il s’agit d’une fonction de la pente et de la zone de contribution en amont:

TWI =

Où:

a = zone de contribution en amont (m2)
b = pente en radians.

Dans ce tutoriel, nous allons implémenter cette équation dans QGIS. La donnée d’entrée est un MDT.

Après ce tutoriel, vous serez en mesure :
  1. Interpoler les vides/lacunes dans un MDT.
  2. Calculer la pente en degrés.
  3. Convertir la pente de degrés en radians.
  4. Calculer TWI.

2. Interpoler les lacunes dans le MNT

Le MDT fourni est dérivé des données LIDAR des Pays-Bas (AHN3). Les données de la projection néerlandaise d’Amersfoort (EPSG : 28992). Le MDT, cependant, a des vides. Les vides sont des pixels sans données. Ces vides doivent être interpolés avant que nous puissions continuer.

1. Lance QGIS.

2. Ajoutez le MDT fourni (dtm.tif) à un projet vide.

Vous pouvez clairement voir les vides/lacunes.

voids

3. Accédez à la boîte à outils traitement : dans le menu principal, choisissez Processing | Toolbox.

4. Dans la boîte à outils de traitement, sélectionnez GDAL | Fill nodata.

5. Dans la boîte de dialogue ‘Fill nodata’ choisir dtm comme couche d’entré ‘Input layer’. Conservez les autres paramètres par défaut (vous pouvez jouer avec ceux-ci si les vides sont plus grands). Enregistrez le résultat en tant que ‘dtm_voidfilled.tif’.

6. Cliquez sur Exécuter ‘run’. Fermez le dialogue après traitement.

7. Retirez la couche dtm du panneau des couches. 

Dans le résultat, vous pouvez clairement voir que les vides ont été interpolés.

voids interpolated

8. Stylez le DTM avec un pseudocolore à bande unique. Vous pouvez utiliser une rampe de topographie du catalogue cpt-city. Mélangez le DTM avec le ‘hillshade’ (voir d’autres tutoriels sur la façon de le faire).

Styled DTM

Maintenant que le DTM est prêt, nous pouvons procéder au calcul de la pente.

Notez que dans QGIS 3.14 vous devez définir la projection de la couche raster après chaque calcul. Ce n’est pas le cas dans les versions LTR 3.10.


3. Calculer la pente en degré et en radian

Pour calculer TWI nous avons besoin de la pente dans les radians. Avec les outils d’analyse raster, nous pouvons calculer la pente en degrés et en pourcentages. Donc, nous allons d’abord calculer la pente en degrés et ensuite convertir les valeurs en radians.

1. Dans le menu principal, allez à Raster | Analysis | Slope...


2. Dans la boîte de dialogue Pente, choisi dtm_voidfilled comme couche d’entrée ‘Input layer’. Cocher la case d’utilisation de la formule ZevenbergenThorne au lieu de celui de Horn. La formule ZevenbergenThorne marche mieux avec moins de relief, lors que nous utilisons celle de Horn pour les zones raides. Choisir slopedegrees.tif comme résultat de la couche de Pente.



3. Cliquez sur Exécuter. Fermez le dialogue après traitement.

Parce que plus tard, nous allons calculer le logarithme naturel (ln) nous devons vérifier s’il ya des pixels avec une pente de 0 degrés.

4. Cliquez à droite sur slopedegrees et choisissez Propriétés.

5. Dans la boîte de dialogue Propriétés de la couche, accédez à l’onglet Informations et vérifiez la valeur minimum (STATISTICS_MINIMUM).



Dans ce cas, la pente minimale est supérieure à 0, donc nous pouvons procéder. Toutefois, si la pente du minumum est de 0, vous devez effectuer les étapes suivantes, que vous pouvez également appliquer à ces données, car il n’aura aucun effet.

6. Cliquez sur Annuler pour fermer la boîte de dialogue Propriétés de la couche.

7. Dans le menu principal, aller à ‘Raster | Raster Calculator’.



8. Dans la calculatrice Raster, créez l’équation suivante:

( "slopedegrees@1" <= 0 ) * 1 + ( "slopedegrees@1" > 0 ) * "slopedegrees@1"

Cela signifie: si les cellules ont une pente inférieure ou égale à 0, alors donner la valeur 1 (True) et multiplier par 1. Pour les cellules qui sont plus grandes que 0, assigner la valeur 1 (True) et multiplier avec sa valeur de pente. Le résultat sera 1 pour les cellules <= 0 et la pente d’origine pour les cellules 0.



9. Enregistrez le résultat sur slopedegrees_modified.tif et cliquez sur OK.

Pour passer des degrés aux radians, nous devons multiplier les valeurs de pixels par 0.01745

 

10. Retournez à la calculatrice Raster et calculez l’équation suivante :

"slopedegrees_modified@1" * 0.01745


11. Enregistrez le résultat sur sloperadians.tif et cliquez sur OK.

12. Styler le résultat du jaune au rouge.

13. Enlever les autres couches de pente.


Slope in radians

L’étape suivante consiste à calculer la zone de pente amont contributive pour chaque pixel du MDT.

4. Calculer la pente amont qui contribue

Dans cette section, nous allons calculer la zone de la pente amont contributive pour chaque pixel dans le MDT. 

1. Accédez à la boîte à outils de traitement et choisissez SAGA | Simulation | Accumulation de flux (qm of esp).


2. Dans le dialogue choisissez dtm_voidfilled comme MNT, conservez le prétraitement par défaut pour remplir temporairement (fill sinks temporarily) les puits et enregistrez la zone contributive sous contributing_upslope_area.sdat.

3. Cliquez sur Exécuter. Fermez le dialogue après traitement. 

Le résultat est pour chaque pixel la quantité de pixels se trouvant en amont.

Contributing Upslope Area

Maintenant, nous sommes prêts à calculer TWI.

5. Calculer l’indice d’humidité topographique TWI

1. Allez à la calculatrice Raster. 

2. Écrivez l’équation suivante :

ln ( ( "contributing_upslope_area@1" * 5 * 5 ) / tan ( "sloperadians@1" ) )

Notez que nous multiplions la zone de la pente ascendante (quantité de pixels) avec la surface du pixel (5 m x 5 m).


3. Enregistrez le résultat sous twi.tif.

Notez que si le résultat CRS n’est pas connu, vous devez le définir pour la couche d’entrée.

4. Cliquer OK.

5. Styler le résultat avec une rampe de couleur du rouge au jaune et au bleu et utiliser le mélange avec la couche hillshade.


TWI

Nous pouvons également visualiser le résultat en 3D.

6. Visualiser l’indice d’humidité topographique dans une vue en 3D

Nous pouvons inspecter davantage les résultats dans la vue 3D. 

1. Dans le menu principal choisir: Voir | Nouvelle vue de carte 3D (View | New 3D Map View) 

2. Configurer la scène en cliquant sur

3. Définissez le type de terrain (Terrain Type) à DEM (Raster Layer) et choisissez la couche dtm_voidfilled pour Elevation. Utilisez une échelle verticale de 3 pour exagérer le terrain.

4. Cliquez sur OK et inspectez le résultat.