Projections (2)

Projections inverses et transformations d'images

On peut avoir l'impression que la transformation qui mène de $(\lambda,\varphi)$ vers les coordonnées planaires $(x,y)$ sur notre carte, c'est la seule chose nous aurons besoin. Cependant, réflechissons comment effectuer une projection non-triviale, et afficher une image bitmap, et non pas des données vectorielles.

En principe on peut considérer que chaque pixel de l'image source constitue un élément primitif à afficher, et effectuer les transformations qui s'imposent. Ceci est très risqué, et donne des résultats pollués par des artefacts, à cause du fait que le bitmap (image - raster) est un objet discret; les coordonnées des pixels sont entières, avec la précision qui dépend de la taille du raster.

On peut avoir des pixels différents (voisins ou non) qui sont mappés sur les même pixel-cible, pixels qui vont "nulle part" (seront invisibles) ou des "trous" non-remplis. Ceci complique le rendu considérablement, en général. Plusieurs procédures d'interpolation peuvent aider à éliminer quelques difficultés, mais en général, ceci n'est pas une bonne approche, et ceci est une vérité générale qui conditionne les techniques du traitement d'images. La solution du problème est basée sur la transformation de la cible vers la source :

La première étape est la même que pour la transformation des données vectorielles. On établit les bornes de l'image-cible ; il faut savoir dans quelles limites se placent les coordonnées $(x,y)$ finales. On alloue le buffer pour l'image-cible.

On parcourt pixel par pixel l'image-cible. Pour chaque point dans la boucle – ...

... – on applique la transformation inverse par rapport à la projection qui transforme les données vectorielles : $(x,y) \to (\lambda,\varphi)$. On récupère le point de l'image-source ou on signale que ce pixel ne correspond pas à un point légal !.
On remplit le pixel de l'image-cible avec la couleur appropriée.

La transformation inverse parfois est plus facile, mais parfois considérablement plus difficile que la directe, et au lieu d'avoir une formula mathématique, il faut appliquer un algorithme itératif. Ici la partie conceptuellement triviale, mais techniquement indispensable : la transformation affine entre les coordonnées "géographiques" et les coordonnées "pixel", se fait deux fois.

Les bornes de l'image-cible correspondent aux limites de l'espace de visualisation, comme avec les données vectorielles. Ici rien ne change, sauf qu'il faut appliquer la transformation à l'envers. Donc :

$x = x_0 + IX \displaystyle{\frac{x_1 - x_0}{W}}\,, \qquad y = y_1 - IY \displaystyle{\frac{y_1 - y_0}{H}}\,.$
si on renverse les axes ! Mais pour les images géographiques la convention visuelle, naturelle est typiquement respectée : le nord se trouve en haut, donc la formule pour $y$ est plutôt similaire à celle de $x$.

Après avoir récupéré par la projection inverse les coordonnées géographiques $(\lambda,\varphi)$ du pixel-source, on dispose de valeurs en radians ($\lambda$ entre 0 et $2\pi$, ou entre $-\pi$ et $\pi$) ; éventuellement en degrés, mais l'image-source (la plate carrée ou autre) a ses dimensions en pixels, donc ici il faut appliquer la transformation affine "directe" non pas pour tracer, mais pour récupérer la couleur du pixel.

Pour les exemples montrés ici c'est tout. Cependant il faut vraiment faire attention à cause de la géométrie discrète. Si une zone d'image-source subit un agrandissement considérable, on pourra voir les effets de pixélisation, la création des "plages". Il faut alors obligatoirement appliquer les algorithmes d'interpolation (bi-linéaire, bicubique, ou autre). Faute de temps, nous n'aborderons pas ce sujet.

Exemples de projections

Voici un échantillon de quelques projections considérées traditionnelles, populaires, utiles, etc. (Au moins utiles pour la connaissance générale des outils de projections, de la terminologie et des besoins des cartographes).

Projection de Mercator

Ceci est la plus célèbre projection cylindrique, conforme, formalisée par Gerardus Mercator en 1569.

La projection de Mercator peut être obtenue de manière similaire à celle montrée. Les rayons projetés du centre de la sphère "mappent" les points de la surface sphérique sur le cylindre. Il est évident que la region des pôles est singulière, la distortion verticale des zones polaires est grande. La correspondance horizontale est bonne.

Si on utilisait les mathématiques de ce dessin, on aurait $y=\tan(\varphi)$, ce qui n'est pas exact. Voici les formules mathématiques concernées (il est utile d'apprendre leur dérivation, mais pas forcément mémoriser). Nous avons la liberté de choisir le méridien central de la carte, disons, que cela correspond à la latitude $\lambda_0$ . Si le résultat final de la projection est le tuple $$(x,y)$, il est évident que $x = \lambda - \lambda_0$, comme pour toute projection cylindrique (R=1). Avec la coordonnée verticale, c'est un peu plus complexe.

La formule $y = R\cdot\tan(\varphi)$ ne fait pas les loxodromies droites, et ne préserve pas des angles. Les vraies formules sont ci-dessous.

