Python à finalité géospatiale: pourquoi débuter avec le module ogr alors qu'il existe des alternatives plus simples et plus didactiques ?

	Niveau	Débutant
	Logiciels utilisés	Python
	Plateforme	Windows | Mac | Linux | FreeBSD

La plupart des débutants dans l'utilisation des modules Python à finalité géospatiale libres commencent par le module GDAL/OGR, nommé osgeo dans la suite, du fait que l'importation du module se fait par la commande (et non open ogr).

from osgeo import ogr

Le module est certes très puissant lorsqu'on le connait bien (voir Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG , par exemple) mais son problème, pour les débutants, est qu'il est relativement difficile à appréhender, complexe à utiliser et relativement peu "Pythonesque" du fait de sa conception (interfaçage de la librairie C++ originale). Ceci n'est pas une critique, loin de là, mais je trouve qu'il est peu adapté pour débuter alors qu'il existe d'autres modules plus faciles à comprendre.

Il suffit de se rendre sur des sites comme GIS StackExchange ou le ForumSig pour s'en rendre compte: la prise en main du module est loin d'être évidente lorsqu'on débute...

De la complexité des choses...

Le très beau site Python GDAL/OGR Cookbook est une des références pour l'utilisation du module, hormis les livres. Il illustre la plupart des traitements qui peuvent être effectués. Prenons le cas de la création d'un simple Polygone qui se fait de la manière suivante (Create a Polygon):

from osgeo import ogr
# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(1179091.1646903288, 712782.8838459781)
ring.AddPoint(1161053.0218226474, 667456.2684348812)
ring.AddPoint(1214704.933941905, 641092.8288590391)
ring.AddPoint(1228580.428455506, 682719.3123998424)
ring.AddPoint(1218405.0658121984, 721108.1805541387)
ring.AddPoint(1179091.1646903288, 712782.8838459781)
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

Il faut

1. créer une géométrie ogr (LinearRing), c'est-à-dire, l'enveloppe linéaire extérieure du polygone;
2. la « peupler » avec des points au format ogr, qui vont la définir;
3. créer une géométrie ogr Polygon;
4. et enfin le construire à partir de son enveloppe extérieure.

Il existe d'autres modules qui permettent de le créer de manière beaucoup plus simple (comme geojson, Shapely, PyGeoif, Shapy Karta, GeoDjango ou GeoPandas)

# Create polygon
poly = {
    'type': 'Polygon',
    'coordinates': [[
        (1179091.1646903288, 712782.8838459781),
        (1161053.0218226474, 667456.2684348812),
        (1214704.933941905, 641092.8288590391),
        (1228580.428455506, 682719.3123998424),
        (1218405.0658121984, 721108.1805541387),
        (1179091.1646903288, 712782.8838459781) ]]}

ou

Polygon([(1179091.1646903288, 712782.8838459781),(1161053.0218226474, 667456.2684348812),(1214704.933941905, 641092.8288590391),(1228580.428455506, 682719.3123998424),(1218405.0658121984, 721108.1805541387),(1179091.1646903288, 712782.8838459781)])

Notons que rien ne vous empêche à ce stade d'utiliser le dictionnaire précédent (poly) pour créer une géométrie ogr

import json
poly = json.dumps(poly)
polygon = ogr.CreateGeometryFromJson(poly)

De la simplification des choses...

Le but de ces nouveaux modules est donc clairement de simplifier la vie des utilisateurs en n'utilisant que des simples objets Python standards, listes, dictionnaires et arrays NumPy pour traiter les données sans faire appel à des commandes externes comme osgeo, PyQGIS ou ArcPy, à titre d'exemple.

1. Leur particularité est qu'ils disposent tous du protocole Geo_interface (voir Les modules Python à finalités géospatiales: quid, quando, ubi ?) où, en simplifiant, tous les traitements se font à l'aide de simples dictionnaires Python (voir Python Geo_interface applications)
2. Certains d'entre eux se basent, tout comme osgeo, sur la librairie GDAL/OGR (simplification des procédures avec un autre interfaçage de la librairie C++ originale), d'autres sur la librairie C++ GEOS, d'autres sont écrits en pur Python et enfin certains combinent le tout avec le module Pandas, module de plus en plus utilisé dans le monde scientifique (équivalent aux dataframes  de R avec une kyrielle d'applications, voir Sam & Max: Le pandas c’est bon, mangez-en, par exemple);
3. Certains fournissent des lignes de commandes comme ogr_info, etc.
4. Ils sont disponibles pour les versions 2.7.x et 3.x de python

