Python : les bases de données géospatiales - 1) traitement classique, principes et problèmes

Accéder avec Python aux bases de données, relationnelles ou non, est facilité par l'existence de la spécification Python Database API (DB API).

Cet API formalise les manières de se connecter à une base de données et d'effectuer tous les traitements: créations, requêtes, mises à jour ou suppression. Théoriquement, car cet API est plus ou moins bien respecté suivant les modules/bibliothèques.

Démarche classique: principes et exemples

Les principes de l'API sont relativement simples si l'on assimile bien la notion de "cursor" (en anglais, curseur, pointeur en français). Toute requête SQL est susceptible de fournir plus d'une ligne de résultats et plus d'un élément par ligne ou un tableau de valeurs. L'ensemble forme un curseur, qui est donc un objet informatique qui permet d'effectuer un traitement, d'en récupérer les résultats, généralement sous forme de tableaux ou de listes et de les traiter à l'aide de pointeurs.

Examinons maintenant la procédure en Python. Toutes les requêtes sont envoyées en SQL.

 

Appliquons ce principe pour effectuer des requêtes à deux bases de données spatiales, Postgis et SpatiaLite:

import psycopg2 #module le plus courant pour la connexion à PostgreSQL
conn = psycopg2.connect("dbname='testpostgis'host='localhost' user='moi'")
cur = conn.cursor()
cur.execute("""SELECT astext(the_geom) from testpoly;""")
listpoly = cur.fetchall()
for i in listpoly:
   print i
('POLYGON((0.09094470046083 0.807557603686636,0.141635944700461 0.80147465437788,0.360622119815668 0.702119815668203,0.240990783410138 0.548018433179723,0.086889400921659 0.584516129032258,0.054447004608295 0.742672811059908,0.09094470046083 0.807557603686636))',)
('POLYGON((0.275460829493088 0.793364055299539,0.537027649769585 0.813640552995392,0.417396313364055 0.592626728110599,0.226797235023042 0.67778801843318,0.226797235023042 0.744700460829493,0.275460829493088 0.793364055299539))',)
cur.close()
cur = conn.cursor()
query = "SELECT nom,  AsText(the_geom), Srid(the_geom) FROM testpoly where nom='%s'" % 'test1'
cur.execute(query)
poly = cur.fetchone()
print poly
('test1', 'POLYGON((0.09094470046083 0.807557603686636,0.141635944700461 0.80147465437788,0.360622119815668 0.702119815668203,0.240990783410138 0.548018433179723,0.086889400921659 0.584516129032258,0.054447004608295 0.742672811059908,0.09094470046083 0.807557603686636))', 4326)

Les résultats fournis par le curseur sont ici sous la forme d'une liste Python. Pour accéder aux divers éléments, il faudra donc utiliser des indices (pointeurs du curseur):

géométrie = poly[0]
srid = poly[1]

