Grass Gis et Paraview: la modélisation géologique 3D pour les démunis, ou les enthousiastes...

À l'heure où il existe des logiciels spécialisés et performants dans la modélisation 3D en géologie, comme Gocad, que peut faire un géologue sans-le-sou, ou qui ne travaille pas dans un organisme (universitaire ou non) qui dispose de licences pour de tels produits ou amateur de logiciels Open Source ?

Les résultats obtenus sont tels que l'on pourrait se dire, arrêtons ! Et pourtant, depuis longtemps, des géologues utilisent Grass pour essayer de le faire, mais les trouver, et surtout obtenir et comprendre leur démarche, s'apparente un peu à la recherche du Graal par les Monty Python ("Sacré Graal")... (Cette démarche a été initiée en souvenir de Louis Franssen (†), géologue, ancien responsable et moteur du programme "Carte géologique de Wallonie" qui m'a fait partager sa passion pour la modélisation 3D en géologie, lui l'entamant dès 1997 avec son cher Pov-Ray )

Principes

Les démarches possibles peuvent être résumées en quelques points:

1) modélisation de surfaces (couche géologique, faille, etc.): Grass ne permet pas de créer directement des TIN à partir de vecteurs, comme d'autres logiciels SIG, mais il a des fonctions beaucoup plus sophistiquées pour modéliser des surfaces sous forme de raster 3D avec des modules comme v.surf.rst. Couplé avec le module r.mapcalc, l'ensemble devient un outil de modélisation extrêmement puissant qui permet, par exemple, de créer de véritables "rubans". Il permet aussi d'importer des fichiers 3D dans divers formats comme le dxf:

                                       

image provenant de grass.itc.it/grass63/screenshots/vector.php

2) modélisations de volumes: puisque Grass sait créer et traiter des volumes (voxels 3D) élaborés à partir de plusieurs surfaces ou par calcul (grass.itc.it/grid3d) 

images  provenant de www-pool.math.tu-berlin.de/~soeren/grass/modules/screenshots/  ou grass.osgeo.org/wiki/Help_with_3D

3) Grass sait exporter des objets au format VTK ("VTK est une bibliothèque portable de visualisation de données scientifiques écrite en C++ et orientée objet. Elle permet de faire des traitements sur ces données en créant très simplement une chaîne d'algorithmes pour produire au final une image 2D/3D", fr.wikipedia.org/wiki/VTK). Il est donc possible d'utiliser des logiciels comme Paraview ou Mayavi en Python (Open Source) pour effectuer tous les traitements et réimporter les résultats dans Grass (grass.osgeo.org/wiki/GRASS_and_Paraview)

image provenant de www-pool.math.tu-berlin.de/~soeren/grass/modules/screenshots/LausanneDemoDataset3dParaview_11.png 

  

Geolmodel: image provenant de la liste Grass-dev : n2.nabble.com/Searching-Docs-about-3D-geological-modelisation-td4269207.html#a4269207

 

Une première solution, pourtant inaccessible...

Les modélisations les plus avancées sont l'oeuvre d'universitaires japonais. Leurs résultats ont été publiés dans des revues scientifiques, mais aussi dans des colloques sur Grass, et donc disponibles sur Internet:

• Masumoto, S. , Raghavan V., Nemoto, T., Shiono K.,(2002), Construction and Visualization of Three Dimensional Geologic Model Using GRASS GIS. In M. Ciolli and P. Zatelli, editors, Open Source Free Software GIS - GRASS users conference 2002 -> www.ing.unitn.it/~grass/conferences/GRASS2002/proceedings/proceedings/pdfs/Masumoto_Shinji.pdf
• Yonezawa , G., Nemoto T., Masumoto S., and Shiono K., (2002) 3-D Geologic Modeling and Visualization of Faulted Structures. In M. Ciolli and P. Zatelli, editors,
  Open Source Free Software GIS - GRASS users conference 2002 -> www.ing.unitn.it/~grass/conferences/GRASS2002/proceedings/proceedings/pdfs/Yonezawa_Go.pdf
• Kajiyama , A., Ikawa, N., Masumoto S., Shiono K., Raghavan V., (2004), Three-dimensional Geological Modeling  by FOSS GRASS GIS –Using Some Field Survey Data, Proceedings of the FOSS/GRASS Users Conference - Bangkok, Thailand, 12-14 September 2004 -> www.gisdevelopment.net/proceedings/fossgrass/papers/Three-dimensional%20Geological%20Modeling.pdf
• Masumoto, S., Nonogaki,S., Nemoto, T, Shiono, K.,Raghavan,(2006), Development of Visualization Tool for Geologic Information on GRASS GIS, GeoInformatics for Spatial-Infrastructure Development in Earth and Allied Sciences (GIS-IDEAS) 2006, Vietnam ->wgrass.media.osaka-cu.ac.jp/gisideas06/viewpaper.php

