Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG

	Niveau	Intermédiaire
	Logiciels utilisés	Python
	Plateforme	Windows | Mac | Linux | FreeBSD

En pratique, le titre est un peu fallacieux, car on a toujours besoin d'un logiciel SIG pour créer les couches à traiter ou pour visualiser les résultats (quoique...), mais je vais présenter ici une introduction à l'utilisation des objets géomatiques (couches vectorielles et matricielles géoréférencées)  dans des scripts Python à finalité géologique, puisque c'est mon métier.

J'en ai écrit un grand nombre basés sur les algorithmes présentés, entre autres, par R. Groshong dans 3-D Structural Geology (Springer, 2008,) ou R. Allmendinger et N. Cardozo  dans Structural Geology Algorithms Vectors and Tensors (Cambridge University Press, 2011), mais les auteurs ne tiennent pas grand compte des possibilités offertes par l'utilisation des objets géomatiques.  Je vais donc essayer de vous expliquer ici comment les traiter et/ou les combiner avec plusieurs modules Python dont le principal est osgeo (gdal et ogr), le tout sans utiliser un logiciel SIG.

En pratique, je vais  refaire en Python ce que j'avais présenté dans Grass Gis et R: superposition des couleurs d'un raster quelconque sur un profil topographique élaboré à partir d'un MNT mais en allant plus loin (en y ajoutant des valeurs structurales de stratification, par exemple). Cette première partie peut intéresser tous les Sigistes.

