L’hydrologie avec un SIG, pour les nuls (que nous sommes): calcul de l’écoulement (4)

Dans les trois derniers articles, nous avons vu un exemple d’application de chacune des méthodes principales de calcul du flux d’écoulement. Bien qu’ils soient assez parlants d’eux mêmes, il faut bien passer par un peu de théorie pour compléter le sujet.
Tout d’abord, voyons ce qui se passe avec la méthode la plus répandue, utilisée en particulier pour tous les calculs hydrologiques d’ArcGis, la méthode D8.

Analyse des problèmes de la méthode D8

La méthode D8 permet à l’écoulement de chaque pixel de se faire en direction d’un seul voisin de réception. Cela revient à traiter les flux qui ont pour origine un pixel (une surface à deux dimensions, X et Y) comme si c’était une source ponctuelle (sans dimension) et de générer une ligne d’écoulement descendante (avec une dimension, la longueur) au lieu d’un chenal d’écoulement (avec deux dimensions, largeur et longueur). Bref, on mélange allégrement des torchons et des serviettes. Le voisin de réception unique impose également des restrictions sur les configurations du chemin d’écoulement possible, parce que le flux peut se faire seulement selon une direction cardinale (haut,bas,gauche,droite) ou sur une diagonale.
Les erreurs résultantes de ces limites sont différentes pour les terrains où le débit est divergent (où la largeur du chemin d’écoulement augmente en descendant la pente), convergent (la largeur du chemin d’écoulement diminue vers l’aval), ou est parallèle (la largeur du chemin d’écoulement est constante tout au long de la pente).

Les pixels amont d’une cellule sont appelés la Surface Contributive Spécifique (SCA) de cette cellule (le « bassin versant » de cette cellule).
La méthode D8 calcule la SCA correctement pour un écoulement parallèle (la largeur de la surface d’écoulement ne varie pas d’amont en aval) que lorsque le flux est dans la direction x ou y. Lorsque le flux forme un angle par rapport à l’orientation de la grille principale, deux types d’erreurs se présentent:

1-les erreurs qui affectent la direction du trajet d’écoulement et,

2- une erreur de sous-estimation de la SCA pour un trajet d’écoulement donnée.

La première source d’erreur se produit lorsque le flux forme un angle différent d’un multiple de 45 °. Par exemple, si le flux forme un angle de 30 ° par rapport à l’Est(mesurés dans le sens antihoraire), la direction de plus grande pente donnée par la méthode D8 sera vers le NE, et chaque pixel se déversera vers son voisin du NE.
Cela se produit par l’approximation de l’angle de l’écoulement à 45 °; le flux modélisé sera détournée de son vrai chemin de 15 °.
La seconde source d’erreur résulte de la projection unidimensionnelle de l’écoulement. Considérez un écoulement parallèle avec un angle de 45°( ou un multiple de 45 °).

Dans l’image, le pixel A est le sommet de la pente.
La surface amont d’écoulement sera calculée par une diagonale de pixels, se touchant par les coins. Pour le pixel situé en haut et à droite, il recevra, d’après le calcul la pluie tombée sur les trois pixels situes en diagonale, jusqu’au sommet.
En réalité elle devrait comprendre aussi une partie (la moitié) de chaque pixel gauche et droit de ces pixels amont. La surface amont d’écoulement est donc sous-estimée par un facteur de 2 lorsque vous utilisez la méthode D8.
Lors de l’utilisation D8, les surfaces amont calculées pour un écoulement parallèle sont correctes seulement pour les écoulement qui se produisent précisément dans une direction cardinale (0°, 90°, 180°, ou 270°), sont sous-estimées par un facteur de deux pour les écoulements qui se font précisément dans une direction diagonale (45°, l35°, 225° ou 3l5°), et sont inférieures à la valeur correcte d’un facteur compris entre 1 et 2 pour le écoulements qui se produisent entre les directions cardinales et les diagonales.

Pour un écoulement divergent comme celui que nous avons sur le MNT semi-sphérique des articles précédente, les choses ne se passent pas de la même manière. Voyons d’abord comment travaille la méthode D8 dans ce cas.

Dans l’image, on considère que le pixel A est le sommet de la montagne, sphérique ou conique. Comme il faut bien sortir de ce pixel vers le bas, supposons que la méthode décide que l’écoulement se fait vers le pixel situé au NE.
Les flèches indiquent les directions définies pour les autres pixels. Considérons maintenant le pixel B. Aucun des pixels amont (pixels Ouest, Sud Ouest et Sud contigus) ne se déverse vers lui. La méthode considère donc qu’il n’y a pas de SCA, de bassin versant, pour ce pixel. C’est seulement la pluie tombée sur lui même qui sera transmise vers son pixel aval.
En réalité, son SCA est dessinée sur l’image suivante:

Le pixel doit recevoir en réalité toute la pluie tombée dan le triangle formé par le sommet de la montagne et ses deux coins « latéraux ». Au lieu de ça, la méthode D8 lui affecte une SCA = 0.

Regardons ce qui se passe avec le pixel situé à gauche de B.

Pour une SCA réelle sensiblement égale à celle du pixel B, on va retrouver une SCA calculée correspondante aux deux pixels situés vers le Sud Ouest, donc SCA=2. Ceci est dû au fait que ce pixel se situe sur une direction multiple de 45° par rapport au sommet.

Ce phénomène va se reproduire tout au long de la pente. Nous aurons ainsi, des pixels tout en bas de la pente qui, comme le pixel B, n’auront pratiquement aucune SCA, et d’autres qui auront cumulé tout au long de la pente. Dans le cas d’école de la montagne sphérique ou conique, tous les pixels en bas de pente devraient avoir exactement la même SCA et le même flux cumulé.

