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)

	Niveau	Intermédiaire
	Logiciels utilisés	QGIS 2.0 matplotlib (module Python)
	Plateforme	Windows | Mac | Linux | FreeBSD

Pour les géologues, créer une coupe géologique à partir d'un logiciel SIG classique est un problème récurrent, surtout lorsque il ne dispose pas d'une application qui gère la 3D (comme GRASS GIS ou ArcGIS avec des extensions très chères). Il s'agit de:

1) créer un profil topographique en fonction d'une ligne de coupe définie et d'un MNT.

• il y a des extensions/plugins QGIS pour le faire comme Profile Tool (qui utilise le module Python matplotlib, comme d'autres extensions);
• il y a des commandes avec GRASS GIS ou SAGA GIS seuls ou dans la Toolbox du menu « Traitements » (ex extension Sextante, maintenant intégrée dans QGIS);
• il y a bien entendu moyen de combiner divers logiciels;
• cela n'aura jamais la souplesse désirée alors qu'il y a moyen de faire exactement ce que l'on veut en Python (voir plus bas).

2) colorier le profil en fonction des couleurs d'une carte et/ou adjoindre les points d'intersection entre les limites géologiques et cette ligne.

Cela a été présenté sur le Portail dans:

• R dans Grass Gis et R: superposition des couleurs d'un raster quelconque sur un profil topographique élaboré à partir d'un MNT;
• Python seul et les modules géospatiaux dans Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG et Python: traitement des couches vectorielles dans une perspective géologique, lecture et enregistrement des couches sous forme de dictionnaires avec le module Fiona;
• et repris sur GIS StackExchange dans Python Script for getting elevation difference between two points.

3) reporter les valeurs de pendages apparents des couches géologiques sur l'axe de la coupe (calculés en fonction de l'angle entre la direction de la coupe  et la direction de la couche, et éventuellement de l'exagération des hauteurs).

• cela a été présenté sur le Portail en Python dans les articles cités

4) dessiner la coupe par interprétation géologique.

• il est toujours plus rapide de le faire à la main;
• cela est presque impossible à mettre en oeuvre sans un logiciel spécialisé;
• j'essaie de le faire par la méthode des kinks (Construct a Fold Cross-Section Using the Kink Method) mais c'est une espèce de défi personnel...
  
  

  

5) utiliser la méthode de « down-plunge profile »

figure reprise de Steve Dutch dans Construct a Down-Plunge Profile

• cela a été présenté sur le Portail en Python dans les articles cités.

Il m'a souvent été demandé comment porter les scripts Python pour le faire directement dans QGIS (collègues ou sur le ForumSIG dans [QGIS 2.x] Coupe geologique avec qgis2). J'avais utilisé les modules Python osgeo (GDAL/OGR), Shapely et Fiona dans les scripts originaux du fait des « lacunes » de la précédente version de PyQGIS et d'un manque de souplesse. Mais depuis le nouvel API (Qgis 2 : changement de l'API python) la tâche a été grandement facilitée puisque certaines fonctions ont été ajoutées et les traitement simplifiés (sans la « pollution » de PyQt4).

Si vous n'aimez pas la console Python et n'avez pas le courage de vous atteler à la confection d'une extension, le nouveau menu « Traitements »  permet de créer très facilement des scripts dans la Toolbox avec une interface simple et de les exécuter avec un clic:

Principes

Tout d'abord le script de traitement des couches complet, très court, et qui sera détaillé dans la suite, sans:

• la manière de sélectionner les données (4 couches suffisent, une ligne de coupe, les limites des couches géologiques, un MNT et un fichier raster: une carte géologique), qu'il suffit d'ajouter en fonction d'un traitement dans la console ou de la création d'un script pour la Toolbox;
• la création des figures, qu'il suffit d'ajouter en fonction de vos besoins.

