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 : »

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 :

import numpy as np
m = 12609
n = 8810
 
# matrice des pixels:
fp = np.matrix([[1,m,m,1],[1,1,n,n]])
# ajout des cordonnées homogènes
newligne = [1,1,1,1]
fp  = np.vstack([fp,newligne])
print fp
matrix([[    1, 12609, 12609,     1],
        [    1,     1,  8810,  8810],
        [    1,     1,     1,     1]])
# matrice des coordonnées
tp = np.matrix([[1616411.91,1617482.52,1617503.82,1616433.22],[9173672.33,9173641.84,9174389.85,9174420.34]])
 
# solution = fp x inverse(tp)
M = tp * fp.I
print M
matrix([[  8.49147367e-02,   2.41854921e-03,   1.61641183e+06],
        [ -2.41830584e-03,   8.49142922e-02,   9.17367225e+06]])
 
# paramètres
A = M[:, 0][0]
B = M[:, 1][0]
C = M[:, 2][0]
D = M[:, 0][1]
E = M[:, 1][1]
F = M[:, 2][1]
 
# dans l'ordre d'un worldfile
>>> print A,D,B,E,C,F
[[ 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 global
print M * fp
matrix([[ 1616411.9125,1617482.5175,1617503.8225,1616433.2175],
        [ 9173672.32999999,9173641.84,9174389.85,9174420.34]])
 
# contrôle par point
xfp = np.matrix([[1],[8810],[1]])
xtp  = M * x0
print xtp
matrix([[ 1616433.2175],
        [ 9174420.34  ]])
 
# nouveau point
xfp = np.matrix([[0],[8810],[1]])
xtp = M * xfp
print xtp
matrix([[ 1616433.13258526],
        [ 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.»

import transformations
m = 12609
n = 8810
fp = [[1,m,m,1],[1,1,n,n]]
tp = [[1616411.91,1617482.52,1617503.82,1616433.22],[9173672.33,9173641.84,9174389.85,9174420.34]]
M = transformations.affine_matrix_from_points(fp, tp)
print M
array([[  8.49147367e-02,   2.41854921e-03,   1.61641183e+06],
       [ -2.41830584e-03,   8.49142922e-02,   9.17367225e+06],
       [  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

from osgeo import gdal
from numpy import *
m = 12609
n= 8810
fp = [[1,m,m,1],[1,1,n,n]]
tp = [[1616411.91,1617482.52,1617503.82,1616433.22],[9173672.33,9173641.84,9174389.85,9174420.34]]
pix = zip(fp[0],fp[1])
coor = zip(tp[0],tp[1])
print pix
([(1, 1), (12609, 1), (12609, 8810), (1, 8810)]
print coor
[(1616411.9099999999, 9173672.3300000001), (1617482.52, 9173641.8399999999), (1617503.8200000001, 9174389.8499999996), (1616433.22, 9174420.3399999999)
])
 
gcps = []
# calcule des paramètres avec gdal.GCP (points d'appui) 
for index, txt in enumerate(pix):
    gcps[index].GCPPixel = pix[index][0]
    gcps[index].GCPLine = 8810-int(pix[index][1])
    gcps[index].GCPX = coor[index][0]
    gcps[index].GCPY = coor[index][1]
    gcps.append(gdal.GCP())
 
geotransform = gdal.GCPsToGeoTransform( gcps )
print geotransform
(1616433.1325852636, 0.084914736675147789, -0.0024185492110491138, 9174420.3424183074, -0.0024183058378076717, -0.084914292200888988)
>>> M = np.matrix([[geotransform[1],geotransform[4],geotransform[0]],[geotransform[2],geotransform[5],geotransform[3]]])
print M
matrix([[  8.49147367e-02,  -2.41830584e-03,   1.61643313e+06],
        [ -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.

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