Pour un écoulement convergent, comme ce serait le cas d’un cratère conique, on retrouve le même type de problème mais avec des écarts entre les directions privilégiées (0,45,90,…) et les autres, moins marqués.

A la lumière de ce qu’on vient de dire, vous pouvez maintenant mieux comprendre le résultat que nous avons obtenu dans nos autres articles pour l’écoulement avec la méthode D8.

Les « traits » blancs sur les directions privilégiées sont les pixels qui ont récupéré les SCA de leurs pixels amont, les zones noires en bord de l’image correspondent aux pixels qui se trouvent être comme notre pixel B de l’exemple.

Analyse des problèmes de la méthode Rho8

La méthode Rho8 tente de résoudre l’un des problèmes de méthode D8, la déviation du sens d’écoulement modélisées vers une directions cardinale ou diagonale. Ceci résulte de la structure arbitraire de la grille où on ne peut choisir que l’une de ces 8 directions.

La méthode utilisée se base sur le hasard et les probabilités, ce qui n’est jamais simple à expliquer. On va essayer d’abord la véritable explication:
La méthode Rho8 introduit une composante stochastique dans la méthode D8, ce qui permet de mieux respecter le sens d’écoulement suivant une pente et direction données. Comme dans la méthode D8. chaque pixel débouche dans l’un de ses huit voisins. Le choix du pixel de réception parmi les voisins est faite de manière stochastique. A l’un des voisins il lui est attribuée une probabilité p d’être choisi et à un autre voisin lui est attribuée une probabilité l – p.
Le système d’attribution de probabilités et l’objectif de cette méthode sont illustrés dans l’exemple suivant : Considérons une surface plane avec une pente orientée de 30 ° (mesurée dans le sens antihoraire à partir de l’Est). Selon la méthode D8, l’écoulement se fera vers le pixel voisin situé au NE. La direction d’écoulement aura donc une différence de 15° avec la véritable direction d’écoulement. La méthode Rho8 attribue une probabilité p pour l’écoulement vers son voisin du NE et une probabilité de 1′-p pour que l’écoulement se fasse vers l’autre possible voisin, le pixel situé à l’Est.
Au fur et à mesure que l’on descend la pente, certains pixels se déverseront vers leur voisin du NE, et le reste se déversera sur le voisin de l’Est. Si le nombre de pixels qui se déversent vers le NE par rapport au voisin de l’Est respecte les bonnes proportions (la proportion attendue est p / (1 – p)), les lignes de flux résultantes auront une direction globale de 30 °. La valeur de p doit être telle que la valeur attendue de la direction de trajet d’écoulement est égal à l’angle de la pente (« aspect »).

En des termes plus caricaturaux, si les conditions se prêtent (longueur de l’écoulement,etc.) le hasard fait que l’orientation est autant modifiée dans un sens comme dans l’autre et qu’on retrouve à l’arrivée l’orientation d’origine.

Bien que cette méthode fournit (en espérance mathématique) la direction de trajet d’écoulement appropriée, tous les autres problèmes identifiés pour la méthode D8 persistent. Mais en plus, la méthode Rho8 introduit ses propres problèmes: le hasard ne garantit pas des résultats reproductibles, et dans des endroits avec flux parallèles, les chemins d’écoulement adjacents ne sont pas parallèles, mais ils sont déviés au hasard, ce qui provoque souvent une convergence de flux qui devraient être parallèles. Ceci arrive en particulier sur des portions de surface planes, où les chemins de flux devraient rester parallèles et où on retrouve, à la place, des pixels avec des fortes accumulations.

Une fois que deux voies d’écoulement ont fusionné en raison de leur déviation aléatoire, il n’existe aucun mécanisme qui puisse les faire diverger à nouveau. Au fur et à mesure qu’on se déplace vers le bas de la pente, les erreurs augmentent à mesure que le flux devient de plus en plus concentré.

En définitive, certains pixels auront les valeurs de SCA surestimées, tandis que d’autres qui ont été manquées par les flux d’écoulement déviés au hasard auront des valeurs sous-estimées.

Analyse des problèmes de la méthode MFD

Les méthodes multi-directionnelles tentent de résoudre la limitation majeure de la méthode D8, la représentation unidimensionnelle de l’écoulement, en distribuant des flux à partir d’un pixel vers tous les pixels voisins qui se situent plus bas que lui. C’est la règle qu’on utilise pour distribuer les écoulements respectifs vers ces voisins qui change, selon la méthode utilisée.

Nous ne rentrerons pas dans les détails de la méthode MFD, basées sur la pente directionnelle des pixels, mais dans tous les cas, chaque pixel va s’écouler vers trois ou quatre autres pixels en fonction de la pente.

Lui-même, va recevoir de l’écoulement de trois ou quatre pixels amont voisins, écoulement de chaque pixel qui ne sera pas complet (comme dans la méthode D8), mais partielle.
Sur un plan incliné, un plan conique, ou circulaire, si l’on dessiner la zone amont de chaque pixel, on peut observer que celle-ci forme un triangle inversé, avec le sommet situé sur le pixel considéré et la base sur la limite amont du plan.

La valeur de SCA calculée pour un pixel (la surface contributive ou bassin amont) est correcte, si et seulement si, elle se trouve assez loin de bords de la zone couverte par le MNT. Si c’est le cas, la zone triangulaire est comprise dans le MNT, sinon la région triangulaire est incomplète et la SCA est sous-estimée.

