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

	Niveau	Intermédiaire
	Logiciels utilisés	Quantum GIS (QGIS) Shapely (module Python) matplotlib (module Python) visvis (module Python)
	Plateforme	Windows | Mac | Linux | FreeBSD

Comment utiliser avec profit dans Quantum GIS, des objets 3D comme les shapefiles dits 3D (PointZ, PolylineZ, PolygonZ,... ) ou les objets 2D qui disposent d'un attribut z. C'est la question qui a été reposée récemment sur la liste QGIS-Developer (osgeo-org.1560.n6.nabble.com/Re-Qgis-user-Using-visualizing-3D-data-Z-values-td5004272.html) et à laquelle j'ai répondu à osgeo-org.1560.n6.nabble.com/visualizing-3D-data-Z-values-or-data-with-z-attribute-a-solution-td5005360.html à l'aide d'autres modules Python que PyQGIS.

La représentation 3D dans QGIS est en effet une question récurrente sur les Forums et les listes de diffusion.  Les solutions proposées sont jusqu'à présent indirectes comme celles de passer par GRASS GIS, SpatiaLite ou PostGIS (www.faunalia.com/content/transfer-3d-shapefiles-z-values-table-attributes) en attendant l'extension Globe, encore en développement.

• Le cas des couches avec un attribut d'altitude est simple à traiter. L'extension Interpolation permet de les traiter et de générer des surfaces au format de grille ASCII ArcInfo (.asc) mais elle n'offre aucun moyen de visualiser le résultat en 3D.
• Le problème des shapefiles dits 3D, où la valeur z est intégrée dans la géométrie, provient du fait que QGIS ou PyQGIS ne les gèrent pas: il sont considérés comme des shapefiles 2D.

Pourtant, il y a moyen, dès à présent, de créer des représentations 3D à l'aide de la console Python et des modules Python qui permettent une telle représentation:

1. Pour le traitement des shapefiles dits 3D, il existe d'autres modules Python qui sont capables d'extraire la valeur z de la géométrie, comme ogr et Shapely, directement depuis la couche QGIS, ou Pyshp/shapefile à partir du fichier shapefile original (code.google.com/p/pyshp/wiki/CreateElevationValues). Ces modules ont déjà été détaillés sur le Portail et peuvent être utilisés dans la Console Python sans problème.
2. Pour la visualisation en 3 dimensions, nous allons utiliser des modules Python qui fonctionnent très bien dans la console Python:
   • le classique module Python Matplotlib;
   • un nouveau module, Visvis. Il est écrit en pur Python et se base sur Numpy et PyOpenGL;
   • nous laissons de côté le très puissant module Mayavi, déjà illustré sur ce Portail, trop complexe à utiliser pour ce que l'on veut faire.

Le problème des fichiers shapefiles dits 3D dans QGIS

Examinons tout d'abord le traitement de ces géométries dans la console Python avec une couche shapefile 3D:

•  traitement avec le module PyQGIS

Le point sélectionné est bien reconnu comme un Point 3D (WKBPoint25D) mais il n'y aucune fonction de la classe Point de PyQGIS qui permet d'extraire la coordonnée z.

• avec le module ogr:

Transformons ce point en objet ogr (à l'aide du format WKB):

Ici la classe Point de ogr permet de récupérer la coordonnée z

• avec le module Shapely

Même démarche avec le module Shapely:

et même conclusion pour la classe Point de Shapely

La solution est donc d'utiliser un de ces modules, en complément de PyQGIS, pour extraire la coordonnée z.

Traitements dans la console Python

Les shapefiles dits 3D

Les utilisateurs de QGIS qui veulent se lancer dans les traitements en Python s'appuient généralement sur les tutoriels de http://www.qgis.org/pyqgis-cookbook/index.html ou de www.qgisworkshop.org/html/index.html. Le problème est que ces deux tutoriels, sans vouloir vexer les auteurs, se révèlent fort complexes et peu  « Pythonesques » à utiliser dans la console comme l'ont souligné diverses personnes dont Sean Gillies (l'auteur de Shapely et Fiona), surtout dans les processus d'itération.

Prenons par exemple, "iterating over a Vector Layer" dans le premier tutoriel qui peut être effectué en quelques lignes de code:

• La première chose est de sélectionner tous les éléments (géométrie et attributs d'une couche). Une simple fonction suffit (sans feat = QgsFeature(), allAttrs = provider.attributeIndexes(), provider.select(allAttrs) et autres éléments...):

def select_tout(couche):
    couche.select([])
    couche.setSelectedFeatures([obj.id() for obj in couche])

• ensuite  la couche vectorielle choisie est activée et la fonction précedente est utilisée pour sélectionner tous les éléments:

macouche = qgis.utils.iface.activeLayer()
select_tout(macouche)

• et toutes les géométries et les attributs de la couche se retrouvent dans la variable macouche. il est alors possible d' utiliser les boucles for, les compréhensions de liste, etc.

Les coordonnées x, y, z sont ensuite extraites grâce à la fonction loads de Shapely:

from shapely.wkb import loads
# lecture de toutes les géométries de la couche choisie et extraction des coordonnées x,y,z
x=[]
y=[]
z=[]
for elem in macouche.selectedFeatures():
           # extraction de la géométrie au format wkb
           geom= elem.geometry() 
           wkb = geom.asWkb()
           # utilisation de Shapely pour extraire les coordonnées x,y,z 
           x.append(loads(wkb).x)
           y.append(loads(wkb).y)
           z.append(loads(wkb).z)

Les couches 2D avec un attribut

Ici, plus besoin d'utiliser un autre module puisque les coordonnées x et y peuvent être directement obtenues par PyQGIS.

Dans la boucle présentée, les attributs sont directement accessibles avec la fonction elem.attribute.map()

• comme dans le cas présent l'attribut z (elev) est en deuxième position, il sera donc obtenu par elem.attribute.map()[1]. Il y aurait moyen d'améliorer le processus en permettant de sélectionner l'attribut par son nom (elev), mais ce n'est pas le sujet ici;
• la partie .toFloat()[0] n'est que la transformation d'un objet Qvariant de PyQt4  en une variable standard de Python.

Le traitement devient donc :

# lecture des coordonnées x,y  et extraction des valeurs de l'attribut z
x=[]
y=[]
z=[]
for elem in macouche.selectedFeatures():
           geom= elem.geometry() 
           x.append(geom.asPoint()[0])
           y.append(geom.asPoint()[1])
           z.append(elem.attributeMap()[1].toFloat()[0])

Les résultats sont, dans les deux cas, 3 listes Python avec les coordonnées x,y et z.

Représentations 3D (points, surfaces à partir des points, grilles .asc)

1) les points

Il est dès lors très facile d'utiliser ces listes avec  Matplotlib dans la console pour représenter les points en 3D: 

from mpl_toolkits.mplot3d.axes3d import *
import matplotlib.pyplot as plt
from matplotlib import cm
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(x,y,z,c=z,cmap=plt.cm.jet)
plt.show()

avec Matplotlib:

avec Visvis :

import visvis
# ouverture d'une fenêtre et affichage des points
f = visvis.gca()
m = visvis.plot(x,y,z, lc='k', ls='', mc='g', mw=2, lw=2, ms='.')
f.daspect = 1,1,10 # z x 10

 

Il y a moyen, bien entendu, d'apporter aux figures toutes les fioritures désirées.

2) surfaces à partir des points

Il y a moyen d'aller beaucoup plus loin puisque des modules Python comme Matplotlib, Numpy ou SciPy disposent de toutes les fonctions d'interpolation que l'on désire.

interpolation avec la fonction griddata de matplotlib

C'est la plus simple à utiliser, elle utilise un algorithme de Delaunay ou celui des plus proches voisins (en fonction de l'installation de Matplotlib) pour interpoler une grille en fonction des coordonnées x, y, z

import numpy as np
from matplotlib.mlab import griddata
# création d'une grille 2D
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolation
Z = griddata(x, y, z, xi, yi)

avec Matplotlib: 

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,linewidth=1, antialiased=True)
plt.show()

 figure

avec Visvis:

f = visvis.gca()
m = visvis.grid(xi,yi,Z) # ou m = visvis.grid(yi,xi,Z) si une ancienne version
f.daspect = 1,1,10 # z x 10

grille avec la couche point superposée:

relief ombragé avec une table de couleur élévation (visvis.CM_JET):

m = visvis.surf(xi,yi,Z)
m.colormap = visvis.CM_JET

interpolation avec une des fonctions splines de SciPy

À titre d'exemple des possibilités, nous allons essayer l'un des algorithmes de spline proposé par le module SciPy: l'algorithme de « Spline en plaque mince (thin-plate spline) » (sigomatic.blogspot.be/2009/09/lissage-par-spline-en-plaque-mince.html)

import scipy as sp
import scipy.interpolate
# construction de la grille
spline = sp.interpolate.Rbf(x,y,z,function='thin-plate')
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolation
Z = spline(X,Y)

avec Visvis

f = visvis.gca()
m = visvis.surf(xi,yi,Z)
m.colormap = visvis.CM_JET 
f.axis.visible = False

3) traitement des fichiers .asc créés par l'extension Interpolation