$y = \log\, \tan(\frac{1}{4} \pi + \frac{\varphi}{2}) = \frac{1}{2} \log \frac{\displaystyle 1+\sin(\varphi)}{\displaystyle 1-\sin(\varphi)}$
$\, = \sinh^{-1}(\tan(\varphi)) = \tanh^{-1}(\sin(\varphi)) = \log(\tan(\varphi) + \sec(\varphi))$

(Plusieurs variantes équivalentes peuvent être utilisées). On peut continuer l'exercice de traçage en implémentant cette transformation, ainsi que son inverse, afin de transformer l'image de la Terre en format plate carrée.

D'où viennent ces formules? Nous savons déjà que quand on s'éloigne de l'équateur, la distorsion (agrandissement) verticale s'accentue. Il faut donc que la carte subisse l'agrandissement vertical correspondant, et les "formules de Cauchy géographiques" nous disent comment.

varphi} \cdot \cos(\varphi)\), et à gauche nous avons une constante, 1, et $y$ évidemment ne dépend pas de $\lambda$, alors $y = \int{\frac{1}{\cos(\varphi)}d\varphi}$. Le calcul intégral nous donne $\log\, \tan(\varphi/2)$, modulo les constantes d'intégration. Pour $\varphi=0$ la valeur est $y=0$.

Quand l'angle $\varphi$ s'approche à 90 degrés, $y$ s'enfuit vers l'infini, même si moins rapidement que $\tan(\varphi)$, ce qui est évident, car l'échelle horizontale explose également. Ainsi la carte de Mercator est inutile dans les zones polaires, ce qui n'a pas dérangé le Flamand Gerhard Kremer vel Mercator, bien dans son XVIème siècle.

On peut aller facilement jusqu'à 85 degrés, ou un peu plus, et rendre la totalité de Groenland, mais cet exemple prouve que les projections cartographiques peuvent cacher plusieurs problèmes mathématiques et du rendu.

Les transformations Mercator inverses sont également simples, la projection est séparable.

$\phi = 2 \tan^{-1}(e^y) -1/2 \pi = \tan^{-1}(\sinh y)$ ,
$\lambda = x + \lambda_0$.

Nous verrons que les zones polaires peuvent être très bien rendu par l'UTM, le Mercator transversal, mais d'abord construisons quelques autres projections, connues de l'école, des livres de géographie, etc.

Si la distance horizontale pas loin des poles retrécit jusqu'à zéro, le modification verticale peut rester fini, donc l'"explosion" d'une carte conforme n'est aucune fatalite.

Le site de Carlos Furuti (Brésil) montre quelques autres, avec des images d'excellente qualité.

À gauche nous avons la carte conforme de Lagrange, peu fréquente, mais intéressante théoriquement. La transformation de la surface sphérique en disque circulaire est un bon point de départ pour d'autres transformations.

Messieurs Eisenlohr et August ont conçu d'autres cartes conformes.

Si vous disposez déjà d'une transformation et de coordonnées $(x,y)$, une fonction analytique complexe qui vous donne $x'+iy' = f(x+iy)$ produit également une transformation conforme.

On peut engendrer des transformation (approximativement) conformes sans aucune formule mathématique, par des processus algorithmiques, dont on parlera encore.

Projection sinusoïdale

On peut la considérer comme la plus simple parmi les projections "sérieuses". Les formules qui la définissent sont : $y=\phi\,;\quad x=\lambda\cos(\phi)$. Le Jacobien de la transformation vaut

$$ \displaystyle{\frac{\partial x}{\partial \lambda}\cdot \frac{\partial y}{\partial \phi} - \frac{\partial x}{\partial \phi} \cdot \frac{\partial y}{\partial \lambda}} = \cos(\phi) $$

ce qui démontre l'équivalence de cette transformation. Les distorsions angulaires polaires sont importantes, l'angle total 360° est rétréci à 115°. (Exercice. Comment trouver ce nombre?).

Les transformations inverses sont évidentes : $\phi=y\,; \quad \lambda = x/\cos(y)$. La gestion de la singularité aux pôles n'est pas difficile. De toute façon, cette expression $0/0$ dans le contexte du mapping des images peut être résolue de manière arbitraire (toute la ligne sur la plate carrée a la même couleur).

La projection est connue sous le nom de Sanson-Flamsteed, mais un de premiers à l'avoir utilisé, c'est Jehan Cossin de Dieppe, autour de 1570.

Cette projection est un cas particulier (équatorial) de la projection équivalente de Bonne, présentée plus tard dans le texte. Elle est souvent un simple point de départ pour la présentation des projections discontinues (découpées), qui seront également brièvement présentées un peu plus tard.

Projection de Mollweide