Un paramètre de la méthode (la convergence) permet de compenser dans une certaine mesure, ce problème.Pour une montagne en forme de cône circulaire, une valeur de 1,1 a donné de très bons résultats. Le souci c’est que cette « calibration » se base sur la symétrie géométrique de la surface et va souffrir dans la mesure la surface naturelle du terrain s’écarte des surfaces géométriques symétriques utilisées pour la calibration de la méthode.

Conclusion (s’il y en a une…)

On va conclure ici cette série d’articles, non parce que le sujet peut être considéré comme clos, mais bien au contraire, parce que le sujet peut être considéré définitivement ouvert.

Le but de ces articles n’est pas d’apporter les réponses, mais de poser les questions. Une fois posées, ces questions, nous aborderons avec plus de modestie la modélisation hydrologique d’un terrain. Et à chacun de poursuivre sa route en essayant de répondre aux nouvelles questions posées.

28 février 2015

L’hydrologie avec un SIG, pour les nuls (que nous sommes): calcul de l’écoulement (3)

Dans les deux articles précédents (L’hydrologie avec un SIG, pour les nuls (que nous sommes): calcul de l’écoulement (1) et L’hydrologie avec un SIG, pour les nuls (que nous sommes): calcul de l’écoulement (2)) on a abordé trois méthodes de calcul du flux d’écoulement: deux méthodes unidimensionnelles (D8 et Rho8) et une méthode bidimensionnelle (MFD).

Maintenant on va voir une variation de la méthode la plus répandue (D8) et que nous avons abordé en premier.

Quatrième méthode : Déterministe infinie (D∞)

La différence de cette méthode avec la méthode D8 est que le flux est réparti dans les deux cellules adjacentes les plus basses et pas seulement dans la cellule la plus basse. C’est une méthode bidimensionnelle.

Dans QGis vous utiliserez, dans la fenêtre de traitements-> SAGA -> Terrain analysis-Hydrology -> Catchment area (parallel)

Le résultat du calcul de l’écoulement avec cette méthode, sur le MNT constitué d’une semi-sphère, est le suivant:

Si nous le comparons avec le résultat de la méthode D8:

On peut observer que, si bien le défaut de la surreprésentation dans les 8 directions (0,45,90,…) est toujours très marqué, les différences sont estompées.

Le rapport entre les cellules avec moins de surface par rapport à celles qui ont beaucoup de surface est de 1 à 700 dans le cas de la méthode D8, et seulement de 1 à 22 dans le cas de la méthode Déterministe infinie.

Cinquième méthode : algorithme de routage cinématique (KRA)

On peut citer ici la méthode KRA: l’algorithme de routage cinématique. Cet algorithme d’écoulement de flux est unidimensionnel. L’écoulement se comporte comme une balle qui roulerait sur le MNT, sans restreindre sa position au centre des cellules. Le résultat de cette méthode est intermédiaire par rapport aux deux précédentes. Le rapport entre les cellules avec moins de surface par rapport à celles qui ont beaucoup de surface est de 1 à 45.

On verra rapidement les méthodes restantes, car elles n’apportent pas beaucoup par rapport à celles déjà abordées.

Sixième méthode:Braunschweiger Digitales Reliefmodell

C’est une méthode multidimensionnelle. Elle s’apparente à la méthode MFD.
Elle est disponible dans QGis dans Fenêtre de traitements-> SAGA -> Terrain analysis-Hydrology -> Catchment area (parallel)
Le résultat du calcul est le suivant:

Celui-ci est à comparer avec le résultat de la méthode MFD:

Le résultat apparaît nettement meilleur que ceux produits par les méthodes unidimensionnelles, mais par rapport à la méthode MFD on voit bien un certain biais radial.

Septième méthode : DEMON (Digital Elevation Model Network)

C’est une méthode bidimensionnelle, mais assez compliquée. Elle se base sur la création d’une matrice de bassin versant amont et aval pour chaque pixel. Impossible de l’expliquer ici. Si vous êtes intéressé, mieux vaut aller à la source (Costa-Cabral, M. C.; Burges,S. J. Digital elevation model networks (DEMON): A model of flow over hillslopes for computation of contributing and dispersal areas. Wat. Resour. Res. 30: 1681–92. (1994))
Retenez une chose, il lui faut beaucoup de temps pour calculer. Sur notre exemple, qui est très simple, il lui faut environ 8 minutes pour calculer le résultat. Alors avec des vrais MNT …
Dans QGis elle est disponible dans QGis dans Fenêtre de traitements-> SAGA -> Terrain analysis-Hydrology -> Catchment area (flow tracing)
Voici le résultat obtenu:

Ce résultat est très proche de la méthode Déterministe infinie (voir en début de cet article), ce qui est normal car les deux sont des méthodes bidimensionnelles.

Nous avons fait le tour des principales méthodes de calcul de l’écoulement. Reste à voir un peu plus en détail les problèmes posés par les principales méthodes. C’est ce que nous aborderons dans le prochain et dernier article sur le sujet

27 février 201530 mars 2017

L’hydrologie avec un SIG, pour les nuls (que nous sommes): calcul de l’écoulement (2)

Cet article fait suite directe de l’article L’hydrologie avec un SIG, pour les nuls (que nous sommes): calcul de l’écoulement(1)

Nous allons voir maintenant le résultat de l’application d’une autre méthode de calcul d’écoulement.

Deuxième méthode : Rho8

