Répondre au commentaire

Python: projections cartographiques, définitions et transformations de coordonnées

sam 06-02-2010 Martin Laloux

Avant de commencer, un petit rappel s'impose. Tout système de coordonnées spatiales est défini par plusieurs paramètres qui sont, en simplifiant à l'extrême:

un ellipsoïde de référence, forme mathématique simplifiée de la terre ;
un datum, adaptation de cet ellipsoïde à la surface étudiée ;
une projection qui transforme les cartes sur la surface de l'ellipsoïde en représentations cartographiques planes;
un point d'origine (0,0,(0)), des axes de coordonnées et l'unité de mesure (degré, mètres...).

image provenant de cfis.savagexi.com/2006/05/02/coordinate-systems-putting-everything-together

La formalisation de ces systèmes (CRS, coordinate reference system ou SRS, spatial referencing system, en anglais) peut se faire suivant différentes syntaxes, WKT (format standard de l'OGC), EPSG, Proj, ESRI et autres (ESRI a développé son propre système WKT, légèrement différent du standard de l'OGC, avec quelques projections différentes. C'est ce que l'on trouve dans les fichiers .prj.). Des formulations de toutes les syntaxes disponibles pour une projection donnée sont fournies par le site www.spatialreference.org/ , sous-projet de l' OSGeo (wiki.osgeo.org/wiki/MetaCRS). Notons que le site est réalisé et géré entièrement en Python avec le framework Django-Geodjango.

Toutes ces notions sont bien expliquées dans:

Notions de base sur les systèmes de coordonnées et le géoréférencement: www.frichtiweb.com/files/12008668211_systemes_de_coordonnees_2008.pdf
Synthèse présentée au Foss4G 2007 (congrès du Free and Open Source Software for Geospatial) : 2007.foss4g.org/labs/L-03/FOSS4G2007_APSG_VictoriaBC.pdf

La librairie Open Source la plus ancienne pour traiter ces systèmes est libproj ou Proj, développée en C depuis 1990. C'est la plus utilisée et les 2 versions actuelles sont Proj.4 de l' OSGeo et libproj4, moins connue (elle a été créée par le programmeur originel suite à un différend expliqué sur son site). La syntaxe d'encodage des CRS peut paraitre rebutante à certains, mais, à ma connaissance, la référence 2 en constitue la meilleure introduction.

Il n'est donc pas étonnant qu'elle serve de base aux principaux modules/bibliothèques Python qui traitent du sujet. Ces modules ne peuvent pas être installés simplement suivant la procédure traditionnelle (easy_install ou setup.py) hormis si Proj.4 et autres librairies C sont préalablement bien installées (ils nécessitent la compilation de certains modules avec ces librairies, ce sont les modules se terminant par .so sur Mac ou Linux et .dll sur Windows).

le gis en Python, Thuban, a développé son propre module projection. Il a été proposé tout un temps de manière autonome sous le nom de py-projection par Howard Butler. Il n'est apparemment plus disponible (Thuban oui);
pyproj est une simple interface aux librairies C de Proj.4 générée par Pyrex.
GDAL/OGR est devenu, en Python, le module osgeo. Son interface pour les projections est osgeo.osr (osr ici). La version Python de GDAL est l'une des premières bibliothèques dérivées et reste la plus utilisée, certaines commandes de Gdal étant d'ailleurs toujours écrites en Python (gdal_merge.py etc.). De très beaux exemples de l'utilisation de la bibliothèque sont proposés dans les cours de www.gis.usu.edu/~chrisg/python/2009/, en particulier sur les projections "Creating Geometries and Handling Projections with OGR" (www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides2.pdf);
Basemap, le module cartographique de Matplotlib a sa propre implémentation.

