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

Les transformations affines (avec NumPy) ou la signification géométrique d'un fichier worldfile (.tfw, .jgw,...) pour le géoréférencement


Quelle est, in fine, la signification réelle d'un fichier worldfile. Le sujet a été traité, en partie, dans   « Python: géoréférencement et worldfiles (tfw, jgw,pngw,...) », mais d'une manière pratique et simplifiée, en particulier dans la définition des paramètres.

Les définitions données dans Wikipedia français ou anglais sont maintenant quasi les mêmes (ce qui n'a pas toujours été le cas) :

«  Ces paramètres représentent une transformation affine selon l’équation suivante : \begin{bmatrix} x\prime \\<br />
y\prime \end{bmatrix}<br />
= \begin{bmatrix} A & B \\<br />
D & E \end{bmatrix}<br />
\begin{bmatrix} x \\<br />
y \end{bmatrix} +<br />
\begin{bmatrix} C \\<br />
F \end{bmatrix} »

ou x' = Ax + Bx + C et y'= Dx + Ey + F

Les variables A, D, B, E, C et F constituent, dans l'ordre, les lignes du fichier worldfile.

Quelle est la signification de cette matrice et comment la calculer sans SIG ni module Python spécialisé (comme GDAL/OGR) ?  C'est ce que nous allons essayer de faire afin de mieux comprendre le processus mathématique de l'opération. On parle de fichier « World File » mais je préfère utiliser le terme fichier worldfile puisque c'est un fichier.

Les transformations géométriques dans le plan

Les transformations géométriques qui peuvent être appliquées à un objet sont :

Ces transformations peuvent être modélisées par des matrices de transformation. Elles permettent de résoudre aisément les problèmes.

La transformation d'un point (x,y) en point (x',y') par une des transformations (hormis la translation, voir plus bas) peut en effet s'écrire sous la forme :

soit x' = Ax + Dy et y'=Bx + Ey

La transformation va dépendre des valeurs données aux variables A, B, D et E. Ce peut être une rotation, un cisaillement, un changement d'échelle ou une combinaison de plusieurs de ces opérations.

Pour les détailler et les illustrer, je ne peux que vous diriger vers le magnifique tutoriel d'Olivier Tarasse, « Les matrices de transformation affine du plan (Matrix en Actionscript) » avec des animations en Actionscript (Flash) interactives. C'est clair et explicite pour tout un chacun.

copie d'écran de l'animation Flash de la page ressources.mediabox.fr/tutoriaux/flashplatform/math_physique/matrix/4-comprendre-transformation

Le problème est que cette matrice ne permettra jamais les translations :

quellle que soit la matrice, le point (0,0) restera le même

Cette difficulté peut être levée par l'introduction d'une troisième composante, neutre, aux vecteurs [x,y] et [x',y']. C'est ce qu'on appelle les coordonnées homogènes. Elles permettent alors de combiner les rotations, les cisaillements, les changements d'échelle et les translations dans une seule matrice, nommée matrice de transformation affine :

c'est-à-dire x' = Ax + By + C, y' = Dx + Ey + F et 1 (0x + 0y + 1)  

qui sont les formules définissant les 6 paramètres d'un fichier worldfile. Cette matrice peut donc être divisée en parties qui ont chacune un rôle précis :

(on peut aussi démontrer que le terme dans la case noire - 1 ici du fait des coordonnées homogènes - produit un changement d'échelle général, c'est-à-dire uniforme selon les deux axes)

Pour la résoudre, il est nécessaire de trouver les 6 paramètres. Il faut 6 équations à 6 inconnues et donc, au minimum,  3 paires de points.  C'est pourquoi on parle de transformation à 6 paramètres. En pratique elle réalise  2 translations, 2 rotations et 2 changements d'échelle, le tout en x et en y. Elle préserve le parallélisme, les rapports de surface, les rapports de longueurs sur une droite  et les coordonnées barycentriques.

En pratique

Lorsqu'un certain nombre de points d'appui sont connus, la transformation est donc:

Le rapport entre les 2 matrices de points donnera la matrice de transformation M comme résultat.

matrice des x', y'
matrice des x, y

Ce qui signifie en calcul matriciel:

(matrice des x',y')  x  inverse(matrice des x,y)

Or, seules les matrices carrées peuvent être inversées. La solution est, encore une fois, les coordonnées homogènes.

T=abcdef

Application en Python avec le module NumPy

En Python, pour traiter les matrices, les 2 modules les plus utilisés et performants sont NumPy et SciPy. Nous le ferons ici avec NumPy et avec son module numpy.matrix.

Nous allons traiter un fichier GEO_*.txt  accompagnant un fichier PCI Image. Il contient les points d'appui suivants qui correspondent à un fichier TIFF de taille l = 12609 et h = 8810 :

(00001,00001) (1616411.91,9173672.33)
(00001,08810) (1616433.22,9174420.34)
(12609,08810) (1617503.82,9174389.85)
(12609,00001) (1617482.52,9173641.84)

Ces points (les 4 premiers d'un fichier GEO_....txt)  correspondent aux 4 coins de l'image (voir www.forumsig.org/showthread.php et georezo.net/forum/viewtopic.php). La solution est donc la suivante :

traitement avec numpy

  1. import numpy as np
  2. m = 12609
  3. n = 8810
  4.  
  5. # matrice des pixels:
  6. fp = np.matrix([[1,m,m,1],[1,1,n,n]])
  7. # ajout des cordonnées homogènes
  8. newligne = [1,1,1,1]
  9. fp = np.vstack([fp,newligne])
  10. print fp
  11. matrix([[ 1, 12609, 12609, 1],
  12. [ 1, 1, 8810, 8810],
  13. [ 1, 1, 1, 1]])
  14. # matrice des coordonnées
  15. tp = np.matrix([[1616411.91,1617482.52,1617503.82,1616433.22],[9173672.33,9173641.84,9174389.85,9174420.34]])
  16.  
  17. # solution = fp x inverse(tp)
  18. M = tp * fp.I
  19. print M
  20. matrix([[ 8.49147367e-02, 2.41854921e-03, 1.61641183e+06],
  21. [ -2.41830584e-03, 8.49142922e-02, 9.17367225e+06]])
  22.  
  23. # paramètres
  24. A = M[:, 0][0]
  25. B = M[:, 1][0]
  26. C = M[:, 2][0]
  27. D = M[:, 0][1]
  28. E = M[:, 1][1]
  29. F = M[:, 2][1]
  30.  
  31. # dans l'ordre d'un worldfile
  32. >>> print A,D,B,E,C,F
  33. [[ 0.08491474]] [[-0.00241831]] [[ 0.00241855]] [[ 0.08491429]] [[ 1616411.82516671]] [[ 9173672.24750401]]

Ce résultat peut être contrôlé et il possible de calculer des points comme  (0, 8810)  :

contrôle

  1. # contrôle global
  2. print M * fp
  3. matrix([[ 1616411.9125,1617482.5175,1617503.8225,1616433.2175],
  4. [ 9173672.32999999,9173641.84,9174389.85,9174420.34]])
  5.  
  6. # contrôle par point
  7. xfp = np.matrix([[1],[8810],[1]])
  8. xtp = M * x0
  9. print xtp
  10. matrix([[ 1616433.2175],
  11. [ 9174420.34 ]])
  12.  
  13. # nouveau point
  14. xfp = np.matrix([[0],[8810],[1]])
  15. xtp = M * xfp
  16. print xtp
  17. matrix([[ 1616433.13258526],
  18. [ 9174420.34241831]])
 

La moyenne des erreurs résiduelles (RMS) ou la somme des erreurs quadratiques moyennes pourraient être directement calculées avec le module numpy.linalg mais ce n'est pas le sujet ici.

figure librement inspirée de celle de www.geo.utexas.edu/courses/371c/Lectures/2010_real_fall/Georeferencing_fall10.pdf

En pratique ces résultats sont équivalents à un ajustement par moindres carrés. Mais il y a d'autres ajustements possibles basés sur d'autres algorithmes (voir www.comp.nus.edu.sg/~cs4243/lecture/register.pdf ou le livre « The mathematics of GIS » de Wolfgang Kainz). En Python, Le module le plus complet est celui de Christoph Gohlke (www.lfd.uci.edu/~gohlke/code/transformations.py.html) , énorme fichier Python qui offre toutes les transformations possibles. Il est basé sur NumPy.

« A library for calculating 4x4 matrices for translating, rotating, reflecting,
scaling, shearing, projecting, orthogonalizing, and superimposing arrays of
3D homogeneous coordinates as well as for converting between rotation matrices,
Euler angles, and quaternions. Also includes an Arcball control object and
functions to decompose transformation matrices.»

avec transformations.py

  1. import transformations
  2. m = 12609
  3. n = 8810
  4. fp = [[1,m,m,1],[1,1,n,n]]
  5. tp = [[1616411.91,1617482.52,1617503.82,1616433.22],[9173672.33,9173641.84,9174389.85,9174420.34]]
  6. M = transformations.affine_matrix_from_points(fp, tp)
  7. print M
  8. array([[ 8.49147367e-02, 2.41854921e-03, 1.61641183e+06],
  9. [ -2.41830584e-03, 8.49142922e-02, 9.17367225e+06],
  10. [ 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])

D'autres algorithmes sont proposés par www.janeriksolem.net/2009/06/affine-transformations-and-warping.html (avec SciPy), code.google.com/p/oldmapsonline/source/browse/trunk/gcps2wld/python-no-dep/gcps2wld.py (avec NumPy) ou elonen.iki.fi/code/misc-notes/affine-fit/ (sans NumPy ni Scipy). Les résultats sont équivalents.

Contrôle avec le module osgeo.gdal

traitement avec osgeo.gdal

  1. from osgeo import gdal
  2. from numpy import *
  3. m = 12609
  4. n= 8810
  5. fp = [[1,m,m,1],[1,1,n,n]]
  6. tp = [[1616411.91,1617482.52,1617503.82,1616433.22],[9173672.33,9173641.84,9174389.85,9174420.34]]
  7. pix = zip(fp[0],fp[1])
  8. coor = zip(tp[0],tp[1])
  9. print pix
  10. ([(1, 1), (12609, 1), (12609, 8810), (1, 8810)]
  11. print coor
  12. [(1616411.9099999999, 9173672.3300000001), (1617482.52, 9173641.8399999999), (1617503.8200000001, 9174389.8499999996), (1616433.22, 9174420.3399999999)
  13. ])
  14.  
  15. gcps = []
  16. # calcule des paramètres avec gdal.GCP (points d'appui)
  17. for index, txt in enumerate(pix):
  18. gcps[index].GCPPixel = pix[index][0]
  19. gcps[index].GCPLine = 8810-int(pix[index][1])
  20. gcps[index].GCPX = coor[index][0]
  21. gcps[index].GCPY = coor[index][1]
  22. gcps.append(gdal.GCP())
  23.  
  24. geotransform = gdal.GCPsToGeoTransform( gcps )
  25. print geotransform
  26. (1616433.1325852636, 0.084914736675147789, -0.0024185492110491138, 9174420.3424183074, -0.0024183058378076717, -0.084914292200888988)
  27. >>> M = np.matrix([[geotransform[1],geotransform[4],geotransform[0]],[geotransform[2],geotransform[5],geotransform[3]]])
  28. print M
  29. matrix([[ 8.49147367e-02, -2.41830584e-03, 1.61643313e+06],
  30. [ -2.41854921e-03, -8.49142922e-02, 9.17442034e+06]])

Les paramètres de translation ne sont pas les mêmes, de même que le signe de B, mais il ne faut pas oublier que :

Il faut donc comparer ce qui est comparable et reprendre la fin du script « contrôle » (plus haut):

xfp = np.matrix([[0],[8810],[1]]) # point 0,0 de gdal
xtp = M * xfp
print xtp
matrix([[ 1616433.13258526],   # paramètre C
        [ 9174420.34241831]])      # paramètre F

et comparer avec les valeurs de translation de la transformation de gdal qui sont

        print M[:, 2][0]                  # paramètre C
        [[ 1616433.13]]
        print M[:, 2][1]                  # paramètre F
        [[ 9174420.34]]

Il faut ensuite ajuster l'ensemble pour ramener le point 0,0 à celui d'un worldfile (voir « Python: géoréférencement et worldfiles (tfw, jgw,pngw,...) »).

Conséquences

Les paramètres d'un fichier worldfile sont donc ceux d'une matrice de transformation affine :

et l'affirmation communément admise que :

  • la ligne 2 (D) correspond au paramètre de rotation autour de l'axe y, généralement 0 pour une translation simple ;
  • la ligne 3 (B) correspond paramètre de rotation autour de l'axe x, généralement 0 pour une translation simple;

n'est pas exacte, hormis dans les cas d'une translation simple.

D'autres transformations peuvent aussi être codées avec NumPy, SciPy comme

  • les tranformations à 4 paramètres ou transformations de Helmert  (très faciles à coder en Python). Elles réalisent 2 translations (en x et en y), 1 rotation et 1 changement d'échelle.  2 points d'appui sont donc au minimum, nécessaires. Elle ne permettent pas le cisaillement car elles préservent les angles et rapports de longueurs ;
  • les transformations polynomiales de premier ordre (= transformation affine), de second ordre ou de troisième ordre  (plus complexes à implémenter en Python avec de 7 à 22 points d'appui nécessaires). Elles permettent, entre autres, de transformer les lignes droites en courbes.
  • les transformations projectives à 8 paramètres (utilisée notamment par le module Pyproj)

Comme la figure précédente le montre, dans les cas de certaines transformations,  notamment les transformations polynomiales, l'image d'origine devra être modifiée. Ceci se fait alors par rééchantillonage de la matrice correspondante avec des algorithmes comme ceux du plus proche voisin, des interpolations linéaires et bilinéaires et/ou des convolutions cubiques. L'aptitude de NumPy et/ou SciPy à traiter les images sous forme matricielle permet de le faire sans problème. Ce sont des opérations très courantes dans le cas général des traitements d'image. Heureusement dans le cas des SIGs,  GDAL/OGR l'a fait pour nous.

Projections cartographiques ?

Vous avez peut-être remarqué que le terme projection géographique n'a jamais été cité. Pour tous les traitements d'image, un simple repère cartésien suffit. La preuve en est que nous avons utilisé le module NumPy qui est un module de traitement mathématique général, non spécifique aux SIGs.

Pour les sigistes, la particularité est que les points géoréférencés sont munis d'une projection. Celle-ci n'intervient donc que d'une manière indirecte. Un fichier worldfile seul ne permet pas de la retrouver, c'est pourquoi il vaut toujours mieux fournir le fichier de projection qui l'accompagne (.prj) ou l'indiquer dans le fichier GeoTiff sinon...

De par leurs caractéristiques, ces transformations servent évidemment aussi à effectuer les changements de projection (modules Pyproj, osgeo.gdal ou osgeo.ogr, par exemple) des rasters, mais aussi des vecteurs (shapefiles).

Conclusions

Les explications et les traitements mathématiques ont été ici simplifiés. On aurait pu discuter sur les propriétés géométriques des transformations (conservation des propriétés initiales d'une image ou non), aller plus loin dans les traitements mathématiques ou les traitements avec les modules NumPy et SciPy ou comparer les diverses transformations, mais ce n'était pas le but. Celui-ci était de montrer la véritable signification des paramètres d'un fichier worldfile et d'en examiner les conséquences, notamment pour les projections cartographiques des fichiers géoréférencés.

  • les paramètres d'un fichier worlfile sont le résultat de transformations affines classiques de traitement d'images où les projections n'interviennent que de manière indirecte ;
  • le système de projection d'un raster géoréférencé doit donc être précisé d'une autre manière ;
  • le fait de reprojeter un raster géoréférencé dans un SIG peut induire de nouvelles transformations géométriques (généralement effectuées de manière transparente par le logiciel)

Et pour finir, une figure tirée du monumental ouvrage de  D'Arcy Wentworth Thompson (1860-1948)    

  « On Growth and Form » (1917)

qui fut un pionnier dans la biologie mathématique. Il a montré que dans plusieurs classes d'organismes, notamment les poissons, les morphologies d'espèces apparentées pouvaient être générées par de simples transformations affines (sans translation, ni projections.... ) . Voir aussi demonstrations.wolfram.com/DArcyThompsonsAffineFishTransformations/ ou www.youtube.com/watch

 

Tous les traitements ont été effectués sur Mac OS X avec Python 2.6.1., NumPy version 2.0.0.dev et osgeo.gdal version 1.8.


Site officiel : Python: géoréférencement et worldfiles (tfw, jgw,pngw,..)
Site officiel : World File (Wikipedia fr.)
Site officiel : World File (Wikipedia angl.)
Site officiel : Les matrices de transformation affine du plan (Matrix en Actionscript)
Site officiel : Coordonnées homogènes
Site officiel : NumPy
Site officiel : Scipy
Site officiel : numpy.matrix
Site officiel : The mathematics of GIS (téléchargeable)
Autres Liens : D'Arcy Wentworth Thompson (Wikipedia)
Autres Liens : On growth and form (téléchargeable)


Creative Commons License
licence Creative Commons Paternité-Pas d'Utilisation Commerciale-Partage des Conditions Initiales à l'Identique Commerciale 2.0 France

Commentaires

Fichier WorldFile, transformation affine

Bonjour,
Bien d'avoir fait cet article.
A mon gout, il est un peu trop orienté dans le sens de l'utilisation par des logiciels particuliers, mais c'est déjà ça.
J'ai regardé l'article de Wikipedia français, il donne toujours comme exemple un cas particulier, en l'occurrence une rotation et une homothétie (+ translation).
Je n'ai pas retrouvé l'article anglais qui était complètement faux.
L'article écrit pas l'auteur du format (ESRI) est parfaitement clair et précis.
Je regrette un peu que vous n'ayez pas donné les avantages de cette transformation et surtout la raison pour laquelle on l'utilise, en particulier dans le domaine de la cartographie.
Pour mémoire la transformation Helmert n'existe pas, il s'agit d'une méthode de calcul, utilisable avec des outils simples (calculatrice avec ou sans mémoire), pour trouver la meilleure transformation (translation + homothétie + rotation) pour effectuer un changement de base à partir de plus de 2 points de calage, en général au moins 6, cette transformation conserve les angles et les rapports de distance. Elle donne de moins bons résultats que la transformation affine pour le calage de deux ou plusieurs images.

Poster un nouveau commentaire

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