Ils sont donc particulièrement adaptés aux débutants comme nous allons le voir. Nous nous focaliserons ici sur les couches vectorielles, laissant les couches rasters pour plus tard.

Lire des couches vectorielles

Nous utiliserons ici Fiona et GeoPandas qui sont les plus connus de cette « nouvelle vague » pour les explications.

Lecture

Je vous laisse découvrir les principes de lecture d'un fichier shapefile avec le module osgeo dans Cookbook: Vector Layers. Vous constaterez qu'il  faut beaucoup de lignes et de commandes à retenir pour simplement extraire le schéma d'un fichier shapefile, sa projection et les géométries et données de la table attributaire....

Avec Fiona

Comme déjà dit auparavant, c'est un nouvel interfaçage de la librairie C++ GDAL/OGR et donc comparable à osgeo.ogr. Tout se fait avec des dictionnaires Python.

import fiona
# j'ouvre le fichier dans un itérateur Python
couche = fiona.open('strati_or.shp')
# le schéma de couche
couche.schema
{'geometry': 'Point', 'properties': OrderedDict([(u'id', 'int:10'), (u'dip', 'int:2'), (u'dip_dir', 'int:3'), (u'dir', 'int:9'), (u'type', 'str:10'), (u'x', 'int:10'), (u'y', 'int:10')])}
# la projection de la couche
couche.crs # au format proj4
{u'lon_0': 4.367486666666666, u'ellps': u'intl', u'y_0': 5400088.438, u'no_defs': True, u'proj': u'lcc', u'x_0': 150000.013, u'units': u'm', u'lat_2': 49.8333339, u'lat_1': 51.16666723333333, u'lat_0': 90}
# premier élément
couche.next()
{'geometry': {'type': 'Point', 'coordinates': (272070.600041, 155389.38792)}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 0), (u'dip', 30), (u'dip_dir', 130), (u'dir', 40),(u'type', u'incliné'), (u'x', 272071), (u'y', 155389)])}
# d'où
elements = [elem for elem in couche]
elements[0]['geometry']
{'type': 'Point', 'coordinates': (272070.600041, 155389.38792)}
elements[0]['properties']
OrderedDict([(u'id', 0), (u'dip', 30), (u'dip_dir', 130), (u'dir', 40), (u'type', u'inclin\xe9'), (u'x', 272071), (u'y', 155389)])
elements[0]['properties']['dip_dir']
130

