Skip to Content

Nous adressons toutes nos pensées à la famille de notre ami Jérôme !

http://www.forumsig.org/showthread.php/43488-Disparition-de-Phoenix

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


python

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.

taille en pixels d'un raster

  1. import sys
  2. from PIL import Image
  3. #ouverture de l'image et gestion des erreurs éventuelles
  4. try:
  5. raster = Image.open("montif.tif")
  6. except IOError:
  7. print "impossible d'ouvrir le fichier"
  8. sys.exit()
  9. #informations diverses
  10. print raster.format, raster.size, raster.mode, raster.info
  11. TIFF (7996, 9995) 1 {'compression': 'tiff_ccitt'}
  12. #c'est donc un fichier TIFF 1 bande noir et blanc (sinon serait précisé RGB...), compression ccitt
  13. # La taille de l'image en pixels est donc:
  14. raster.size
  15. (7996, 9995)
  16. rasterx = raster.size[0]
  17. rastery = raster.size[1]
  18. print rasterx,rastery
  19. 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)

taille d'un pixel

  1. ppx = (x2-x1)/rasterx
  2. ppy = (y2-y1)/rastery
  3. print ppx,ppy
  4. 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:  

coordonnées du centre du pixel

  1. xcentre = x1 + (ppx * .5)
  2. ycentre = y1 + (ppy * .5) #puisque ppy est négatif
  3. print xcentre, ycentre
  4. 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:

création du fichier worldfile

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

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

traitement de tous les fichiers dans un dossier

  1. import sys, os, glob
  2. dossier = 'dossier d'images à traiter'
  3. from PIL import Image
  4. for dir, subdir, files in os.walk(dossier):
  5. for file in files:
  6. if glob.fnmatch.fnmatch(file,"*.tif"):
  7. fich = dir+"/"+file #nom complet du fichier
  8. try:
  9. f = Image.open(fich)
  10. except IOError:
  11. print "erreur de lecture du fichier:", file
  12. sys.exit()
  13. print f.info
  14. etc. (traitement du fichier)

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

lecture et traitement du worldfile

  1. from __future__ import with_statement #nécessaire pour les versions < 2.6
  2. from PIL import Image
  3. tfw = "montif.tfw"
  4. raster = Image.open("montif.tif")
  5.  
  6. # dimension du raster
  7. rasterx = raster.size[0]
  8. rastery = raster.size[1]
  9.  
  10. #paramètres du fichier tfw
  11. #lecture du fichier
  12. with open(tfw) as f:
  13. data = [line.strip('\n').strip() for line in f]
  14. #paramètres
  15. taillepixelx = float(data[0])
  16. taillepixely = float(data[3])
  17. topcentregauchex = float(data[4])
  18. topcentregauchey = float(data[5])
  19.  
  20. #calculs
  21. xsupgauche = topcentregauchex - (taillepixelx * .5)
  22. ysupgauche = topcentregauchey - (taillepixely * .5)
  23. print xsupgauche,ysupgauche
  24. 162012.677881 108172.869385
  25. #taille du raster en m
  26. lx = rasterx * taillepixelx
  27. print "dimension en x: ", lx
  28. 8002.3097846
  29. ly = rastery * taillepixely
  30. print "dimension en y: ", ly
  31. -10002.8872308
  32. #calcul des coordonnées du point inférieur droit
  33. xinfdroit = xsupgauche + lx
  34. yinfdroit = ysupgauche + ly
  35. print xinfdroit,yinfdroit
  36. 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.

traitement

  1. from osgeo import gdal
  2. gdal.AllRegister() #tous les formats d'images de Gdal sont ainsi reconnus
  3. raster = gdal.Open("montif.tif")
  4. #dimension du raster
  5. rasterx = raster.RasterXSize
  6. rastery = raster.RasterYSize
  7. bands = raster.RasterCount
  8. print rasterx,rastery,bands
  9. 7996 9995 1
  10. #paramètres du géoréférencement
  11. geotransform = raster.GetGeoTransform()
  12. #paramètres
  13. xsupgauche = geotransform[0]
  14. ysupgauche = geotransform[3]
  15. taillepixelx = geotransform[1]
  16. taillepixely = geotransform[5]
  17. print xsupgauche,ysupgauche, taillepixelx,taillepixely
  18. 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.


Site officiel : Geotiff
Site officiel : Gdal osgeo
Site officiel : Python Imaging Library (PIL)
Site officiel : numpy
Site officiel : gdalinfo
Autres Liens : worldfile
Autres Liens : Python: projections cartographiques, définitions et transformations de coordonnées

Commentaires

Petite erreur

Bonjour,
Tout d'abord un grand merci pour vos articles qui sont vraiment une précieuse mine d'informations.

Il me semble qu'il y a une petite coquille dans la partie "coordonnées du centre du pixel" : il faut enlever en x et ajouter en y et dans le code il faut faire une soustraction:
xcentre = x1 - (ppx * .5)
ycentre = y1 - (ppy * .5) #puisque ppy est négatif

Par la suite le code est d'ailleurs correct sous "lecture et traitement du worldfile":
xsupgauche = topcentregauchex - (taillepixelx * .5)
ysupgauche = topcentregauchey - (taillepixely * .5)

Poster un nouveau commentaire

Le contenu de ce champ sera maintenu privé et ne sera pas affiché publiquement.