Elle est une modification de la méthode précédente. C’est toujours un flux unidimensionnel: toute l’eau de la cellule passe vers une et une seule autre cellule. Si dans la méthode D8, le flux est calculé selon la pente la plus importante entre la cellule considérée et ses 8 cellules contiguës, dans la méthode Rho8 on va introduire un facteur aléatoire, calculé à partir de l’orientation de la cellule centrale et de la direction des deux cellules contiguës qui se situent dans cette direction.

Dans QGis vous utiliserez, dans la fenêtre de traitements-> SAGA -> Terrain analysis-Hydrology -> Catchment area (parallel)

Le résultat de cette méthode sur le MNT en forme semi-sphérique, pour la couche « catchment area » est le suivant:

Si vous comparez maintenant avec le résultat de la méthode D8:

Vous pouvez observer que la concentration sur les 8 directions principales (0,45,90,135,…) est nettement moins marquée. Par contre, étant donné que nous sommes sur un sphère parfaite, l’eau devrait s’écouler de manière homogène et non sous forme de « fils ».

Pour mieux comprendre ce commentaire, nous allons passer à une autre méthode de calcul.

Troisième méthode : Multiple Flow Direction (MFD) aussi dénommée FD8

Contrairement aux deux méthodes précédentes, l’écoulement à partir d’une cellule ne se fait pas en direction d’une et une seule cellule. Selon la pente du terrain, toutes les cellules situés plus bas que la cellule concernée recevront une partie de l’écoulement.

La figure suivante est tirée de la publication à l’origine de cette méthode.

Vous pouvez consulter cette publication directement :Quinn, P.F.; Beven, K.J.; Chevallier, P.; Planchon, O.; The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models, Hydrological Processes, 5: 59–79. 1991

Pour mettre en oeuvre cette méthode, avec QGis, vous utiliserez, la fenêtre de traitements-> SAGA -> Terrain analysis-Hydrology -> Catchment area (parallel)

Le résultat de cette méthode sur le MNT en forme semi-sphérique, pour la couche « catchment area » est le suivant:

Ici nous avons un écoulement globalement uniforme, sur une semi-sphère uniforme. Si on zoome, on s’aperçoit que la structure en 8 imprime un certain tramage. En effet, même si les flux sont multiples, il sont calculés de centre de la cellule centrale aux centres des cellules contiguës. La trame 0,45,90,… apparaît ici, mais quand on regarde au niveau du détail et disparait quand on regarde le phénomène de plus haut.

Dans le prochain article on abordera les autres méthodes de calculs restantes.

26 février 201530 novembre 2018

L’hydrologie avec un SIG, pour les nuls (que nous sommes): calcul de l’écoulement(1)

Une citation pour débuter l’article:
« De la même manière que le manuel d’un traitement de texte ne vous apprend pas à écrire une nouvelle ou un poème ou qu’un tutoriel de CAO ne vous montre pas comment calculer la taille d’une poutre pour un bâtiment, ce guide ne vous apprendra rien sur l’analyse spatiale. A la place, il vous montrera comment utiliser l’environnement des Traitements QGIS qui est un outil puissant pour réaliser des analyses spatiales. Il est de votre responsabilité d’apprendre les concepts qui seront indispensables à la compréhension de ce type d’analyses. Sans eux, vous n’arriverez pas à utiliser l’environnement et ses algorithmes même si vous pouvez être tenté de le faire. » (Documentation QGis)

Dans cette série d’articles, on ne va pas essayer de vous apprendre l’hydrologie. Nous allons juste ouvrir quelques portes qui restent, le plus souvent fermées. Si cela ne vous apprend pas l’hydrologie, au moins nous espérons glisser le doute sur ce que vous faites, doute qui vous permettra d’apprendre (si vous le souhaitez).

Dans les articles consacrés au calcul des bassins versants (ArcHydro : détermination des bassins versants d’un territoire (1), ArcHydro : détermination des bassins versants d’un territoire (2)) vous pourrez remarquer que la base de pratiquement tous les calculs que nous avons fait est la couche des directions d’écoulement (Flow direction)

Cette couche constitue l’élément clé de l’analyse hydrologique. De nombreuses alternatives existent, chacune d’elles ayant ses avantages et ses inconvénients. Les expliquer en détail prendrait trop de temps, et nous allons simplement présenter les idées qui sont intéressantes d’un point de vue pratique, afin que vous puissiez avoir une idée de la différence d’utiliser une méthode ou une autre.

Fondamentalement, les méthodes peuvent être divisés en deux groupes: celles qui considèrent le flux comme un déplacement entre les centres des cellules et celles dans lesquelles les flux se déplacent «librement» dans le MNT (algorithmes « Flow tracing »). Celles du premier groupe sont en relation avec la méthode D8 (la plus ancienne et la seule que vous trouverez dans ArcGis), tandis que celles du deuxième sont plus complexe et son utilisation est plutôt restreinte.

Une autre classification peut être faite séparant celles qui considèrent un flux unidimensionnel (communément appelés algorithmes de direction d’écoulement unique) et ceux qui considèrent un écoulement bidimensionnel (algorithmes de direction d’écoulement multiple). En d’autres termes, les premier considèrent que l’écoulement d’une cellule se fait sur une, et une seule, autre cellule, tandis que les deuxièmes considèrent que l’écoulement peut se faire sur plus d’une cellule contiguë de la cellule considérée.

Pour rester sur des logiciels connus, nous allons voir les différentes options proposées par QGis, sachant que sous ArcGis il n’y a pas de choix possible. D’autres méthodes peuvent être utilisées par des logiciels spécifiques au calcul hydrologique.

Avec QGis, vous avez les choix entre 8 méthodes de calcul de la direction d’écoulement:

Déterministe 8 (D8): La méthode classique, implémentée dans ArcGis. Le flux va du centre d’une cellule jusqu’au centre d’une (et seulement une) des cellules environnantes. Les directions de flux sont donc limitées à des multiples de 45 °, ce qui est la raison principale des inconvénients de la méthode. (O’Callaghan, J. F. y Mark D.M. The extraction of drainage networks from digital elevation data. Computer
Vision, Graphics and Image Processing 28: 323–44. 1984).

Rho8: Comme la précédente, mais avec une composante stochastique qui devrait améliorer les résultats. Le sens de circulation est déterminée par un argument aléatoire qui dépend de la différence entre l’orientation (« aspect ») et la direction des deux cellules voisines adjacentes. (Fairfield, J.; Leymarie P. Drainage networks from grid digital elevation models. Wat. Resour. Res. 27(5):709–717, 1991).

Déterministe infinie (D∞): Le flux passe d’une cellule à deux cellules environnantes contiguës, constituant ainsi un écoulement bidimensionnel, ce qui permet de résoudre les inconvénients de la méthode D8. (Tarboton, D.G.; Shankar, U. (1998), The Identification and Mapping of Flow Networks from DigitalElevation Data, Invited Presentation at AGU Fall Meeting, San Francisco, 1998).

Braunschweiger Digitales Reliefmodell: Un autre algorithme de direction de flux multiple. L’écoulement est divisé entre la cellule environnante dont l’orientation est la plus proche de l’orientation (« aspect ») de la cellule centrale et de ses deux cellules adjacentes. (Bauer, J.; Rohdenburg, H.; Bork, H.-R. Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluess, Landschaftsgenese und Landschaftsoekologie, H.10, Parameteraufbereitung fuer deterministische Gebiets-Wassermodelle, Grundlagenarbeiten zu Analyse von
Agrar-Oekosystemen, (Eds.: Bork, H.-R.; Rohdenburg, H.), p.1–15, 1985).

FD8 : un algorithme de calcul d’itinéraire d’écoulement bidimensionnel dérivé du D8. (Quinn, P.F.; Beven, K.J.; Chevallier, P.; Planchon, O.; The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models, Hydrological Processes, 5: 59–79. 1991).

algorithme de routage cinématique (KRA). Un algorithme d’écoulement de flux unidimensionnel. L’écoulement se comporte comme une balle qui roulerait sur le MNT, sans restreindre sa position au centre des cellules. (Lea, N. L. An aspect driven kinematic routing algorithm. En Parsons, A. J.; Abrahams, A. D. Overland Flow: Hydraulics and Erosion Mechanics, New York, Chapman & Hill. 1992).

Réseau Modèle Numérique de Terrain (Network DEMON): Le plus complexe. Un algorithme de traçage de flux bidimensionnel. Très gourmand en temps de calcul (Costa-Cabral, M. C.; Burges,S. J. Digital elevation model networks (DEMON): A model of flow over hillslopes for computation of contributing and dispersal areas. Wat. Resour. Res. 30: 1681–92. (1994)).

Le meilleur moyen de les comparer c’est de les appliquer tour à tour sur le même terrain et voir le résultat du calcul du flux d’écoulement.

Pour cela, nous allons utiliser QGis et son module de traitements. Pour beaucoup d’utlilisateurs de QGis, ce module est inconnu, car dans les versions précédentes il fallait l’activer pour qu’il soit accessible. A partir de la version 2.6 il est installé par défaut. Voilà comment y accéder.

Et voilà comment ça se présente;

Cette fenêtre permet de lancer des traitement extérieurs à QGis. Nous allons utiliser les traitements disponibles dans SAGA pour illustrer le calcul du flux d’écoulement.

Pour mieux les comparer, au lieu de les appliquer à un terrain plus ou moins complexe, nous allons les appliquer à un MNT fictif correspondant à une semi-sphère.

Nous allons donc calculer les directions d’écoulement sur ce terrain. Afficher les directions de chaque cellule ne nous apporterait pas grand chose, par contre afficher le résultat de l’accumulation de flux est beaucoup plus parlant. Cette couche découle du calcul des directions d’écoulement et montre le nombre de cellules amont pour chaque cellule du raster.

Si vous souhaitez faire vous-même les test, vous pouvez télécharger le mnt de test en cliquant ici.

Première méthode : D8

C’est la méthode utilisée par ArcGis et ArcHydro. Le flux est calculé selon la pente la plus importante entre la cellule considérée et ses 8 cellules contiguës.
C’est donc un flux unidimensionnel: toute l’eau de la cellule passe vers une et une seule autre cellule.

Dans QGis vous utiliserez, dans la fenêtre de traitements-> SAGA -> Terrain analysis-Hydrology -> Catchment area (parallel)

Parmi les différentes couches résultat, celle qui nous intéresse c’est la couche « catchment area ». Voici le résultat pour cette méthode:

Etant donné que nous sommes sur une semi-sphère, le flux d’écoulement devrait être homogène. Par contre on voit bien que le fait de n’avoir que des options séparées par 45° (d’une cellule on ne peut que aller vers une des 8 directions correspondantes aux 8 cellules contiguës) introduit un biais en créant des zones d’accumulation forte dans ces directions (0°,45°,90°,…)

Dans le prochain article nous verrons les résultats pour les autres méthodes de calcul.

25 février 2015

ArcHydro : détermination des bassins versants d’un territoire (2)

Dans l’article précédent, ArcHydro : détermination des bassins versants d’un territoire (1), nous nous sommes arrêtés à la définition de notre réseau hydrographique. Nous allons terminer e processus en abordant les trois derniers points:

Segmentation du réseau d’écoulement

Il ne faut pas oublier que nous sommes en train de travailler en mode raster. Chaque pixel (cellule) peut contenir une ou plusieurs valeurs, mais elles sont toutes « indépendantes ». Contrairement au réseau hydrographique de type vecteur, où nous avons des tronçons identifiés, ici rien ne nous permet de dire que deux pixels, mêmes contigus, appartiennent à un même tronçon hydrographique.

L’étape de segmentation à laquelle nous allons procéder répond à ce besoin: la commande va chercher les pixels contigus entre deux nœuds hydrographiques.

Les nœuds considérés sont: les sources, les jonctions, les exutoires.
Les pixels contigus entre deux nœuds auront le même identifiant de tronçon. Et chaque tronçon sera numéroté différemment.

Pour créer ce réseau, nous utilisons la commande « Terrain Preprocessing » -> « Stream Segmentation »

Visuellement, le raster résultant est identique au précédent. Pour voir la différence vous pouvez changer la symbologie en « Valeurs uniques ». Vous verrez alors les pixels appartenant à chaque tronçon avec la même couleur.

Détermination du bassin versant de chaque tronçon du réseau d’écoulement

Maintenant que nous avons construit les tronçons du réseau hydrographique nous allons déterminer les bassins versants de chaque tronçon, c’est à dire les cellule qui se déversent sur chacun des tronçons. Contrairement à notre notion classique de bassin versant qui prend en comte le ruissellement des berges di tronçon ainsi que les apports des tronçons en amont de celui qui est considéré, dans cette étape nous allons seulement prendre en compte le ruissellement direct sur chaque tronçon, sans tenir compte des apports qui peuvent arriver par les biais des affluents du tronçon.

Pour cela nous utiliserons la commande « Terrain Preprocessing » -> « Catchment Grid delineation »

Les entrées sont le raster avec les directions d’écoulement et le raster avec les tronçons hydrographiques.
La sortie sont les bassins versants de chaque tronçon.

Maintenant, nous souhaitons voir les bassins versants tels que nous les connaissons, pas tronçon par tronçon, mais bien par cours d’eau principal. C’est ce que nous allons faire dans la prochaine étape.

Fusion des bassins versants de chaque élément pour constituer des bassins versant de la taille souhaitée

Pour continuer dans cette étape, nous allons quitter le monde raster pour rejoindre le monde vecteur. Même s’il y a moyen de faire ce que nous souhaitons avec des rasters, il est beaucoup plus simple de travailler en mode vecteur.
Nous allons donc convertir en vecteur les bassins versants des tronçons et le réseau de tronçons.
Pour convertir les bassins versants des tronçons, nous utilisons la commande « Terrain Preprocessing » -> « Catchment Polygon Processing »

Pour convertir le réseau de tronçons en polylignes, nous utilisons la commande « Terrain Preprocessing » -> « Drainage Lines Processing »

Vous pouvez voir que les propriétés de chaque tronçon permettent de savoir les nœuds qui le limitent ainsi que le tronçon aval suivant.

Nous allons maintenant calculer les bassins contigus avec la fonction « Adjoint Catchment Processing ».

Cette fonction génère le bassin versant agrégé, cumulé et contigu en amont de chaque bassin de tronçon hydrographique. Pour chaque bassin versant généré, le polygone construit délimite toute la zone contigüe amont du tronçon. Si Deux tronçons confluent vers le tronçon en question, le polygone entoure les trois bassins de tronçons. Puis, on cherche quels bassins de tronçons sont contigus à ce polygone et on en génère un nouveau bassin versant contigu. On effectue ce processus jusqu’à ne plus avoir de bassins de tronçons qui s’écoulent vers le bassin versant généré.

Pour exécuter cette commande on utilise « Terrain Preprocessing » -> « Adjoint Catchment Processing ».

Les entrés sont les lignes d’écoulement (« DrainageLine ») et la couche de polygones bassins versants (« Catchment”). Le résultat est stocké dans une nouvelle couche vectorielle de type polygone (« AdjointCatchment ») que vous pouvez renommer.

Cette couche nous permettra de réaliser toute une série de calculs intéressants, mais la commande a aussi ajouté une champ attribut dans la table de « Catchment » avec l’identifiant du bassin aval et un autre champ attributaire dans la table de DrainageLine avec l’identifiant du bassin versant du tronçon.

24 février 2015

ArcHydro : détermination des bassins versants d’un territoire (1)

Nous allons utiliser ArcHydro pour calculer les différents bassins versants d’un territoire. Il est important de connaître les différentes étapes pour arriver au résultat, mais aussi de savoir sur quelles hypothèses on se base. En effet, les hypothèses utilisées par ArcHydro ne sont pas les seules existantes pour faire ce travail, par contre elles sont les seules disponibles dans le cadre d’ArcGis.

La détermination des bassins versants se fait en suivant les étapes:

1- Détermination de la direction d’écoulement

Chaque cellule du MNT est entourée par 8 autres cellules.

Le principe d’ArcHydro est d’établir vers quelle cellule s’écoule l’eau, à partir de la cellule centrale. Pour cela il calcule la pente entre la cellule centrale et les 8 cellules environnantes. Il considère que l’eau s’écoule vers la cellule qui a la pente la plus forte.
Cette hypothèse est une des bases d’ArcHydro, mais elle n’est pas universelle. D’autres algorithmes (non disponibles avec ArcGis) considèrent que, bien qu’il y ait un écoulement plus important vers la cellule avec plus forte pente, il y aura un écoulement vers toutes les cellules plus basses que la cellule centrale.

