QGIS et géologie structurale: création automatique de rosaces directionnelles et de canevas stéréographiques avec PyQGIS

	Niveau	Expert
	Logiciels utilisés	Qgis
	Plateforme	Windows | Mac | Linux | FreeBSD

 

La géologie structurale traite des relations entre des objets spatiaux en trois dimensions:

• une carte géologique n'est que la représentation 2D de l'intersection de ces structures géologiques et d'un plan horizontal (mais un géologue la voit toujours en 3D). C'est ce que l'on peut visualiser dans un SIG;

• cette 3D est approchée par les coupes géologiques (intersection avec un plan vertical ou incliné, comme les falaises ou les parois de montagne) ou les modélisations géologiques 3D. Il faut alors utiliser des logiciels de modélisation géologique spécialisés ou éventuellement un SIG comme dans Grass Gis et Paraview: la modélisation géologique 3D pour les démunis, ou les enthousiastes... ou PyQGIS (QGIS 2): les coupes géologiques (colorisation d'un profil topographique élaboré à partir d'un MNT en fonction des couleurs d'un raster et placement des points d'intersection avec les limites des couches géologiques) sur le Portail.

figure tirée de Martin-Luther-Univestät Halle-Wittenberg: 3D Modelling Geology and Mining: 3D geomodeller

Mais il est aussi important d'analyser les rapports d'orientation de ces éléments dans l'espace, indépendamment de leur position géographique (modèle d'orientation préférentielle des couches, des fractures, ou de n'importe quel objet linéaire sur une carte pour en extraire des classes, des tendances, des relations).

Dans ce système, chaque objet est défini par:

1. sa direction (ou azimut)
2. son pendage (ou plongement)

Suivant ces critères, il est  possible d'utiliser:

• les histogrammes et les diagrammes de direction ou rosaces directionnelles (« rose diagram ») pour analyser les caractéristiques d'une variable (pendage ou direction). Dans certains cas comme l'analyse des orientations des linéaments sur photo aérienne, c'est le seul moyen puisque l'on ne connait que leurs directions;

• les canevas stéreographiques (« stereonets ») qui permettent d'analyser les deux variables en même temps (un objet est défini par son pendage et sa direction).

Comme le but ici n'est pas un cours de géologie, je laisse à d'autres le soin d'expliquer ces méthodes pour me focaliser sur leurs implémentations directes dans QGIS avec Python.

Objectif

Il existe déjà beaucoup d'applications, libres, gratuites ou payantes, pour effectuer ces traitements. Il en existe même certaines qui permettent de placer directement les points de mesure sur un fond Google Map, comme Stereonet de Rick Allmendinger, ou d'extraire les valeurs structurales à partir d'une carte géologique sous format image, comme GeolMapDataExtractor, du même auteur (programmes gratuits pour une utilisation non commerciale, versions Windows et Mac OS X).

Mais le problème de toutes ces solutions est que pour les utiliser à partir d'un SIG, il est nécessaire de passer par un fichier intermédiaire (du SIG vers le logiciel). Il existe néanmoins quelques solutions natives pour ArcGIS, sous forme d'extension, comme Orientation Analysis Tools for GIS (ArcGIS).

Le but ici est donc d'obtenir les résultats/graphiques directement dans QGIS, et ce,  sans logiciel tiers.

Je n'ai pas le temps de réinventer la roue et il existe heureusement des solutions en Python, prêtes à l'emploi:

• Rose diagram code pour le traitement des diagrammes en rosaces géologiques;
• mplstereonet pour le traitement des projections stéréographiques (stereonet). Il dispose nativement de beaucoup de traitements (voir les exemples).

Elles utilisent toutes les deux le module Python matplotlib pour les figures. Il y a donc moyen de sauver les figures obtenues dans divers formats dont le format vectoriel SVG.

Dans les deux cas, je me suis contenté d'adapter quelque peu les scripts originaux et de rajouter quelques traitements en utilisant mes versions en Python des scripts originaux pour Matlab publiés dans Structural Geology Algorithms Vectors and Tensors de R.W. Allmendinger, N.Cardozo, et D.M. Fisher (2012) et disponibles dans la partie ressources du site (versions Python réalisées avec leur accord).

Résultats

Mes premières réalisations ont été directement effectuées dans la console Python de QGIS ou avec un script externe à l'aide du plugin ScriptRunner (il y a des exemples des résultats obtenus dans QGIS : lancer des scripts Python ou des commandes Shell depuis la Console Python ou avec Script Runner (sans créer d'extension) et Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG.

Le défaut principal de ces solutions est qu'elles ne sont pas interactives. Je voulais:

1. pouvoir traiter directement des éléments sur la carte, en en sélectionnant certains ou dans leur totalité;
2. cliquer sur un simple élément pour obtenir les résultats
3. recommencer autant de fois que je le désire avec d'autres paramètres (en modifiant rapidement les scripts, si nécessaire).

La solution idéale était donc de créer des scripts Python dans la Boîte à outil de Traitements (Processing): il suffit de cliquer dessus, et l'interface est créée automatiquement (voir plus bas, similaire à la création de modules dans GRASS GIS).

1) rosaces géologiques

Interface obtenue avec le module processing: (angle = intervalle angulaire pour créer le diagramme)

