Google Git
Sign in
apache / sis-site / f218c657452b9ce9135e434595966157297c7259 / . / book / fr / referencing / CoordinateOperations.html
blob: 425745c87e8b1f79641a1a6dbf3d713b47f57abd [file] [log] [blame]
<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE html>
<!--
Licensed to the Apache Software Foundation (ASF) under one
or more contributor license agreements. See the NOTICE file
distributed with this work for additional information
regarding copyright ownership. The ASF licenses this file
to you under the Apache License, Version 2.0 (the
"License"); you may not use this file except in compliance
with the License. You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing,
software distributed under the License is distributed on an
"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
KIND, either express or implied. See the License for the
specific language governing permissions and limitations
under the License.
-->
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="fr">
<head>
<title>CoordinateOperations</title>
<meta charset="UTF-8"/>
<link rel="stylesheet" type="text/css" href="../../../content/book/book.css"/>
</head>
<body>
<!--
Content below this point is copied in "(…)/content/book/fr/developer-guide.html" file
by the 'org.apache.sis.internal.book.Assembler' class in 'sis-build-helper' module.
-->
<section>
<header>
<h2 id="CoordinateOperations">Opérations sur les coordonnées</h2>
</header>
<p>
Étant donné un système de référence des coordonnées (<abbr>CRS</abbr>) <em>source</em> selon lequel sont exprimés des coordonnées existantes
et un système de référence des coordonnées <em>destination</em> selon lequel les coordonnées sont désirées,
Apache <abbr>SIS</abbr> peut fournir une <em>opération sur les coordonnées</em> 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 <abbr>CRS</abbr>
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 <abbr>SIS</abbr> qu’on s’intéresse à une région plus petite
peut lui permettre de sélectionner une opération plus précise.
</p>
<div class="example"><p><b>Exemple:</b>
la base de données géodésiques <abbr>EPSG</abbr> (dans sa version 7.9) définit 77 opérations sur les coordonnées
allant du système géographique <cite>North American Datum 1927</cite> (EPSG:4267)
vers le système <cite>World Geodetic System 1984</cite> (EPSG:4326).
Il y a une opération valide seulement pour la transformation de coordonnées au Québec,
une autre opération valide pour la transformation de coordonnées au Texas mais à l’ouest de 100°W,
une autre opération pour le même état mais à l’est de 100°W, <i>etc</i>.
Si l’utilisateur ne spécifie pas la région géographique qui l’intéresse,
alors le comportement par défaut de Apache <abbr>SIS</abbr> est de sélectionner l’opération valide dans la plus grande région géographique.
Dans cet exemple, ce critère entraîne la sélection d’une opération valide pour le Canada, mais qui n’est pas valide pour les États-Unis.</p>
</div>
<h3 id="CRS.findOperation">Obtention d’une opération sur les coordonnées</h3>
<p>
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é <code>org.apache.sis.referencing.CRS</code>:
</p>
<pre><code>CoordinateOperation cop = <code class="SIS">CRS.findOperation</code>(sourceCRS, targetCRS, areaOfInterest);</code></pre>
<p>
Parmi les information fournies par l’objet <code>CoordinateOperation</code> obtenu, on note en particulier:
</p>
<ul>
<li>Le <cite>domaine de validité</cite>, soit comme une description textuelle telle que « Canada – onshore and offshore »
ou comme les coordonnées géographiques d’une boîte englobante.</li>
<li>La <cite>précision</cite>, qui peut être n’importe quoi entre 1 centimètre et quelques kilomètres.</li>
<li>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:
<ul class="verbose">
<li>
Les <strong>conversions</strong> 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.
</li><li>
Les <strong>transformations</strong> 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
<cite>Nouvelle Triangulation Française</cite> (<abbr>NTF</abbr>) vers le
<cite>Réseau Géodésique Français 1993</cite> (<abbr>RGF93</abbr>) entrent dans cette catégorie.
</li>
</ul>
</li>
</ul>
<p>
Lorsque l’opération sur les coordonnées est une instance de <code>Transformation</code>,
il est possible que l’instance choisie par <abbr>SIS</abbr> 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 <code>Conversion</code>.
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.
</p>
<h3 id="MathTransform">Exécution de opérations</h3>
<p>
L’objet <code>CoordinateOperation</code> introduit dans la section précédente fournit des informations de haut-niveau
(<abbr>CRS</abbr> source et destination, zone de validité, précision, paramètres de l’opération, <i>etc</i>).
Le travail mathématique réel est effectué par un objet séparé, obtenu par un appel à <code>CoordinateOperation.getMathTransform()</code>.
Contrairement aux instances de <code>CoordinateOperation</code>, les instances de <code>MathTransform</code> ne contiennent aucune méta-données.
Elles sont une sorte de boîte noire qui ignore tout des <abbr>CRS</abbr> source et destination
(en fait la même instance de <code>MathTransform</code> peut être utilisée pour différentes paires de <abbr>CRS</abbr> si le travail mathématique est le même).
En outre une instance de <code>MathTransform</code> peut être implémentée d’une manière très différente à ce que <code>CoordinateOperation</code> 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, <i>etc.</i>)
sont implémentées par <code>MathTransform</code> comme des <a href="#AffineTransform">transformations affines</a>
et concaténées pour des raisons d’efficacité, même si <code>CoordinateOperation</code> les affiche comme une chaîne d’opérations telles que la projection de Mercator.
La section « <a href="#CoordinateOperationSteps">chaîne d’opération conceptuelle versus réelle</a> » explique plus en détails les différences.
</p>
<p>
Le code Java suivant effectue une projection cartographique à partir de coordonnées géographiques selon le référentiel
<cite>World Geodetic System 1984</cite> (<abbr>WGS84</abbr>) vers des coordonnées selon le système <cite>WGS 84 / UTM zone 33N</cite>.
Afin de rendre l’exemple un peu plus simple, ce code utilise des constantes pré-définies dans la classe de commodité <code>CommonCRS</code>.
Mais des applications plus avancées voudront souvent utiliser des codes <abbr>EPSG</abbr> plutôt.
Notons que toutes les coordonnées géographiques dans ce code ont la latitude avant la longitude.
</p>
<pre><code>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 = <code class="SIS">CRS.findOperation</code>(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);
}
}</code></pre>
<h3 id="TransformDerivative">Dérivées partielles des opérations</h3>
<p>
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,
<a href="https://www.ogc.org/standards/ct">OGC 01-009</a>, aujourd’hui un peu oubliée mais pourtant encore utile.
Appelons <var>P</var> une projection cartographique qui convertit une latitude et longitude (<var>φ</var>, <var>λ</var>) en degrés
vers une coordonnée projetée (<var>x</var>, <var>y</var>) 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):
</p>
<div class="row-of-boxes">
<div style="min-width:350px; padding-right:60px">
<div class="caption">Équation</div>
<math xmlns="http://www.w3.org/1998/Math/MathML" display="block" alttext="MathML capable browser required">
<mi>P</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo>
<mo>=</mo>
<mfenced open="[" close="]">
<mtable>
<mtr><mtd><mi>x</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo></mtd></mtr>
<mtr><mtd><mi>y</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo></mtd></mtr>
</mtable>
</mfenced>
</math>
</div>
<div style="min-width:500px; padding-left:60px">
<div class="caption">Code Java</div>
<pre style="margin:0"><code>DirectPosition geographic = new DirectPosition2D(<var>φ</var>, <var>λ</var>);
DirectPosition projected = <var><b>P</b></var>.transform(geographic, null);
double <var>x</var> = projected.getOrdinate(0);
double <var>y</var> = projected.getOrdinate(1);</code></pre>
</div>
</div>
<p>La dérivée de la projection cartographique en ce même point peut se représenter par une matrice Jacobienne:</p>
<div class="row-of-boxes">
<div style="min-width:350px; padding-right:60px">
<div class="caption">Équation</div>
<math xmlns="http://www.w3.org/1998/Math/MathML" display="block" alttext="MathML capable browser required">
<msup><mi>P</mi><mo>′</mo></msup><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo>
<mo>=</mo>
<msub><mi>JAC</mi><mrow><mi>P</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo></mrow></msub>
<mo>=</mo>
<mfenced open="[" close="]">
<mtable>
<mtr>
<mtd><mfrac><mrow><mo>∂</mo><mi>x</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo></mrow><mrow><mo>∂</mo><mi>φ</mi></mrow></mfrac></mtd>
<mtd><mfrac><mrow><mo>∂</mo><mi>x</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo></mrow><mrow><mo>∂</mo><mi>λ</mi></mrow></mfrac></mtd>
</mtr>
<mtr>
<mtd><mfrac><mrow><mo>∂</mo><mi>y</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo></mrow><mrow><mo>∂</mo><mi>φ</mi></mrow></mfrac></mtd>
<mtd><mfrac><mrow><mo>∂</mo><mi>y</mi><mo>(</mo><mi>φ</mi><mo>,</mo><mi>λ</mi><mo>)</mo></mrow><mrow><mo>∂</mo><mi>λ</mi></mrow></mfrac></mtd>
</mtr>
</mtable>
</mfenced>
</math>
</div>
<div style="min-width:500px; padding-left:60px">
<div class="caption">Code Java</div>
<pre style="margin:0"><code>DirectPosition geographic = new DirectPosition2D(<var>φ</var>, <var>λ</var>);
Matrix jacobian = <var><b>P</b></var>.derivative(geographic);
double dx_dλ = jacobian.getElement(0,1);
double dy_dφ = jacobian.getElement(1,0);</code></pre>
</div>
</div>
<p>
Les équations suivantes dans cette section abrégeront
∂<var>x</var>(<var>φ</var>, <var>λ</var>) par ∂<var>x</var> ainsi que
∂<var>y</var>(<var>φ</var>, <var>λ</var>) par ∂<var>y</var>,
mais il faut garder à l’esprit que chacune de ces valeurs dépendent de la coordonnée (<var>φ</var>, <var>λ</var>) 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 ∂<var>φ</var> degrés de latitude
à partir de la position (<var>φ</var>, <var>λ</var>)
— c’est-à-dire si on se déplace à la position geographique (<var>φ</var> + ∂<var>φ</var>, <var>λ</var>) —
alors la coordonnée projetée subira un déplacement de (∂<var>x</var>, ∂<var>λ</var>) metres
— c’est-à-dire qu’elle deviendra (<var>x</var> + ∂<var>x</var>, <var>y</var> + ∂<var>λ</var>).
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 ∂<var>λ</var> 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, <var>P</var><sub>1</sub> et <var>P</var><sub>2</sub>,
pour mieux illustrer le fait que le résultat varie en chaque point.
Dans cette figure, les vecteurs <var>U</var> et <var>V</var> désignent respectivement
la première et deuxième colonne des matrices Jacobiennes.
</p>
<div class="row-of-boxes">
<div>
<img style="border: solid 1px" src="../../../content/book/images/Derivatives.png" alt="Exemple de dérivées d’une projection cartographique"/>
</div><div>
<p>où les vecteurs sont reliés à la matrice par:</p>
<math xmlns="http://www.w3.org/1998/Math/MathML" display="block" alttext="MathML capable browser required">
<mtable><mtr>
<mtd>
<mover><mi>U</mi><mo>→</mo></mover><mo>=</mo>
<mfenced open="[" close="]">
<mtable>
<mtr>
<mtd><mfrac><mrow><mo>∂</mo><mi>x</mi></mrow><mrow><mo>∂</mo><mi>φ</mi></mrow></mfrac></mtd>
</mtr>
<mtr>
<mtd><mfrac><mrow><mo>∂</mo><mi>y</mi></mrow><mrow><mo>∂</mo><mi>φ</mi></mrow></mfrac></mtd>
</mtr>
</mtable>
</mfenced>
</mtd>
<mtd><mtext>et</mtext></mtd>
<mtd>
<mover><mi>V</mi><mo>→</mo></mover><mo>=</mo>
<mfenced open="[" close="]">
<mtable>
<mtr>
<mtd><mfrac><mrow><mo>∂</mo><mi>x</mi></mrow><mrow><mo>∂</mo><mi>λ</mi></mrow></mfrac></mtd>
</mtr>
<mtr>
<mtd><mfrac><mrow><mo>∂</mo><mi>y</mi></mrow><mrow><mo>∂</mo><mi>λ</mi></mrow></mfrac></mtd>
</mtr>
</mtable>
</mfenced>
</mtd>
</mtr></mtable>
</math>
</div>
</div>
<p>
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à.
</p>
<h4 id="DerivativeAndEnvelope">Utilité des dérivées pour la reprojection d’enveloppes</h4>
<p>
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.
</p>
<div class="row-of-boxes">
<div>
<div class="caption">Enveloppe avant la projection</div>
<img style="border: solid 1px; margin-right: 15px" src="../../../content/book/images/GeographicArea.png" alt="Enveloppe dans un CRS géographique"/>
</div><div>
<div class="caption">Forme géométrique après la projection</div>
<img style="border: solid 1px; margin-left: 15px" src="../../../content/book/images/ConicArea.png" alt="Forme géométrique dans un CRS projeté"/>
</div>
</div>
<p>
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 <abbr>SIG</abbr> 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.
</p><p>
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 <i>f(<var>x</var>)</i> = <var>C</var>₀ + <var>C</var>₁<var>x</var> + <var>C</var>₂<var>x</var>² + <var>C</var>₃<var>x</var>³),
dont on peut calculer les coefficients <var>C</var>.
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 <var>C</var>.
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.
</p>
<div class="row-of-boxes">
<div>
<img src="../../../content/book/images/Approximation.png" alt="Approximation cubique d’un bord d’une enveloppe"/>
</div><div>
Légende:
<ul>
<li><b>En bleue:</b> la forme géométrique correspondant à la projection de l’enveloppe.
C’est la forme dont on souhaite avoir le rectangle englobant.</li>
<li><b>En rouge</b> (sous les hachures): L’approximation
<var>y</var> = <var>C</var>₀ + <var>C</var>₁λ + <var>C</var>₂λ² + <var>C</var>₃λ³.</li>
<li><b>En vert</b> (pointillés): La position λ<sub>m</sub> du minimum de l’approximation, trouvée en résolvant
0 = <var>C</var>₁ + 2<var>C</var>₂λ<sub>m</sub> + 3<var>C</var>₃λ<sub>m</sub>².
Il peut y avoir jusqu’à deux minimums pour une même courbe cubique.</li>
</ul>
</div>
</div>
<p>
Dans la pratique Apache <abbr>SIS</abbr> 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 <abbr>SIS</abbr> 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).
</p><p>
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 à <var>n</var> 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.
</p><p>
L’approche décrite dans cette section est implémentée dans Apache <abbr>SIS</abbr>
par la méthode statique <code>Envelopes.transform(CoordinateOperation, Envelope)</code>.
Une méthode <code>Envelopes.transform(MathTransform, Envelope)</code> existe aussi comme alternative,
mais cette dernière ne devrait être utilisée que si on ne connaît pas l’objet <code>CoordinateOperation</code> utilisé.
La raison est que les objets de type <code>MathTransform</code> ne contiennent pas d’information sur le système de coordonnées sous-jasent,
ce qui empêche la méthode <code>Envelopes.transform(…)</code> de savoir comment gérer les points aux pôles.
</p>
<h4 id="DerivativeAndRaster">Utilité des dérivées pour la reprojection d’images</h4>
<p>
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 <em>destination</em>, on obtient la coordonnées
du pixel correspondant dans l’image source en utilisant <em>l’inverse</em> 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.
</p>
<div style="display:flex; flex-direction:column; align-items:center">
<div style="display:flex; justify-content:space-between; width:564px">
<div class="caption">Image source</div>
<div class="caption">Image destination</div>
</div>
<img src="../../../content/book/images/RasterProjection.png" alt="Reprojection d’images"/>
</div>
<p>
Toutefois, calculer la projection inverse pour chacun des pixels peut être relativement lent.
Afin d’accélérer les calculs, on utilise parfois une <cite>grille d’interpolation</cite>
dans laquelle on a pré-calculé les <em>coordonnées</em> 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 <cite>transformations affines</cite>.
Cette approche est implémentée par exemple dans la bibliothèque <cite>Java Advanced Imaging</cite> avec l’objet <code>WarpGrid</code>.
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.
</p><p>
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:
</p>
<div class="row-of-boxes">
<div>
<img src="../../../content/book/images/WarpGrid.png" alt="Warp grid"/>
</div><div>
Légende:
<ul>
<li><b>Points bleus:</b> première itération (9 points).</li>
<li><b>Points verts:</b> seconde itération (25 points, dont 16 nouveaux).</li>
<li><b>Points rouges:</b> troisième itération (81 points, dont 56 nouveaux).</li>
</ul>
Si l’on continue…
<ul>
<li>Quatrième itération: 289 points, dont 208 nouveaux.</li>
<li>Cinquième itération: 1089 points, dont 800 nouveaux.</li>
<li>Sixième itération: 4225 points, dont 3136 nouveaux.</li>
<li>…</li>
</ul>
</div>
</div>
<p>
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 <strong>après</strong> 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 <em>somme</em> du nombre de nouveaux points de <em>toutes</em> les itérations précédentes.
</p><p>
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 <strong>avant</strong> 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).
</p>
<p style="text-align:center"><img style="border: solid 1px" src="../../../content/book/images/IntersectionOfDerivatives.png" alt="Intersection of derivatives"/></p>
<p>
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).
</p>
<h4 id="GetDerivative">Obtention de la dérivée en un point</h4>
<p>
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 <abbr>SIS</abbr>.
Les implémentations des projections cartographiques dans Apache <abbr>SIS</abbr> 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 à <abbr>SIS</abbr> de réutiliser un maximum de termes communs.
Exemple:</p>
<pre><code>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.<code class="SIS">transform</code>(sourcePoint, 0, targetPoint, 0, true);</code></pre>
<p>
Si seule la matrice Jacobienne est désirée (sans la projection du point), alors la méthode
<code>MathTransform.derivative(DirectPosition)</code> offre une alternative plus lisible.
</p><p>
Apache <abbr>SIS</abbr> 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, <i>etc.</i>
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:
</p>
<pre><code>@Override
public Matrix derivative(DirectPosition p) throws TransformException {
Matrix jac = inverse().derivative(transform(p));
return Matrices.inverse(jac);
}</code></pre>
</section>
</body>
</html>