Insérer des données se fait de la même façon. Le processus sera appliqué à SpatiaLite, module spatial de SQLite (attention, il faut une version de SQLite qui accepte de charger des extensions, ce qui n'est pas le cas avec celle livrée avec Python > 2.5 pour le module sqlite3. La solution est donc de compiler une version de pysqlite2 à partir d'un SQLite adéquat. C'est facile avec Linux ou Mac OS X, pour Windows il faudra les télécharger sur le site de SpatiaLite (il y a aussi des applications complètes à télécharger: spatialite-gui et spatialite-gis).

from pysqlite2 import dbapi2 as sqlite
conn = sqlite.connect('test.db3')
#chargement de l'extension spatialite
conn.enable_load_extension(True);
conn.execute('SELECT load_extension("/Users/Shared/source/spatialite/lib/libspatialite.dylib")') #ici sur Mac Os X, so sur Linux, dll sur Windows)
cur = con.cursor()
query='SELECT InitSpatialMetaData()'
cur.execute(query)
query="INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, ref_sys_name, proj4text) VALUES (4326, 'epsg', 4326, 'WGS 84', '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')"
cur.execute(query)
conn.commit() # obligatoire ici puisque l'on veut modifier la table en insérant des données
conn.close()

Démarche classique: problèmes et solutions

Les procédures, bien que conformes à l'API, sont légèrement différentes. En pratique, des tables Postgis et SpatiaLite ne pourraient pas être traitées tout à fait de la même manière.

De plus pour accéder à chaque système de bases de données, il peut exister plusieurs modules (8 pour Postgresql, 2 pour Oracle...), chacun avec ses particularités propres et chacun respectant plus ou moins bien les règles de l'API.

Un autre problème est le manque d'unité du monde des bases de données. Bien qu'un standard SQL existe, les différences sont plus ou moins marquées entre les dialectes SQL, chacun y allant de ses ajouts ou de ses modifications pour se démarquer, avec les meilleurs arguments du monde. 

La démarche classique se voit donc pénalisée:

1. avant de commencer, il faut connaitre les caractéristiques du module de connexion choisi et les multiples avatars des dialectes SQL. Ainsi, le traitement de bases Oracle et MySQL peut être fort différent et le processus n'est pas portable;
2. le principe du curseur est lourd et non "pythonesque". Il ne permet pas la pagination, c'est-à-dire qu'il ne permet pas l'accès à un élément particulier avant d'avoir obtenu tous les résultats en bloc.  Cette technique n'est pas limitée au domaine des bases de données. Le module d'ESRI, argisscripting l'utilise et cela lui est souvent reproché (michalisavraam.org/2010/02/7-wishes-for-the-new-geoprocessor/ ou gissolved.blogspot.com/2009/04/using-for-loops-for-cursors.html );
3. toute opération externe sur la base de données pourrait affecter tous les traitements SQL que les scripts pourraient réaliser.

La visibilité et le traitement des résultats peuvent être améliorés par l'utilisation de dictionnaires comme résultats:

import psycopg2
conn = psycopg2.connect("dbname='testpostgis'host='localhost' user='moi'")
import psycopg2.extras
dict_cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
dict_cur.execute("""SELECT nom,the_geom from testpoly;""")
rec = dict_cur.fetchone()
print rec['nom']
'test1'
print rec['the_geom']
'01030000000100000007000000B64939E22648B73F210C2F0B83D7E93F0399366B2021C23F5B07A42CAEA5E93FA3A3B4CC6E14D73F3164C5F9C377E63FAB25B136C9D8CE3F179603F45D89E13FD2DA533F623EB63FB4B2452B5BB4E23FA1C950507AE0AB3FA12E64C5F9C3E73FB64939E22648B73F210C2F0B83D7E93F'
dict_cur.execute("""SELECT nom,astext(the_geom) from testpoly;""")
rec2 = dict_cur.fetchone()
print rec['astext'] #puisque l'on demande la fonction Postgis AsText, on pourrait faire la même chose avec askml etc.
'POLYGON((0.09094470046083 0.807557603686636,0.141635944700461 0.80147465437788,0.360622119815668 0.702119815668203,0.240990783410138 0.548018433179723,0.086889400921659 0.584516129032258,0.054447004608295 0.742672811059908,0.09094470046083 0.807557603686636))'

Dans le premier cas (the_geom), c'est le résultat brut au format binaire utilisé par Postgis (format WKB de l' OGC) qui est obtenu, dans le second, c'est un format texte qui est obtenu par la fonction AsText de Postgis (format WKT de l'OGC).

Tous les modules de connexion aux bases de données ne permettent cependant pas de faire ce genre de traitement. On peut alors utiliser la fonction toDict de code.activestate.com/recipes/528939-converting-dbi-results-to-a-list-of-dictionaries/, plus universelle:

## fonction toDict: Recipe 528939 - Converting DBI results to a list of dictionaries 
def toDict(curs):
    cols = [column[0] for column in curs.description]
    return [dict(zip(cols, row)) for row in curs.fetchall()]
 
 
cur = conn.cursor()
query = "SELECT nom,  AsText(the_geom), Srid(the_geom) FROM testpoly where nom='%s'" % 'test1'
cur.execute(query)
print toDict(cur)
[{'srid': 4326, 'nom': 'test1', 'astext': 'POLYGON((0.09094470046083 0.807557603686636,0.141635944700461 0.80147465437788,0.360622119815668 0.702119815668203,0.240990783410138 0.548018433179723,0.086889400921659 0.584516129032258,0.054447004608295 0.742672811059908,0.09094470046083 0.807557603686636))'}]

ou un peu plus compliqué comme, code.activestate.com/recipes/534114-list-wrapper-for-generators-indexable-subscriptabl/ qui permet d'accéder à n'importe quelle partie du résultat sans devoir en charger la totalité.

On le voit donc, l'idéal serait un procédé permettant d'accéder de la même manière à n'importe quelle base et de ne traiter que des objets Python, sans se préoccuper des particularités SQL ou autres. Des modules permettant de tels niveaux d'abstraction existent, depuis le simple parseur qui transforme les résultats d'un curseur (en dictionnaire comme l'exemple précédent), jusqu'à ceux qui libèrent le programmeur de toute commande SQL  (wiki.python.org/moin/HigherLevelDatabaseProgramming). Parmi ces derniers, ceux qui utilisent le mapping objet-relationnel (ORM) sont les plus aboutis (modules SQLAlchemy, SQLObject ou GeoAlchemy, frameworks web comme Django, TurboGears, web2py ou MapFish).

Ce dernier thème est développé dans la suite :

"Python : les bases de données géospatiales - 2) mapping objet-relationnel (ORM, SQLAlchemy, SQLObject, GeoAlchemy, Django-GeoDjango, TurboGears ou MapFish)" www.portailsig.org/content/python-les-bases-de-donnees-geospatiales-2-mapping-objet-relationnel-orm-sqlalchemy-sqlobjec

Tous les traitements ont été effectués sur Mac OS X avec Python 2.5.4 et Inkscape pour la figure.