Mais, en soi, créer un profil topographique ne sert pas à grand-chose, surtout qu'il faut le compléter pour en faire une coupe géologique. Je présenterai en final un des résultats de mes traitements plus géologiques avec la création d'un profil aval-pendage (« down-plunge profile » ou coupe suivant le plongement de l'axe du pli, voir www.uwgb.edu/dutchs/structge/SL165DownPlunge.HTM),  à partir du script originel en MATLAB présenté dans Structural Geology Algorithms Vectors and Tensors, pp. 51-55, converti en Python par mes soins.

Nous verrons donc dans la suite:

• Comment lire les données de type vectoriel ou matriciel (depuis un fichier shapefile ou un fichier raster, depuis une geodatabase fichier de GRASS GIS et depuis un service WMS);
• Comment extraire les valeurs d'un pixel d'un fichier matriciel;
• Comment croiser un fichier de type vectoriel et un fichier de type matriciel pour en extraire les données (altitude, couleurs RGB);
• Comment créer un profil topographique en combinant une couche vectorielle, un MNT et un fichier raster ;
• Comment ajouter les points d'intersection entre 2 fichiers vectoriels (profil et limites des objets géologiques);
• Comment ajouter les valeurs structurales de stratification (pendage et direction) et calculer le pendage apparent sur la ligne de coupe;
• Exemple d'application géologique.

Les résultats seront visualisés ici avec le module matplotlib (et le module visvis pour une figure), mais tous les résultats peuvent être exportés en taille réelle dans divers formats dont le vectoriel (PDF ou SVG) pour être utilisés dans d'autres logiciels.

Comment lire les données

Les fichiers de type vectoriel

Je ne reviendrai pas ici sur le cas des fichiers vectoriels de type shapefile, déjà traité sur le Portail.  Les formats que le module ogr dans osgeo peut lire sont obtenus par la commande:

ogrinfo --formats

Cela peut varier suivant les installations, mais, il s'agit de la majorité des formats classiques avec, depuis peu, ceux qui sont présents dans les geodatabases fichiers d'ESRI (FileGDB).

Ouvrir un fichier vecteur dans la geodatabase fichier de GRASS GIS n'est qu'une question de chemin:

from osgeo import ogr
# ouverture d'une couche vectorielle de GRASS GIS
ds = ogr.Open('/Users/martinlaloux/grassdata/geol/MNT/vector/ligne_de_coupe/head')
couche = ds.GetLayer(0)
detail = couche.GetFeature(0)
geom = detail.GetGeometryRef()
# noeuds de la ligne
for i in range(geom.GetPointCount()):
     xy = geom.GetPoint(i)
     print xy
(206643.21517600652, 125181.18058575876)
(201007.33432923188, 121517.85552059766)

Les fichiers de type matriciel (raster, MNT,...)

De la même manière, les formats que le module gdal dans osgeo permet de lire sont indiqués par la commande:

gdalinfo --formats

Je vais présenter ici quelques exemples:

Fichier raster géoréférencé:

Si je dispose d'un classique raster géoréférencé (vu ici dans QGIS):

La lecture du fichier se fait de la manière suivante:

from osgeo import gdal
fichier = '/Users/Shared/cartes_geol/hen_raer.jpg' 
couche = gdal.Open(fichier) 
gt =couche.GetGeoTransform()
bandes = couche.RasterCount 
# ce qui donne les paramètres de géoréférencement
print gt 
(258012.37107330866, 2.11668210080698, 0.0, 163176.63853988209, 0.0, -2.1168501270110074) 
# et le nombre de bandes de couleurs print bandes
3

La transformation de ce fichier raster en vecteurs (array) numpy se fait avec la commande gdal ReadAsArray (ici pour la première bande de couleur):

bande1 = couche.GetRasterBand(1)
donnee1 = bande1.ReadAsArray()
etc...

Une fois que les bandes sont transformées en vecteurs, il est facile d'extraire le morceau de carte compris dans la zone limitée par un rectangle rouge et de visualiser le résultat avec la commande imshow du module matplotlib:

Couche MNT dans GRASS GIS:

C'est la même démarche avec une couche matricielle présente dans une geodatabase fichier de GRASS GIS, il suffit d'ajuster le chemin d'accès:

fichier = '/Users/martinlaloux/grassdata/geol/MNT/cellhd/dem10m'
couche = gdal.Open(fichier)
gt =couche.GetGeoTransform()
bandes = couche.RasterCount
 # ce qui donne les paramètres de géoréférencement
print gt
(41995.0, 10.0, 0.0, 172185.0, 0.0, -10.0)
 # et le nombre de bandes de couleurs
print bandes
1

(le MNT est illustré ici sous forme d'une surface ombragée dans GRASS GIS)

Service WMS

L'organisme où je travaille offre aussi un service WMS pour les cartes géologiques (geoservices.wallonie.be/arcgis/rest/services/SOL_SOUS_SOL/CARTE_GEOLOGIQUE_SIMPLE/MapServer, image vue ici avec OpenLayers). GDAL sait accéder aux services WMS via un fichier de configuration XML dont les détails figurent dans www.gdal.org/frmt_wms.html. Ce fichier XML décrit habituellement l'URL, la projection, le type d'image et la (les) couche(s) requise(s), mais il est aussi possible de sélectionner l'étendue désirée (et donc pas besoin de numpy ici).

En ouvrant le fichier XML, le résultat est:

fichier = '/Users/martinlaloux/wms/carte_geol.xml'
couche = gdal.Open(fichier)

service WMS visualisé avec OpenLayers

Comment extraire la valeur d'un pixel

Une fois ces fichiers ouverts, comment obtenir les valeurs des rasters situés sous un objet vectoriel comme le point bleu ici (sur une carte géologique et un MNT):

Commençons par le MNT (dans GRASS ici, mais ça peut être un fichier de type ASCII Grid d'ESRI, .asc ou  GeoTIFF):

fichier = '/Users/martinlaloux/grassdata/geol/MNT/cellhd/dem10m'
couche = gdal.Open(fichier)
gt =couche.GetGeoTransform()
bandes = couche.RasterCount
print bandes
1
print gt
(263104.72544800001, 10.002079999999999, 0.0, 155223.647811, 0.0, -10.002079999999999)

Les valeurs de géoréférencement obtenues ont été expliquées dans Python: géoréférencement et worldfiles (tfw, jgw,pngw,...) et je ne reviendrai pas dessus . Dès lors, passer des coordonnées xy du point à celles du raster est une simple règle de trois:

x,y  = (263220.5,155110.6)
rasterx = int((x - gt[0]) / gt[1])
rastery = int((y - gt[3]) / gt[5])
print rasterx,rastery
11,11

(Attention, je ne tiens pas compte ici de la particularité de GDAL pour le géoréférencement, voir Python: géoréférencement et worldfiles (tfw, jgw,pngw,...))

)

Puisqu'il n'y a qu'une bande, la valeur du pixel sous le point est:

print couche.GetRasterBand(1).ReadAsArray(px,py, 1, 1)
array([[222]], dtype=int32)

Dans le cas présent (MNT), cette valeur correspond à une altitude (222 m).

Dans le cas de la couche géologique, nous avons vu qu'il y a 3 bandes. Après avoir recalculé les valeurs px et py en fonction du raster,  les résultats sont:

bande1 = couche.GetRasterBand(1)
bande2 = couche.GetRasterBand(2)
bande3 = couche.GetRasterBand(3)
print bande1.ReadAsArray(px2,py2, 1, 1)
[[253]]
print bande2.ReadAsArray(px2,py2, 1, 1)
[[215]]
print bande3.ReadAsArray(px2,py2, 1, 1)
[[118]] 

qui constituent, en pratique, les valeurs RGB de la couleur du pixel.

C'est la même démarche avec le service WMS:

fichier = 'carte_geol.xml'
couche = gdal.Open(fichier)
gt =couche.GetGeoTransform()
bande1 = couche.GetRasterBand(1)
print bande1.ReadAsArray(px,py, 1, 1)
[[247]]
etc.

Quelque soit le type de raster utilisé, il est donc possible d'obtenir la valeur du pixel situé sous un point de coordonnées x, y. On peut alors créer une fonction universelle qui permet d'extraire la valeur d'un pixel sous forme de liste, quel que soit le nombre de bandes. Contrairement à Python ou la valeur des indices commence à 0 dans les boucles (j ici, si trois bandes, les valeurs sont 0,1,2), celle de GetRasterBand commence à 1 (valeurs 1, 2, 3). Il n'y a donc pas de GetRasterBand(0) d'où le j+1 dans la boucle:

def Val_raster(x,y,src_ds,bandes,gt)
      col=[]
      px = int((x - gt[0]) / gt[1])
      py =int((y - gt[3]) / gt[5])
      for j in range(bandes):
         bande = src_ds.GetRasterBand(j+1)
         data = band.ReadAsArray(px,py, 1, 1)
         col.append(data[0][0])
      return col

Application:

# MNT
px1 = int((x - gt1[0]) / gt1[1])
py1 = int((y - gt1[3]) / gt1[5])
# couche géologique
px2 = int((x - gt2[0]) / gt2[1])
py2 = int((y - gt2[3]) / gt2[5])
# résultats
print Val_raster(x,y,couche1, bandes1,gt1)
[222] # altitude 
print Val_raster(x,y,couche2, bandes2,gt2)
[253, 215, 118] # valeurs RGB 

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

Ligne de profil et fichiers rasters

Il est donc possible de croiser n'importe quel fichier vectoriel (en le transformant en point) avec des fichiers rasters pour en extraire les valeurs. 

Le profil topographique (en bleu) est ici figuré sur la carte géologique et sur le MNT (représenté ici par un relief ombragé provenant de GRASS GIS.

La première chose à faire est d'ouvrir le fichier vecteur (ligne bleue, couche GRASS ou fichier shapefile) avec le module ogr, puis de générer des points équidistants sur la ligne. Ogr ne permet pas encore de le faire, mais le module Shapely, oui, avec la fonction interpolate. Comme cette ligne de profil n'a pas besoin d'être droite (il est possible de traiter une ligne avec plusieurs segments), nous utiliserons la fonction Union de ogr pour n'obtenir qu'un seul objet géométrique.

# traitement de la ligne de coupe
from shapely.wkb import loads
# création d'une couche ogr LineString pour regrouper tous les segments
# éventuels d'une couche ligne
profilogr = ogr.Geometry(ogr.wkbLineString)
# ouverture du fichier shapefile et intégration des segments éventuels
source = ogr.Open('profilHC2.shp')
cshp = source.GetLayer()
for element in cshp:
      geom =element.GetGeometryRef()
      profilogr = profilogr.Union(geom)
# transformation en géométrie Shapely
profilshp = loads(profilogr.ExportToWkb())
# création des points équidistants sur la ligne avec un pas de 20m
longueur=profilshp.length
x = []
y = []
z = []
# couleur RGB
couleur = []
# distance pour le profil topographique
dista = []
for distance in range(0,longueur,20):
         # création du point sur la ligne 
         point = profilshp.interpolate(distance)
         xp,yp=point.x, point.y
         x.append(xp)
         y.append(yp)
         # extraction de la valeur altitude à partir du MNT
         z.append(Val_raster(paramètres de la couche MNT)[0])
         # extraction des couleurs RGB à partir de la couche raster
         couleur.append(Val_raster(paramètres de la couche raster))
         dista.append(distance)

Ce qui permet de dessiner le profil en 3D (x,y,z) (modules matplotlib)  ou de le placer sur le MNT à titre de contrôle (module visvis):

et de dessiner le profil 2D classique (distance, z x 10):

et sans exagération des hauteurs:

Comment ajouter les points d'intersection entre 2 fichiers vectoriels (profil et limites des objets géologiques)

Il est aussi possible d'utiliser la couche vectorielle des limites des formations géologiques (en rouge ici):

Le module ogr ne dispose pas encore de la possibilité de calculer directement l'intersection entre 2 couches vectorielles (ce sera disponible dans la version 1.10) mais en transformant les 2 couches en 2 objets géométriques uniques (avec la fonction Union de ogr)  il est possible d'utiliser la fonction Intersection de ogr  (il est aussi possible de le faire avec le module Shapely):

# ouverture de la couche 
limites = ogr.Geometry(ogr.wkbLineString)
source2 = ogr.Open('lim_form_HC2.shp')
lshp = source2.GetLayer()
for element in lshp:
     geom =element.GetGeometryRef()
     limites = limites.Union(geom)
# points d'intersection entre les 2 couches
ptintersec = profilogr.Intersection(limites)

avec comme résultat la série de points d'intersection:

  

et le profil résultant est (z x 10) :

et sans exagération de hauteur:

Comment ajouter des valeurs structurales de stratification (pendage et direction) en calculant le pendage apparent sur la ligne de coupe.

Si, comme sur cette carte géologique, on dispose de mesures de pendage et de direction, il y a moyen d'aller beaucoup plus loin en se basant sur les scripts Python apparentdip.py de Geologize (geo.distortions.net/search/label/apparent%20dip) qui permet de calculer le pendage apparent de la mesure sur la ligne de coupe et les intersections ou lineorientation_07.py de Scientific and Gis programming with Python and Qt (www.malg.eu/lineangle.php) qui permet de calculer les angles entre  des lignes qui se croisent et leurs intersections. Je les ai combinés en utilisant Shapely en plus d'ogr pour arriver à un script qui répond à mes besoins.

Procédure

Le script commence par créer une zone tampon (buffer) à partir de la ligne de coupe pour sélectionner automatiquement les points à retenir (il est possible de choisir la distance de la zone tampon). Je peux aussi choisir les mesures qui me conviennent en me basant sur la connaissance de la structure géologique, sans zone tampon:

Les lignes de direction des couches géologiques sont ensuite interpolées et tracées et les points d'intersection avec la ligne de profil sont calculés:

ainsi que les pendages apparents sur la coupe (en fonction des angles entre les lignes de direction et la ligne de profil):

ce qui permet de placer ces points sur la coupe avec l'inclinaison des couches:

• sans exagération des hauteurs

 

• et avec exagération des hauteurs (x10) et inclinaison du pendage calculée en fonction de l'exagération (= arctan(pendage x 10)):

Exemples d'application géologique:

R. Allmendinger et N. Cardozo (Structural Geology Algorithms Vectors and Tensors) ont proposé un script MATLAB pour créer des « down-plunge profile » (DownPlunge.m). Transformé en script Python, il est ici appliqué à leur exercice 3.5  de la p.65 (anticlinal de Big Elk, Idaho, USA). Une animation du processus est disponible à www.geo.cornell.edu/geology/classes/RWA/EAS_5220/problem_sets/big_elk_anticline_graphics/big_elk_anticline_animation.html

• les données de départ sont des points situés le long des limites des couches (avec des coordonnées x,y, et z). Ils peuvent être automatiquement obtenus avec les procédures décrites:

• ces coordonnées x,y,z  sont transformées dans le système de coordonnées du pli avec l'axe z' parallèle à l'axe du pli (26 vers 125 ici) et l'axe y' perpendiculaire à z'. Le plan de la coupe est  x'y' (changement de système de coordonnées, voir l'animation):

• et le processus est réalisé directement dans la console Python de QGIS ou avec un script Python seul (figure précédente):

Conclusions

Pour faire cela, il est nécessaire de connaître Python, me direz-vous. Oui, certes, mais il me semble qu'il est de plus en plus nécessaire de connaître un langage de programmation au vu de la puissance des ordinateurs. C'est d'ailleurs ce que montrent R. Allmendinger et N. Cardozo (Structural Geology Algorithms Vectors and Tensors) avec leurs nombreux scripts en MATLAB qui permettent de faire beaucoup plus rapidement ce que j'ai appris à faire avec une combinaison de méthodes graphiques et de calculs de trigonométrie plane. Ils sont disponibles dans la partie « Resources » du site du livre, mais il faut l'acheter pour les remettre dans leurs contextes. Ces scripts sont très facilement transformables en scripts Python avec l'avantage de ne pas devoir payer la licence de MATLAB... (il faut aussi signaler que la plupart des scripts fonctionnent sur les clones libres de MATLAB que sont Octave (www.gnu.org/software/octave/) ou Freemat (freemat.sourceforge.net/), hormis quelques fonctions graphiques ou spécialisées)

Il est aussi possible de sélectionner les valeurs à partir d'un fichier vectoriel et d'utiliser d'autres scripts comme ceux de github.com/joferkington/mplstereonet ou de yonggeng.wordpress.com/2010/11/16/rose-diagram-code/ (entre autres) qui utilisent matplotlib.

et de les utiliser directement (ici dans la console Python de QGIS) à partir des valeurs sélectionnées (en rouge):

Tous les traitements ont été faits sur Mac OS X et les figures créées avec matplotlib et visvis hormis les captures d'écran (QGIS, GRASS GIS et OpenLayers).

Site officiel : module osgeo (gdal, ogr)
Site officiel : module Shapely
Site officiel : module numpy
Site officiel : module matplotlib
Site officiel : modile visvis
Autres Liens : 3-D Structural Geology
Autres Liens : Structural Geology Algorithms, Vectors and Tensors
Autres Liens : Grass Gis et R: superposition des couleurs d'un raster quelconque sur un profil topographique élaboré à partir d'un MNT
Autres Liens : Python: géoréférencement et worldfiles (tfw, jgw,pngw,...)

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