Publié le

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)
fenêtre pour la méthode déterministe infinie

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

Si nous le comparons avec le résultat de la méthode D8:
résultat du calcul d'écoulement avec 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.

résultat du calcul de l'écoulement avec la méthode kra

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:
résultat de la méthode Braunschweiger Digitales Reliefmodell
Celui-ci est à comparer avec le résultat de la méthode MFD:
résultat du calcul avec 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:

résultat de la méthode DEMON de calcul de l'écoulement

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

Publié le

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)
méthode Rho8 pour le calcul de l'écoulement
Le résultat de cette méthode sur le MNT en forme semi-sphérique, pour la couche « catchment area » est le suivant:
résultat de l'écoulement Rho8

Si vous comparez maintenant avec le résultat de la méthode D8:
résultat de la méthode de calcul de l'écoulement 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.
écoulement multiple

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)

fenêtre de traitement sage : MFD

Le résultat de cette méthode sur le MNT en forme semi-sphérique, pour la couche « catchment area » est le suivant:
résultat de l'écoulement calculé avec la méthode MFD (FD8)

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.

Publié le

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.
    ouverture de la fenêtre traitements de QGis
    Et voilà comment ça se présente;
    fenêtre traitements de QGis

    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.
    mnt fictif pourt tester les méthodes de calcul de l'écoulement

    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)

    commande de calcul pour la méthode D8

    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:

    flux d'écoulement calculé avec la méthode D8

    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.

    Publié le

    Comment trouver le systèmes de coordonnées d’une couche dans QGis ?

    Normalement, toutes les couches de données géographiques possèdent la description relative au système de coordonnées. Malheureusement, ce n’est pas toujours le cas. En général se problème se double avec celui de ne pas (ou plus) pouvoir contacter le producteur de la donnée.

    Diagnostic du problème

    A- Vous demandez à QGis de charger une couche vecteur et vous avez la fenêtre « Sélectionneur du Système de Coordonnées de référence » qui s’ouvre. Si c’est une couche au format « shape », les données doivent manquer de système de coordonnées défini. Dans le cas de fichiers de formes, il manque le fichier PRJ.

    B- Vous demandez à QGis de charger la couche et, sans avoir de fenêtre qui s’ouvrer, vous ne voyez pas vos nouvelles données, ou elles sont complètement ailleurs de là où elles devraient être.
    Avant de vous dire qu’il y a un problème dans les données,vérifiez que vous avez coché la Projection à la volée dans « Propriétés du Projet »-> « SCR ».
    Si la case est cochée, les données ont un système de coordonnées défini, mais il est faux. Le fichier PRJ existe , mais son contenu est erroné.

    Procédure de recherche

    Dans cet article nous essayons de vous apporter un guide pour déterminer le système de projection d’une couche de données SIG, quand elle est inconnue.

    Nous nous sommes limités aux systèmes couramment utilisés en France par les différents organismes et administrations : les projections Lambert et Lambert93, les projections UTM et les données non projetées (géographiques en latitude/longitude).

    Bien sûr, il en existe beaucoup d’autres, utilisées ponctuellement. Mais il est impossible de les déterminer sans un véritable travail de détective.

    Nous donnons comme acquis que vous savez utiliser QGis et que vous possédez d’autres données de référence qui vous permettront de comparer et de juger des résultats obtenus.

    Vous pouvez télécharger le document pdf de cet article sur le site de NASCA, en cliquant ici.

    Vous devez trouver la plage de valeurs X et Y contenues dans les données.

    1- Démarrez QGis avec un nouveau projet vide
    2- Ajoutez les données dotées du système de coordonnées inconnu. Les données ne doivent pas avoir de système de coordonnées défini. Dans le cas de fichiers de formes, il ne doit pas posséder de fichiers PRJ. S’il y en a un, renommez-le différemment.
    3- Dans la fenêtre « Sélectionneur du Système de Coordonnées de référence »,sélectionne le système WGS84, EPSG:4326

    4- Sur la couche chargée, ouvrez le menu contextuel
    5- Cliquez sur Propriétés pour ouvrir la boîte de dialogue Propriétés de la couche,
    6- Sélectionnez l’onglet « Métadonnées », puis descendez dans le champ « Propriétés » jusqu’à l’item « Emprise ».

    La ligne « xMin,yMin -13146.9,6755689.52 : xMax,yMax 169765.86,6890746.18 » vous donne l' »Etendue » que nous allons utiliser par la suite.
    Vous aurez l’étendue des Y : yMin à yMax (6755689.52 , 6890746.18) et,
    l’étendue des X : xMin à Xmax (-13146.9 , 169765.86)

    Le terme système de coordonnées peut s’appliquer à des données exprimées en degrés décimaux (coordonnées géographiques) ou à un système de coordonnées projetées exprimé en mètres.
    Si les coordonnées trouvées pour l’Étendue sont exprimées en degrés décimaux, elle seront comprises entre -180 et +180 pour les longitudes (étendue des X) et entre -90 et +90 pour les latitudes (étendue des Y). Le système des données (SCR) sera à rechercher parmi les Systèmes de coordonnées géographiques. Il reste à trouver le système géodésique (Datum) des données. (Voir plus loin)

    Si les coordonnées pour l’Étendue sont de l’ordre des centaines de milliers ou des millions, il s’agit de mètres. Le système des données sera à rechercher parmi les systèmes de coordonnées projetées, et il restera aussi à trouver le système géodésique (Datum) des données.

    Il suffit maintenant de reprendre la méthode décrite dans l’article Comment trouver le systèmes de coordonnées d’une couche dans ArcGis 10.X? en utilisant les valeurs trouvés pour l' »Étendue », à partir du point Trouver le système de projection.

    Quand vous aurez trouvé le système de projection, revenez ici, pour finir le travail en déterminant le DATUM de vos données.

    Trouver le système géodésique (DATUM)

    Une fois franchies les deux premières étapes il reste un dernier point à déterminer. Tout système de localisation se réfère obligatoirement à un centre de la Terre. Comme la Terre n’est pas une sphère parfaite, et qu’il faut calculer son centre, il y a plusieurs manières de le calculer et par conséquent, plusieurs « centres » différents. La différence n’est pas énorme et jusqu’à il y a quelques décennies ce n’était qu’une discussion plutôt théorique, la différence de positionnement résultante étant, en général, inférieure à 300m.
    En principe, un système de coordonnées est toujours associé à un système géodésique.
    Les projections Lambert 1, 2, 3,4 et 2 étendue sont toujours associées au système NTF (Nouvelle Triangulation Française)
    Les projections Lambert 93, et CC42 à 50 sont toujours associées au système RGF.
    Donc, si dans l’étape précédente vous êtes arrivés à définir une de ces projections, le travail est fini.
    Pour les projections UTM 30 à 32, en principe le système associé est le système WGS84. Mais elles peuvent aussi être associées au système Europe 50.
    Dans le cas des données géographiques (non projetées) elles sont aussi associées en générale au système WGS84, mais on peut aussi les trouver associées au système NTF ou Europe 50.
    Dans l’étape précédente nous avons choisi, par défaut, le systèmes WGS84. Mais nous devoinbs maintenant vérifier que ce chois est le bon.
    Comment faire pour le savoir?
    • Il faut disposer d’une couche de données de référence, avec le système de projection défini correctement et surtout, ayant une bonne précision (détail).
    • Dans un projet nouveau dans QGis, chargez cette couche de référence.
    • Assurez-vous que l’option « Activer la Projection à la volée » est bien cochée dans « Propriétés du projet »
    • Chargez la couche à déterminer dans QGis, en indiquant la projection trouvée dans l’étape précédente dans la fenêtre « Sélectionneur de Système de coordonnées de référence ».
    Si
    • les données apparaissent au bon endroit et qu’il n’y a pas de léger décalage (100-300m), vous avez fini. La définition que vous avez adoptée est la bonne.
    • vous avez un décalage de toutes vos entités, de l’ordre de 100 à 300m, vos données ne sont pas en WGS84. Le système de coordonnées (UTM ou géographique) est bon, mais vos données ne sont pas en WGS84. Elles doivent être en Europe 50 ou NTF.

    Pour savoir quel est le bon Datum, il faut changer le SRC de la couche chargée :
    1-Ouvrez le menu contextuel en cliquant droit sur la couche,
    2-cliquez sur Définir le SCR d’une couche,
    Si la projection trouvée dans l’étape précédente est
    UTM fuseau 30N : rentrez dans « Filtre » le code IGNF :UTM30
    UTM fuseau 31N : rentrez dans « Filtre » le code IGNF :UTM31
    UTM fuseau 32N : rentrez dans « Filtre » le code IGNF :UTM33
    3-Dans la Liste des SCR mondiaux, sélectionnez cette projection et cliquez sur OK
    Si maintenant vous n’observez plus de décalage général des données, vous avez trouvé le bon système de coordonnées. Si ce n’est pas le cas, vous êtes devant une couche qui sort du cadre de cet article.

    Si le système trouvée dans l’étape précédente est un Système de coordonnées géographiques (Etendue en X entre -180 et 180 et en Y de -90 à 90)
    1-Ouvrez le menu contextuel en cliquant droit sur la couche,
    2-cliquez sur Définir le SCR d’une couche,
    3-Rentrez dans filtre ED50 et sélectionnez le SRC EPSG :4230
    Si maintenant vous n’observez plus de décalage général des données, vous avez trouvé le bon système de coordonnées (EPSG :4230) : vos données sont en géographique Europe 50.

    Si ce n’est pas le cas,

    1-Ouvrez le menu contextuel en cliquant droit sur la couche,
    2-cliquez sur Définir le SCR d’une couche,
    3-Rentrez dans filtre NTF et sélectionnez le SRC EPSG :4275
    Si maintenant vous n’observez plus de décalage général des données, vous avez trouvé le bon système de coordonnées (EPSG :4275) : vos données sont en géographique NTF.

    Si ce n’est pas le cas, vous êtes devant une couche qui sort du cadre de cet article.