Python : Comment convertir les graphiques obtenus avec le module matplotlib en couches shapefiles, application aux champs vectoriels (streamplot()), contourages (countour(),contourf()) et explications

mar 08-03-2016 Martin Laloux

SIG OpenSource


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

Matplotlib est la librairie Python la plus ancienne et la plus complète pour les graphismes scientifiques en Python. Il y a bien sur des alternatives qui ont été proposées depuis sa naissance comme pandas, Seaborn et ggplot, simples aménagements de matplotlib, Bokeh ne dépendant pas de matplotlib et plutôt réservés aux notebooks Jupiter/IPython (du fait de sa conception), Pygal qui concerne le format SVG ou plot.ly, pour les visualisations en ligne, Vispy, basé sur OpenGL, ou Mayavi pour la 3D (inventaire non exhaustif...). Neamoins ce module reste le plus utilisé pour des graphiques complexes, du fait des multiples paramètres disponibles, complexes à souhait pour les débutants (d’ou certaines des alternatives proposées). Il permet même de créer toute sorte de cartes avec le Basemap Toolkit (nombreux exemples sur Internet). Ce thème ne sera pas abordé ici.

Matplotlib est très utilisé dans les plugins QGIS pour créer des diagrammes à partir de diverses données (trop nombreux pour être cités) et son utilisation a été de nombreuses fois illustrées dans des articles sur le Portail (comme dans QGIS et géologie structurale: création automatique de rosaces directionnelles et de canevas stéréographiques avec PyQGIS 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).

Mais savez-vous qu'il est possible de faire l'inverse, c'est-à-dire de récupérer les résultats des graphiques pour en créer des couches vectorielles ?

C'est ce que fait déjà le Contour plugin de QGIS qui utilise le module shapely pour convertir les résultats obtenus avec matplolib en une couche shapefile (d'une manière plus complexe, car il extrait aussi les valeurs des lignes de contour, ce que je ne présente pas ici).

Principes

Les résultats des commandes maplotlib sont des points, des lignes ou des polygones basés sur les données fournies (dans un repère cartésien, sans projection bien entendu). Si l’on connait les principes du module, il est donc possible d’extraire ces éléments géométriques et de les transformer en géométries shapely (Points, LineString, Polygon) et qui dit shapely dit aussi QGIS avec le plugin QuickWKT ou la création d’un shapefile avec le module Fiona.

La solution pratique se trouve dans la manière dont matplotlib produit ces figures: matplotlib.path.

Mais, pour commencer, je vais vous ouvrir l'appétit avec un exemple concret d'application puis je fournirai des explications.

Exemple concret: calcul des lignes tangentes (streamlines) d'un champ scalaire/vectoriel

C’est le genre de concept qui fait mathématiquement peur avec des calculs de dérivées et/ou d’intégrales (mécanique des fluides, par exemple). Théoriquement les notions ont été vues à l’école et sortent du cadre de cet article. Je laisse donc à d’autres, plus qualifiés que moi, le soin de les expliquer. Heureusement, matplotlib dispose de la commande pour faire cela simplement, streamplot(), mais cela ne vous dispense pas du fait de comprendre ce que vous faites avant de l’utiliser, sinon... (quelques liens sur les traitements Python dans le domaine spatial ou des plugins QGIS sont données en fin d'article)

1) Au départ, je dispose d’une grille régulière de points obtenue par traitements géostatistique (effectué avec R) avec des angles de direction sous forme de cosinus directeurs dont je veux obtenir les streamlines (pour simplifier ce sont les lignes de champ ou lignes tangentes à chaque point de la grille).

exemple de streamlines

ma grille de départ

2) Après importation du fichier shapefile et traitements, le résultat de la commande matplotlib quiver() (trace un champ 2D avec des flèches) est :

3) avec la commande streamplot() (calcul des streamlines) en plus:

4) avec la commande streamplot() seule:

5) résultat dans QGIS après exportation des lignes résultantes sous forme de shapefile:

Le script d'exportation après traitement avec matplotlib est simple

from shapely.geometry import LineString, mapping
import fiona
# lignes de résultat de la commande streamplot()
lines = test.lines.get_paths()
from fiona.crs import from_epsg
# définition du schéma du nouveau shapefile
monschema = {'geometry': 'LineString', 'properties': {'id' : 'int:2'}}
# définition du crs du nouveau shapefile
crs = from_epsg(31370)
# création du shapefile résultant
with fiona.open('stramplot_Leon3det05_1.shp','w',driver='ESRI Shapefile', crs=crs,schema= monschema) as sortie:
       for n,l in enumerate(lines):
              # le résultat de l.vertices est un array Numpy
          vertices = [(i[0], i[1]) for i in zip(l.vertices.T[0],l.vertices.T[1])]
          geom = mapping(LineString(vertices))
          prop = {'id': n}
          sortie.write({'geometry':geom, 'properties': prop})

Explication avec un exemple simple

Comment transformer les résultats de la commande contour() en géométries shapely ? Cette question a été posée dans Converting Matplotlib contour objects to Shapely objects sur GIS Stack Exchange et je vais détailler la réponse que j’ai fournie.

Commençons par les lignes (contour seul -> contour())

import matplotlib.pyplot as plt
x = [1,2,3,4]
y = [1,2,3,4]
m = [[15,14,13,12],[14,12,10,8],[13,10,7,4],[12,8,4,0]]
cs = plt.contour(x,y,m)
plt.show()

Si l’on s’en réfère aux paramètres de contour(), le nombre d’objets sur la figure (= collection de lignes) et les coordonnées de la première ligne sont donnés par:

print len(cs.collection)
7
# coordonnées de la première ligne
ligne = cs.collections[0].get_paths()[0]
v = ligne.vertices
x = v[:,0]
y = v[:,1]
print x, y
[ 4.   3.5] [ 3.5  4. ]