Vous comprenez alors aisément que tous les traitements se résument à la simple manipulation de dictionnaires (modifier le schéma d'un shapefile ou modifier un élément de la table attributaire par exemple). Pour traiter ensuite les données, vous avez l'embarras du choix avec Shapely (géométries, ou les autres déjà cités)  NumPy, SciPy, Scikit-learn ou autres modules spécifiques.

Avec GeoPandas

Avec ce  module, vous entrez dans un autre monde encore peu connu des géomaticiens, le module Pandas. Vous créez directement un DataFrame, ou une Serie Pandas (respectivement GeoDataFrame et GeoSeries) à l'aide des modules Shapely et Fiona. Le traitement des DataFrames Pandas est abondamment illustré sur Internet.

import geopandas as gp
# création d'un GeoDataframe à partir du fichier shapefile
couche = gp.GeoDataFrame.from_file('strati_or.shp')
# affichage classique avec Pandas des 5 premiers éléments du DataFrame couche 
couche.head()
   dip  dip_dir  dir                             geometry  id      type     x       y  
0   30      130   40   POINT (272070.600041 155389.38792)   0   incliné   272071  155389 
1   55      145   55  POINT (271066.032148 154475.631377)   1  retourné   271066  154476
2   40      155   65  POINT (273481.498868 153923.492988)   2   incliné   273481  153923
3   80      120   30  POINT (270977.604378 153144.810665)   3  retourné   270978  153145
4   40      130   40  POINT (267162.940066 150990.398109)   4   incliné   267163  150990
# aperçu rapide des propriétés de la couche
couche.describe()
             dip     dip_dir         dir         id              x         y  
count  12.000000   12.000000   12.000000  12.000000      12.000000      12.000000 
mean   59.416667  170.416667   80.416667   5.500000  270592.666667  152742.166667 
std    21.094844   71.147042   71.147042   3.605551    2508.060364    2508.06036 
min    30.000000  120.000000   30.000000   0.000000  267163.000000  2150442.000000  
25%    40.000000  133.750000   43.750000   2.750000  268675.500000   150910.750000
50%    57.500000  145.000000   55.000000   5.500000  270824.500000  152728.000000    
75%    78.500000  155.000000   65.000000   8.250000  272423.500000  154375.000000   
max    90.000000  327.000000  237.000000  11.000000  274188.000000  155389.000000   

Vous pouvez alors utiliser toutes les fonctions de GeoPandas ou de Pandas pour traiter les données géométriques et numériques et utiliser le pandas Ecosystem (Statsmodels, Sklearn-pandas, etc) en plus des modules précédemment cités.

Avec les autres

Ce protocole de lecture des données sous forme de dictionnaire (geo_interface) peut être appliqué avec d'autres modules, voir ogr_geointerface.py,  PyShp_geointerface.py,  PostGIS_geointerface.py, spatialite_geointerface.py, ou mapnik_geointerface.py.

Exemple avec PyShp (shapefile):

def records(filename):  
    # generateur 
    reader = shapefile.Reader(filename)  
    fields = reader.fields[1:]  
    field_names = [field[0] for field in fields]  
    for sr in reader.shapeRecords():  
        geom = sr.shape.__geo_interface__  
        atr = dict(zip(field_names, sr.record))  
        yield dict(geometry=geom,properties=atr)    
import shapefile
elem = records('strati_res.shp')
elem.next()
{'geometry': {'type': 'Point', 'coordinates': (272070.600041, 155389.38792)}, 'properties': {'dip_dir': 130, 'dip': 30, 'cosa': -0.6428, 'sina': -0.6428}}

 Exemple avec osgeo:

def records(shapefile):  
    # generateur
    reader = ogr.Open(shapefile)
    layer = reader.GetLayer(0)
    for i in range(layer.GetFeatureCount()):
        feature = layer.GetFeature(i)
        yield json.loads(feature.ExportToJson())
from osgeo import ogr
elem = records('strati_res.shp')
elem.next()
{'geometry': {'type': 'Point', 'coordinates': (272070.600041, 155389.38792)}, 'properties': {'dip_dir': 130, 'dip': 30, 'cosa': -0.6428, 'sina': -0.6428}}

Sauvegarder des couches vectorielles

Avec Fiona

Vous sauvegardez des dictionnaires Python. Par exemple dans l'exemple suivant, je veux créer un nouveau shapefile basé sur le précédent en créant de nouveaux champs et en éliminant certains:

import fiona
from fiona.crs import from_epsg
 # définition des fonctions de calcul
import math
sind = lambda x: math.cos( math.radians(x))
cosd = lambda x: math.cos( math.radians(x))
# définition du schéma du nouveau shapefile
schéma = {'geometry': 'Point', 'properties': {'dip' : 'int:2', 'dip_dir' :'int:3', 'cosa': 'float:11.4','sina':'float:11.4'}}
# définition du crs du nouveau shapefile
crs = from_epsg(31370) # ou en reprenant simplement le crs du shapefile en entrée, comme dans la suite
# je le remplis avec le dictionnaire
with fiona.open('strati_or.shp') as entree:
   with fiona.open('strati_res.shp','w',driver='ESRI Shapefile', crs=entree.crs,schema= schéma) as sortie:
       for elem in entree:
          # conctruction du dictionnaire et sauvegarde 
          geom = elem['geometry'] # puisque c'est la même géométrie
          prop = {'dip': elem['properties']['dip'],'dip_dir': elem['properties']['dip_dir'], 'cosa': cosd(elem['properties']['dip_dir']), 'sina': sind(elem['properties']['dip_dir']) }
          sortie.write({'geometry':geom, 'properties': prop})

Avec GeoPandas

C'est encore plus simple, car après avoir effectué les changements dans le GeoDataFrame précédent  (couche),  il suffit, pour  sauvegarder le shapefile résultant de:

couche.to_file('strati_res.shp')

Quelques exemples pour vous convaincre

Je reprendrai ici quelques-unes de mes réponses sur GIS StackExchange

• Problem intersecting shapefiles with OGR in Python (général)
• Using Fiona to write a new shapefile from scratch (avec Fiona)
• Can't write decimal data into shape (shp) attribute table with Python (avec Fiona
• intersecting two shapefiles from Python or command line (avec Fiona ou GeoPandas);
• More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (avec Fiona) à comparer avec la version GeoPandas, More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc;
• csv to shp missing file  (fichier .cvs vers shapefile, avec PyShp)
• Python buffer points and clip with line-segment (OS Open Roads in BNG) (fichier .cvs vers shapefile, avec Fiona)
• Rewriting csv file for QGIS (manipuler des données avec Pandas)
• ....

Pour ceux qui sont allergiques au format GeoJSON

Vous pouvez utiliser Shapely ou GeoPandas pour transformer vos géométries au format WKT par exemple, mais il y a une solution plus simple, le module geomet

elem= {'geometry': {'type': 'Point', 'coordinates': (272070.600041, 155389.38792)}, 'properties': {'dip_dir': 130, 'dip': 30, 'cosa': -0.6428, 'sina': -0.6428}}
from geomet import wkt
print wkt.dumps(elem['geometry'], decimals=2)
POINT (272070.60 155389.39)

Pour ceux qui préfèrent les commandes (ogr_info, ogr2ogr, etc.)

Fiona dispose de la commande Fio avec autant de possibilités que les commandes d'OGR (voir Fiona-Rasterio-Shapely Cheat Sheet)

$ fio info strati_or.shp
{"count": 12, "crs": "+ellps=intl +lat_0=90 +lat_1=51.1666672333 +lat_2=49.8333339 +lon_0=4.36748666667 +no_defs +proj=lcc +units=m +x_0=150000.013 +y_0=5400088.438", "driver": "ESRI Shapefile", "bounds": [267162.940066, 150441.970164, 274188.072595, 155389.38792], "crs_wkt": "PROJCS[\"Belge_1972_Belgian_Lambert_72\",GEOGCS[\"GCS_Belge 1972\",DATUM[\"Reseau_National_Belge_1972\",SPHEROID[\"International_1924\",6378388,297]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",51.16666723333333],PARAMETER[\"standard_parallel_2\",49.8333339],PARAMETER[\"latitude_of_origin\",90],PARAMETER[\"central_meridian\",4.367486666666666],PARAMETER[\"false_easting\",150000.013],PARAMETER[\"false_northing\",5400088.438],UNIT[\"Meter\",1]]", "schéma": {"geometry": "Point", "properties": {"id": "int:10", "dip": "int:2", "dip_dir": "int:3", "dir": "int:9", "type": "str:10", "x": "int:10", "y": "int:10"}}}

Conclusions

Alors bien sûr, ces modules ne peuvent pas faire tout ce que fait le module osgeo, car ils ne s'adressent qu'à des fichiers (et non à l'accès des données provenant de SBGD comme PostGIS, MongoDB, CouchDB, SpatiaLite et autres ou le traitement des services WFS, par exemple, c'est d'ailleurs souligné dans les "Rules of Thumb" du Manuel de Fiona par Sean Gillies).  Mais ils peuvent servir à créer des fichiers à partir de ces données comme dans la question  Forum SIG: Python avec osgeo.org, ouvrir un WFS - Débutant;

D'après mon expérience, ils sont beaucoup plus faciles à maîtriser pour les débutants et sont plus rapides que PyQGIS (ou ArcPy d'après mes collègues) pour les mêmes traitements. Une fois compris, il est beaucoup plus facile d'aborder le module osgeo.

Alors pourquoi débuter avec osgeo.ogr ?

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