Python: géoréférencement et worldfiles (tfw, jgw,pngw,...)

Le processus de géoréférencement d'une image (raster) a comme résultat divers paramètres qui sont, soit enregistrés dans un fichier particulier, de type worldfile (fichiers .tfw, .jgw, .pngw etc.), soit insérés comme métadonnées dans l'image elle même (format GeoTIFF).  Ces worldfiles, dont le format a été créé par ESRI, restent d'actualité pour tous les types différents du TIFF.

En règle générale, un  SIG  est utilisé pour ces traitements et le fichier est créé automatiquement. Mais en pratique, comment font ces logiciels ?  Pour le comprendre, nous allons, à partir d'un exemple simple, reconstruire en Python le worldfile obtenu par l'un d'entre eux. La démarche permettra de mieux appréhender les principes et éventuellement d'éviter le phénomène de «boite noire» (les données sont rentrées d'un coté, les résultats sont récupérés de l'autre, sans maitrise sur le procédé utilisé, que se passe-t-il en cas d'erreur, données erronnées ou artefacts du traitement ?)

Dans un premier temps, ce travail sera effectué sans aucun module géospatial comme osgeo de GDAL, précédemment utilisé pour les projections, car il fournit directement les réponses, sans explication. Il sera abordé dans la partie finale.

Paramètres d'un worldfile

L'image utilisée est rectangulaire et au format TIFF noir et blanc. Le worldfile créé par le SIG est:

• 1.00078911763392
• 0.00000000000000
• 0.00000000000000
• -1.00078911763392
• 162013.17827588637000
• 108172.36899085061000

Ces valeurs correspondent aux paramètres suivants:

• ligne 1: largeur de chaque pixel suivant l'unité de mesure choisie (qui n'est pas précisée);
• ligne 2: paramètre de rotation autour de l'axe y, généralement 0 pour une translation simple;
• ligne 3: paramètre de rotation autour de l'axe x, généralement 0 pour une translation simple;
• ligne 4: longueur de chaque pixel suivant l'unité de mesure choisie (qui n'est pas précisée). Elle est généralement négative du fait du point de référence de l'image qui est située dans le coin supérieur gauche;
• ligne 5: coordonnée x du centre du pixel situé dans le coin supérieur gauche, sans unité de mesure;   
• ligne 6: coordonnée y du centre du pixel situé dans le coin supérieur gauche, sans unité de mesure.

La projection n'est pas spécifiée, le travail se fait dans un repère cartésien (x, y en mètres, degrés, etc.). C'est pourquoi certains worldfiles rajoutent

• ligne 7: système de coordonnées;
• Line 8: datum.

Mais tous les SIG ne savent pas utiliser ces paramètres et préfèrent utiliser un fichier projection (.prj).

Traitement sans module géospatial

Les traitements d'images peuvent se faire avec différents modules/bibliothèques, mais ici le plus abouti, Python Imaging Library (PIL) sera utilisé. Il permet de lire, traiter, modifier ou créer tous les types de formats.

import sys
from PIL import Image
#ouverture de l'image et gestion des erreurs éventuelles
try:
    raster = Image.open("montif.tif")
except IOError:
    print "impossible d'ouvrir le fichier"
    sys.exit()
#informations diverses
print raster.format, raster.size, raster.mode, raster.info
TIFF (7996, 9995) 1 {'compression': 'tiff_ccitt'}
#c'est donc un fichier TIFF 1 bande noir et blanc (sinon serait précisé RGB...), compression ccitt
# La taille de l'image en pixels est donc:
raster.size
(7996, 9995)
rasterx = raster.size[0]
rastery = raster.size[1]
print rasterx,rastery
7996,9995

Ce qui donne la largeur et la longueur en pixels de l'image. À ce stade, tous les traitements permis par le module pourraient être effectués ou tout autre traitement mathématique en  transformant l'image en tableau numpy (voir francoislouislaillier.developpez.com/Python/Tutoriel/InitiationNumpy/Tuto1/)

Ce n'est pas le but poursuivi et pour la suite du traitement, les valeurs des coins supérieur gauche et inférieur droit fournies par le SIG (en mètres) seront utilisées, avec le nombre de décimales fourni pour "coller" à la réalité:

x1=162012.677881
y1=108172.869385
x2=170014.987666
y2=98169.9821547

Puisque nous avons choisi des coins opposés, le simple rapport

intervalle entre les points xy en mètres / dimensions de l'image en pixels