Une autre projection équivalente, souvent présente dans les livres de géographie (elle est considérée plus esthétique que la sinusoïdale), cartes scolaires, etc., a été proposée par Karl Mollweide en 1805. (Jacques Babinet l'a popularisée en France sous le nom "homolographique").

Il n'y a pas de formule close pour cette projection. Introduisons un angle auxiliaire $\theta$, défini par : $2\theta + \sin(2\theta) = \pi \sin(\varphi)$. Notez que ceci est une équation pour $\theta$, qui'il faut résoudre, non pas une définition directe. Seulement quand $\varphi=\pm \pi/2$ on voit que $\theta= \pm \pi/2$ également (le cas de zéro est également trivial). On peut utiliser plusieurs méthodes itératives pour résoudre cette équation, celle de Newton marche très bien. La solution de l'équation $f(x)=0$, si on connaît une approximation $x_0$, peut être obtenue par l'itération : $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$, jusqu'à la convergence. Ici :

$$ \theta_0 = \phi\,; \quad \theta_{n+1} = \theta_n -\displaystyle{\frac{2\theta_n + \sin(2\theta_n) - \pi \sin(\phi)} {2 + 2 \cos(2\theta_n)}}\,. $$

(mieux : $\theta_0 = \arcsin(2\phi/\pi)$ ; 5 - 6 itérations suffisent). Avec $\theta$, la projection est : $$ x = \frac{2\sqrt{2}}{\pi} \lambda \cos(\theta)\,; \quad y = \sqrt{2} \sin(\theta)\,. $$ Les paramètres numériques (échelle) viennent uniquement de l'équivalence de la projection, et de la volonté de Mollweide d'avoir une "jolie" ellipse, approximativement deux fois plus large que haute.

Pour la transformation des images, pixel par pixel, des itérations comme ça ralentisseraient le processus considérablement, mais ici la transformation inverse est – comme on le voit – donnée directement :

$$ \theta = \arcsin(y/\sqrt{2})\,; \quad \phi=\arcsin((2\theta+\sin(2\theta))/\pi)\,; \quad \lambda = \displaystyle{\frac{\pi x}{2\sqrt{2} \cos(\theta)}}\,. $$

Il y a 200 ans la connaissance de la géometrie et de méthodes numériques était assez bonne, mais imaginez le travail manuel nécessaire pour tracer une carte avec une précision suffisante ! Remarque générale. Les images et les formules utilisent un répère particulier, le méridian central de la carte est $\lambda=0$, et verticalement les cartes sont centrées sur l'équateur. Pour les projections (pseudo)cylindriques, le choix d'un autre méridien est trivial, il suffit de translater horizontalement/cycliquement les données ou les images. Verticalement ceci est plus délicat et complique les formules, mais parfois c'est une bonne option : on peut ainsi déplacer les singularités sur des zones moins intéressantes, et rendre mieux les zones polaires.

Exercice. Effectuer une rotation des données standard ("plate carrée") autour de l'axe $y$, perpendiculaire au méridien zéro sur l'équateur. Regardez l'applet sur la section précédente afin de voir la distorsion de la graticule. Voici un échantillon.

Mercator Transverse ; UTM

La Transverse Mercator (connue en Allemagne comme Gauss-Krüger), est une projection cylindrique conforme importante. Le cylindre de projection est "horizontal", tangent à un méridien, les distorsions augmentent horizontalement, et la projection globale, est inadaptée aux cartes mondiales.

L'UTM ("Universelle") développé par le Corps des Ingénieurs de l'armée des ÉtatsUnis, ajoute des précisions à Mercator. Dans plusieurs pays elle est officielle. UTM constitue un système de coordonnées $(E,N)$, utilisé pour le géo-réferencement, comme alternative à $(\lambda,\varphi)$.

Le schéma officiel découpe la surface en 60 zones verticales de 6 degrés (parfois modifiées), et considère que les zones polaires, au dessus de 85° restent dehors du système (même si on a éliminé les vraies singularités de projection).

Les formules de Mercator transversal sont basées sur le Mercator d'origine, adaptées. Nous présentons des formules simplifiées, uniquement pour référence. Les vraies, elliptiques, utilisent les excentricités et sont peu lisibles. Cherchez le Web... Ci-dessous le facteur de correction d'échelle $k_0$ est très proche de 1, et le facteur d'échelle global $a$ est quelconque.

$$ x = \frac{1}{2}k_0 a \log\left(\displaystyle{\frac{1+\sin\lambda \cos\phi} {1-\sin\lambda\cos\phi}}\right)\,; \quad y = k_0 a \arctan(\tan\phi \sec\lambda)\,. $$

Et les inverses :

$$ \lambda = \arctan\left(\sinh\displaystyle{\frac{x}{k_0 a}} \sec\displaystyle{\frac{y}{k_0 a}}\right)\,; \quad \phi = \arcsin\left(\mathrm{sech}\displaystyle{\frac{x}{k_0 a}} \sin\displaystyle{\frac{y}{k_0 a}}\right)\,. $$

La construction des coordonnées $(E,N)$ est basée sur le développement de cela. Les formules sont longues et dépassent le cadre de ce cours. Mais elles ont la précision meilleure qu'un centimètre !