Étant donné un système de référence des coordonnées (CRS) source selon lequel sont exprimés des coordonnées existantes et un système de référence des coordonnées destination selon lequel les coordonnées sont désirées, Apache SIS peut fournir une opération sur les coordonnées qui effectuera le travail de conversion ou de transformation. La recherche d’une opération peut utiliser un troisième argument, optionnel mais recommandé: la région géographique des données à transformer. Ce dernier argument est recommandé parce que les opérations sur les coordonnées sont souvent valides seulement dans une région géographique (typiquement un pays ou une province particulière), et plusieurs transformations peuvent exister pour la même paire de CRS source et destination mais avec des domaines de validité différents. Il peut aussi y avoir des différentes transformations qui sont différents compromis entre la précision et le domaine de validité, de sorte que spécifier à Apache SIS qu’on s’intéresse à une région plus petite peut lui permettre de sélectionner une opération plus précise.

Obtention d’une opération sur les coordonnées

La façon la plus facile d’obtenir une opération sur les coordonnées à partir des informations présentées ci-dessus est d’utiliser la classe de commodité org.apache.sis.referencing.CRS:

CoordinateOperation cop = CRS.findOperation(sourceCRS, targetCRS, areaOfInterest);

Parmi les information fournies par l’objet CoordinateOperation obtenu, on note en particulier:

• Le domaine de validité, soit comme une description textuelle telle que « Canada – onshore and offshore » ou comme les coordonnées géographiques d’une boîte englobante.
• La précision, qui peut être n’importe quoi entre 1 centimètre et quelques kilomètres.
• Le sous-type de l’opération sur les coordonnées. Parmi les sous-types possibles, deux offrent les mêmes fonctionnalités mais avec une différence conceptuelle notable:
  • Les conversions de coordonnées sont entièrement définies par une formule mathématique. Les conversions s’effectueraient avec une précision infinie s’il n’y avait pas les inévitables erreurs d’arrondissements inhérents aux calculs sur des nombres réels. Les projections cartographiques entrent dans cette catégorie.
  • Les transformations de coordonnées sont définies de manière empirique. Elles ont souvent des erreurs de quelques mètres qui ne sont pas dues à des limites de la précision des ordinateurs. Ces erreurs découlent du fait que la transformation utilisée n’est qu’une approximation d’une réalité plus complexe. Les changements de référentiels tel que le passage de la Nouvelle Triangulation Française (NTF) vers le Réseau Géodésique Français 1993 (RGF93) entrent dans cette catégorie.

Lorsque l’opération sur les coordonnées est une instance de Transformation, il est possible que l’instance choisie par SIS ne soit qu’une parmi plusieurs possibilités en fonction de la région d’intérêt. En outre, sa précision sera certainement moindre que la précision centimétrique que l’on peut attendre d’une Conversion. Vérifier la zone de validité ainsi que la précision déclarées dans les méta-données de la transformation prend alors une importance particulière.

Exécution de opérations

L’objet CoordinateOperation introduit dans la section précédente fournit des informations de haut-niveau (CRS source et destination, zone de validité, précision, paramètres de l’opération, etc). Le travail mathématique réel est effectué par un objet séparé, obtenu par un appel à CoordinateOperation.getMathTransform(). Contrairement aux instances de CoordinateOperation, les instances de MathTransform ne contiennent aucune méta-données. Elles sont une sorte de boîte noire qui ignore tout des CRS source et destination (en fait la même instance de MathTransform peut être utilisée pour différentes paires de CRS si le travail mathématique est le même). En outre une instance de MathTransform peut être implémentée d’une manière très différente à ce que CoordinateOperation dit. En particulier, plusieurs opérations conceptuellement différentes (par exemple rotations de la longitude, changements d’unités de mesure, conversions entre deux projections de Mercator qui utilisent le même référentiel, etc.) sont implémentées par MathTransform comme des transformations affines et concaténées pour des raisons d’efficacité, même si CoordinateOperation les affiche comme une chaîne d’opérations telles que la projection de Mercator. La section « chaîne d’opération conceptuelle versus réelle » explique plus en détails les différences.

Le code Java suivant effectue une projection cartographique à partir de coordonnées géographiques selon le référentiel World Geodetic System 1984 (WGS84) vers des coordonnées selon le système WGS 84 / UTM zone 33N. Afin de rendre l’exemple un peu plus simple, ce code utilise des constantes pré-définies dans la classe de commodité CommonCRS. Mais des applications plus avancées voudront souvent utiliser des codes EPSG plutôt. Notons que toutes les coordonnées géographiques dans ce code ont la latitude avant la longitude.

import org.opengis.geometry.DirectPosition;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.operation.TransformException;
import org.opengis.util.FactoryException;
import org.apache.sis.referencing.CRS;
import org.apache.sis.referencing.CommonCRS;
import org.apache.sis.geometry.DirectPosition2D;