Diagramme des directions de toutes les failles représentées sur une carte:

Résultats avec quelques failles sélectionnées:

Diagramme des directions de quelques mesures de schistosité sélectionnées:

2) canevas stéréographiques

interface obtenue avec le module processing:

Projections stéréographiques des pôles des plans de stratification d'un ensemble de points choisis sur la carte (canevas de Schmidt, hémisphère inférieur, il y a aussi moyen d'utiliser le canevas de Wullf):

Pôles et grands cercles de quelques mesures de stratification sélectionnées sur la carte:

Contour de quelques mesures de stratification sélectionnées sur la carte (plusieurs algorithmes sont disponibles, Khamb etc.) :

 

Estimation de l'axe d'un pli  :

J'utilise ici ma version Bingham.py du programme Bingham.m pour Matlab (Structural Geology Algorithms Vectors and Tensors, pp;93-97) .

 

Réalisation:

Création de l'interface:

Le simple en-tête d'un script:

#Interface
#==================================
##[Mes scripts GEOL]=group
##entree=vector
##entree2=raster
##dip_dir=field entree
##dip=field entree
##angle=number 10
##contour=number 0
##sortie=output vector

fournira l'interface suivante:

Il y a moyen de définir (voir Using processing algorithms)

• une ou plusieurs couches en entrées: vector1=vector, vector2=vector, raster1 = raster;
• une ou plusieurs couches en sorties: vector3 = vector output, raster3 = raster output;
• plusieurs champs à sélectionner/choisir;
• des valeurs numériques, alphanumériques, etc.

Obtention des données à partir des couches vectorielles:

Les données sont de deux types:

• les valeurs d'un champ (ou de plusieurs champs) qui doivent être extraites de la table attributaire (pendages et les directions, par exemple);
• des calculs à partir des coordonnées géographiques (les valeurs de direction/azimut, par exemple).

De plus il faut pouvoir choisir entre:

• traiter toutes les données;
• ne traiter que les valeurs sélectionnées.

1) obtenir l'azimut de tous les segments d'une couche polyligne

Ce calcul a été présenté dans PyQGIS: des géométries, des vecteurs, de l'algèbre vectorielle ou des cosinus directeurs, exemples d'applications  et GIS Stack Exchange: How do I find vector line bearing in QGIS or GRASS?. Il suffit d'utililiser la fonction Point1.azimuth(Point2) de PyQGIS.

Les résultats (angles entre 0 et 360°) se retrouvent sous forme de liste dans une variable.

2) extraction des valeurs des champs de la table choisis: exemple pour 2 champs, choix entre utiliser les valeurs sélectionnées ou la totalité des valeurs:

# Interface
##[Mes scripts GEOL]=group
##entree=vector
##dip_dir=field entree
##dip=field entree
# Traitements
from qgis.core import *
layer = processing.getObject(entree)
dip_dir_index = layer.fieldNameIndex(dip_dir)
dip_index = layer.fieldNameIndex(dip)
if layer.selectedFeatureCount():
    T,P = zip(*((elem.attributes()[dip_dir_index],elem.attributes()[dip_index]) for elem in layer.selectedFeatures()))
else:
    T,P = zip(*((elem.attributes()[dip_dir_index],elem.attributes()[dip_index]) for elem in layer.getFeatures()))

Les valeurs de direction se retrouvent dans la variable T, celles de pendage dans la variable P.

Traitements proprement dits

Il suffit d'importer et d'utiliser les modules cités (après en avoir compris les principes en se basant sur les exemples).

1) importation du premier module, NorthPolarAxes (rosaces géologiques): et traitements

import numpy as np
import matplotlib.pyplot as pl
from matplotlib.projections import register_projection
from NorthPolarAxes import NorthPolarAxes
register_projection(NorthPolarAxes)
# traitement de la variable T suivant les scripts d'exemples

En pratique, j'ai inséré la classe NorthPolarAxes dans mon script (au lieu de l'importer).

2) importation du deuxième module, mplstereonet et traitements

import matplotlib.pyplot as plt
import mplstereonet
fig = plt.figure()
ax = fig.add_subplot(111, projection='stereonet')
ax.pole(T, P, 'bo', markersize=5)
ax.grid()
plt.show()

Conclusions

Mon objectif de géologue utilsant un SIG dans une perspective géologique est atteint. Á partir de QGIS et de scripts dans la Boîte à outils de Traitements, je peux maintenant intéractivement:

• créer des modèles 3D: QGIS, représentation 3D des couches vectorielles (shapefiles dits 3D ou shapefiles avec attributs z) avec les modules Python Matplotlib ou Visvis à partir de la console Python
• entamer une coupe géologique: PyQGIS (QGIS 2): les coupes géologiques (colorisation d'un profil topographique élaboré à partir d'un MNT en fonction des couleurs d'un raster et placement des points d'intersection avec les limites des couches géologiques)
• créer des rosaces géologique ou des canevas stéréographiques à partir d'éléments choisis

Tous les traitements ont été effectués sur Mac OS X ici, mais aussi Linux et Windows avec plusieurs versions de QGIS 2.x

Site officiel : Visible Geology: Stereonet (en ligne)

licence Creative Commons Paternité-Pas d'Utilisation Commerciale-Pas de Modification 2.0 France