Ils ont développé une technique et des modules spécialisés avec Grass seul pour des résultats assez extraordinaires, mais, malheureusement, ils ne désirent pas encore les mettre dans le domaine public (discussions sur la liste Grass-dev). Il ne nous reste donc malheureusement plus que nos yeux pour admirer...

Une deuxième solution, accessible

Une autre solution complète est développée par le géologue polonais Maciej Tomaszczyk sur son site opengeology.pl/ (en polonais), particulièrement dans la partie "Geologia 3D". Il combine Grass avec VTK et Paraview (opengeology.pl/index.php) et travaille tant avec les surfaces que les volumes. À l'inverse des précédents, sa démarche est accessible: je lui ai demandé de l'aide et il l'a fait spontanément pour mon plus grand bonheur. Il m'a aussi donné l'autorisation de reproduire ses figures, et je l'en remercie pour tout  (avec de nombreux courriels...).

 images obtenues avec Paraview à partir de données Grass, tirées de opengeology.pl/index.php

Les résultats précédents sont obtenus à partir de forages qui permettent d'obtenir les coordonnées x,y,z de limites de couches. Hormis l'interprétation volumétrique, de tels résultats pourraient être obtenus avec de nombreux autres SIG.

Mais comment faire si l'on n'en dispose pas ?  Maciej Tomaszczyk utilise alors des méthodes élaborées par le géologue canadien Eric de Kemp qui a publié dans la revue "Computers & Geosciences" entre, 1998 et 2000, trois articles qui  expliquent comment réaliser des modélisations géologiques 3D à partir de données de surface (affleurements et/ou cartes géologiques), les forages éventuels ne servant qu'au contrôle ou à une éventuelle amélioration postérieure du processus.

• Three-dimensional projection of curvilinear geological features through direction cosine interpolation of structural field observations, Computers & Geosciences, Volume 24, n° 3 (1998), pages 269-284 (www.sciencedirect.com/science);
  
• Visualization of complex geological structures using 3-D Bézier construction tools, Computers & Geosciences, Volume 25, n° 5 (1999), pages 581-597 (www.sciencedirect.com/science);
  
• 3-D visualization of structural field data: examples from the Archean Caopatina Formation Abitibi greenstone belt, Québec, Canada, Computers & Geosciences,Volume 26, n° 5, (2000), pages 509-530 (www.sciencedirect.com/science).
  

Ces articles et une synthèse sont aussi présentés dans sa thèse de Doctorat défendue à l'UQAC (Université du Québec à Chicoutimi) en 2000:

• "Intégration et visualisation 3-D de données structurales de terrain: outils pour la cartographie géologique régionale"(en anglais)
  

et disponible sur le site de l'Université du Québec à Chicoutimi à bibvir.uqac.ca/theses/12127333/12127333.pdf (196 pages).