fournit les paramètres 1 et 3 du worldfile (taille d'un pixel)

ppx = (x2-x1)/rasterx
ppy = (y2-y1)/rastery
print ppx,ppy
1.00078911768 -1.00078911759

Comment obtenir maintenant les paramètres 5 et 6, dont les coordonnées ne correspondent visiblement pas à celles du coin supérieur gauche ? La solution est dans la définition des paramètres: "coordonnées du centre du pixel". 

image inspirée par celle de www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides4.pdf

Il faut donc ajouter (x) ou enlever (y) la valeur d'un 1/2 pixel (ppx, ppy) aux coordonnées du point supérieur gauche pour avoir la bonne valeur du centre du pixel:  

xcentre = x1 + (ppx * .5)
ycentre = y1 + (ppy * .5) #puisque ppy est négatif
print xcentre, ycentre
162013.178276 108172.36899

Ces résultats sont équivalents à ceux du worldfile. En Python, la commande print arrondit les valeurs réelles qui sont (1.0007891176838408, -1.000789117588794) et (162013.17827555886, 108172.3689904412) mais il faut se méfier de la précision obtenue sans module adéquat. L’utilisation du module standard decimal donnerait des résultats plus fiables. 

Il est donc possible, de créer un script complet de géoréférencement d'images simples (translation) dont les coins xy sont connus:

worldfile = open('montif.tfw', "w")   # le nom du fichier peut être obtenu automatiquement à partir de nom du raster ('montif.tif'.split('.'))
worldfile.write(str(ppx)+"\n")
worldfile.write(str(rotaty)+"\n") # = 0
worldfile.write(str(rotatx)+"\n") # = 0
worldfile.write(str(ppy)+"\n")
worldfile.write(str(centrex)+"\n")
worldfile.write(str(centrey)+"\n")
worldfile.close()

et  éventuellement de traiter, par exemple, un ensemble de fichiers dans un dossier:

import sys, os, glob
dossier = 'dossier d'images à traiter'
from PIL import Image
for dir, subdir, files in os.walk(dossier):
   for file in files:
       if glob.fnmatch.fnmatch(file,"*.tif"):
           fich = dir+"/"+file #nom complet du fichier
           try:
               f = Image.open(fich)
           except IOError:
               print "erreur de lecture du fichier:", file
               sys.exit()
           print f.info
           etc. (traitement du fichier)

À l'inverse, on peut aussi lire et interpréter tout fichier worldfile:

from __future__ import with_statement #nécessaire pour les versions < 2.6
from PIL import Image
tfw = "montif.tfw"
raster = Image.open("montif.tif")
 
# dimension du raster
rasterx = raster.size[0]
rastery = raster.size[1]
 
#paramètres du fichier tfw
#lecture du fichier
with open(tfw) as f:
     data = [line.strip('\n').strip() for line in f]
#paramètres
taillepixelx = float(data[0])
taillepixely = float(data[3])
topcentregauchex  = float(data[4])
topcentregauchey  = float(data[5])
 
#calculs
xsupgauche = topcentregauchex - (taillepixelx * .5)
ysupgauche = topcentregauchey - (taillepixely * .5)
print xsupgauche,ysupgauche
162012.677881 108172.869385
#taille du raster en m
lx = rasterx * taillepixelx
print "dimension en x: ", lx
8002.3097846
ly = rastery * taillepixely
print "dimension en y: ", ly
-10002.8872308
#calcul des coordonnées du point inférieur droit
xinfdroit = xsupgauche + lx
yinfdroit = ysupgauche + ly
print xinfdroit,yinfdroit
170014.987666 98169.9821547

Contrôle avec gdalinfo

et maintenant, à titre de comparaison, le résultat de gdalinfo sur le fichier est:

$ gdalinfo montif.tif
Driver: GTiff/GeoTIFF
Files:montif.tif
Size is 7996, 9995
Coordinate System is `'
Origin = (162012.677881327545037,108172.869385409416282)
Pixel Size = (1.000789117633920,-1.000789117633920)
Metadata:
  TIFFTAG_SOFTWARE=Arc/Info 6.0.1
  TIFFTAG_XRESOLUTION=0.99900001
  TIFFTAG_YRESOLUTION=0.99900001
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
  COMPRESSION=CCITTRLE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  162012.678,  108172.869)
Lower Left  (  162012.678,   98169.982)
Upper Right (  170014.988,  108172.869)
Lower Right (  170014.988,   98169.982)
Center      (  166013.833,  103171.426)
Band 1 Block=7996x9 Type=Byte, ColorInterp=Palette
  Image Structure Metadata:
    NBITS=1
  Color Table (RGB with 2 entries)
    0: 0,0,0,255
    1: 255,255,255,255
 

Certaines  métadonnées du fichier TIFF auraient pu être obtenues avec PIL mais elles ne sont pas très intéressantes pour notre démarche.

Traitement avec le module osgeo.gdal

Il est donc possible de créer et d'interpréter les fichiers worldfiles en Python sans bibliothèque/module géospatial. Mais avec le module osgeo.gdal, tout devient beaucoup plus simple: pas besoin de traiter séparément les fichiers images et worldfiles, ce qui permet de traiter aussi les GeoTIFF.

from osgeo import gdal
gdal.AllRegister() #tous les formats d'images de Gdal sont ainsi reconnus
raster = gdal.Open("montif.tif")
#dimension du raster
rasterx = raster.RasterXSize
rastery = raster.RasterYSize
bands = raster.RasterCount
print rasterx,rastery,bands
7996 9995 1
#paramètres du géoréférencement
geotransform = raster.GetGeoTransform()
#paramètres
xsupgauche = geotransform[0]
ysupgauche = geotransform[3]
taillepixelx = geotransform[1]
taillepixely = geotransform[5]
print xsupgauche,ysupgauche, taillepixelx,taillepixely
162012.677881 108172.869385 1.00078911763 -1.00078911763

Une remarque s'impose, le point d'origine donné par GDAL ou osgeo est  le point supérieur gauche du pixel et non pas le centre du pixel comme dans un fichier worldfile. C'est ce qui est indiqué sur la figure.

Pour la suite des traitements équivalents  avec osgeo, je vous renvoie au cours de www.gis.usu.edu/~chrisg/python/2009/ et notamment à "Reading raster data with GDAL" (www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides4.pdf et aux corrections de travaux demandés aux étudiants www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw4a.py, www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw4b.py).

Conclusions

Il est évident que cette démarche ne peut être appliquée qu'à des images simples dont les coordonnées xy des coins sont connues et qui ne nécessitent qu'une simple translation. Dans tous les autres cas de figure (images aux contours irréguliers, points de contrôle situés ailleurs, nécessitant des rotations et donc des transformations...), c'est une autre histoire. C'est possible avec le module osgeo, qui permet directement de lire, de traiter ou de créer des rasters en les géoréférençant (www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides5.pdf) mais le but ici était simplement de comprendre le processus.

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