L’utilité des projections cartographiques est discutée dans de nombreux livres et sites web.
Il s’agit de représenter une surface courbe (la Terre) sur une surface plane (une carte ou un écran d’ordinateur) en contrôlant les déformations: on peut préserver les angles ou les surfaces, mais pas les deux à la fois. Les propriétés géométriques à conserver dépendent de l’objet d’étude et du travail à effectuer. Par exemple les pays plutôt allongés dans le sens Est-Ouest utilisent souvent une projection de Lambert, alors que les pays plutôt allongés dans le sens Nord-Sud préfèrent une projection de Mercator Transverse.
Les bibliothèques de logiciels cartographiques offrent le nécessaire pour calculer les coordonnées d’un point géographique dans une projection au choix. Mais il existe une autre opération moins connue, qui consiste à calculer non pas la coordonnée projetée d’un point, mais plutôt la dérivée de la fonction de projection cartographique en ce point. Cette opération était définie dans une ancienne spécification du consortium Open Geospatial, OGC 01-009,
aujourd’hui un peu oubliée mais pourtant encore utile. Apache SIS est l’une des rares bibliothèques libres offrant cette opération.
Appelons P une projection cartographique qui convertit une longitude et latitude (λ,φ) en degrés
vers une coordonnée projetée (x,y) en mètres. Soit P (λ, φ) composée par 2 fonctions η(λ, φ)
pour la longitude et ξ(λ, φ) pour la latitude tel que : P (λ, φ) = (η(λ, φ), ξ(λ, φ)). Dans l’expression ci-dessous, nous représentons le résultat de la projection cartographique sous forme d’une matrice colonne (la raison sera plus claire bientôt) :
La dérivée de la projection cartographique en ce même point peut se représenter par la matrice Jacobienne définie tel que :
Le premier élément de la matrice (∂x/∂λ) nous indique à quel déplacement vers l’Est (x en mètres) correspond un déplacement de un degré de longitude (λ). De même, le dernier élément de la matrice (∂y/∂φ) nous indique à quel déplacement vers le Nord (y en mètres) correspond un déplacement de un degré de latitude (φ). Les autres éléments (∂x/∂φ et ∂y/∂λ) sont des termes croisés (par exemple à quel déplacement en mètres vers le Nord correspond un déplacement de un degré de longitude). Ces valeurs ne sont généralement valides qu’à la position géographique (λ,φ). Si on se déplace un peu, ces valeurs changent légèrement.
On peut se représenter visuellement cette matrice comme ci-dessous. Cette figure représente la dérivée en deux points, P1 et P2, pour mieux illustrer le fait que le résultat varie en chaque point. Dans cette figure, les vecteurs U et V désignent respectivement la première et deuxième colonne des matrices de dérivées.
où les vecteurs sont reliés à la matrice par:
Cette figure nous montre déjà une utilisation possible des dérivées: elles donnent la direction des parallèles et des méridiens à une position donnée dans une projection cartographique. Par extension, on peut aussi s’en servir pour déterminer si des axes sont interchangés, ou si la direction d’un axe est renversée. Mais l’intérêt des dérivées ne s’arrête pas là.
Projection d’enveloppes
Les systèmes d’information géographiques ont très fréquemment besoin de projeter une enveloppe. Mais l’approche naïve, qui consisterait à projeter chacun des 4 coins du rectangle, ne suffit pas. La figure ci-dessous montre une enveloppe avant le projection, et la forme géométrique que l’on obtiendrait si on projetait finement l’enveloppe (pas seulement les 4 coins). Cette forme géométrique est plus complexe qu’un simple rectangle à cause des courbures induites par la projection cartographique. Construire une enveloppe rectangulaire qui engloberait les 4 coins de cette forme géométrique ne suffit pas, car la surface en bas de la forme est plus basse que les 2 coins du bas. Cette surface serait donc en dehors du rectangle.
| Enveloppe avant la projection | Forme géométrique après la projection |
|---|---|
 |
 |