Pour suivre le processus nous allons utiliser le MNT corrigé lors des articles précédents (ArcHydro : 2- Préparer un MNT corrigé pour l’hydrologie – 1ère partie et ArcHydro : Préparer un MNT corrigé pour l’hydrologie – 2ème partie).

Pour calculer la direction d’écoulement cliquez sur « Terrain Preprocessing »-> « Flow direction »

Le résultat montre les différentes directions d’écoulement calculées à partir du MNT.

2- Détermination de l’accumulation d’écoulement

La deuxième étape consiste à calculer combien de cellules s’écoulent, en amont, vers chacune des cellules du MNT.

Quand aucune cellule s’écoule vers la cellule considérée, la valeur est de 0. Quand plusieurs cellules s’écoulent vers la cellule considérée, la valeur est égale à la somme de valeurs d’écoulement de ces cellules.

Pour calculer la direction d’écoulement cliquez sur « Terrain Preprocessing »-> « Flow accumulation »
Les seuls champs à renseigner sont le raster en entrée (Flow direction) et le nom du raster en sortie.

Le résultat se présente comme suit:

3- Détermination du réseau d’écoulement

Pourquoi réseau d’écoulement et pas réseau hydrographique? Le réseau hydrographique répertorie les « cours d’eau ». La qualification de cours d’eau repose essentiellement sur les deux critères suivants :

la présence et la permanence d’un lit naturel à l’origine, distinguant ainsi un cours d’eau d’un canal ou d’un fossé creusé par la main de l’homme;

la permanence d’un débit suffisant une majeure partie de l’année.

Dans notre cas, on ne tiendra pas compte de ces critères. Nous allons construire un réseau seulement à partir des données d’accumulation. On construit un réseau à partir d’un seuil d’accumulation, par exemple 150. Ceci revient à dire que, du moment où l’accumulation atteint 150 cellules, on considère que l’on est sur un tronçon de notre réseau hydrographique. Ceci est complètement indépendant du fait que cette cellule se trouve ou non sur le tracé d’un cours d’eau « officiel ». Il est clair qu’à partir d’un certain seuil, on se trouvera forcément sur les cours d’eau « officiels ».

La définition du seuil ne répond à aucune règle absolue. Il dépend de l’objectif de l’étude, de la taille de votre zone d’étude, du type de terrain, bref, de beaucoup de paramètres. Ceci dit, ça ne vous avance pas beaucoup. Alors, comment trouver la bonne valeur?

Dites vous que vous allez, à la fin du traitement, obtenir une photographie de l’ensemble de votre territoire. La valeur que vous allez choisir ici, pour créer votre réseau hydrographique, s’apparente au grain de votre photo. C’est la taille la plus petite pour laquelle vous pourrez avoir une information. En fonction de votre objectif, et sachant que plus vous ajoutez de la définition à votre photo, plus long et compliqué est le processus, vous devez trouver un bon compromis. Il est inutile d’ajouter de la définition très fine si vous ne recherchez qu’une vue globale, et il serait malvenue de choisir un gros grain si vos objectifs sont des études de détail de votre territoire.

Pour calculer le réseau d’écoulement cliquez sur « Terrain Preprocessing »-> « Stream definition ». La fenêtre de définition s’ouvre:

Vous pouvez rentrer soit le nombre de cellules d’accumulation, soit la surface moyenne de bassin versant à partir desquelles on considère qu’il y a « cours d’eau ». Ce qui est très utile c’est que, dès que vous rentrez une des deux valeurs, ArcHydro calcule automatiquement l’autre.

Pour notre exemple, nous rentrons 150 cellules, ce qui correspond à un bassin versant d’environ 0,84 km².

Le résultat de la commande est le suivant:

Dans le prochain article on utilisera ce réseau calculé pour obtenir les bassins versants.

22 février 2015

ArcHydro : Préparer un MNT corrigé pour l’hydrologie – 2ème partie

Nous avons vu dans l’article précédent,ArcHydro : Préparer un MNT corrigé pour l’hydrologie – 1ère partie, comment éliminer les dômes provoqués essentiellement par la végétation.
Nous allons maintenant effectuer deux opérations:

prendre en compte des surfaces d’eau (lacs, étangs)

imprimer le réseau hydrologique connu sur le MNT

Puis, en dernière opération de préparation du MNT nous allons combler les cuvettes existantes dans notre MNT.

Prise en compte des surfaces d’eau

Le but de cette étape est de permettre au MNT d’avoir des surfaces lisses pour représenter les surfaces occupées par les lacs et les étangs. Pour cela il nous faut une couche vectorielle de polygones correspondants aux lacs. ArcHydro croisera cette couche vectorielle avec le MNT et, pour les cellules correspondantes, il affectera une valeur constante d’élévation.
Quelle valeur?
Deux possibilités:
– vous indiquez la valeur d’élévation dans un champ attribut de la table
– ArcHydro cherche la valeur minimale d’élévation sur les cellules qui constituent le bord du lac et affecte cette valeur (éventuellement avec un offset)

Pour lancer le processus, cliquez sur « Terrain Preprocessing »-> »DEM Manipulation » -> « Level DEM »

La fenêtre permet de configurer la sortie et la façon d’affecter une valeur unique à la surface de chaque lac.

Impression du réseau hydrologique sur le MNT

Bien sûr, on peut utiliser le MNT pour calculer les cours d’eau. Mais si on dispose d’une couche vectorielle des cours d’eau, on peut corriger les cellules proches de ceux-ci. Les corrections comprennent le biais de pente introduit par la végétation des rives, mais aussi assurent que le point le plus bas des deux berges se retrouve sur la cellule correspondante au cours d’eau.