Il y a d'autres modules comme Python Cartographic Library (PCL), Mapnik, des frameworks cartographiques complets ("web mapping") comme Django-GeoDjango ou Mapfish, ou même la possibilité d'utiliser directement les fonctions de Postgis en Python. Nous nous limiterons ici à pyproj et osr et les données de type x,y ou géométriques (pas les rasters qui nécessitent d'autres modules et traitements).

Définition d'une projection

Dans les deux modules, il est nécessaire, en premier lieu, de créer une projection:

import pyproj
a = pyproj.Proj("paramètres de la projection") # création de la projection

et

from osgeo import osr ou import osgeo.osr as osr
maref = osr.SpatialReference() # définition d'une projection, vide ici
maref.ImportFrom... # création de la projection proprement dite

Les paramètres de pyproj sont ceux de Proj. Voici quelques exemples:

a = pyproj.Proj('+proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62')
b = pyproj.Proj(proj='utm',zone=10,ellps='WGS84')
c = pyproj.Proj(init='epsg:26915')
d = pyproj.Proj('+init=IGNF:LAMBE')
e = pyproj.Proj('+init=IGNF:LAMBE +wktext') #suivant les recommandations de l' IGN

Ceux de osr sont plus riches et variés:

a = maref.ImportFromWkt('PROJCS["Belge 1972 / Belgian Lambert 72",GEOGCS["Belge 1972",DATUM["Reseau_National_Belge_1972",SPHEROID["International 1924",6378388.0,297.0,
AUTHORITY["EPSG","7022"]],TOWGS84[-99.1,53.3,-112.5,0.419,-0.83,1.885,-1.0],
AUTHORITY["EPSG","6313"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.017453292519943295],AXIS["Lon",EAST],AXIS["Lat",NORTH],
AUTHORITY["EPSG","4313"]],PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["central_meridian",4.367486666666666],PARAMETER["latitude_of_origin",90.0],
PARAMETER["standard_parallel_1",51.166667233333335],PARAMETER["false_easting",150000.013],
PARAMETER["false_northing",5400088.438],PARAMETER["standard_parallel_2",49.83333389999999],
UNIT["m",1.0],AXIS["x",EAST],AXIS["y",NORTH],AUTHORITY["EPSG","31370"]]')

b = maref.ImportFromEPSG(4326)
c = maref.ImportFromProj4('+init=IGNF:LAMB93')
d = maref.ImportFromESRI('projection au format ESRI')

D'autres formats sont possibles comme ImportFromPCI(<params>), ImportFromUSGS(<params>) ou importFromXML(<xml>).

Avec osr, les erreurs éventuelles peuvent être contrôlées par la réponse donnée par le module après définition d'une projection:

maref = osr.SpatialReference()
maref.ImportFromEPSG(4326)
0

0 indique que maref est une projection valide, > 0 indique qu'il y a un problème (projection vide car non valide, par exemple). Avec pyproj, Python signalera simplement qu'il y a une erreur sans plus de détail permettant de la résoudre.

osr fournit aussi la traduction entre tous les formats définis, sous la forme ExportTo...:

maref = osr.SpatialReference()
maref.ImportFromEPSG(4326)
a = maref.ExportToProj4()
print a
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

Ceci permet, si on le veut, à titre de contrôle, d'utiliser simultanément les 2 bibliothèques.

b = pyproj.Proj(maref.ExportToProj4())
b.srs
'+units=m +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '

Transformation de coordonnées x, y

La suite est alors aisée:

on crée les 2 systèmes de projection, source et destination;
on crée la transformation voulue (de source à destination);
on applique cette transformation à n'importe quel objet spatial.

Prenons l'exemple de l'IGN français lambert93.ign.fr/index.php puisqu'il soutient l'OSGeo (www.osgeo.org/sponsors, wiki.osgeo.org/wiki/Accord_cadre_OSGeo-fr_/_IGN) et participe, entre autres, aux développements de GDAL.

L'outil permettant de réaliser des conversions de coordonnées est cs2cs.
Pour convertir une liste de coordonnées rentrées interactivement,
du Lambert 2 étendu au Lambert 93, écrire:
cs2cs -f "%.3f" +init=IGNF:LAMBE +to +init=IGNF:LAMB93 (return)
565767.9060 2669005.7300 (return)

On obtient alors le résultat suivant:
619119.460      7102502.980 0.000

Appliquons-le avec osr:

from osgeo import osr
LAMBE = osr.SpatialReference()
LAMBE.ImportFromProj4("+init=IGNF:LAMBE +wktext")
LAMB93 = osr.SpatialReference()
LAMB93.ImportFromProj4("+init=IGNF:LAMB93")
x = 565767.9060
y = 2669005.7300
transformation1 = osr.CoordinateTransformation(LAMBE,LAMB93)
res1 = transformation1.TransformPoint(x, y)
print res1[0],res1[1]
619119.46051729936, 7102502.9797217501
print "%.3f, %.3f" % (res1[0],res1[1]) #formatage similaire à celui de cs2cs
619119.461 7102502.98

# "reprojection" inverse

transformation2 = osr.CoordinateTransformation(LAMB93,LAMBE)
res2 = transformation2.TransformPoint(res[0],res[1])
print res2[0],res[1]
565767.90599998576, 2669005.729997803
print "%.3f, %.3f" % (res2[0],res2[1])
565767.906 2669005.73