Une façon simple d’atténuer le problème est d’échantillonner un plus grand nombre de points sur chacun des bords de la forme géométrique. On trouve ainsi des bibliothèques de SIG qui vont par exemple échantillonner 40 points sur chaque bord, et construire un rectangle qui englobe tout ces points. Mais même avec 40 points, les échantillons les plus proches peuvent encore être légèrement à côté du point le plus bas de la figure. Une petite portion de la forme géométrique peut donc toujours se trouver en dehors du rectangle. Il est tentant de considérer cette légère erreur comme négligeable, mais ces quelques pixels manquants entraînent divers artefacts comme une apparence de quadrillage lors de l’affichage d’images tuilées, ou une “pointe de tarte” manquante lors de la projection d’images sur un pôle. Augmenter artificiellement d’un certain pourcentage la taille de l’enveloppe projetée peut éliminer ces artefacts dans certains cas. Mais un pourcentage trop élevé fera traiter plus de données que nécessaire (en particulier lorsque cela entraîne le chargement de nouvelles tuiles d’images), alors qu’un pourcentage trop faible laissera quelques artefacts.
Les dérivées des projections cartographiques permettent de résoudre ce problème d’une manière plus efficace que la force brute. La figure ci-dessous reprend la forme projetée en exagérant des déformations. L’approche consiste à calculer la projection cartographiques des 4 coins comme dans l’approche naïve, mais en récupérant aussi les dérivées de la projection de ces 4 coins. Entre deux coins et avec leurs dérivées, on peut faire passer une et une seule courbe cubique (de la forme f(x) = C₀ + C₁x + C₂x² + C₃x³), dont on peut calculer les coefficients C.
Cette approximation (représentée en rouge ci-dessous) ne correspond pas tout-à-fait à la courbe désirée (en bleue) mais s’en rapproche. Ce qui nous intéresse n’est pas vraiment les valeurs de l’approximation, mais plutôt la position de son minimum, en particulier la longitude λ où se trouve ce minimum dans notre exemple (ligne pointillée verte). L’avantage est que la position du minimum d’une courbe cubique est facile à calculer lorsque l’on connaît les valeurs de C. En supposant que la longitude du minimum de la courbe cubique est proche de la longitude du minimum de la courbe réelle, il suffit de calculer la projection cartographique d’un point à cette longitude plutôt que d’échantillonner 40 points sur le bord de l’enveloppe.
Légende:
- En bleue: la forme géométrique correspondant à la projection de l’enveloppe.
C’est la forme dont on souhaite avoir le rectangle englobant. - En rouge (sous les hachures): L’approximation y = C₀ + C₁λ + C₂λ² + C₃λ³.
- En vert (pointillés): La position λm du minimum de l’approximation, trouvée en résolvant
0 = C₁ + 2C₂λm + 3C₃λm².
Il peut y avoir jusqu’à deux minimums pour une même courbe cubique.
Dans la pratique la mise en oeuvre dans le projet MapFaces utilise 8 points, soit les 4 coins plus un point au centre de chaque bord du rectangle à projeter, afin de réduire le risque d’erreur qu’induirait une courbe trop tordue entre deux points. Selon nos tests, l’utilisation de ces seuls 8 points avec leurs dérivées comme décrit ci-haut donne un résultat plus précis que l’approche “force brute” utilisant un échantillonnage de 160 points sur les 4 bords du rectangle. La précision du code pourrait être encore améliorée en répétant le processus à partir du minimum trouvée (une ou deux itérations suffiraient peut-être). Une économie de 150 points n’est pas énorme vu les performances des ordinateurs d’aujourd’hui. Mais toute la discussion précédente utilisait une forme géométrique à deux dimensions en guise d’exemple, alors que l’algorithme est applicable dans un espace à n dimensions. Et de fait, l’implémentation fonctionne pour un nombre arbitraire de dimensions. Les économies apportées par cet algorithme par rapport à la force brute augmentent de manière exponentielle avec le nombre de dimensions.
Projection d’images
La projection cartographique d’une image s’effectue en préparant une image initialement vide qui contiendra le résultat de l’opération, puis à remplir cette image en itérant sur tous les pixels. Pour chaque pixel de l’image destination, on obtient la coordonnées du pixel correspondant dans l’image source en utilisant l’inverse de la projection cartographique que l’on souhaite appliquer. La position obtenue ne sera pas nécessairement au centre du pixel de l’image source, ce qui implique qu’une interpolation de la valeur (ou de la couleur dans l’image ci-dessous) peut être nécessaire.
| Image source | Image destination |
|---|---|
 |