# traitement de la ligne de coupe et intégration des segments éventuels (=UNION)
geomcoupe = QgsGeometry.fromWkt('GEOMETRYCOLLECTION EMPTY')
for elem in coupeori.getFeatures():
      geomcoupe = geomcoupe.combine(elem.geometry())
# longueur de la ligne de coupe
longueur=geomcoupe.length()
# traitement des limites des couches géologiques (polylines) 
limites = QgsGeometry.fromWkt('GEOMETRYCOLLECTION EMPTY')
for elem in lim.getFeatures():
    limites = limites.combine(elem.geometry())
 
# valeur d'un pixel sous un point
def Val_raster(point,raster):
      return raster.dataProvider().identify(point, QgsRaster.IdentifyFormatValue).results().values()
 
# calcul des des valeurs sur la ligne de coupe
x,y,z,couleur,dista = map(list,zip(*[(geomcoupe.interpolate(distance).asPoint().x(), geomcoupe.interpolate(distance).asPoint().y(),Val_raster(geomcoupe.interpolate(distance).asPoint(),MNT)[0], Val_raster(geomcoupe.interpolate(distance).asPoint(),carte_geol),distance)  for distance in xrange(0,longueur,20)]))
 
# intersection entre la la ligne de coupe et les limites des coupes
ptintersec = geomcoupe.intersection(limites)
# calcul des valeurs des points d'intersection sur la ligne de coupe
origine = geomcoupe.interpolate(0) 
xd,yd,zd,distd = map(list,zip(*[(point.x(),point.y(),Val_raster(point,MNT)[0],QgsGeometry.fromPoint(point).distance(origine))  for point in ptintersec.asMultiPoint()]))

Pour ce faire, nous reprendrons le plan général de la version en Python seul, de la manière la plus simple possible et sans toutes les fioritures qui complexifient la tâche dans les tutoriels sur PyQGIS. Comme vous pouvez le voir, il y a moyen de le faire en un minimum de lignes de code avec les fonctions de programmation fonctionnelle comme map() ou les compréhensions de liste (voir La compréhension de liste en Python, une syntaxe moderne pour map() et filter()). Quelques exemples de traitements, utilisés dans le script présenté, seront signalés en italique.

Comme je n'ai pas développé l'aspect graphique (avec le module matplotlib) dans les scripts originaux, il sera, cette fois-ci, détaillé.

Accéder aux couches vectorielles ou matricielles

• dans la console Python

Le parcours des éléments (géométrie et attributs) est direct, sans sélection préalable:

macouche = qgis.utils.iface.activeLayer()
for elem in macouche.getFeatures():
     geom= elem.geometry()
     attr =elem.attributes()
      (traitements)

Mais il est pénible de devoir sélectionner à chaque fois une couche active pour y accéder dans le script,  il est beaucoup plus pratique d'utiliser un dictionnaire des couches activées sur le canevas:

canvas= qgis.utils.iface.mapCanvas()
# sélection des couches affichées
couches = dict((k.name(),i) for (i, k) in enumerate(canvas.layers()))
print couches
{u'ligne_coupe': 1, u'lim_form': 0, u'MNT': 2, u'carte_geol': 3}

Ceci permet de sélectionner directement une couche par son nom:

coupeori = canvas.layer(couches['ligne_coupe'])
# ou/au lieu de
coupeori = canvas.layer(1)

La grande différence réside dans le fait que dans le cas de l'API 1.x, le résultat est un itérable (liste, par exemple) alors que dans l'API 2.x, c'est un générateur/itérateur (comme dans Fiona) Quelle est la différence ?