Il y a 2 méthodes pour lire les fichiers .asc créés par l'extension Interpolation dans la console Python:

• avec le module Python osgeo/gdal
• depuis très longtemps, David Finlayson a proposé diverses versions de ses fonctions pour extraire les points x, y et z d'une grille ASCII ArcInfo (ou d'une grille Golden Surfer), en pur Python. La dernière version disponible se trouve à bitbucket.org/david_finlayson/pyraster/src.

Puisque c'est un grille régulière, il a des fonctions Numpy/SciPy pour les traiter directement.

résultat de l'extension Interpolation sur les points (par l'algorithme de Delaunay de l'extension)

résultat du traitement dans la Console Python (animation réalisée avec le module visvis):

L'extension QGIS gSurf, qui permet de calculer l'intersection d'un MNT et d'une surface géologique et de la figurer en 3D, utilise cette technique avec Matplotlib (www.malg.eu/gsurf_en.php).

Autres traitements

D'autres traitements peuvent être aussi effectués comme la représentation des lignes de contours en 3 dimensions (alors que l'extension Contour, elle aussi basée sur Matplotlib, ne le fait qu'en 2 dimensions)

résultat de l'extension contour sur les points

représentation 3D des lignes de contour (Matplotlib)

 

En pratique, toutes les lignes 3D ou les polygones 3D peuvent être représentés (Visvis):

Comme je suis géologue, il est possible de représenter des sondages (à partir d'un shapefile avec la cote des formations géologiques). Il est aussi possible de générer des surfaces pour chaque limite de formation.

Pour aller plus loin, l'exportation des surfaces

Tout ça c'est bien beau me direz-vous, mais ce n'est que de la visualisation, et on ne peut pas exploiter ces résultats. Et bien non, il y a très facilement moyen de sauver les grilles résultantes en grilles ASCII en utilisant les fonctions proposéés par David Finlayson ou en utilisant Numpy/SciPy.

Comparaison avec les résultats dans GRASS GIS

À aucun moment dans ces traitements en Python, n'interviennent les projections, hormis par l'entremise des coordonnées elles-mêmes (référentiel purement cartésien). Il est donc intéressant de comparer les résultats obtenus avec ceux obtenus par GRASS GIS avec ces mêmes données, dans un système de projection:

résultat avec Visvis (z x 1)                résultat avec GRASS GIS (z x 1)

résultat avec Visvis (z x 10)                                      GRASS GIS (z x 10)

Les différences principales ne résultent pas du fait que les données soient projetées ou non (pour autant que le système de projection soit en x,y) mais de l'algorithme d'interpolation, du coefficient d'exagération des hauteurs utilisé et du type de perspective. L'utilisation du module Basemap permettrait sans doute de traiter le problème des projections mais je n'en ai pas besoin pour le moment.

Conclusions

Il est clair que ces méthodes ne remplaceront jamais l'utilisation de logiciels comme GRASS GIS ou l'extension Globe toujours en développement:

figure reprise de www.faunalia.com/quantumgis

Mais elles me sont très utiles en première approche. Elles me permettent de visualiser ce que je veux directement à partir de points sélectionnés  et/ou de combiner les diverses sélections sans devoir passer par GRASS GIS.

Leur utilisation a déjà amené à des corrections dans les codes du module visvis (groups.google.com/forum/). Elles pourraient, à mon avis, facilement être transposées à ArcGIS avec ArcPy.

En faire une extension me semble difficile, car:

• tous les cas doivent être pris en considération avant de commencer, ce qui va demander beaucoup de temps, que je n'ai pas;
• de plus, je suis très content avec la console Python et si il le faut, il est possible de créer des scripts spécifiques avec l'extension Script Runner (www.portailsig.org/content/qgis-lancer-des-scripts-python-ou-des-commandes-shell-depuis-la-console-python-ou-avec-scrip).

Tous les traitements ont été effectués sur Mac OS X  avec QGIS, version "master branch"  de Larry Shaffer disponible à qgis.dakotacarto.com/ et QGIS, version 1.8 de  William Kyngesburye (www.kyngchaos.com/)

Site officiel : QGIS (Quantum GIS)
Autres Liens : Python et les shapefiles (suite) ou le module le plus simple à utiliser (un seul fichier, sans dépendance)
Autres Liens : Python: le module Shapely, géométries, prédicats spatiaux, analyse spatiale, matrices de Clementini (DE-9IM), traitements de shapefiles ou autres, et représentation avec matplotlib

licence Creative Commons Paternité-Pas d'Utilisation Commerciale-Partage des Conditions Initiales à l'Identique Commerciale 2.0 France