public class MyApp {
    public static void main(String[] args) throws FactoryException, TransformException {
        CoordinateReferenceSystem sourceCRS = CommonCRS.WGS84.geographic();
        CoordinateReferenceSystem targetCRS = CommonCRS.WGS84.universal(40, 14);  // Obtient la zone valide pour 14°E.
        CoordinateOperation operation = CRS.findOperation(sourceCRS, targetCRS, null);

        // Les lignes précédentes sont coûteuses et ne devraient être exécutées qu’une seule fois avant
        // de transformer plusieurs points.  Dans cet exemple, l’opération que nous obtenons est valide
        // pour des coordonnées dans la région géographique allant de 12°E à 18°E (zone 33) et 0°N à 84°N.

        DirectPosition ptSrc = new DirectPosition2D(40, 14);           // 40°N 14°E
        DirectPosition ptDst = operation.getMathTransform().transform(ptSrc, null);

        System.out.println("Source: " + ptSrc);
        System.out.println("Target: " + ptDst);
    }
}

Dérivées partielles des opérations

La section précédente indiquait comment projeter les coordonnées d’un système de référence vers un autre. 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. Appelons P une projection cartographique qui convertit une latitude et longitude (φ, λ) en degrés vers une coordonnée projetée (x, y) en mètres. 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):

DirectPosition geographic = new DirectPosition2D(φ, λ);
DirectPosition projected = P.transform(geographic, null);
double x = projected.getOrdinate(0);
double y = projected.getOrdinate(1);

La dérivée de la projection cartographique en ce même point peut se représenter par une matrice Jacobienne:

DirectPosition geographic = new DirectPosition2D(φ, λ);
Matrix jacobian = P.derivative(geographic);
double dx_dλ = jacobian.getElement(0,1);
double dy_dφ = jacobian.getElement(1,0);

Les équations suivantes dans cette section abrégeront ∂x(φ, λ) par ∂x ainsi que ∂y(φ, λ) par ∂y, mais il faut garder à l’esprit que chacune de ces valeurs dépendent de la coordonnée (φ, λ) donnée au moment du calcul de la matrice Jacobienne. La première colonne de la matrice nous dit que si l’on effectue un petit déplacement de ∂φ degrés de latitude à partir de la position (φ, λ) — c’est-à-dire si on se déplace à la position geographique (φ + ∂φ, λ) — alors la coordonnée projetée subira un déplacement de (∂x, ∂λ) metres — c’est-à-dire qu’elle deviendra (x + ∂x, y + ∂λ). De même la dernière colonne de la matrice nous indique quel sera le déplacement que subira la coordonnée projetée si on effectue un petit déplacement de ∂λ degrés de longitude de la coordonnée géographique source. On peut se représenter visuellement ces déplacements comme ci-dessous. Cette figure représente la dérivée en deux points, P₁ et P₂, 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 Jacobiennes.

où les vecteurs sont reliés à la matrice par:

$$\begin{array}{ccc}\overrightarrow{U}=\left[\begin{array}{c}\frac{\partial x}{\partial \phi }\\ \frac{\partial y}{\partial \phi }\end{array}\right]& \text{et}& \overrightarrow{V}=\left[\begin{array}{c}\frac{\partial x}{\partial \lambda }\\ \frac{\partial y}{\partial \lambda }\end{array}\right]\end{array}$$

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à.

Utilité des dérivées pour la reprojection 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 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 Apache SIS 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 de SIS 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 de Apache SIS 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.

L’approche décrite dans cette section est implémentée dans Apache SIS par la méthode statique Envelopes.transform(CoordinateOperation, Envelope). Une méthode Envelopes.transform(MathTransform, Envelope) existe aussi comme alternative, mais cette dernière ne devrait être utilisée que si on ne connaît pas l’objet CoordinateOperation utilisé. La raison est que les objets de type MathTransform ne contiennent pas d’information sur le système de coordonnées sous-jasent, ce qui empêche la méthode Envelopes.transform(…) de savoir comment gérer les points aux pôles.

Utilité des dérivées pour la reprojection 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).

Dans l’approche sans dérivées, l’itération s’arrête lorsque la distance entre la ligne pointillée (positions interpolées) et la ligne rouge (positions projetées) est inférieure au seuil de tolérance, ce qui implique de calculer la position projetée. Dans l’approche avec dérivées, on remplace la position projetée par l’intersection des deux tangentes (carré bleu foncé). Si la courbe n’est pas trop tordue – ce qui ne devrait pas être le cas entre deux points suffisamment proches – la courbe réelle passera à quelque part entre la droite pointillée et l’intersection. On s’évite ainsi des projections cartographiques, en apparence une seule dans cette illustration, mais en fait beaucoup plus dans une grille de transformation d’image (3× la somme des itérations précédentes).

Obtention de la dérivée en un point

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 aux 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. 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);

Si seule la matrice Jacobienne est désirée (sans la projection du point), alors la méthode MathTransform.derivative(DirectPosition) offre une alternative plus lisible.

Apache SIS est capable 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. Les opérations inverses (des systèmes projetés vers géographiques) sont souvent beaucoup plus compliquées à implémenter que les opérations originales (des systèmes géographiques vers projetés), mais par chance la matrice Jacobienne d’une fonction inverse est simplement l’inverse de la matrice Jacobienne de la fonction originale. Une fonction inverse peut donc implémenter le calcul de sa dérivée comme suit:

@Override
public Matrix derivative(DirectPosition p) throws TransformException {
    Matrix jac = inverse().derivative(transform(p));
    return Matrices.inverse(jac);
}