Eric de Kemp propose diverses stratégies (et algorithmes) suivant l'échelle d'étude choisie (d'une zone sur une carte géologique à l'ensemble d'un bassin géologique) et les applique à des exemples variés. Je me suis donc lancé, avec les explications envoyées par Maciej Tomaszczyk, dans une application "simple" à partir d'un extrait de carte géologique avec le premier algorithme proposé pour des petites surfaces. Je ne rentrerai pas ici dans le détail des nombreux problèmes géologiques posés par la démarche pour me focaliser uniquement sur son application avec Grass et Paraview, à coup d'essais et d'erreurs (nombreuses !) car cela m'a permis de mieux comprendre les divers principes de la modélisation 3D. Pour toutes les explications géologiques et la résolution des problèmes posés, je renvoie le lecteur aux articles ou à la thèse d'Eric de Kemp. Ce n'est en effet pas le lieu sur ce portail qui n'est pas spécialisé en géologie.

Quelques petits rappels s'imposent néanmoins. Pour un géologue, toute couche géologique peut être assimilée, à l'endroit où sont prises les mesures, à un plan dont on connait la direction et la pente par rapport à l'horizontale (le pendage).

Les mesures effectuées sont ponctuelles, pendage, direction et les coordonnées x, y et z, si l'on dispose d'un outil fiable de mesure. La méthode des cosinus directeurs, proposée par Eric de Kemp,  permet de générer des points supplémentaires avec des x, y et z calculés en les projetant suivant la ligne de pente avec une distance et des intervalles choisis:

 

Le processus est répété pour chaque point examiné:

Une surface 3D de la couche est ensuite générée (à partir de tous ces points).  Il est évident que le résultat final sera déterminé par le nombre de points de mesure. Mais au vu des structures géologiques, il est aussi évident que l'on ne peut pas interpoler ces lignes de pente à l'infini et choisir les points de mesure au hasard (tous ces thèmes sont traités dans les travaux d'Eric de Kemp, qui en propose diverses améliorations, dont l'utilisation de courbes de Bézier programmées, ou un algorithme de recherche automatique de points situés au voisinage d'une limite géologique pour pouvoir les utiliser, par exemple).

Un exemple d'application: démarche et problèmes rencontrés

Nous allons appliquer cette méthode sur un extrait d'une carte géologique à 1/25000. Le MNT utilisé sera le SRTM 3 de la Nasa. Les résultats obtenus sont bien loin de ceux obtenus par Maciej Tomaszczyk, mais je progresse peu à peu et ils permettent de bien expliquer comment appliquer la méthode.

Extrait de la carte géologique, la barre plus noire au NW (coin supérieur gauche) indique l'échelle: elle représente 3 km

Il s'agit d'une structure plissée avec deux failles inverses au SE et une discordance au NW. Cet extrait est publié avec l'autorisation de la Région wallonne et du programme "Carte géologique de Wallonie" auquel je participe. La carte complète peut être visualisée en ligne à carto1.wallonie.be/geologie/viewer.htm et le site est environnement.wallonie.be/cartosig/cartegeologique/index.htm. Le module Nviz est utilisé pour la visualisation en 3D (et Paraview pour trois figures).

Carte géologique drapée sur le MNT (exagération verticale x 3)

Nous allons nous intéresser à la formation géologique marquée par un fort relief. Les deux images précédentes donnent une idée de l'échelle générale de la zone étudiée. Nous n'en spécifierons plus dans la suite, car  les meilleures illustrations possible pour illustrer le propos ont été choisies, sans trop se préoccuper de l'échelle à ce stade (l'échelle sera donnée dans la figure finale).

La première démarche est de "récolter" les points de mesure à la limite de la couche. Leurs altitudes (z) sont ici obtenues automatiquement à partir du MNT (commande v.drape, valeurs z interpolées à partir d'un MNT). Au vu de la structure géologique (pendages de 65 à 90°) et des coupes géologiques, les points seront projetés sur les lignes de pente sur une distance raisonnable de 300 m.

points projetés à partir des mesures sur la carte géologique

 

points  (vecteur 3D) utilisés pour créer la surface

La démarche suivante est la plus complexe, car Grass permet de jouer avec de nombreux paramètres difficiles à appréhender de prime abord. Il faut véritablement les tester avec de nombreuses tentatives avant de trouver une solution acceptable. Les modules utilisés sont v.surf.rst pour générer la surface puis r.mapcalc pour restreindre la surface aux points. 

   résultat brut de v.surf.rst                   suppression de la partie supérieure          suppression de la partie inférieure

La surface brute obtenue déborde largement les limites des points, mais ce serait la même chose si l'on appliquait un algorithme comme celui de Delaunay (ici avec Paraview).

Dans le cas qui nous concerne, les paramètres ont été fixés pour maximiser la surface limitée aux points en donnant deux fois plus d'importance à l'axe d'élongation (directions d'ennoyage des plis). Cela veut dire que toute la surface située en dehors des points est virtuelle et ne peut pas être utilisée. Celle-ci doit donc être "nettoyée" à l'aide du module r.mapcalc, en deux étapes (suppression des parties supérieure et inférieure). J'ai laissé des "dentelures" dans la partie supérieure de la couche pour montrer l'importance du rôle des paramètres utilisés. La partie inférieure de la couche est plus "correcte" puisque le traitement a éliminé les dentelures. Le résultat final est un véritable raster 3D (et non une surface vectorielle) qui peut donc être soumis à tous les traitements de rasters qui sont une des spécialités de Grass.

Ces résultats peuvent aussi être exportés au format VTK (v.out.vtk) et traités avec Paraview. Le traitement de ce type de données se fait avec des enchainements d'algorithmes. Les résultats sont représentés dans les figures suivantes (couche seule puis traitement de deux limites de couches et celle d'une des failles (en rouge):

Ces résultats peuvent être réimportés dans Grass et visualisés et traités:

 couche, faille et points utilisés, placés sur le MNT (légèrement transparent)

avec la carte géologique elle-même

et enfin avec les limites des couches et une échelle, sans exagération verticale ( x 1)

Ces résultats peuvent déjà paraître très intéressants, mais je suis encore loin du compte. En effet, Maciej Tomaszczyk va encore plus loin en créant des volumes (voxels 3D) à partir des diverses surfaces créées et là, j'avoue que je coince un peu, même si je comprends le principe (volume entre les couches jaune et verte, par exemple).

 

volumes visualisés avec Paraview (image tirée de opengeology.pl/index.php)

volumes visualisés avec Nviz (image tirée de opengeology.pl/index.php )

Pour un géologue, tant Paraview que Grass permettent de sectionner les volumes obtenus suivant diverses directions et divers angles et obtenir de véritables coupes géologiques. Il est aussi possible de placer des forages éventuels (en 3D) pour contrôler la validité du modèle (limites des couches géologiques).

Conclusions: suite ?

Et non, la quête ne s'arrête pas ici, même si les résultats déjà obtenus sont plus qu'intéressants, car pour l'instant ce n'est que de la visualisation 3D de structures géologiques et non de volumes comme ceux obtenus précédemment ou par les géologues japonais dans la première solution.

C'est la démarche qui me parait la plus passionnante, car elle permet de "pousser" Grass dans ces derniers retranchements et d'apprendre toutes ses subtilités. En pratique, il est vrai, en râlant beaucoup (sur moi et sur Grass ou Paraview), et en travaillant peu à peu, en soirée et le WE, puisque c'est en plus de mon travail et pour mon plaisir (tant Grass que Paraview sont "scriptables" en Python).

Quand le virus est là, on ne peut plus s'arrêter et on peut ainsi reprendre d'anciens travaux comme la modélisation des failles en forme de "coques de bateau" de la structure géologique connue en Belgique sous le nom de "Massif de la Tombe» (voir www2.ulg.ac.be/geolsed/geolwal/geolwal.htm), à partir de leurs isohypses publiées sur la feuille "Fontaine-l’Evêque - Charleroi”, feuille 46/7-8 de la nouvelle carte géologique de Wallonie (DELCAMBRE, B. & PINGOT, J.L., 2000).C'est avec ce projet que tout à commencé.

Modélisation Grass-Paraview à partir des données résultantes du travail de Louis Franssen. Les éléments représentés sont les diverses failles de la structure, la rouge étant la "Faille du Midi"

Il y aurait encore tant de choses à faire, comme essayer d'appliquer les autres algorithmes d'Eric de Kemp ou arriver au niveau de Maciej Tomaszczyk, mais enfin, pour paraphraser un autre film des Monty Python : "Always look on the bright side of life".

Une remarque, il est nécessaire de suivre régulièrement ou de fouiller la liste Grass-dev pour trouver ce que l'on veut puisque c'est là que ça se passe. Je fournis ici les liens vers quelques posts intéressants qui m'ont beaucoup servi (il y en a beaucoup d'autres, en recherchant "3D geology"):