Pour imprimer le réseau hydrologique nous utilisons la commande « Terrain Preprocesing » -> « DEM Manipulation » -> « DEM Reconditioning ».

En plus des trois premiers champs (mnt en entrée, couche linéaire de cours d’eau, mnt en sortie) il y a trois paramètres que vous devez configurer. Pour cela vous devez avoir une connaissance de la zone.
Le premier « Stream Buffer » et la largeur de la zone à « creuser » dans le mnt, de part et d’autre de la ligne de cours d’eau. Selon la taille de cellule du MNT, vous devez définir l’empreinte à faire pour chaque cours d’eau. Dans notre exemple, la maille du MNT est de 75m. En entrant une valeur de 2 nous provoquons une empreinte de 150m de part et d’autre du cours d’eau (la terrasse alluviale).
Le deuxième « Smooth raise/drop » correspond au dénivellement, en mètres, de la terrasse alluviale ; différence de hauteur entre la cellule immédiatement à l’extérieur du « Stream Buffer » et le bord du cours d’eau. Dans notre exemple nous indiquons 5 mètres.
Le dernier « Sharp raise/drop » est la profondeur moyenne du cours d’eau, c’est à dire la cellule directement en dessous de la ligne de la couche « cours d’eau ». Nous indiquons ici 2mètres.
Une fois exécutée la commande pouvez utiliser le bouton « profil » de la barre d’outils 3D Analyste pour voir l’effet produit.

Sur le profil de gauche vous voyez la couper AVANT la commande. Sur le profil de droite vous voyez la même coupe APRES l’exécution de la commande.

Remplissage des cuvettes

La dernière étape pour obtenir un MNT hydrologiquement correct est de remplir les cuvettes existantes, soit dans le MNT d’origine, soit créées pendant les différentes étapes du pré-traitement. Les cuvettes sont des zones entourées par des cellules plus hautes et qui vont bloquer les directions d’écoulement.

Pour les remplir nous utilisons la commande « Fill sinks » (Terrain Preprocessing -> DEM Manipulation -> Fill Sinks). Le seul élément à bien observer c’est de cocher la case Fill All, pour que toutes les cuvettes soient repérées et remplies.

L’image suivante montre l’effet du remplissage sur une cuvette du MNT.

Le résultat c’est un MNT Hydroloigiquemenbt correct:

Nous allons utiliser ce MNT pour voir les différents calculs possibles sur le réseau hydrologique d’un territoire.

Pour commencer, dans le prochain article, nous verrons comment calculer les bassins versants des cours d’eau.

19 février 2015

ArcHydro : modeliser les ressources hydriques d’un territoire : 1-Installation

Présentation d’Arc Hydro

Arc Hydro est un ensemble de modèles de données et d’outils qui fonctionnent sous ArcGIS pour prendre en charge les analyses de données surv les ressources hydriques superficielles.
Vous pouvez utilisez Arc Hydro pour délimiter et caractériser les lignes de partage des eaux, définir et analyser le réseau hydrogéométrique, gérer les données chronologiques et exporter les données vers des modèles numériques.

Arc Hydro se compose de deux éléments clés :

le modèle de données Arc Hydro ;

les outils Arc Hydro.

Ces deux composants structurent une base de données et un ensemble d’outils qui facilitent les analyses effectuées dans le domaine des ressources hydrologiques. Arc Hydro est destiné à fournir les fonctionnalités initiales qui peuvent ensuite être étendues par l’ajout de structures de base de données et de fonctions requises par une tâche ou une application spécifique.

Les outils Arc Hydro fonctionnent dans l’environnement ArcGIS. Certaines des fonctions nécessitent l’extension Spatial Analyst. La majorité des outils sont accessibles via la barre d’outils Arc Hydro, où ils sont regroupés par fonctionnalité en six menus et neuf outils. D’autres outils ont été développés dans l’environnement de géotraitement et sont disponibles dans la boîte à outils Arc Hydro utilisables à la fois dans ArcMap et dans ArcCatalog.

Les outils Arc Hydro ont deux objectifs principaux. Le premier est de manipuler (attribuer) des attributs clés dans le modèle de données Arc Hydro. Ces attributs constituent la base des analyses ultérieures. Ils comprennent les identifiants clés (tels que HydroID, DrainID, NextDownID, etc.) et les attributs de mesure (tels que LengthDown). Le second est de fournir certaines fonctionnalités de base souvent utilisées dans les applications de ressources hydrologiques. Cela inclut la délimitation des bassins versants basés sur le MNT, la génération de réseaux et le traçage basé sur les attributs.

Pour télécharger la version 2.0 des modèles de données et des outils Arc Hydro, pour ArcGis 10.X
http://downloads.esri.com/ArcHydro/ArcHydro/Setup/
Le site propose les versions d’ArcHydro disponibles pour les versions 10.0,10.1,10.2.1 et 10.2.2, et même pour la 10.3.

En cliquant sur la version d’ArcGis qui vous correspond, vous trouverez deux téléchargements possibles:

le premier est la version 32 bits

le deuxième avec x64 dans le nom est la version 64bits

Une fois le fichier téléchargé il suffit de double-cliquer sur celui-ci pour installer l’extension.
L’extension comporte trois éléments:

une barre d’outils et

deux boîtes à outils

Nous allons aborder dans les prochains articles comment générer, pas à pas, les bassins versants d’un territoire avec ArcHydro.