L'IGN précise que l'IGNF:LAMB93 peut être remplacé par l'EPSG:2154, testons-le:

LAMB93_2 = osr.SpatialReference()
LAMB93_2.ImportFromEPSG(2154)
transformation3 = osr.CoordinateTransformation(LAMBE,LAMB93_2)
res3 = transformation3.TransformPoint(x, y)
print res3[0],res3[1]
619119.46051729936, 7102502.9797217492
print "%.3f, %.3f" % (res3[0],res3[1])
619119.461, 7102502.980

Et effectivement, les résultats de la transformation sont équivalents.

Avec pyproj, les principes sont à peu près les mêmes, mais le traitement de ces dernières données provoque chez moi une erreur Python sans explication claire.

Reprenons donc les résultats de la géolocalisation de "Place de l'Université, 1348 Louvain-la-Neuve, Belgique" développée dans www.portailsig.org/content/python-geocodage-geolocalisation, avec comme résultat:

lat: 50.6701607 long:4.6151754

import pyproj
wgs84 = pyproj.Proj('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
# avec les valeurs de spatialreference.org/ref/sr-org/epsg31370-correct-datum-shift/
be31370 =pyproj.Proj('+proj=lcc +lat_1=51.16666723333334 +lat_2=49.83333389999999 +lat_0=90 +lon_0=4.367486666666666
+x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.1,53.3,-112.5,0.419,-0.83,1.885,-1.0 +units=m +no_defs')
lat,long=wgs84(50.6701607,4.6151754)
x2, y2 = pyproj.transform(wgs84,be31370,lat,lon,radians=True)
print "%.2f, %.2f" % (x2,y2)
167419.91, 151090.66

Je ne peux pas, ici, placer des extraits de carte de l' I.G.N. belge sans autorisation, mais ces coordonnées Lambert belge correspondent bien au lieu indiqué sur Google Maps, ce qui n'est pas le cas avec l'EPSG:31370 officiel (spatialreference.org/ref/epsg/31370/, elles diffèrent au niveau du paramètre +towgs84 d'où le nom "correct datum shift").

Traitement des objets géométriques (shapefiles par exemple)

Les traitements précédents ont été limités à des coordonnées x,y mais osr permet de traiter tout objet géométrique (point, ligne, polygone, multi...) encodé suivant les standards de l'OGC (voir www.portailsig.org/content/python-visualiser-et-traiter-des-donnees-spatiales-de-type-xyz-shapefiles-ou-mnt-srtm). Pour mémoire, un polygone, par exemple, sera encodé de la manière suivante:

'POLYGON ((0.0909447004608295 0.8075576036866359, 
0.1416359447004608 0.8014746543778801,...))'
même syntaxe pour les points, lignes, multi...

Les transformations seront effectuées avec la commande:

geom.Transform(transformation)

N'importe quel shapefile sera alors traité de la manière suivante (osgeo.ogr est utilisé ici, mais nous aurions pu utiliser n'importe quel module analysé dans www.portailsig.org/content/python-et-les-shapefiles) pour lire un shapefile:

from osgeo import ogr, osr
# lecture des shapefiles et de ses éléments
driver = ogr.GetDriverByName('ESRI Shapefile')
leshape = driver.Open('mon.shp')
couche = leshape.GetLayer()

# si le shapefile a une projection sous forme de fichier .prj
spatialrefor = couche.GetSpatialRef()
# sinon vous savez maintenant créer une projection

element = couche.GetNextFeature()
# premier élément géométrique
geom = element.GetGeometryRef()

# Transformation
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
transformation = osr.CoordinateTransformation(spatialrefor, wgs84)
geom.Transform(transformation)

et ainsi de suite pour tous les autres élements/géométries du shapefile.

Un exemple complet est donné à www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw2b.py qui est la correction d'un travail demandé aux étudiants. Le script lit un shapefile, dbf compris, en change la projection et enregistre le résultat sous la forme d'un nouveau shapefile.

Un autre bel exemple est fourni par www.perrygeo.net/wordpress/ qui montre les déformations induites par une projection sous la forme de l'indicatrice de Tissot (transformation d’un cercle par la projection).

Les exemples précédents sont limités aux shapefiles mais tout système produisant des géométries conformes peut être utilisé (module shapely, résultats provenenant de Postgis etc.)

Tous les traitements ont été effectués sur Mac OS X avec Python 2.5.4