Python : visualiser et traiter des données spatiales de type xyz (shapefiles ou MNT SRTM de la Nasa)

La bibliothèque graphique la plus utilisée dans le monde Python est Matplotlib. Elle permet de créer pratiquement tous les graphiques 2D (et 3D avec une bibliothèque supplémentaire). Elle possède son propre module cartographique, Basemap qui autorise un grand nombre de traitements géospatiaux. Son utilisation est un peu plus complexe et nous limiterons ici au module Matplotlib seul.

Shapefiles

En utilisant l'une des nombreuses bibliothèques qui permettent de les lire (www.portailsig.org/content/python-et-les-shapefiles), il est parfaitement possible de les traiter et de les visualiser puisque ce sont des points, lignes ou polygones avec des coordonnées xy(z). Par exemple, dans un shapefile de type polygone, chaque élément est représenté par les coordonnées xy de ses sommets sous la forme :
 

POLYGON:
////Propriétés////
        numPoints=7
        numParts=1
////shapeparts///
        GC:[0.090944700460829503, 0.80755760368663587]
        GC:[0.14163594470046084, 0.80147465437788012]
        GC:[0.36062211981566822, 0.7021198156682027]
        GC:[0.24099078341013827, 0.54801843317972343]
        GC:[0.086889400921658977, 0.58451612903225802]
        GC:[0.054447004608294942, 0.74267281105990768]
        GC:[0.090944700460829503, 0.80755760368663587]
ou
x=  [0.090944700460829503, 0.14163594470046084, 0.36062211981566822, 0.24099078341013827, 0.086889400921658977, 0.054447004608294942, 0.090944700460829503]
y=  [0.80755760368663587, 0.80147465437788012, 0.7021198156682027, 0.54801843317972343, 0.58451612903225802, 0.74267281105990768, 0.80755760368663587]

Il est alors très facile de le dessiner avec la commande

pylab.plot(x,y)

exemple de résultat avec un shapefile polygone contenant 2 éléments

La bibliothèque Shapely permet de manipuler et de traiter les géométries spatiales. En la couplant avec le traitement précédent, nous disposons des commandes géométriques d'un gis classique en Python. Le même polygone transformé en un format Shapely aura la forme suivante :

'POLYGON ((0.0909447004608295 0.8075576036866359, 0.1416359447004608 0.8014746543778801, 0.3606221198156682 0.7021198156682027, 0.2409907834101383 0.5480184331797234, 0.0868894009216590 0.5845161290322580, 0.0544470046082949 0.7426728110599077, 0.0909447004608295 0.8075576036866359))'

Cette bibliothèque dispose de nombreuses commandes comme :

buffered=polygon.buffer(1.0) -> buffer

polygon1.intersection(polygon2) -> intersection de 2 polygones
etc. voir gispython.org/shapely/manual.html

Après toutes les manipulations, il est possible de recréer et enregistrer une nouvelle shapefile ou de visualiser les résultats (sgillies.net/blog/500/rendering-shapely-geometries-in-matplotlib/)

Résultats de l'intersection des 2 polygones précédents

Les mêmes procédures peuvent être effectuées avec les bibliothèques Ogr ( jgomezdans.googlepages.com/ogr%2Cpythonymatplotlib) ou Basemap (stevendkay.wordpress.com/2009/10/12/scatter-plots-with-basemap-and-matplotlib/) mais de manière un peu moins simple (en particulier pour les polygones).

Il est encore possible d'aller plus loin dans le choix des couleurs des polygones ou en les rendant "cliquables" (par exemple pour afficher les données du fichier dbf associé). On peut aussi jouer avec les projections (www.perrygeo.net/wordpress/).

Modèles numériques de terrain SRTM de la Nasa

Les fichiers SRTM de la Nasa (.hgt) peuvent aussi être traités de la même manière. Comment, des fichiers rasters ?

Et non, ce sont des fichiers de type xyz, en fait des images matricielles binaires. Un fichier hgt est formé d'une grille de 1201 x 1201 cellules (pour la résolution de trois secondes d'arc SRTM3, 3601x3601 pour la résolution d'une seconde d'arc, SRTM 1), chaque cellule représentant une coordonnée z codée en format big-endian avec 2 bytes par cellule. Chaque dalle (fichier hgt) couvre 1° (projection WGS84, EPSG:4326).

Le module STRUCT de python (déjà utilisé pour les shapefiles, www.portailsig.org/content/python-et-les-shapefiles) permet de les décoder d'une manière très simple :

 

fichier = open("fichier.hgt","rb")

contenu=fichier.read()

z= struct.unpack(">1442401H", contenu) ---- [> = codage big-endian, 1442401 = 1201*1201 valeurs, H=entier]

On obtient alors un tableau de 1442401 valeurs de z qu'il est possible de transférer vers des fichiers CSV ou texte ou de l'utiliser en Python avec diverses bibliothèques mathématiques. Il est même possible de l'intégrer à Postgis avec la bibliothèque srtm2postgis.

 

stevendkay.wordpress.com/2009/09/05/beginning-digital-elevation-model-work-with-python/ montre la démarche utilisée pour visualiser une dalle SRTM avec les bibliothèques Matplotlib et Numpy .

SRTM 3 N55W004.hgt

Il va encore plus loin en effectuant divers traitements mathématiques pour transformer ces données : stevendkay.wordpress.com/2009/09/18/mapping-water-drainage-with-srtm-elevation-data-in-python/

 

"mapping water drainage" : image tirée du site précédemment cité

Ces mêmes données peuvent aussi être visualisées par la bibliothèque Mayavi2, qui autorise le vrai 3D : code.enthought.com/projects/mayavi/docs/development/html/mayavi/auto/example_canyon_decimation.html

Grand canyon du Colorado avec le SRTM 1 N36W113.hgt

Détail dans le canyon

La démarche suivante serait de "draper" un shapefile sur un SRTM. C'est théoriquement possible, mais c'est une autre histoire...

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

Site officiel : Python
Site officiel : Matplotlib
Site officiel : Basemap
Site officiel : Shapely
Site officiel : Numpy
Site officiel : Mayavi2
Site officiel : srtm2postgis
Autres Liens : SRTM 3, version 2.1: téléchargement
Autres Liens : SRTM 1, version 2.1: téléchargement