• l'entièreté de la liste doit être chargée en mémoire avant de la parcourir avec une boucle for, par exemple;
• dans un générateur, les valeurs sont générées une par une (c'est la seule chose en mémoire, une valeur) et une fois affichée, cette valeur n'existe plus. C'est ce qui permet de traiter les grandes masses de données (voir par exemple Comment utiliser yield et les générateurs en Python ?).

Il est donc parfaitement possible de faire:

coupeori = canvas.layer(couches['ligne_coupe'])
elem = coupeori.getFeatures().next()
geomcoupe= elem.geometry()
....

et d'appliquer la geo_interface, voir dans Les modules Python à finalités géospatiales: quid, quando, ubi ?

• dans un script de Traitements (processing, ex Sextante)

Il y a ici moyen de créer une interface simple qui permet le choix des couches à traiter et de les récupérer dans le script:

# interface
##ligne_coupe=vector
##lim_form= vector
##MNT=raster
##carte_geol= raster
#Traitements
from qgis.core import *
coupeori = processing.getobject(ligne_coupe)
lim = processing.getobject(lim_form)
MNT = processing.getobject(MNT)
carte_geol = processing.getobject(carte_geol) 

(Attention, dans les versions master 2.1, il s'agit de .getObject)

Ceci permet, lorsqu'on clique sur le bouton d'obtenir le dialogue suivant:

Comment extraire la valeur d'un pixel d'un raster sous un point

PyQGIS offre maintenant une solution beaucoup plus simple que celle vue avec les autres modules pour extraire la valeur d'un pixel sous un point:

# nombre de bandes
MNT.bandCount()
1L
carte_geol.geol.bandCount()
3L
# valeur d'un pixel sous un point
MNT.dataProvider().identify(QgsPoint(229774,111171), QgsRaster.IdentifyFormatValue).results()
{1: 221.0}
carte_geol.dataProvider().identify(QgsPoint(229774,111171), QgsRaster.IdentifyFormatValue).results()
{1: 180.0, 2: 142.0, 3: 125.0}

Le résultat obtenu est un dictionnaire avec comme clé, la bande, et comme valeur, l'élement recherché (altitude dans le cas d'un MNT, couleur RGB dans le cas d'un couche multibandes). Il est donc très facile de réécrire la fonction  Val_raster() établie dans les scripts originaux.

def Val_raster(point,raster):
      return raster.dataProvider().identify(point, QgsRaster.IdentifyFormatValue).results().values()

Comment croiser un fichier de type vectoriel et un fichier de type matriciel pour en extraire les données.

Nous reprenons le script présenté dans les articles cités. Le nouvel API de PyQGIS dispose maintenant, lui aussi, de la fonction interpolate pour générer des points équidistants sur une ligne (il n'y donc plus besoin du module Shapely pour cet aspect). Plus le nombre de points est élévé, plus la représentation sera précise.

# traitement de la ligne de coupe
# création d'une couche vide pour regrouper tous les segments
# éventuels d'une couche ligne (= UNION)
geomcoupe = QgsGeometry.fromWkt('GEOMETRYCOLLECTION EMPTY')
# intégration des segments éventuels (=UNION)
for elem in coupeori.getFeatures():
      geomcoupe = geomcoupe.combine(elem.geometry())
#-------------------------------------
# puisque c'est un générateur, et dans le cas d'une ligne simple,sans segment,
# il n'y a pas besoin de boucle
elem = coupe.getFeatures().next()
geomcoupe= elem.geometry()
longueur = geomcoupe.length()
#-------------------------------------
# création des points équidistants sur la ligne avec un pas de 20 m
longueur=geomcoupe.length()
x = []
y = []
z = []
# couleur RGB
couleur = []
# distance pour le profil topographique
dista = []
for distance in xrange(0,longueur,20):
         # création du point sur la ligne 
         point = geomcoupe.interpolate(distance)
         point = point.asPoint() 
         xp,yp=point.x(), point.y()
         x.append(xp)
         y.append(yp)
         # extraction de la valeur altitude à partir du MNT
         z.append(Val_raster(point,MNT)[0])
         # extraction des couleurs RGB à partir de la couche raster
         couleur.append(Val_raster(point,carte_geol))
         dista.append(distance)

Ce qui en programmation fonctionnelle peut se faire en une seule ligne:

x, y,z, couleur,dista = map(list,zip(*[(geomcoupe.interpolate(distance).asPoint().x(), geomcoupe.interpolate(distance).asPoint().y(),Val_raster(geomcoupe.interpolate(distance).asPoint(),MNT)[0], Val_raster(geomcoupe.interpolate(distance).asPoint(),carte_geol),distance)  for distance in xrange(0,longueur,20)]))

Comment obtenir les points d'intersection entre 2 fichiers vectoriels (profil et limites des objets géologiques) et les placer sur la coupe géologique.

L'API 2.x ne dispose pas de la possibilité de calculer directement l'intersection entre 2 couches vectorielles, mais comme dans le script original, en transformant les 2 couches en 2 objets géométriques uniques, il est possible d'utiliser la fonction intersection.

# ouverture de la couche 
limites = QgsGeometry.fromWkt('GEOMETRYCOLLECTION EMPTY')
for elem in lim.getFeatures():
     limites = limites.combine(elem.geometry())
# points d'intersection entre les 2 couches
ptintersec = geomcoupe.intersection(limites)

Il faut maintenant les placer sur la coupe géologique (distance, z). Cela se fait en calculant la distance des points à partir du point d'origine de la coupe  avec la fonction point.distance(point) (puisqu'ils sont sur la ligne de coupe).

# origine de la coupe 
origine = geomcoupe.interpolate(0) 
# x,y si représentation 3D souhaitée 
xd = [] 
yd = [] 
zd = []    
distd = []       
for point in ptintersec.asMultiPoint():    
         xd.append(point.x())    
         yd.append(point.y())    
         zd.append(Val_raster(point,MNT)[0])    
         distd.append(QgsGeometry.fromPoint(point).distance(origine))

Ce qui en programmation fonctionnelle peut se faire en une seule ligne:

xd, yd,zd,distd = map(list,zip(*[(point.x(),point.y(),Val_raster(point,MNT)[0],QgsGeometry.fromPoint(point).distance(origine))  for point in ptintersec.asMultiPoint()]))

Comment représenter les coupes avec le module matplotlib

C'est ici, apparemment, que se posent le plus de problèmes pour ceux qui ont essayé le portage. Il est vrai que matplotlib est un monde en soi qui permet de tout faire, mais qui n'a rien à voir avec les SIGs. C'est de la représentation graphique pure et il faut bien connaître Python.

Pour la représentation de la ligne de coupe en 3D:

C'est le plus direct, car le script fournit les valeurs x, y et z des points générés le long de la coupe:

from mpl_toolkits.mplot3d.axes3d import * 
import matplotlib.pyplot as plt 
fig = plt.figure() 
ax = Axes3D(fig) 
# ligne reliant les points, de couleur rouge  
ax.plot3D(x,y,z, c='r') 
plt.show() 

Résultat (il y a moyen de faire tourner le diagramme et de le sauver sous divers formats, dont le SVG, vectoriel):

Si l'on veut utiliser les couleurs de la carte géologique, matplotlib permet de le faire à partir de tuples avec les valeurs R, G, B de chaque point. C'est ce qui a été obtenu dans la liste couleurs, mais matplotlib utilise des valeurs de 0 à 1 et donc:

color = [[j/255.0 for j  in i] for i in couleur] 
fig = plt.figure() 
ax = Axes3D(fig) 
points 3D ici 
ax.scatter3D(x,y,z, c = tuple(color)) 
plt.show() 

Résultat:

Avec matplotlib, il n'est pas encore possible de bien contrôler le rapport entre l'échelle des x et celle des y en 3D. Il faut alors passer par un module comme visvis.

import visvis 
f = visvis.gca() 
m = visvis.plot(x,y,z, lc='k', ls='', mc='g', mw=2, lw=2, ms='.') 
# z x 1 
f.daspect = 1,1,1 
# z x 10
f.daspect = 1,1,10

Pour la coupe 2D:

Ici c'est le couple (distance, z) qui est pris en compte. Pour utiliser les couleurs du raster dans la figuration d'une ligne, je vais donner une couleur à chaque segment basée sur la couleur d'un de ses points limite.

Je le ferai ici avec toutes les enjolivures que je vous laisse découvrir, c'est comme ça que j'ai appris,  à coup d'essais et d'erreurs.

fig = plt.figure() 
ax = fig.gca() 
# dessin des segments entre chaque point 
for i in range(len(dista)-1): 
      ax.plot([dista[i],dista[i+1]],[z[i],z[i+1]], c=tuple(color[i+1]),linewidth=2) 
   
# enjolivures 
# échelle de l'axe y x 10 
ax.set_aspect(10) 
xrange = [0, 7500] 
yrange = [200, 400] 
ax.set_xlim(*xrange) 
ax.set_ylim(*yrange) 
ax.set_yticks([200,300,400]) 
ax.set_xlabel('Distance (m)') 
ax.set_ylabel('Altitude x 10 (m)') 
ax.set_title('N                                                                    S') 
plt.show() 

Résultat:

Pour placer les intersections avec les limites des formations géologiques

Il suffit de rajouter au script précédent la ligne:

ax.scatter(distd,zd,s=200, marker='|') 

qui permettra de dessiner une ligne verticale à la place du point.

Pour placer les angles des pendages

C'est un peu plus complexe et je ne présenterai ici que les principes:

• tout comme dans les scripts originaux, il faut aller chercher l'angle du pendage dans la table d'attributs;
• le pendage apparent doit être calculé en fonction de l'angle entre la direction de la couche (table d'attributs aussi) et l'azimut de la ligne de coupe (calcul déjà présenté sur le Portail) et éventuellement de l'exagération des hauteurs;

Ensuite on utilise le principe de représentation des angles de matplotlib:

angles = [0,10,20,30,40,50,60,70,80,90]
ax = plt.gca()
ax.plot(x,y)
for i in range(len(x)):
      ax.scatter(x[i],y[i],s=300, marker=(2,0,(angles[i]))) 

Résultat final:

Conclusions

Vous savez maintenant:

• ajouter au début du script traitement_couches.py, la sélection des couches nécessaires, en fonction d'une utilisation dans la console ou de la création d'un script pour la Toolbox;
• ajouter à la fin du script, la figuration souhaitée  en fonction de vos besoins.

Quels sont les avantages de cette démarche ?

• le traitement peut être effectué directement dans QGIS, soit dans la console, soit avec un script placé dans la Toolbox du menu « Traitements ».

Quels sont les désavantages de cette démarche  par rapport à un script en Python seul comme ceux présentés auparavant ?

• la limitation du traitement aux couches activées sur le canevas;
• les couches provenant d'un WMS ou de GRASS GIS par l'extension GRASS ne sont pas prises en compte lorqu'il s'agit d'extraire une valeur (aucune valeur dans mon cas);
• vous ne pouvez donc pas aller « piocher » où vous voulez, sans tenir compte de l'origine, comme avec les scripts en Python seul.

Que faudrait-il améliorer/ajouter (surtout dans le cas d'un script dans la Toolbox) ?

• la gestion des erreurs comme lorsqu'une mauvaise couche a éte choisie, ce qui provoquera une erreur;
• la possibilité de choisir divers paramètres et le type de représentation désiré;
• et donc une standardisation de la table d'attributs (pour les valeurs de pendage et de direction);
• ....

Tous les traitements ont été faits sur Mac OS X avec la version QGIS Dufour de KyngChaos et testés avec succès sur:

• la versions QGIS Dufour standalone pour Windows;
• la version QGIS Dufour pour Ubuntu;
• la version QGIS 2.1 master (attention, il s'agit de .getObject);
• le module matplotlib fait « planter » QGIS dans la version script Toolbox sur Mac OS X, pas dans la version console (dans mon cas).

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