Il y a donc moyen de transformer très facilement ces géométries en LineString shapely :

from shapely.geometry import LineString
for i in range(len(cs.collections)):
    p = cs.collections[i].get_paths()[0]
    v = p.vertices
    x = v[:,0]
    y = v[:,1]
    line = LineString([(j[0], j[1]) for j in zip(x,y)])
    print line.wkt
LINESTRING (4 3.5, 3.5 4)
LINESTRING (4 3, 3 4)
LINESTRING (4 2.5, 3.333333333333333 3, 3 3.333333333333333, 2.5 4)
LINESTRING (4 2, 3 2.666666666666667, 2.666666666666667 3, 2 4)
LINESTRING (4 1.5, 3 2, 2 3, 1.5 4)
LINESTRING (4 1, 3 1.333333333333333, 2 2, 1.333333333333333 3, 1 4)
LINESTRING (2 1, 1 2)

Résultat dans QGIS en utilisant le plugin QuickWKT :

Avec les polygones (contourf())

import matplotlib.pyplot as plt
x = [1,2,3,4]
y = [1,2,3,4]
m = [[15,14,13,12],[14,12,10,8],[13,10,7,4],[12,8,4,0]]
cs = plt.contourf(x,y,m)
plt.show()

from shapely.geometry import polygon
for i in range(len(cs.collections)):
    p = cs.collections[i].get_paths()[0]
    v = p.vertices
    x = v[:,0]
    y = v[:,1]
    poly = Polygon([(i[0], i[1]) for i in zip(x,y)])
    print poly.wkt
POLYGON ((4 3.5, 4 4, 3.5 4, 4 3.5))
POLYGON ((4 3, 4 3.5, 3.5 4, 3 4, 4 3, 4 3))
POLYGON ((4 2.5, 4 3, 3 4, 2.5 4, 3 3.333333333333333, 3.333333333333333 3, 4 2.5))
POLYGON ((4 2, 4 2.5, 3.333333333333333 3, 3 3.333333333333333, 2.5 4, 2 4, 2.666666666666667 3, 3 2.666666666666667, 4 2, 4 2))
POLYGON ((3 2, 4 1.5, 4 2, 3 2.666666666666667, 2.666666666666667 3, 2 4, 1.5 4, 2 3, 3 2, 3 2))
POLYGON ((2 2, 3 1.333333333333333, 4 1, 4 1.5, 3 2, 2 3, 1.5 4, 1 4, 1.333333333333333 3, 2 2, 2 2))
POLYGON ((2 1, 3 1, 4 1, 3 1.333333333333333, 2 2, 1.333333333333333 3, 1 4, 1 3, 1 2, 2 1, 2 1))
POLYGON ((1 2, 1 1, 2 1, 1 2))

Résultat dans QGIS en utilisant le plugin QuickWKT :

Explications

Les géométries matplotlib sont toujours encodées de la manière suivante :

des vertices (points) sous forme d’array Numpy;
des codes qui indiquent comment matplotlib traite les relations géométriques:
- 0 = STOP: (fin du chemin), non requis;
- 1 = MOVETO: déplace le « crayon » vers un point donné;
- 2 = LINETO : dessine une ligne depuis un point donné vers un autre;
- 79 = CLOSEPOLY: dessine une ligne depuis un point jusqu’à une ligne existante (fermeture de polygone).

Et d’autres qui nous intéressent moins ici comme CURVE3 et CURVE4 pour dessiner des courbes de Bézier.

Il y a moyen de visualiser ces paramètres en utilisant la manière préconisée pour travailler dans matplotlib avec les contours linéaires précédemment présentés (avec iter.segments()):

for i in range(len(cs.collections)):
    p = cs.collections[i].get_paths()[0]
    for vert, code in p.iter_segments():
            print code, vert
1 [ 4.   3.5]
2 [ 3.5  4. ]
1 [ 4.  3.]
2 [ 3.  4.]
1 [ 4.   2.5]
2 [ 3.33333333  3.        ]
2 [ 3.          3.33333333]
2 [ 2.5  4. ]
1 [ 4.  2.]
2 [ 3.          2.66666667]
2 [ 2.66666667  3.        ]
2 [ 2.  4.]
1 [ 4.   1.5]
2 [ 3.  2.]
2 [ 2.  3.]
2 [ 1.5  4. ]
1 [ 4.  1.]
2 [ 3.          1.33333333]
2 [ 2.  2.]
2 [ 1.33333333  3.        ]
2 [ 1.  4.]
1 [ 2.  1.]
2 [ 1.  2.]

Une polyligne est donc définie par

1:[ 4. 2.] (départ) suivi de 2:[ 3. 2.66666667] (déplacement), de 2:[ 2.66666667 3] et de 2:[ 2. 4.] (fin de la ligne)

Ce qui qui correspond bien à la LINESTRING (4 2, 3 2.666666666666667, 2.666666666666667 3, 2 4) précédemment trouvée.

Des exemples sont donnés dans Path tutorial

Notons aussi que Matplotlib dispose aussi de certains prédicats spatiaux comme:

from matplotlib.path import Path
p = Path([[25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]]) 
p1=[27.254629577800088, -76.728515625]
p2=[27.254629577800088, -74.928515625]
p.contains_point(p1)
 
p.contains_point(p2)
1
# à comparer avec 
from shapely.geometry import Point, Polygons
poly = Polygon([[25.774252, -80.190262], [18.466465, -66.118292], [32.321384, -64.75737]]) 
p1=Point(27.254629577800088, -76.728515625)
p2=Point(27.254629577800088, -74.928515625)
poly.contains(p1)
False
poly.contains(p2)
True