Python: traitement des couches vectorielles dans une perspective géologique, lecture et enregistrement des couches sous forme de dictionnaires avec le module Fiona

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

Lorsque j'essaie d'expliquer à mes collègues que Python est simple pour traiter des données géospatiales, ils me regardent d'un drôle d'œil, même s'ils connaissent le langage pour leurs traitements géologiques. Il est vrai que travailler avec des modules comme osgeo.ogr, PyQGIS ou ArcPy n'est pas une sinécure avec des while provider.nextFeature(feat), des feat.GetGeometryRef() ou des arcpy.UpdateCursor(points) et autres joyeusetés qui obligent à avoir en permanence  un aide-mémoire à côté de soi. Cela vient du fait que ces modules veulent tout faire en même temps, lire et écrire les données et les traiter, sans parler de la création des interfaces comme dans les plugins de Quantum GIS, le tout dans un joyeux mélange.... Seul Pyshp échappe en partie à cette règle, mais demande toujours de jongler entre les géométries et les tables d'attributs. De plus, le module  rencontre de sérieux problèmes si, par mégarde, un des champs renseignés est vide (aucune donnée).

Longtemps j'ai rêvé que les traitements de données spatiales en Python soient aussi simples que pour les autres types de données. Sauvegarder sous forme de shapefile les résultats d'un traitement de données géologiques (avec des angles...) m'a souvent demandé plus de travail  que le traitement lui-même....

Pourtant, il existe une solution simple qui permet de manipuler ces données avec des structures classiques comme des dictionnaires ou des listes. Si comme moi, vous avez des scripts tout prêts et que vous voulez les utiliser tels quels, sans devoir tout modifier pour créer des résultats sous forme de fichier shapefile, par exemple, alors la solution est le module Fiona de Sean Gillies.  Son seul but est simplement de lire et d'écrire les données à traiter dans un format géospatial, rien de plus, rien de moins.

Puisqu'un petit exemple vaut mieux qu'un beau discours :

import fiona
couche = fiona.open('testligne.shp')
rec = couche.next()
print rec

dont le résultat est:

{'geometry': {'type': 'LineString', 'coordinates': [(269884.20917418826, 151805.1917153612), (270409.89083992655, 153146.21637285672), (272298.05355768028, 154047.38494269375), (272941.74539327814, 155484.96337552898), (272169.31519056071, 156117.92701386689)]}, 'id': '1',  'properties': {'faille': u'de Salinas', 'type': u'normale'}}

et on reconnait un simple dictionnaire Python dans lequel tous les éléments d'un objet spatial (géométrie et attributs) sont représentés. Plus de notions de couches et d'attributs, pas de curseurs, plus d'opérations géométriques ni de transformations entre systèmes de coordonnées à gérer, avec Fiona, lire et écrire des données spatiales se résument à lire ou à placer des éléments dans un dictionnaire Python. Toutes les autres opérations sont laissées à d'autres modules Python comme Shapely,  Pyproj, Numpy, Scipy ou Math. En pratique, Fiona est particulièrement adapté à Shapely, auteur oblige...

Tout cela est rendu possible par l'utilisation du format GeoJSON car Fiona n'est en fait qu'un habillage du module ogr et ici, les utilisateurs de Windows n'ont plus d'excuse, car le module prêt à l'emploi est disponible chez Christoph Gohlke (www.lfd.uci.edu/~gohlke/pythonlibs/#fiona).

J'aborderai l'utilisation de Fiona avec des exemples d'applications géologiques:

• Principes de Fiona, lire et écrire des fichiers;
• Premier exemple: les azimuts ou directions d'objets géologiques linéaires;
• Deuxième exemple: les intersections entre une ligne de coupe et les directions des couches géologiques;
• Troisième exemple: les transformations affines, translations, rotations, mises à l'échelle et cisaillements.

Principes de Fiona

Fiona travaille avec des itérateurs et des générateurs. Comprendre comment ceux-ci fonctionnent est très important si l'on veut utiliser Python de manière productive. Il y a une multitude d'explications disponibles sur Internet et je ne vais donc pas en proposer une énième fois. Néanmoins, en simplifiant:

• un itérateur permet de parcourir une liste d'objet en implémentant la fonction next() (ce qui a été fait dans couche.next()) qui va retourner la prochaine valeur (et l'exception StopIteration quand il n'y a plus d'élément). Vous avez déjà tous utilisé les itérateurs lorsque vous parcourez un objet « itérable » comme dans le cas des boucles for élément in sequence (voir sametmax.com/dis-papa-dis-papa-dis-moi-dis-moi-comment-cest-fait-dans-une-boucle-for/).
• un générateur est une fonction qui va retourner les valeurs sous forme d'itérateur à l'aide de la fonction yield() (qui est elle même un itérateur, disposant de la méthode next()). Les compréhensions de listes sont en pratique des générateurs (voir sametmax.com/comment-utiliser-yield-et-les-generateurs-en-python/)