|
Toutefois, calculer la projection inverse pour chacun des pixels peut être relativement lent.
Afin d’accélérer les calculs, on utilise parfois une grille d’interpolation dans laquelle on a pré-calculé les coordonnées de la projection inverse de seulement quelques points. Les coordonnées des autres points se calculent alors par des interpolations bilinéaires entre les points pré-calculés, calculs qui pourraient éventuellement tirer parti d’accélérations matérielles sous forme de transformations affines. Cette approche est implémentée par exemple dans la bibliothèque Java Advanced Imaging avec l’objet WarpGrid. Elle offre en outre l’avantage de permettre de réutiliser la grille autant de fois que l’on veut si on a plusieurs images de même
taille à projeter aux mêmes coordonnées géographiques.
Mais une difficulté de cette approche est de déterminer combien de points il faut pré-calculer pour que l’erreur (la différence entre une position interpolée et la position réelle) ne dépasse pas un certain seuil (par exemple ¼ de pixel). On peut procéder en commençant par une grille de taille 3×3, puis en augmentant le nombre de points de manière itérative:
Légende:
- Points bleus: première itération (9 points).
- Points verts: seconde itération (25 points, dont 16 nouveaux).
- Points rouges: troisième itération (81 points, dont 56 nouveaux).
Si l’on continue…
- Quatrième itération: 289 points, dont 208 nouveaux.
- Cinquième itération: 1089 points, dont 800 nouveaux.
- Sixième itération: 4225 points, dont 3136 nouveaux.
- …
L’itération s’arrête lorsque, après avoir calculé de nouveaux points, on a vérifié que la différence
entre les coordonnées projetées et les coordonnées interpolées de ces nouveaux points est inférieure au seuil qu’on s’est fixé. Malheureusement cette approche nous permet seulement de déterminer après avoir calculé de nouveaux points… que ce n’était pas la peine de les calculer. C’est un peu dommage vu que le nombre de nouveaux points requis par chaque itération est environ 3 fois la somme du nombre de nouveaux points de toutes les itérations précédentes.
Les dérivées des projections cartographiques nous permettent d’améliorer cette situation en estimant si c’est la peine d’effectuer une nouvelle itération avant de la faire.
L’idée de base est de vérifier si les dérivées de deux points voisins sont presque pareilles,
auquel cas on présumera que la transformation entre ces deux points est pratiquement linéaire. Pour quantifier “presque pareil”, on procède en calculant l’intersection entre les tangentes aux deux points (une information fournie par les dérivées), et en calculant la distance entre cette intersection et la droite qui relie les deux points (la ligne pointillée dans la figure ci-dessous).
Implémentation dans Apache SIS
Cette discussion n’aurait pas un grand intérêt si le coût du calcul des dérivées des projections cartographiques était élevé par rapport au coût de la projection des points. Mais lorsque l’on dérive analytiquement les équations des projections, on constate que les calculs des positions et de leurs dérivées ont souvent plusieurs termes en commun. En outre, le calcul des dérivées est simplifié lorsque le code Java effectuant les projections ne se concentre que sur le “noyau” non-linéaire, après s’être déchargé des parties linéaires en les déléguant aux transformations affines comme le fait SIS (ce sujet est brièvement mentionné dans le blog précédent).
Les implémentations des projections cartographiques dans Apache SIS tirent parti de ces propriétés en ne calculant les dérivées que si elles sont demandées, et en offrant une méthode qui permet de projeter un point et obtenir sa dérivée en une seule opération afin de permettre à SIS de réutiliser un maximum de termes communs.
Exemple:
AbstractMathTransform projection = …; // Une projection cartographique de Apache SIS.
double[] sourcePoint = {longitude, latitude}; // La coordonnée géographique que l’on veut projeter.
double[] targetPoint = new double[2]; // Là où on mémorisera le résultat de la projection.
Matrix derivative = projection.transform(sourcePoint, 0, targetPoint, 0, true);
Apache SIS est capable de combiner les dérivées des projections cartographiques de la même façon que pour les projections de coordonnées: concaténation d’une chaîne de transformations, inversion, opérer sur un sous-ensemble des dimensions, etc.
L’approche décrite dans la section “projection d’enveloppes” est implémentée par la méthode statique Envelopes.transform(CoordinateOperation, Envelope) de SIS. L’approche décrite dans la section “projection d’images” n’est pour l’instant disponible que dans le projet Geotk (la source du code que nous migrons vers SIS), mais sera proposée à terme au projet Apache SIS.
Commentaire sur “A quoi servent les dérivées des projections cartographiques”