Copied: sis/site/trunk/book/math/PixelToGeographicSameAxisDirections.html (from r1745833, sis/site/trunk/book/fr/referencing.html)
URL: http://svn.apache.org/viewvc/sis/site/trunk/book/math/PixelToGeographicSameAxisDirections.html?p2=sis/site/trunk/book/math/PixelToGeographicSameAxisDirections.html&p1=sis/site/trunk/book/fr/referencing.html&r1=1745833&r2=1746518&rev=1746518&view=diff
==============================================================================
--- sis/site/trunk/book/fr/referencing.html (original)
+++ sis/site/trunk/book/math/PixelToGeographicSameAxisDirections.html Thu Jun 2 04:57:06 2016
@@ -20,886 +20,46 @@
under the License.
-->
-<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="fr"
- xmlns:xi = "http://www.w3.org/2001/XInclude">
+<html xmlns="http://www.w3.org/1999/xhtml">
<head>
- <title>Systèmes de références spatiales par coordonnées</title>
+ <title>Pixel indices to geographic coordinates</title>
<link rel="stylesheet" type="text/css" href="../book.css"/>
</head>
<body>
- <header>
- <h1 id="Referencing">Systèmes de références spatiales</h1>
- </header>
- <p>
- Pour donner une position sur la Terre on peut utiliser des noms tels que celui dâune ville ou une adresse postale
- â on parle alors de références spatiales par <cite>identifiants</cite> â
- ou on peut donner des valeurs numériques valides dans un système de coordonnées donné
- â on parle alors de références spatiales par <cite>coordonnées</cite>.
- Chaque système implique des approximations telles que:
- </p>
- <ul>
- <li>Le choix de la forme géométrique (géoïde, ellipsoïde, <i>etc.</i>) utilisée comme approximation de la forme de la Terre.</li>
- <li>Le choix des propriétés géométriques (angles, distances, <i>etc.</i>) à préserver lors de la représentation dâune carte sur une surface plane.</li>
- <li>Les pertes de précision lorsque lâon doit transformer des positions exprimées selon un système vers des positions exprimées selon un autre système.</li>
- </ul>
- <p>
- En lâabsence dâindication contraire, la précision recherchée pour les coordonnées sur la Terre est de 1 centimètre.
- Mais la maîtrise de cette précision nécessite le respect de certaines conditions:
- </p>
- <ul>
- <li>Rester dans la zone de validité du système telle que donnée par <code>ReferenceSystem.getDomainOfValidity()</code>.</li>
- <li>Savoir que les mesures de distances dans une projection cartographique donnée ne sont vraies quâà certains endroits,
- appelés par exemple « parallèles standards ».</li>
- <li>Vérifier la précision des transformations de coordonnées, par exemple avec
- <code>CoordinateOperation.getCoordinateOperationAccuracy()</code>.</li>
- </ul>
- <p>
- Le module <code>sis-referencing</code> de Apache <abbr>SIS</abbr> fournit un ensemble de classes implémentant
- les différentes spécialisations de lâinterface <code>ReferenceSystem</code> ainsi que leurs composantes.
- Ces implémentations permettent de stocker une description des systèmes de références spatiales
- ainsi que leurs méta-données telles que la zone de validité.
- Toutefois ces objets nâeffectuent aucune opération sur les coordonnées.
- Ces opérations sont le travail dâune autre famille de classes, dont la racine est lâinterface <code>CoordinateOperation</code>.
- Ces classes seront discutées dans <a href="#CoordinateOperation">une autre section</a>,
- mais nous mentionnons ici deux spécialisations en rapport avec le sujet de la précision des coordonnées:
- </p>
- <ul>
- <li>
- <p>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.</p>
- <div class="example"><p><b>Exemple:</b> les projections cartographiques.</p></div>
- </li><li>
- <p>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.</p>
-
- <div class="example"><p><b>Exemple:</b> 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>),
- même lorsque la méthode de projection cartographique (Lambert conique conforme) ne change pas.</p></div>
- </li>
- </ul>
-
-
-
- <h2 id="CRS_Components">Composantes dâun système de références par coordonnées</h2>
- <p>
- Les systèmes de références spatiales par coordonnées fournissent les informations nécessaires pour faire
- correspondre des coordonnées numériques à des positions dans le monde réel. Dans Apache <abbr>SIS</abbr>,
- ils sont pratiquement tous représentés par des classes dont le nom se termine en <abbr>CRS</abbr>
- (lâabréviation de <i>Coordinate Reference System</i> en anglais). Ces objets contiennent:
- </p>
- <ul>
- <li>Un référentiel (<i>datum</i> en anglais),
- qui indique entre autres quel ellipsoïde utiliser comme approximation de la forme de la terre.</li>
- <li>Une description de chaque axe (nom, direction, unité de mesure, limites).</li>
- <li>Parfois une liste de paramètres permettant de convertir les coordonnées dâun autre système.
- Câest le cas notamment lorsquâon utilise des projections cartographiques.</li>
- </ul>
- <p>
- Ces systèmes sont décrits par la norme <abbr>ISO</abbr> 19111 (<i>Referencing by Coordinates</i>),
- qui remplace en grande partie une norme plus ancienne mais encore utilisée pour certains aspects,
- <abbr>OGC 01-009</abbr> (<i>Coordinate Transformation Services</i>).
- Ces normes sont complétées par deux autres standards définissant des formats dâéchanges:
- <abbr>ISO</abbr> 19136 et 19162 pour respectivement
- le <cite>Geographic Markup Language</cite> (<abbr>GML</abbr>) â un format <abbr>XML</abbr> précis mais verbeux â
- et le <cite>Well-Known Text</cite> (<abbr>WKT</abbr>) â un format texte plus facile à lire par les humains.
- </p>
-
- <h3 id="Ellipsoid">Géoïde et ellipsoïde</h3>
- <p>
- La surface topographique réelle étant difficile à représenter mathématiquement, elle nâest pas utilisée directement en cartographie.
- Une autre surface un peu plus facilement utilisable est le géoïde,
- une surface sur laquelle la force gravitationnelle a partout la même valeur (surface équipotentielle du champ de gravité terrestre).
- Cette surface est en tout point perpendiculaire à la direction indiquée par un fil à plomb (verticale du lieu).
- Le géoïde coïnciderait avec le niveau moyen des mers sâil nây avait ni vent ni courants marins permanents comme le Gulf Stream.
- </p><p>
- Tout en étant nettement plus lisse que la surface topographique,
- le géoïde présente des creux et des bosses liés à lâinégale distribution des masses de la Terre.
- Pour une utilisation mathématiquement plus aisée, le géoïde est donc approximé par un ellipsoïde.
- Cette « figure de la Terre » est représentée dans GeoAPI par lâinterface <code>Ellipsoid</code>,
- qui constitue un élément fondamental des systèmes de références de type <code>GeographicCRS</code> et <code>ProjectedCRS</code>.
- Plusieurs dizaines dâellipsoïdes sont couramment employés pour la définition de référentiels.
- Certains offrent une excellente approximation pour une région précise
- au détriment des régions pour lesquelles le référentiel nâa pas été conçu,
- et dâautres offrant un compromis pour lâensemble de la planète.
- </p>
- <div class="example"><p><b>Exemple:</b>
- au début du XX<sup>e</sup> siècle aux Ãtats-Unis, lâétat du Michigan utilisait pour ses cartes un ellipsoïde basé
- sur lâellipsoïde « Clarke 1866 » mais auquel la longueur des axes a été allongée de 800 pieds.
- Cette modification visait à tenir compte du niveau moyen de lâétat au dessus du niveau de la mer.</p>
- </div>
-
- <h3 id="GeodeticDatum">Référentiel géodésique</h3>
- <p>
- Pour définir un système géodésique dans un pays, lâétat met en place un ellipsoïde de référence
- qui épouse au mieux sur lâensemble du pays la forme locale du géoïde.
- Lâécart entre cet ellipsoïde de référence et les creux et les bosses du géoïde reste généralement inférieur à 100 mètres.
- Les paramètres qui permettent de lier un <code>Ellipsoid</code> à la surface de la Terre (par exemple la position de son centre)
- sont représentées par un objet de type <code>GeodeticDatum</code>, que lâon traduit en français par « référentiel géodésique ».
- Plusieurs <code>GeodeticDatum</code> peuvent utiliser le même <code>Ellipsoid</code>, mais centré ou orienté différemment.
- </p><p>
- Avant lâavènement des satellites, les mesures géodésiques se déroulaient exclusivement à la surface de la terre.
- En conséquence, deux îles ou continents qui ne sont pas à portée visuelle lâun de lâautre nâétaient pas rattachés géodésiquement entre eux.
- Ainsi les référentiels <cite>North American Datum 1983</cite> (<abbr>NAD83</abbr>) et
- <cite>European Datum 1950</cite> (<abbr>ED50</abbr>) sont indépendants lâun de lâautre:
- leurs ellipsoïdes de référence ont des centres distincts et des dimensions différentes.
- Une même coordonnée géographique correspondra à des positions différentes dans le monde réel
- selon que la coordonnée se réfère à lâun ou lâautre de ces référentiels.
- </p><p>
- Lâinvention du <abbr title="Global Positioning System">GPS</abbr> a précipité la création dâun système géodésique mondial,
- nommé <abbr title="World Geodetic System 1984">WGS84</abbr>.
- Lâellipsoïde de référence est alors unique et centré au centre de gravité de la terre.
- Les <abbr>GPS</abbr> donnent à tout moment la position absolue du récepteur rapportée à ce système géodésique.
- Mais <abbr>WGS84</abbr> étant un système mondial, il peut diverger significativement des systèmes locaux.
- Par exemple lâécart entre <abbr>WGS84</abbr> et le système européen <abbr>ED50</abbr> est de lâordre de 150 mètres,
- et lâécart moyen par rapport au système de lâîle de la Réunion 1947 est de 1,5 kilomètres.
- Il ne faut donc pas rapporter aveuglement des positions <abbr>GPS</abbr> sur une carte.
- Des correspondances avec les systèmes régionaux peuvent être nécessaires
- et sont représentées dans GeoAPI sous forme dâobjets de type <code>Transformation</code>
- (une classe dâopérations mentionnée dans lâ<a href="#Referencing">introduction de ce chapitre</a>).
- </p><p>
- Les généralisation de lâusage du système <abbr>WGS84</abbr> tend à réduire le besoin dâutiliser
- les objets <code>Transformation</code> pour les données récentes, mais ne lâélimine pas complètement.
- La Terre bouge sous lâeffet de la tectonique des plaques et de nouveaux systèmes sont définis chaque année pour en tenir compte.
- Par exemple bien que le référentiel <abbr>NAD83</abbr> a été défini à lâorigine comme pratiquement équivalent à <abbr>WGS84</abbr>,
- il existe maintenant (en 2016) un écart dâenviron 1.5 mètres entre ces deux systèmes.
- Le référentiel <cite>Japanese Geodetic Datum 2000</cite> était aussi défini comme pratiquement équivalent à <abbr>WGS84</abbr>,
- mais le nouveau référentiel <cite>Japanese Geodetic Datum 2011</cite> sâen écarte.
- Le référentiel <abbr>WGS84</abbr> lui-même, sensé correspondre à une définition à un instant donné,
- a subit des révisions dues notamment à lâamélioration de la précision des instruments.
- Ainsi il existe aujourdâhui au moins six versions de <abbr>WGS84</abbr>.
- En outre beaucoups de bordures ont été définies légalement dans des référentiels plus anciens, par exemple <abbr>NAD27</abbr> aux Ãtats-Unis.
- Mettre à jour dans un nouveau référentiel peut obliger à transformer des lignes droites ou des formes géométriques simples en des formes plus irrégulières
- si on ne veut pas que des parcelles de terrain changent de propriétaire.
- </p>
- <article>
- <header>
- <h1>Bibliothèques de type « early binding » versus « late binding »</h1>
- </header>
- <p>
- Le caractère universel du système <abbr>WGS84</abbr> rend tentante lâidée de lâutiliser comme système pivot,
- afin de simplifier lâimplémentation dâune bibliothèque de transformation de coordonnées.
- La transformation dâune coordonnée dâun référentiel <var>A</var> vers un référentiel <var>B</var>
- pourrait se faire en transformant dâabord de <var>A</var> vers <abbr>WGS84</abbr>, puis de <abbr>WGS84</abbr> vers <var>B</var>.
- Il suffirait ainsi de stocker dans chaque objet <code>GeodeticDatum</code> les informations nécessaires à la transformation vers <abbr>WGS84</abbr>.
- Cette approche était encouragée dans la version 1 du format <abbr>WKT</abbr>, qui définissait un élément <code>TOWGS84[â¦]</code> remplissant ce rôle.
- Cette approche est désignée par <abbr>EPSG</abbr> sous le nom de « early binding »
- car elle associe des informations sur la transformations de coordonnées très tôt dans la définition des objets géodésiques,
- souvent directement au moment de la construction dâun object <code>GeographicCRS</code>.
- Bien que <abbr>EPSG</abbr> reconnaisse que cette approche soit couramment employée, elle nâest pas recommandée pour plusieurs raisons:
- </p>
- <ul>
- <li>Il existe parfois plusieurs transformations allant dâun référentiel <var>A</var> vers <var>B</var>,
- chacune étant plus précise pour une région géographique donnée.</li>
- <li>Certaines opérations sont conçues spécifiquement pour transformer de <var>A</var> vers <var>B</var>
- et nâont pas la même précision quâaurait une autre transformation faisant un détour par <abbr>WGS84</abbr>.</li>
- <li><abbr>WGS84</abbr> lui-même subit parfois des révisions, ce qui en fait une cible mouvante (bien que très lentement)
- pour les bibliothèques de transformations de coordonnées.</li>
- <li>Il existe dâautres systèmes globaux qui pourraient servir de pivot, par exemple le <cite>Galileo Reference Frame</cite> (<abbr>GTRF</abbr>)
- mis en place par le concurrent européen du <abbr>GPS</abbr>.</li>
- </ul>
- <div class="example"><p><b>Exemple:</b>
- la base de données géodésiques <abbr>EPSG</abbr> définie une cinquantaine de transformations de <abbr>NAD27</abbr> vers <abbr>NAD83</abbr>.
- Dans une approche de type « early binding », le même système de référence « <abbr>NAD27</abbr> » représenté dans le format <abbr>WKT</abbr> 1
- aurait besoin dâêtre défini avec un élément <code>TOWGS84[-8, 160, 176]</code> pour des coordonnées aux Ãtats-Unis,
- ou avec un élément <code>TOWGS84[-10, 158, 187]</code> pour coordonnées aux Canada.
- Différents paramètres existent aussi pour dâautres régions telles que Cuba.
- Il nâest donc pas possible de représenter une telle diversité en associant un seul élément <code>TOWGS84[â¦]</code> à un système de référence des coordonnées.
- Même en restreignant lâusage dâun référenciel au domaine de validité de son élément <code>TOWGS84[â¦]</code>,
- ces transformations resteraient approximatives avec une précision de 10 mètres dans le cas des Ãtats-Unis.
- Des transformations plus précises existent sous la forme des grilles de changements de référentiel <abbr>NADCON</abbr>,
- mais ces grilles sont pour des transformations de <abbr>NAD27</abbr> vers <abbr>NAD83</abbr>
- (qui se déplacent ensemble sur la même plaque continentale) et non vers <abbr>WGS84</abbr> (qui se déplace indépendamment).
- Cette différence était souvent ignorée lorsque lâon considérait que <abbr>NAD83</abbr> et <abbr>WGS84</abbr>
- étaient pratiquement équivalents, mais cette supposition est aujourdâhui à prendre avec plus de précautions.
- </p></div>
- <p>
- <abbr>EPSG</abbr> recommande plutôt dâutiliser une approche dite « late binding »,
- selon laquelle les méthodes et paramètres nécessaires aux transformations de coordonnées sont définis pour des paires
- de référentiels « <var>A</var> vers <var>B</var> » (éventuellement complétées par leurs domaines de validité)
- plutôt quâassociés à des référentiels pris isolément.
- Apache <abbr>SIS</abbr> est une implémentation de type « late binding »,
- bien quâune réminiscence de lâapproche « early binding » existe toujours
- sous la forme de la propriété <code>DefaultGeodeticDatum.getBursaWolfParameters()</code>.
- <abbr>SIS</abbr> nâutilise cette dernière que comme solution de dernier recours
- sâil ne peut pas appliquer lâapproche « late binding » avec les systèmes de références donnés.
- </p>
- </article>
-
- <h3 id="CoordinateSystem">Systèmes de coordonnées</h3>
- <p style="color: red">TODO</p>
-
- <h3 id="GeographicCRS">Systèmes géographiques</h3>
- <p style="color: red">TODO</p>
-
- <h4 id="GeographicWKT">Format <i>Well-Known Text</i></h4>
- <p style="color: red">TODO</p>
-
- <h3 id="ProjectedCRS">Projections cartographiques</h3>
- <p>
- Les projections cartographiques représentent 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.
- </p>
- <p style="color: red">TODO</p>
-
- <h4 id="ProjectedWKT">Format <i>Well-Known Text</i></h4>
- <p style="color: red">TODO</p>
-
- <h3 id="CompoundCRS">Dimensions verticales et temporelles</h3>
- <p style="color: red">TODO</p>
-
- <h4 id="CompoundWKT">Format <i>Well-Known Text</i></h4>
- <p style="color: red">TODO</p>
-
-
-
- <h2 id="GetCRS">Obtention dâun système de référence spatial</h2>
- <p style="color: red">TODO</p>
-
- <h3 id="CRSAuthorityFactory">Systèmes prédéfinis par des autorités</h3>
- <p style="color: red">TODO</p>
-
- <h3 id="CRSParsing">Lecture dâune définition au format GML ou WKT</h3>
- <p style="color: red">TODO</p>
-
- <h3 id="CRSFactory">Construction programmatique explicite</h3>
- <p style="color: red">TODO</p>
-
- <h3 id="CRS_UserCode">Ajout de définitions</h3>
- <p style="color: red">TODO</p>
-
-
-
- <h2 id="CoordinateOperation">Opérations sur les coordonnées</h2>
- <p style="color: red">TODO</p>
-
- <h3 id="MathTransform">Exécution de opérations</h3>
- <p style="color: red">TODO</p>
-
- <h4 id="AffineTransform">Les transformations affines</h4>
- <p>
- Parmi les sortes dâopérations quâun <abbr>SIG</abbr> doit effectuer sur les coordonnées spatiales, il en est une à la fois simple et très fréquente.
- Ce sont les opérations linéaires, constituées uniquement dâune combinaison dâadditions et de certaines multiplications.
- Ces opérations nâeffectuent pas de projections cartographiques, plus complexes, mais couvrent de nombreux autres cas:
- </p>
- <ul>
- <li>Changer lâordre des axes, par exemple de (<var>latitude</var>, <var>longitude</var>) vers (<var>longitude</var>, <var>latitude</var>).</li>
- <li>Changer la direction des axes (par exemple lâaxe des <var>y</var> des images pointe souvent vers le bas).</li>
- <li>Changer de méridien dâorigine (par exemple de <cite>Paris</cite> vers <cite>Greenwich</cite>).</li>
- <li>Changer le nombre de dimensions (par exemple passer des coordonnées 3D vers 2D).</li>
- <li>Convertir des unités de mesures (par exemple convertir des pieds en mètres).</li>
- <li>Convertir les coordonnées pixels dâune image en coordonnées géographiques
- (par exemple la conversion exprimée dans les fichiers <code>.tfw</code> qui accompagnent certaines images <abbr>TIFF</abbr>).</li>
- <li>Prendre en charge une petite partie des projections cartographiques
- (par exemple les paramètres <cite>False Easting</cite>, <cite>False Northing</cite> et <cite>Scale factor</cite>).</li>
- <li>Appliquer des rotations, translations, échelles ou cisaillements (des transformations dites <cite>affines</cite>).</li>
- </ul>
- <p>
- Les opérations linéaires peuvent se combiner efficacement:
- peu importe le nombre dâopérations linéaires que lâon enchaîne, le résultat sera exprimable par une seule opération linéaire.
- Cette propriété est plus facilement visible lorsque les opérations linéaires sont exprimées sous forme de matrices:
- pour combiner les opérations, il suffit de multiplier les matrices.
- </p>
- <div class="example"><p><b>Example:</b>
- supposons que nous disposons dâune image dont les coordonnées des pixels sont représentées par (<var>i</var>,<var>j</var>).
- Supposons que la taille de chaque pixel correspond à un nombre fixe de degrées de longitude et de latitude
- dans un système géographique donné et quâil nây a pas de rotation.
- La conversion des coordonnées pixels (<var>i</var>,<var>j</var>) vers les coordonnées géographiques (<var>λ</var>,<var>Ï</var>)
- est alors linéaire et peut être représentée par la matrice suivante:</p>
-
- <table class="hidden"><tr><td>
- <xi:include href="../math/PixelToGeographic.html"/>
- </td><td style="vertical-align: middle; padding-left: 30px">
- où
- </td><td style="vertical-align: middle">
- <ul>
- <li><var>S</var> est un facteur dâéchelle (<cite>Scale</cite>) correspondant dans cet exemple à la taille des pixels en degrés.</li>
- <li><var>H</var> est un terme de cisaillement (<cite>Shear</cite>), habituellement zéro sauf si lâimage a une rotation.</li>
- <li><var>T</var> est une translation (<cite>Translation</cite>) correspondant dans cet exemple à la coordonnée géographique dâun coin de lâimage.</li>
- </ul>
- </td></tr></table>
- <p>
- Concentrons notre attention sur la matrice du milieu dans lâéquation ci-dessus.
- Si nous nâinterchangeons ni nâinversons la direction dâaucun axe, alors une conversion des coordonnées pixels vers les coordonnées géographiques
- pourrait sâexprimer par la matrice « conversion originale » ci-dessous.
- Mais si lâon veut en outre inverser la direction de lâaxe des <var>j</var> pour se conformer à la convention la plus courante appliquée aux images
- (« changement 1 ») et interchanger lâordre des axes pour exprimer la latitude avant la longitude (« changement 2 »),
- alors on peut exprimer ces modifications par des multiplications matricielles comme suit
- (lâordre dans laquelle les opérations sont effectuées sur les coordonnées se lit de droite à gauche):
- </p>
- <table class="hidden"><tr>
- <th>Changement 2</th><th></th>
- <th>Changement 1</th><th></th>
- <th>Conversion originale</th><th></th>
- <th>Conversion modifiée</th>
- </tr><tr>
- <td style="vertical-align: middle">
- <math xmlns="http://www.w3.org/1998/Math/MathML" display="block" alttext="MathML capable browser required">
- <mfenced open="[" close="]">
- <mtable>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>1</mn></mtd>
- <mtd><mn>0</mn></mtd>
- </mtr>
- <mtr>
- <mtd><mn>1</mn></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>0</mn></mtd>
- </mtr>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>1</mn></mtd>
- </mtr>
- </mtable>
- </mfenced>
- </math>
- </td>
- <td style="vertical-align: middle; padding-left: 15px; padding-right: 15px">Ã</td>
- <td style="vertical-align: middle">
- <math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
- <mfenced open="[" close="]">
- <mtable>
- <mtr>
- <mtd><mn>1</mn></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>0</mn></mtd>
- </mtr>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>-1</mn></mtd>
- <mtd><mo>(</mo><msub><mi>j</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>j</mi><mrow>min</mrow></msub><mo>)</mo></mtd>
- </mtr>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>1</mn></mtd>
- </mtr>
- </mtable>
- </mfenced>
- </math>
- </td>
- <td style="vertical-align: middle; padding-left: 15px; padding-right: 15px">Ã</td>
- <td style="vertical-align: middle">
- <math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
- <mfenced open="[" close="]">
- <mtable>
- <mtr>
- <mtd><mfrac>
- <mrow>
- <msub><mi>λ</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>λ</mi><mrow>min</mrow></msub>
- </mrow><mrow>
- <msub><mi>i</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>i</mi><mrow>min</mrow></msub>
- </mrow>
- </mfrac></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><msub><mi>λ</mi><mrow>min</mrow></msub></mtd>
- </mtr>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mfrac>
- <mrow>
- <msub><mi>Ï</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>Ï</mi><mrow>min</mrow></msub>
- </mrow><mrow>
- <msub><mi>j</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>j</mi><mrow>min</mrow></msub>
- </mrow>
- </mfrac></mtd>
- <mtd><msub><mi>Ï</mi><mrow>min</mrow></msub></mtd>
- </mtr>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>1</mn></mtd>
- </mtr>
- </mtable>
- </mfenced>
- </math>
- </td>
- <td style="vertical-align: middle; padding-left: 15px; padding-right: 15px">=</td>
- <td style="vertical-align: middle">
- <math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
- <mfenced open="[" close="]">
- <mtable>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mo>-</mo><mfrac>
- <mrow>
- <msub><mi>Ï</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>Ï</mi><mrow>min</mrow></msub>
- </mrow><mrow>
- <msub><mi>j</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>j</mi><mrow>min</mrow></msub>
- </mrow>
- </mfrac></mtd>
- <mtd><msub><mi>Ï</mi><mrow>max</mrow></msub></mtd>
- </mtr>
- <mtr>
- <mtd><mfrac>
- <mrow>
- <msub><mi>λ</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>λ</mi><mrow>min</mrow></msub>
- </mrow><mrow>
- <msub><mi>i</mi><mrow>max</mrow></msub><mo>-</mo>
- <msub><mi>i</mi><mrow>min</mrow></msub>
- </mrow>
- </mfrac></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><msub><mi>λ</mi><mrow>min</mrow></msub></mtd>
- </mtr>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>1</mn></mtd>
- </mtr>
- </mtable>
- </mfenced>
- </math>
- </td>
- </tr></table>
- <p>
- Lâidée principale est quâil nây a pas besoin dâécrire un code dédié à lâinversion des axes.
- Cette opération, et bien dâautres, est prise en compte naturellement par lâalgèbre matricielle.
- On y gagne en généricité du code et en performance.
- </p>
- </div>
-
- <p style="color: red">TODO</p>
-
- <article>
- <header>
- <h1>Particularités dâune bibliothèque de calculs matriciels pour un <abbr>SIG</abbr></h1>
- </header>
- <p>
- Les <abbr>SIG</abbr> font un usage intensif de matrices afin dâafficher leurs cartes ou transformer des coordonnées.
- On pourrait croire que le marché est suffisamment bien pourvu en excellentes bibliothèques de calculs matriciels, libres ou commerciales.
- Pourtant, les <abbr>SIG</abbr> ont des besoins spécifiques qui divergent un peu des objectifs de plusieurs bibliothèques existantes.
- Des manipulations de matrices comme lâexemple précédent interviennent dans quasiment toutes les opérations
- appliquées par Apache <abbr>SIS</abbr> sur des coordonnées.
- Mais lâanalyse de ces opérations révèle quelques patterns:
- </p>
- <ul>
- <li><p>Ces matrices sont presque toujours de petites tailles, dépassant rarement 5 lignes par 5 colonnes.</p></li>
- <li><p>Les opérations matricielles « lourdes » (multiplications ou inversions de matrices) ne surviennent pas dans des endroits où la performance est importante.
- Dans la quasi-totalité des cas, elles ne sont effectuées quâune fois pour toute, à la lecture dâun fichier,
- ou lors des étapes de préparation avant de convertir des coordonnées.
- Elles ne surviennent quasiment jamais dans la boucle convertissant chacune des coordonnées.</p></li>
- <li><p>Dans une succession de multiplications et dâinversions de matrices, les erreurs dâarrondissement sâaccumulent et grandissent rapidement
- au point de se confondre avec certaines opérations légitimes, notamment les changements de référentiel.
- Ces dernières sâexpriment souvent par un changement de la taille, position et orientation de lâellipsoïde
- choisi comme approximation de la forme de la Terre. Les changements de la taille sâexpriment en parties par million et
- les rotations en arc-secondes. Retranscrites dans une matrice, ces valeurs sont donc assez petites.</p></li>
- <li><p>Il arrive fréquemment que des matrices sâannulent en tout ou en partie,
- câest-à -dire que leurs multiplications ramènent des facteurs dâéchelles à 1 et des translations à 0.
- Toutefois les erreurs dâarrondissements font que les valeurs obtenues sont rarement exactes,
- mais plutôt des valeurs sâen rapprochant telles que 0,9999â¦97 à la place de 1.
- Malheureusement, les erreurs dâarrondissement sont parfois telles quâil est difficile de savoir
- si certains coefficients de la matrices sont des artefacts ou proviennent dâun réel changement de référentiel.</p></li>
- </ul>
- <p>
- Ces points font que, pour un <abbr>SIG</abbr>, la précision dâune bibliothèque de calculs matriciels
- est plus importante que la performance. Paradoxalement, un bon moyen de gagner en performance est justement dâinvestir davantage de temps de CPU
- pour effectuer des opérations matricielles plus précises, car on augmente ainsi les chances de détecter correctement quelles opérations sâannulent.
- Lâeffort investit dans cette détection permet de sauver du temps là où ça compte: quand viendra le moment de boucler sur des millions de coordonnées à transformer.
- </p><p>
- Mais les bibliothèques dédiées aux calculs matriciels sont souvent conçues pour opérer de manière très performante
- sur des matrices de grandes tailles, ayant par exemple des milliers de lignes et colonnes.
- Elles sont ainsi conçues pour être capable de résoudre efficacement des systèmes dâéquations linéaires comportant des centaines dâinconnues.
- Les problèmes quâelles résolvent sont certes difficiles, mais assez différents de ceux qui intéressent Apache <abbr>SIS</abbr>.
- Pour cette raison, et aussi à cause dâun autre besoin spécifique détaillé dans les paragraphes suivants,
- Apache <abbr>SIS</abbr> utilise ses propres fonctions de calculs matriciels.
- Ces fonctions tentent de résoudre le problème de précision en utilisant lâarithmétique « double-double »
- (une technique permettant de simuler une précision dâenviron 120 bits)
- au prix de la performance dans une partie du code où elle nâest pas jugée critique.
- </p>
- <h2>Que faire des matrices qui ne sont pas carrées (et pourquoi)</h2>
- <p>
- Apache <abbr>SIS</abbr> a très souvent besoin dâinverser des matrices,
- afin dâobtenir une conversion de coordonnées qui fasse le contraire de la conversion originale.
- Mais on nâinverse habituellement que des matrices carrées.
- Or, Apache <abbr>SIS</abbr> a besoin dâeffectuer des inversions de matrices non-carrées.
- Selon que lâon ait plus de lignes ou plus de colonnes:
- </p>
- <ul>
- <li>Pour <abbr>SIS</abbr>, une matrice non-carrée est une conversion qui ajoute ou supprime une dimension aux coordonnées.</li>
- <li>Pour les bibliothèques dâalgèbre linéaire, une matrice non-carrée est un système dâéquations sous-déterminé ou surdéterminé.</li>
- </ul>
- <p>
- Pour mieux comprendre les difficultés que causerait une transposition trop directe des bibliothèques dâalgèbre linéaire aux <abbr>SIG</abbr>,
- imaginons une conversion qui projetterait les points dâun espace 3D vers une surface 2D:
- </p>
- <table class="hidden">
- <tr>
- <td>(λâ, Ïâ, <var>h</var>) â (λâ, Ïâ)</td>
- <td style="padding-left: 30px">où</td>
- <td><ul style="margin-top: 0">
- <li>λ est la longitude.</li>
- <li>Ï est la latitude.</li>
- <li>(λâ, Ïâ) nâégale pas forcement (λâ, Ïâ) si la hauteur <var>h</var> nâest pas perpendiculaire à la surface.</li>
- </ul></td>
- </tr>
- </table>
- <p>
- Pour des bibliothèques dâalgèbre linéaire, la matrice représentant cette conversion serait un système dâéquations sous-déterminé, et donc insoluble.
- Câest-à -dire quâon ne peut pas inverser cette conversion pour obtenir (λâ, Ïâ) â (λâ, Ïâ, <var>h</var>) puisquâon ne sait pas quelle valeur donner à <var>h</var>,
- ce qui implique quâon ne peut pas trouver (λâ, Ïâ) non-plus car ces valeurs dépendent peut-être de <var>h</var>.
- Toutefois, dans le cas des <abbr>SIG</abbr>, lâaxe des <var>h</var> est très souvent perpendiculaire à la surface sur laquelle sont exprimées les coordonnées (<var>λ</var>,<var>Ï</var>).
- Cette perpendicularité rend λâ et Ïâ indépendants de <var>h</var>. Dans ce cas particulier, et ce cas seulement, on peut encore sauver les meubles.
- </p><p>
- Apache <abbr>SIS</abbr> procède en vérifiant si les coordonnées <var>h</var> sont indépendantes des coordonnées λ et Ï.
- Nous reconnaissons ce cas en vérifiant quels coefficients de la matrice ont la valeur zéro.
- Si <abbr>SIS</abbr> arrive à identifier des dimensions indépendantes,
- il peut les exclure temporairement de manière à inverser sans ambiguïté la conversion dans les dimensions restantes.
- Sâil ne trouve pas de dimension indépendante, alors une exception est levée.
- </p><p>
- Si une inversion a été possible, alors il reste à décider du sort des dimensions que <abbr>SIS</abbr> avait temporairement exclues.
- Dans notre exemple, <abbr>SIS</abbr> assignera la valeur <code>NaN</code> (<cite>Not-a-Number</cite>) aux valeurs de <var>h</var> dans la conversion (λâ, Ïâ) â (λâ, Ïâ, <var>h</var>).
- Là encore, le choix du coefficient à mettre à <code>NaN</code> dans la matrice est basé sur la présomption quâelle représente une conversion de coordonnées.
- </p><p>
- Le traitement particulier fait par <abbr>SIS</abbr> permet donc dâinverser des matrices que lâon rencontre couramment dans les <abbr>SIG</abbr>,
- même si en principe le système est sous-déterminé.
- Dans notre exemple la coordonnée <var>h</var> reste inconnue â nous ne faisons pas surgir de lâinformation du néant â mais au moins les coordonnées (<var>λ</var>,<var>Ï</var>) ont pu être récupérées.
- </p><p>
- Le problème inverse, celui des systèmes surdéterminés, est plus subtil.
- Une approche classique des bibliothèques dâalgèbre linéaire est de résoudre les systèmes surdéterminés par la méthode des moindres carrées.
- Transposée à notre exemple, cette approche proposerait une conversion (λâ, Ïâ, <var>h</var>) â (λâ, Ïâ)
- qui semble le meilleur compromis pour diverses valeurs de λâ, Ïâ et <var>h</var>, tout en nâétant (sauf cas particuliers) une solution exacte pour personne.
- De plus, les éventuelles combinaisons linéaires entre ces trois variables sont délicates compte tenu de lâhétérogénéité des unités de mesures,
- où les <var>h</var> sont en mètres et (<var>λ</var>,<var>Ï</var>) en degrés.
- Apache <abbr>SIS</abbr> procède plutôt comme pour les systèmes sous-déterminés: en exigeant que certaines dimensions soient indépendantes des autres,
- faute de quoi la matrice sera considérée non-inversible.
- Dans le cas des systèmes surdéterminés <abbr>SIS</abbr> refusera donc dâeffectuer certaines opérations que les bibliothèques dâalgèbre linéaire auraient faite,
- mais garantira que les conversions obtenues seront exactes (aux erreurs dâarrondissement prêts).
- </p>
- <p>
- En résumé, les besoins qui ont amené Apache <abbr>SIS</abbr> à fournir ses propres fonctions de calculs matriciels sont:
- </p>
- <ul>
- <li>Structure légère pour les petites matrices, particulièrement celles de taille 3Ã3.</li>
- <li>Précision accrue avec lâarithmétique « double-double », quitte à sacrifier un peu de performance dans des endroits où elle nâest pas critique.</li>
- <li>Traitement particulier de lâinversion des matrices non-carrées pour des conversions de coordonnées.</li>
- </ul>
- </article>
-
- <h3 id="TransformDerivative">Dérivées partielles des opérations</h3>
- <p>
- La section précédente indiquait comment 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 <em>coordonnées projetée</em> dâun point,
- mais plutôt la <em>dérivée de la fonction de projection cartographique</em> en ce point.
- Cette opération était définie dans une ancienne spécification du consortium Open Geospatial,
- <a href="http://www.opengeospatial.org/standards/ct">OGC 01-009</a>, aujourdâhui un peu oubliée mais pourtant encore utile.
- </p>
-
- <p>
- Appelons <var>P</var> une projection cartographique qui convertit une longitude et latitude (<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>
-
- <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>
-
- <p>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:</p>
-
- <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>
+ <math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
<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>
+ <mtd><mfrac>
+ <mrow>
+ <msub><mi>λ</mi><mrow>max</mrow></msub><mo>-</mo>
+ <msub><mi>λ</mi><mrow>min</mrow></msub>
+ </mrow><mrow>
+ <msub><mi>N</mi><mrow><mi>x</mi></mrow></msub>
+ </mrow>
+ </mfrac></mtd>
+ <mtd><mn>0</mn></mtd>
+ <mtd><msub><mi>λ</mi><mrow>min</mrow></msub></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>
+ <mtd><mn>0</mn></mtd>
+ <mtd><mfrac>
+ <mrow>
+ <msub><mi>Ï</mi><mrow>max</mrow></msub><mo>-</mo>
+ <msub><mi>Ï</mi><mrow>min</mrow></msub>
+ </mrow><mrow>
+ <msub><mi>N</mi><mrow><mi>y</mi></mrow></msub>
+ </mrow>
+ </mfrac></mtd>
+ <mtd><msub><mi>Ï</mi><mrow>min</mrow></msub></mtd>
+ </mtr>
+ <mtr>
+ <mtd><mn>0</mn></mtd>
+ <mtd><mn>0</mn></mtd>
+ <mtd><mn>1</mn></mtd>
</mtr>
</mtable>
</mfenced>
</math>
-
- <p>
- Dans la suite de ce texte nous abrégerons â<var>x</var>(<var>λ</var>,<var>Ï</var>) par â<var>x</var> et de même pour â<var>y</var>,
- mais il faut garder à lâesprit que chacune de ces valeurs dépendent de la coordonnée (<var>λ</var>,<var>Ï</var>) originale.
- Le premier élément de la matrice (â<var>x</var>/â<var>λ</var>) nous indique à quel déplacement vers lâEst
- (<var>x</var> en mètres) correspond un déplacement de un degré de longitude (<var>λ</var>).
- De même, le dernier élément de la matrice (â<var>y</var>/â<var>Ï</var>) nous indique à quel déplacement vers le Nord
- (<var>y</var> en mètres) correspond un déplacement de un degré de latitude (<var>Ï</var>).
- Les autres éléments (â<var>x</var>/â<var>Ï</var> et â<var>y</var>/â<var>λ</var>) sont des termes croisés (par exemple à quel déplacement
- en mètres vers le <em>Nord</em> correspond un déplacement de un degré de <em>longitude</em>).
- Ces valeurs ne sont généralement valides quâà la position géographique (<var>λ</var>,<var>Ï</var>) donnée.
- Si on se déplace un peu, ces valeurs changent légèrement.
- Cette matrice nous donne toutefois une bonne idée du comportement de la projection dans le voisinage du point projeté.
- </p>
-
- <p>
- On peut se représenter visuellement cette matrice 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 de dérivées.
- </p>
-
- <table class="hidden"><tr>
- <td><img style="border: solid 1px" src="../images/Derivatives.png" alt="Exemple de dérivées dâune projection cartographique"/></td>
- <td style="padding-left: 30px; vertical-align: middle">
- <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>
- </td>
- </tr></table>
-
- <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>
- <table class="hidden">
- <tr>
- <th>Enveloppe avant la projection</th>
- <th>Forme géométrique après la projection</th>
- </tr>
- <tr>
- <td><img style="border: solid 1px; margin-right: 15px" src="../images/GeographicArea.png" alt="Envelope in a geographic CRS"/></td>
- <td><img style="border: solid 1px; margin-left: 15px" src="../images/ConicArea.png" alt="Shape in a projected CRS"/></td>
- </tr>
- </table>
- <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>
- <table class="hidden"><tr><td>
- <img src="../images/Approximation.png" alt="Approximation cubique dâun bord dâune enveloppe"/>
- </td><td style="padding-left: 60px">
- 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>
- </td></tr></table>
- <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>
- <table class="hidden">
- <tr>
- <th style="text-align: left">Image source</th>
- <th style="text-align: right">Image destination</th>
- </tr>
- <tr>
- <td colspan="2"><img src="../images/RasterProjection.png" alt="Intersection of derivatives"/></td>
- </tr>
- </table>
- <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>
- <table class="hidden"><tr>
- <td><img src="../images/WarpGrid.png" alt="Intersection of derivatives"/></td>
- <td style="padding-left: 60px">
- 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>
- </td></tr></table>
- <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="../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>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);</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>@Override
-public Matrix derivative(DirectPosition p) throws TransformException {
- Matrix jac = inverse().derivative(transform(p));
- return Matrices.inverse(jac);
-}
-</pre>
</body>
</html>
Copied: sis/site/trunk/book/math/PixelToGeographicTerms.html (from r1745833, sis/site/trunk/book/math/PixelToGeographic.html)
URL: http://svn.apache.org/viewvc/sis/site/trunk/book/math/PixelToGeographicTerms.html?p2=sis/site/trunk/book/math/PixelToGeographicTerms.html&p1=sis/site/trunk/book/math/PixelToGeographic.html&r1=1745833&r2=1746518&rev=1746518&view=diff
==============================================================================
--- sis/site/trunk/book/math/PixelToGeographic.html (original)
+++ sis/site/trunk/book/math/PixelToGeographicTerms.html Thu Jun 2 04:57:06 2016
@@ -27,41 +27,52 @@
</head>
<body>
<math xmlns="http://www.w3.org/1998/Math/MathML" display="block" alttext="MathML capable browser required">
- <mfenced open="[" close="]">
- <mtable>
- <mtr><mtd><mi>λ</mi></mtd></mtr>
- <mtr><mtd><mi>Ï</mi></mtd></mtr>
- <mtr><mtd><mn>1</mn></mtd></mtr>
- </mtable>
- </mfenced>
- <mo>=</mo>
- <mfenced open="[" close="]">
- <mtable>
- <mtr>
- <mtd><msub><mi>S</mi><mrow>λ</mrow></msub></mtd>
- <mtd><msub><mi>H</mi><mrow>λ</mrow></msub></mtd>
- <mtd><msub><mi>T</mi><mrow>λ</mrow></msub></mtd>
- </mtr>
- <mtr>
- <mtd><msub><mi>H</mi><mrow>Ï</mrow></msub></mtd>
- <mtd><msub><mi>S</mi><mrow>Ï</mrow></msub></mtd>
- <mtd><msub><mi>T</mi><mrow>Ï</mrow></msub></mtd>
- </mtr>
- <mtr>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>0</mn></mtd>
- <mtd><mn>1</mn></mtd>
- </mtr>
- </mtable>
- </mfenced>
- <mo>Ã</mo>
- <mfenced open="[" close="]">
- <mtable>
- <mtr><mtd><mi>i</mi></mtd></mtr>
- <mtr><mtd><mi>j</mi></mtd></mtr>
- <mtr><mtd><mn>1</mn></mtd></mtr>
- </mtable>
- </mfenced>
+ <mtable>
+ <mtr>
+ <mtd><mo>λ</mo></mtd>
+ <mtd><mo>=</mo></mtd>
+ <mtd><msub><mi>S</mi><mrow>λ</mrow></msub><mi>x</mi><mo>+</mo><msub><mi>T</mi><mrow>λ</mrow></msub></mtd>
+ <mtd><mtext>        where        </mtext></mtd>
+ <mtd><msub><mi>S</mi><mrow>λ</mrow></msub></mtd>
+ <mtd><mo>=</mo></mtd>
+ <mtd>
+ <mfrac>
+ <mrow>
+ <msub><mi>λ</mi><mrow>max</mrow></msub><mo>-</mo>
+ <msub><mi>λ</mi><mrow>min</mrow></msub>
+ </mrow><mrow>
+ <msub><mi>N</mi><mrow><mi>x</mi></mrow></msub>
+ </mrow>
+ </mfrac>
+ </mtd>
+ <mtd><mtext>    and    </mtext></mtd>
+ <mtd><msub><mi>T</mi><mrow>λ</mrow></msub></mtd>
+ <mtd><mo>=</mo></mtd>
+ <mtd><msub><mi>λ</mi><mrow>min</mrow></msub></mtd>
+ </mtr>
+ <mtr>
+ <mtd><mo>Ï</mo></mtd>
+ <mtd><mo>=</mo></mtd>
+ <mtd><msub><mi>S</mi><mrow>Ï</mrow></msub><mi>y</mi><mo>+</mo><msub><mi>T</mi><mrow>Ï</mrow></msub></mtd>
+ <mtd><mtext>        where        </mtext></mtd>
+ <mtd><msub><mi>S</mi><mrow>Ï</mrow></msub></mtd>
+ <mtd><mo>=</mo></mtd>
+ <mtd>
+ <mfrac>
+ <mrow>
+ <msub><mi>Ï</mi><mrow>max</mrow></msub><mo>-</mo>
+ <msub><mi>Ï</mi><mrow>min</mrow></msub>
+ </mrow><mrow>
+ <msub><mi>N</mi><mrow><mi>y</mi></mrow></msub>
+ </mrow>
+ </mfrac>
+ </mtd>
+ <mtd><mtext>    and    </mtext></mtd>
+ <mtd><msub><mi>T</mi><mrow>Ï</mrow></msub></mtd>
+ <mtd><mo>=</mo></mtd>
+ <mtd><msub><mi>Ï</mi><mrow>min</mrow></msub></mtd>
+ </mtr>
+ </mtable>
</math>
</body>
</html>
|