L'intérêt des générateurs est multiple, en effet les valeurs ne sont pas conservées en mémoire. Une fois qu'une valeur a été générée et envoyée, elle est perdue, ce qui facilite les traitements de grosses masses de données.

Mais ne vous inquiétez pas, tout cela est transparent dans Fiona !

Lecture des données

pour bien comprendre comment fonctionne Fiona: parcours du résultat de manière brute avec la fonction next():

Je vais ici parcourir l'itérateur renvoyé par Fiona avec la fonction next() pour bien saisir le module: la fonction  va renvoyer les divers élements de la couche shapefile sans les conserver en mémoire:

 # au niveau de la couche (type de géométrie + définition des champs)
print couche.schema
{'geometry': 'LineString', 'properties': {u'faille': 'str:20', u'type': 'str:20', u'id': 'int'}}
print couche.crs # (crs au format Proj4)
{'lon_0': 4.3674866666666663, 'ellps': 'intl', 'y_0': 5400088.4380000001, 'no_defs': True, 'proj': 'lcc', 'x_0': 150000.01300000001, 'units': 'm', 'lat_2': 49.8333339, 'lat_1': 51.166667233333328, 'lat_0': 90}
 print couche.bounds # emprise de la couche
(269884.20917418826, 151805.1917153612, 275291.22059321037, 156117.92701386689)
 # au niveau de l'entrée (valeurs)
print rec['geometry']
{'type': 'LineString', 'coordinates': [(269884.20917418826, 151805.1917153612), (270409.89083992655, 153146.21637285672), (272298.05355768028, 154047.38494269375), (272941.74539327814, 155484.96337552898), (272169.31519056071, 156117.92701386689)]}
print rec['geometry']['coordinates']
[(269884.20917418826, 151805.1917153612), (270409.89083992655, 153146.21637285672), (272298.05355768028, 154047.38494269375), (272941.74539327814, 155484.96337552898), (272169.31519056071, 156117.92701386689)]
 # transformation directe en géométrie Shapely
from shapely.geometry import shape
geom = shape(rec['geometry'])
print geom
LINESTRING (269884.2091741882613860 151805.1917153612012044, 270409.8908399265492335 153146.2163728567247745, 272298.0535576802794822 154047.3849426937522367, 272941.7453932781354524 155484.9633755289833061, 272169.3151905607082881 156117.9270138668944128)
python># valeur d'un attribut
 print rec['properties']['faille']
u'de Salinas'
# entrée suivante
rec= couche.next()
print rec
{'geometry': {'type': 'LineString', 'coordinates': [(270892.65971662494, 151805.1917153612), (272416.06372753991, 152856.55504683769), (274003.8369220146, 153478.79048791563), (274819.17991377192, 154283.40528241298), (275291.22059321037, 154648.16398925177)]}, 'id': '1', 'properties': {'faille': u'Arnao', 'type': u'chevauchement', 'id': None}}
 etc.

Si je continue indéfiniment avec couche.next(), cela va générer une exception lorsqu'il n'y aura plus d'éléments à rendre.

avec un générateur

Pour cela, je pourrai créer ma propre fonction:

def iterateur(c):
    ligne =  c.next()
    while ligne:
       yield ligne
       ligne = c.next()
 # traitement
for ligne in iterateur(couche):
     print ligne

Elle ne permet toujours que de lire les éléments, sans soulever l'exception StopIteration mais si je veux les conserver je dois utiliser une compréhension de liste. Mais comme Fiona fournit son propre générateur nommé collection, je continue avec celui-ci:

iterateur  = fiona.collection('testligne.shp')
# sans conservation en mémoire des éléments
for i in iterateur:
    print i
{'geometry': {'type': 'LineString', 'coordinates': [(269884.20917418826, 151805.1917153612), (270409.89083992655, 153146.21637285672), (272298.05355768028, 154047.38494269375), (272941.74539327814, 155484.96337552898), (272169.31519056071, 156117.92701386689)]}, 'id': '0', 'properties': {'faille': u'de Salinas', 'type': u'normale', 'id': None}}
{'geometry': {'type': 'LineString', 'coordinates': [(270892.65971662494, 151805.1917153612), (272416.06372753991, 152856.55504683769), (274003.8369220146, 153478.79048791563), (274819.17991377192, 154283.40528241298), (275291.22059321037, 154648.16398925177)]}, 'id': '1', 'properties': {'faille': u'd' Arnao', 'type': u'chevauchement', 'id': None}}
# avec conservation en mémoire des éléments dans une liste (par compréhension de liste)
iterateur  = fiona.collection('testligne.shp')
data = [elem for elem in iterateur]
print data
[{'geometry': {'type': 'LineString', 'coordinates': [(269884.20917418826, 151805.1917153612), (270409.89083992655, 153146.21637285672), (272298.05355768028, 154047.38494269375), (272941.74539327814, 155484.96337552898), (272169.31519056071, 156117.92701386689)]}, 'id': '0', 'properties': {'faille': u'de Salinas', 'type': u'normale', 'id': None}}, 
{'geometry': {'type': 'LineString', 'coordinates': [(270892.65971662494, 151805.1917153612), (272416.06372753991, 152856.55504683769), (274003.8369220146, 153478.79048791563), (274819.17991377192, 154283.40528241298), (275291.22059321037, 154648.16398925177)]}, 'id': '1', 'properties': {'faille': u'd'Arnao', 'type': u'chevauchement', 'id': None}}]

Écriture

L'écriture est tout aussi facile, il suffit de créer un dictionnaire de le remplir:

from shapely.geometry import LineString
import fiona
# deux géométries Shapely
lignes = [LineString([(272830.63,155125.73),(273770.32,155467.75)]),LineString([(273536.47,155914.07),(272033.12,152265.71)])]
# définition du schéma général de la couche avec un attribut 
# ou création du dictionnaire
schema = {'geometry': 'LineString','properties': {'titi': 'int'}}
# écriture dans un nouveau shapefile
with fiona.collection('monshp.shp', 'w', 'ESRI Shapefile', schema) as couche:
    for ligne in lignes:
        # remplissage du dictionnaire schema
        elem = {}
        # utilisation de la fonction mapping de Shapely
        elem['geometry'] = mapping(ligne) 
        # valeur de l'attribut (un seul ici) 
        elem['properties'] = {'titi': 145}
        # écriture de l'élément dans le fichier
        couche.write(elem)

Le module ogr sur lequel se base Fiona a de gros problèmes lorsqu'il s'agit de modifier une couche, il vaut mieux alors en faire une copie. Avec Fiona, c'est un jeu d'enfant:

from shapely.geometry import mapping, shape
import fiona
# Lecture du fichier shapefile originel
with fiona.collection('monshp.shp', 'r') as entree:
     # le fichier de sortie a le même schema -> copie du dictionnaire
     schema = entree.schema.copy()
     # écriture du nouveau fichier shapefile
     with fiona.collection('monshp_copie.shp', 'w', 'ESRI Shapefile', schema) as sortie:
         for elem in entree:
             sortie.write({'properties': elem['properties'],'geometry': mapping(shape(elem['geometry']))})

Cela nous sera très utile dans la suite.

Premier exemple: les azimuts ou directions d'objets géologiques linéaires

Trouver les azimuts ou directions d'objets géologiques linéaires est une tâche courante en géologie pour les traitements postérieurs. Or si calculer la direction d'une ligne simple (A) est facile,  ce n'est pas le cas du calcul des azimuts des segments de lignes plus complexes (B):

La démarche est alors la suivante:

1. lire le fichier shapefile originel
2. récupérer les segments de ligne et calculer leurs azimuts
3. les insérer dans un nouveau champ d'un nouveau fichier shapefile contenant aussi les attributs de l'original

Puisque je dispose déjà d'une fonction permettant de calculer l'azimut d'une ligne à partir de deux points (géométrie Shapely) en fonction de l'axe N-S, je vais la réutiliser (je travaille dans un système de coordonnées cartésiennes et je veux toutes mes valeurs entre 0 et 180°):

def azimut(point1, point2):
       '''Retourne l'azimuth de la ligne entre 2 points shapely'''
       angle = math.atan2(point2.x - point1.x, point2.y - point1.y)
       return math.degrees(angle)if angle >0 else math.degrees(angle) + 180

Je dispose aussi d'un générateur qui me permet de traiter les éléments dans une liste par paires (tiré de stackoverflow.com/questions/1257413/iterate-over-pairs-in-a-list-circular-fashion-in-python)

def paires(liste):
     '''parcourt une liste par paires'''
     for i in range(1, len(liste)):
         yield liste[i-1], liste[i]

je vais l'utiliser pour extraire les segments de ligne: 

with fiona.collection('testligne.shp', 'r') as entree:
    # copie du schema de la couche et création d'un nouveau champ 'azimut'
    schema = entree.schema.copy()
    schema['properties']['azimut'] = 'int'
    # création d' une nouvelle couche avec le schéma résultant
    with fiona.collection('testligne_azim.shp', 'w', 'ESRI Shapefile', schema) as sortie:
        for ligne in entree:
        # utilisation de la fonction paire() pour extraire les segments de lignes
            for seg_start, seg_end in paires(ligne['geometry']['coordinates']):
                # création d'une ligne en fonction des points des segments
                line_start = Point(seg_start)
                line_end = Point(seg_end)
                segment = LineString([line_start.coords[0],line_end.coords[0]])
                # copie des attributs d'entrée et ajout de la valeur résultante
                # de la fonction azimut()
                elem = {}                  
                elem['properties'] = ligne['properties'] 
                elem['properties']['azimut'] = azimut(line_start, line_end)
                elem['geometry'] = mapping(segment)
                sortie.write(elem)

et le résultat est:

 

Si l'on veut rajouter la longueur de chaque segment, il suffit d'insérer le résultat de la fonction suivante dans une nouvelle clé du dictionnaire:

def distance_euclid2pt(p1, p2):
 '''distance euclidienne entre 2 points shapely''' 
    vect_x = p2.x - p1.x
    vect_y = p2.y - p1.y 
    return math.sqr(vect_x**2 + vect_y**2))

ou, plus simplement avec la fonction length de Shapely.

 return LineString([p1.coords[0],p2.coords[0]]).length

 

Deuxième exemple: les intersections entre une ligne de coupe et les directions des couches géologiques

Je vais vous montrer ici comment résoudre une partie de la procédure de « Comment ajouter des valeurs structurales de stratification (pendage et direction) en calculant le pendage apparent sur la ligne de coupe » dans Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG sans passer par des couches intermédiaires (lignes en pointillés noirs) et en utilisant une fonction que j'avais déjà, sans la modifier:

def ptor2ptfin(p, angle, dist):
    ''' calcul de la position du point situé à une distance dist du point p (Shapely) dans la direction angle'''
    dist_x, dist_y = (dist * sin(radians(angle)),dist * cos(radians(angle)))
    xfinal, yfinal = (p.x + dist_x, p.y + dist_y)
    return Point(xfinal,yfinal)
   
with fiona.collection('tesligne_coupe.shp', 'r') as coupe:
    # comme il n'y a qu'une seule ligne de coupe, j'utilise next()
    geomcp = shape(coupe.next()['geometry'])
    # ouverture de la couche avec des mesures structurales
    with fiona.collection('struct_point.shp', 'r') as struct:
        schema = struct.schema.copy()
        # création directe des points d'intersection entre 
        # la ligne de coupe et les directions structurales
        with fiona.collection('pt_intersect.shp', 'w', 'ESRI Shapefile', schema) as sortie:
            for pt in struct:
                geompt = shape(pt['geometry'])
                angle = pt['properties']['pend_orienté']
                # calcul du point situé à une distance de 1000 m
                # suivant la direction indiquée par le point (pendage orienté - 90)
                geomptres = ptor2ptfin(geompt, angle-90, 1000) 
                # création de la ligne virtuelle reliant les 2 points
                ligne_virt = LineString([ geomptres.coords[0],ptres.coords[0]])
                # point d'intersection entre cette ligne virtuelle et la ligne de coupe
                pt_inters = ligne_virt.intersection(geomcp)
                ptres = {} 
                # copie des attributs des points d'origine
                ptres['properties'] = pt['properties']
                ptres['geometry'] = mapping(pt_inters)
                sortie.write(ptres)

Obtenir le pendage apparent sur la ligne de coupe ne consiste ensuite qu'à ajouter le résultat du calcul d'une fonction dans le dictionnaire.

Troisième exemple: les transformations affines, translations, rotations, mises à l'échelle et cisaillements

Le module Shapely dispose depuis sa dernière version de tous les éléments permettant d'effectuer les transformations affines en 2D ou en 3D (github.com/Toblerity/Shapely/blob/master/shapely/affinity.py). Cela permet toutes les combinaisons possibles avec les objets Shapely, même en utilisant des matrices de transformation. Pour vous donnez une idée de ce qu'il est possible de faire, je vous conseille d'examiner le script github.com/Toblerity/Shapely/blob/master/shapely/tests/test_affinity.py. Les algorithmes utilisent les principes déjà exposés dans Les transformations affines (avec NumPy) ou la signification géométrique d'un fichier worldfile (.tfw, .jgw,...) pour le géoréférencement 

Combiné avec Fiona, cela devient un vrai plaisir pour un géologue.

exemple de rotation d'un fichier shapefile d'un angle de 35°:

with fiona.collection('testligne.shp', 'r') as entree:
    schema = entree.schema.copy()       
    with fiona.collection('rotation35.shp', 'w', 'ESRI Shapefile', schema) as sortie:
       for ligne in entree:
          geom = shape(ligne['geometry'])
          rotation = affinity.rotate(geom, 35.0)
          sortie.write({'properties': ligne['properties'],'geometry': mapping(rotation)})

Résultat:

 

mais il est aussi possible d'utiliser de vraies matrices de transformations affines 2D et 3D:

matrix2d = (0.422, 0.022,
            0.052,-0.998,
            0.905, 0.0474)
with fiona.collection('testligne.shp', 'r') as entree:
    schema = entree.schema.copy()       
    with fiona.collection('matrix.shp', 'w', 'ESRI Shapefile', schema) as sortie:
        for ligne in entree:
            geom = shape(ligne['geometry'])
            aff = affinity.affine_transform(geom, matrix2d)
            sortie.write({'properties': ligne['properties'],'geometry': mapping(aff)})

Résultat:

 

Je peux donc transformer le « down-plunge profile » présenté dans Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG en fichier shapefile:

Il y a donc, théoriquement, moyen de « géoréférencer » un fichier shapefile par transformation affine.

Conclusions

Sean Gillies dans son tutoriel de Fiona (toblerity.github.com/fiona/manual.html) souligne ses avantages et ses désavantages:

In what cases would you benefit from using Fiona?

• If the features of interest are from or destined for a file in a non-text format like ESRI Shapefiles, Mapinfo TAB files, etc.
• If you’re more interested in the values of many feature properties than in a single property’s value.
• If you’re more interested in all the coordinate values of a feature’s geometry than in a single value.
• If your processing system is distributed or not contained to a single process.

In what cases would you not benefit from using Fiona?

• If your data is in or destined for a JSON document you should use Python’s json or simplejson modules.
• If your data is in a RDBMS like PostGIS, use a Python DB package or ORM like SQLAlchemy or GeoAlchemy. Maybe you’re using GeoDjango already. If so, carry on.
• If your data is served via HTTP from CouchDB or CartoDB, etc, use an HTTP package (httplib2, Requests, etc) or the provider’s Python API.
• If you can use ogr2ogr, do so.

Mais tel quel, il me permet d'utiliser tous mes scripts existants, sans modification et de sauver simplement tous les résultats sans devoir me rappeler toutes les fioritures de ogr ou de PyQGIS. Transformer un fichier de données en fichier shapefile ne demande que quelques lignes de code... Pour voir toutes les autres possibilités de Fiona qui n'ont pas été abordées ici, il suffit de chercher sur le blog de Sean Gillies (sgillies.net/blog).

Ah, et les projections diront certains, il n'y a pas moyen d'utiliser les systèmes de projection, c'est donc rédhibitoire...  Et non, car ce n'est pas le but de Fiona, comme déjà souligné, mais comme il laisse à d'autre modules le soin de faire le travail (pyproj ici), d'une manière tout aussi simple:  sgillies.net/blog/1112/coordinate-reference-systems-for-fiona/ ,  sgillies.net/blog/1125/fiona-pyproj-shapely-in-a-functional-style/ ou github.com/sgillies/Fiona/blob/master/examples/with-pyproj.py.

Tous les traitements ont été faits sur Mac OS X avec la version 0.9.1 de Fiona.

Site officiel : Fiona
Site officiel : pyproj
Autres Liens : Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG
Autres Liens : Les transformations affines (avec NumPy) ou la signification géométrique d'un fichier worldfile (.tfw, .jgw,...) pour le géoréférencement

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