• n2.nabble.com/Searching-Docs-about-3D-geological-modelisation-td4269207.html#a4269207
• n2.nabble.com/Plot-3D-Drillholes-in-Grass-GIS-td2457531.html#a2458123
• n2.nabble.com/GRASS-user-GRASSLIST-1105-tin-to-raster-dem-td1840009.html#a1840013
• n2.nabble.com/GRASS-user-voxel-from-dems-td1839294.html#a1839296

Et enfin, signalons qu'une version spéciale de Paraview a été développée spécialement pour les applications géologiques, ParaviewGeo, Open Source, mais, malheureusement, seule une version Windows est disponible. Elle permet l'import ou l'export des données Gocad et autres (paraviewgeo.mirarco.org/index.php/Supported_Formats#Supported_Readers) et offre des algorithmes de traitement plus adaptés au domaine géologique ou minier.

 

Tous les traitements ont été effectués sur Mac OS X avec Grass 6.4 RC4 de www.kyngchaos.com/ , diverses versions de Paraview, suivant leur mise à jour, le script en awk d'Eric de Kemp "pusch.awk" pour les cosinus directeurs, et Inkscape pour les figures.

Site officiel : Gocad
Site officiel : Paraview
Site officiel : Mayavi
Site officiel : ParaviewGeo
Autres Liens : Courbe de Bézier
Autres Liens : TIN: Triangulated irregular network