Auteur : Paul Passy Licence : |cc_by_nc_sa| .. |cc_by_nc_sa| image:: figures/Cc-by-nc-sa_icone.png :width: 80 px .. _intro-SCR: Les Systèmes de Coordonnées de Références =========================================== Dans cette section nous allons découvrir quelques points théoriques et pratiques sur les systèmes de coordonnées de références (SCR), aussi appelés abusivement *projections*. Nous nous efforcerons de voir les implications de la théorie sur les applications en géomatique. Nous verrons ensuite quelques manipulations pratiques faisant intervenir les SCR sur des données vecteur et raster. .. contents:: Table des matières :local: .. _formes-terre: Les formes de la Terre ------------------------ La représentation de la forme de la Terre a évolué au cours du temps et peut être *choisie* en fonction des applications souhaitées. Définir la forme de la Terre est l'étape primordiale afin de s'y situer et donc, pour nous, d'y faire de l'analyse spatiale basée sur les systèmes d'information géographiques et la télédétection. Le choix d'un système géodésique, puis d'un système de projection a des répercussions sur le croisement des couches au sein d'un SIG et sur les calculs qu'il est possible de faire sur les longueurs et les surfaces des objets étudiés, qu'ils soient de type *vecteur* ou *raster*. .. _terre-sphere: La Terre comme une sphère *************************** Connaître la forme de la Terre est une question qui se pose depuis très longtemps. Dans le monde méditerranéen, les savants grecs furent parmi les premiers à avancer des réponses à cette question. Le premier (du moins le premier documenté) a avoir fait l'hypothèse d'une Terre sphérique serait Pythagore vers 500 avant notre ère. Cette idée viendrait des observations faites par les astronomes antiques lors des éclipses de Soleil. La Terre passe devant le Soleil et y forme alors une ombre en forme de croissant. Suite à cette hypothèse, un autre savant, du nom de Ératosthène (:numref:`eratosthene`), mit au point une méthode pour déterminer la circonférence de cette sphère. Il calcule une circonférence de 250 000 stades, une unité de mesure utilisée alors en Égypte. Selon certains historiens, un stade valait 157.5 mètres, ce qui donne une circonférence de 39 375 km. Cette valeur est impressionnante car diffère seulement de 2 % de la valeur actuellement connue. Cependant, nous ne sommes pas tout à fait certains de la longueur en mètres d'un stade égyptien. Quoi qu'il en soit, pour une époque si reculée, ces résultats sont impressionnants. .. figure:: figures/fig_Eratosthene_enseignant_a_Alexandrie.jpg :width: 35em :align: center :alt: Eratosthène enseignant à Alexandrie :figclass: align-center :name: eratosthene Eratosthène enseignant à Alexandrie, peinture de Bernardo Strozzi, 1635. Selon le géographe antique romain Strabon, Cratès de Mallos aurait été le premier a réalisé un globe terrestre pour représenter physiquement la Terre sous sa forme sphérique. Il y aurait placé le monde méditerranéen connu de son époque, plus des continents inconnus faisant symétrie au monde connu. Selon lui, l'hémisphère nord et sud devaient être séparées par une zone océanique appelée *zone torride*. .. _terre-geoide: La Terre comme un géoïde *************************** Durant le Moyen-Âge, une image de la Terre plate a été répandue par le clergé catholique, même si certains savants pensaient toujours que la Terre était sphérique. Cette sphéricité de la Terre est finalement revenue totalement sur le devant de la scène à la Renaissance, jusqu'à ce que des observations plus fines au cours des 19ème et 20ème siècles précisent mieux la forme de notre planète pour la définir comme un *géoïde*. Un géoïde est une surface équipotentielle du champ de pesanteur. Nous pouvons nous aider d'une image pour visualiser cet objet. Imaginons une planète totalement recouverte par un océan qui ne subirait aucune force de marée. La surface océanique d'une telle planète épouserait alors parfaitement la forme du geoïde de la planète en question. Sur Terre, il s'agit donc du niveau moyen des mers avec sa prolongation imaginaire sous les continents. Le géoïde n'est pas une sphère pour deux raisons. Tout d'abord, du fait de sa rotation, la Terre est légèrement aplatie au niveau des pôles. Ainsi, le rayon terrestre entre le centre de la Terre et l'Équateur est d'à peu près 6378 km et entre le centre de la Terre et les pôles d'environ 6357 km. D'autre part, l'intérieur de la planète n'étant pas homogène, certaines zones sont constituées de roches plus denses que d'autres. Au niveau des zones de roches les plus denses, le géoïde tend à se *creuser*, alors que là où les roches sont moins denses, il tend à se *bomber* (:numref:`geoide-esa`). .. figure:: figures/fig_geoid_esa.jpg :width: 16em :align: center :alt: Illustration du géoïde :figclass: align-center :name: geoide-esa Illustration, très exagérée, du géoïde terrestre (ESA). La forme précise du geoïde peut être définie par des mesures in-situ ou via des missions satellites dédiées à la mesure du champ de pesanteur terrestre, (`comme la mission Grace`_). Un modèle souvent utilisé pour le géoïde terrestre est le EGM96 (*Earth Gravitational Model 96*). Selon ce modèle, le point le plus *bas* du géoïde se situe au sud de l'Inde dans l'Océan Indien et le point le plus *haut* se situe au large de l'Islande dans l'Atlantique Nord. La différence entre ces deux points atteint un peu moins de 200 de mètres. Au sud de l'Inde, le géoïde est à `-106 m de sa hauteur moyenne`_ et au large de l'Islande il est à +85 m de sa hauteur moyenne (:numref:`egm96`). .. figure:: figures/fig_Earth_Gravitational_Model_1996.png :width: 35em :align: center :alt: Géoïde EGM96 :figclass: align-center :name: egm96 Le modèle EGM96 comme représentation du géoïde terrestre. Les ondulations du géoïde sont finalement peu importantes par rapport aux ondulations de la topographie. Le point culminant de la surface terrestre, l'Everest, culmine à 8849 m et le point le plus bas, la Fosse des Mariannes, s'enfonce à près de 10 985 mètres sous le niveau de la mer. La différence topographique est donc d'une vingtaine de kilomètres. Cette différence atteint tout de même près de 9200 m pour les seules surfaces continentales puisque le point le plus bas, la Mer Morte, se situe à -434 m. .. _terre-ellipsoide: La Terre comme un ellipsoïde ****************************** Travailler avec un objet aussi complexe qu'un géoïde s'avère extrêmement compliqué. C'est pour cette raison que nous définissons un nouvel objet conceptuel appelé *ellipsoïde*. Un ellipsoïde est une simplification du géoïde et peut être perçu comme la *moyenne* du géoïde. Finalement, en 3D un ellipsoïde est au géoïde ce que l'ellipse est au cercle en 2D. Étant une forme géométrique idéale, un ellipsoïde peut être définie par une série de paramètres : - le demi grand axe *a* : correspond au rayon entre le centre de l'ellipsoïde et l'Équateur - le demi petit axe *b* : correspond au rayon entre le centre de l'ellipsoïde et un des Pôles - l'inverse de l'aplatissement : :math:`\frac{1}{f}=\frac{a}{a-b}` - la première excentricité : :math:`e=\sqrt{\frac{a^2-b^2}{a^2}}` - le carré de l'excentricité : :math:`e^2=\frac{a^2-b^2}{a^2}` - la deuxième excentricité : :math:`e'=\sqrt{\frac{a^2-b^2}{b^2}}` Par endroits, l'ellipsoïde est plus *haut* que le géoïde et à d'autres, il est plus *bas*. La figure suivante (:numref:`geoide-ellipsoide`) présente une vue schématique de ces deux objets où l'ondulation du géoïde a été largement exagérée pour illustration. Si nous faisons abstraction des marées et des courants, la surface des océans suit la surface du géoïde. .. figure:: figures/fig_geoide_ellipsoide.png :width: 35em :align: center :alt: L'ellipsoïde comme simplification du géoïde :figclass: align-center :name: geoide-ellipsoide L'ellipsoïde comme simplification du géoïde. Il est ensuite possible de définir un **système géodésique** en faisant un choix d'ellipsoïde. La forme de l'ellipsoïde, l'orientation de ses axes et sa position par rapport au centre de la Terre permettent de définir un système géodésique, aussi connu sous le nom de *Datum* en anglais. La figure suivante (:numref:`datums`) présente deux systèmes géodésiques différents de part leurs formes et leurs orientations. .. figure:: figures/fig_datums.png :width: 35em :align: center :alt: Deux systèmes géodésiques différents :figclass: align-center :name: datums Deux systèmes géodésiques (A et B) différents selon leurs formes et leurs orientations. Il existe des ellipsoïdes globaux et locaux. Un ellipsoïde global va moyenner le géoïde à l'échelle de toute la Terre. Tandis qu'un ellipsoïde local va le moyenner à l'échelle d'une région donnée. Le second sera donc localement plus précis que le premier mais ne pourra être utilisé que sur une partie du globe. Par exemple, en Amérique du Nord l'ellipsoïde de *Clarke 1866* a été utilisé au 19ème siècle et début 20ème car il était plus proche du géoïde pour cette partie du monde. Néanmoins, avec l’avènement des systèmes de positionnement globaux par satellite (comme le GPS), la plupart des systèmes géodésiques actuels repose sur le même ellipsoïde global nommé *Geodetic Reference System 1980*, connu simplement comme l'ellipsoïde **GRS80**. Une (très faible) variante de cet ellipsoïde est utilisé dans le système géodésique mondial **WGS84**. La différence entre ces deux ellipsoïdes n'est pas significative dans les applications liées à la géomatique comme l'indique le tableau suivant. .. list-table:: Paramètres de trois ellipsoïdes :widths: 15 25 15 25 :header-rows: 1 * - Ellipsoïde - Demi grand axe (a) - Demi petit axe (b) - Inverse de l'aplatissement * - Clarke 1866 - 6 378 206.4 m - 6 356 583.8 m - 294.978698214 * - GRS80 - 6 378 137.0 m - 6 356 752.314140 m - 298.257222100 * - WGS84 - 6 378 137.0 m - 6 356 752.314245 m - 298.257223563 Entre l'ellipsoïde qui était utilisé en Amérique du Nord (Clarke 1866) et les deux utilisés actuellement à l'échelle globale, la différence est de moins de 100 m sur les deux axes. Quant aux deux ellipsoïdes *GRS80* et *WGS84* ils sont pour ainsi dire identiques, car seul leur petit axe diffère de moins de 1 cm. Une fois le système géodésique défini, il est possible de s'y positionner et d'y faire des traitements spatiaux. .. _coord-geo: Les coordonnées géographiques ****************************** Afin de se positionner sur un ellipsoïde, il faut définir un repère géographique avec deux axes d'origine, un pour les latitudes et un pour les longitudes (:numref:`long-lat`). .. figure:: figures/fig_longitude_latitude.png :width: 27em :align: center :alt: Longitudes et latitudes :figclass: align-center :name: long-lat Coordonnées du point *M* en longitude (λ), latitude (φ) et hauteur par rapport à l'ellipsoïde (h). Sur la figure :numref:`long-lat`, le point M peut se repérer par trois coordonnées : sa longitude, sa latitude et son altitude. La longitude, notée *λ* ou *X*, est la valeur de l'angle entre le méridien d'origine et le méridien de M. La longitude s'exprime en degrés *ouest* et *est* selon que le point se trouve à moins de 180° à l'ouest ou à l'est du méridien d'origine. La latitude, notée *φ* ou *Y* est la valeur entre le parallèle d'origine et la droite coupant perpendiculairement le méridien de M. Elle est exprimée en degrés *nord* ou *sud* par rapport au parallèle d'origine, entre 0° et 90° (correspondant aux Pôles). Enfin, l'altitude de M est l'écart, exprimé en mètres, entre M et le niveau de référence vertical à savoir le géoïde généralement. L'axe d'origine pour les latitudes ne pose aucun problème, il s'agît de l'**Équateur**, ligne imaginaire coupant la Terre en deux hémisphères égaux. L'Équateur est également le parallèle le plus long, tous les autres rétrécissant vers les Pôles. Les parallèles 90° nord ou sud ne sont d'ailleurs plus que des points. Par contre, le méridien d'origine n'es pas aussi évident. En effet, tous les méridiens ont la même longueur et il n'est pas possible de définir des hémisphères est et ouest comme nous définissons les hémisphères nord et sud. Pendant longtemps, chaque pays utilisait son propre méridien d'origine passant généralement par sa capitale. Ainsi en France, le méridien de Paris était utilisé. Mais durant l'Antiquité, le point le plus à l'ouest connu, dans l'Atlantique à l'ouest du Détroit de Gibraltar, était utilisé. Ce n'est que lors de l'`International Meridian Conference de 1884`_ tenue à Washington aux États-Unis, que le **méridien de Greenwich** a été choisi comme le seul méridien d'origine (:numref:`eq-green`). Il était déjà utilisé par les anglais dont la flotte dominait alors les mers du globe. C'est en partie en raison de cette prédominance maritime anglaise que ce méridien a été choisi. Une petite mise à jour a été faite depuis cette date, et le méridien d'origine actuel passe 102 mètres plus à l'ouest que le méridien d'origine historique. Pour en savoir plus sur l’épopée du calcul des longitudes, vous pouvez vous référer au `très bon livre de Dava Sobel`_ intitulé *Longitude*. .. figure:: figures/fig_equateur_greenwich.png :width: 40em :align: center :alt: Équateur et Greenwich :figclass: align-center :name: eq-green L'Équateur faisant office de latitude 0° passant au travers d'un café en Ouganda (A) et matérialisation du méridien de Greenwich de longitude 0° (B). Dans ce repère géographique constitué de méridiens, qui graduent la longitude, et de parallèles, qui graduent la latitude, il est possible de repérer tout point de la planète. L'arc entre deux méridiens est subdivisé en 60 portions égales, figurant 60 minutes. Puis chaque arc de minute est lui-même subdivisé en 60 portions égales, figurant 60 secondes. Il en est de même pour les latitudes. Ainsi, par exemple la `Place de l'Étoile Rouge à Cotonou`_ au Bénin a pour coordonnées 6°22'17.1"N 2°24'35.5"E. Le long d'un méridien, un degré équivaut à peu près à 111.3 km. Une minute vaut alors 1.855 km, unité utilisée pour le mile nautique. Enfin, une seconde vaut à peu près 30 mètres. Pour les parallèles, ces valeurs ne sont valables que pour l'Équateur. La taille des parallèles diminuant à mesure que nous nous éloignons de l'Équateur, la taille de chaque subdivision diminue également. Certains parallèles présentent des propriétés astronomiques qui leur ont valu d'être identifés nommément (:numref:`globe-long-lat`). .. figure:: figures/fig_globe_longitude_latitude.png :width: 30em :align: center :alt: Longitudes et latitudes :figclass: align-center :name: globe-long-lat Les longitudes et les latitudes comme repère géographique et les parallèles notables. D'un point de vue astronomique, l'Équateur est la ligne où la durée du jour, 12 heures, est toujours égale à la durée de la nuit, quel que soit la période de l'année. Les tropiques sont les parallèles extrêmes où le Soleil passe au zénith au moins un jour dans l'année, respectivement au solstice de juin pour le Tropique du Cancer et au solstice de décembre pour le Tropique du Capricorne. Entre les Tropiques et les Pôles, le Soleil ne sera jamais au zénith. Enfin, les cercles polaires sont les latitudes extrêmes où le Soleil disparaît sous l'horizon au moins pendant 24 heures. Cette disparation a lieu lors du solstice de juin côté Antarctique et lors du solstice de décembre côté Arctique. En plus de ces coordonnées en longitude et latitude, la troisième composante du positionnement est l'altitude. L'altitude s'exprime en mètre et correspond à l'écart entre la surface topographique et un niveau de référence. Ce niveau de référence peut être la surface de l'ellipsoïde mais le plus souvent il s'agît de la surface du géoïde, appelée communément *niveau de la mer* (:numref:`altitude-geoide`). .. figure:: figures/fig_altitude_geoide.png :width: 35em :align: center :alt: Altitude et référentiel :figclass: align-center :name: altitude-geoide Altitude et hauteurs par rapport à l'ellipsoïde. D'un point de vue pratique, les coordonnées géographiques s'obtenaient à l'aide du `sextant`_, pour la latitude, du `chronomètre de marine`_, pour la longitude et du `baromètre`_ pour les altitudes. Aujourd'hui ces trois composantes sont fournies par les systèmes globaux de positionnement par satellite, **GNSS** en anglais pour *Global Navigation Satellite System*. Il existe quatre systèmes actuellement, maintenus par quatre puissances géopolitiques différentes : *GPS* maintenu par les États-Unis, *GLONASS (ГЛОНАСС)* développé par l'`Union Soviétique`_ et maintenu aujourd'hui par la Russie, *Galileo* opéré par l'Union Européenne et *Beidu (北斗)* maintenu par `la Chine`_. Ces quatre systèmes sont interopérables et utilisables par la plupart des récepteurs GPS et smartphones pour se localiser. La plupart des récepteurs GPS intègrent un modèle de géoïde ce qui permet d'obtenir l'altitude. Néanmoins, ces coordonnées géographiques sont parfaites pour se positionner sur un globe en trois dimensions, mais ne sont pas idéales pour la cartographie ou la géomatique qui a besoin de se faire sur des surfaces en deux dimensions. Il est donc nécessaire de transformer cet ellipsoïde en un plan muni d'un nouveau repère en deux dimensions. Ce processus s'appelle une **projection**. .. _projections-definition: Les projections, éléments de définition ---------------------------------------- Pour passer de l'objet 3D ellipsoïde à l'objet 2D carte ou écran, nous utilisons une projection. La figure suivante (:numref:`geoide-proj`). récapitule la partie précédente et y place la projection. .. figure:: figures/fig_geoide_to_projection.png :width: 40em :align: center :alt: Du géoïde à la projection :figclass: align-center :name: geoide-proj De la Terre à la carte en passant par la projection. Mettre *à plat* un ellipsoïde passe donc par un processus de **projection** présentant trois composantes : - la nature de la projection : * Conforme * Équivalente * Équidistante * Aphylactique - le type du plan de la projection : * Cylindrique * Conique * Azimutal - la position du plan de projection : * Parallèle à l'Équateur (normale) * Perpendiculaire à l'Équateur (transverse) * Oblique Ces différentes composantes sont détaillées dans la partie qui suit. .. _projections-natures: La nature des projections *************************** Mathématiquement, il est impossible de projeter à plat une forme sphérique sans entraîner des déformations. Pour ce qui intéresse le monde de la géomatique, nous devons choisir entre **conserver les surfaces** et donc utiliser une projection **équivalente** ou **conserver les formes**, ce qui revient à conserver les angles, et donc utiliser une projection **conforme**. La figure suivante (:numref:`proj-conf-eq`) éclaire le propos en faisant varier un cercle en surface ou en forme selon la projection choisie. Sur cette figure, le cercle central est répété à l'identique vers les latitudes supérieures et inférieures. Le changement visible sur le cercle n'est dû qu'à la nature de la projection choisie. .. figure:: figures/fig_cercle_conforme_equivalent.png :width: 40em :align: center :alt: Déformations selon la projection. :figclass: align-center :name: proj-conf-eq La projection conforme conserve la forme du cercle mais pas sa taille (A) et la projection équivalente conserve sa taille mais pas sa forme (B). Sur la projection conforme, le cercle reste un cercle mais il semble de plus en grand à mesure que nous nous éloignons de l'Équateur. Au contraire, sur la projection équivalente, le cercle conserve sa surface initiale mais il se déforme de plus en plus à mesure que nous nous dirigeons vers les Pôles. Pour se rendre compte de ce phénomène sur une carte du monde, il est possible d'utiliser les polygones de Tissot. C'est le même principe que les cercles mentionnés précédemment mais placés sur une carte du monde de projection donnée. Sur la figure suivante (:numref:`tissot-mercator`) nous avons une carte du monde en projection de Mercator, sans doute la projection la plus célèbre. Il s'agît d'une projection *conforme*, sur laquelle les formes et les angles sont conservées mais pas les surfaces comme nous le montrent les polygones de Tissot placés sur cette projection. Sur cette projection, les terres au niveau de l'Équateur présentent une forme et une superficie proches de la réalité. Par contre, vers les Pôles la forme est conservée mais la superficie est extrêmement exagérée. Par exemple, le Groenland semble faire la taille de l'Amérique du Sud, alors qu'en réalité il est un peu plus petit que l'Algérie. Cependant cette projection est pratique pour la navigation car les angles, et donc les azimuts, correspondent aux angles réels. .. figure:: figures/fig_Tissot_Mercator.png :width: 40em :align: center :alt: Projection de Mercator et polygones de Tissot. :figclass: align-center :name: tissot-mercator Polygones de Tissot sur projection conforme de Mercator (à gauche) et le Monde en projection de Mercator (à droite). Inversement, sur une projection équivalente le cercle perd sa forme de cercle mais conserve sa superficie sur toute la surface projetée. À mesure que nous nous éloignons du, ou des parallèles sécants au plan de projection un cercle se déforme de plus en plus mais conserve sa surface. La figure suivante (:numref:`tissot-behrmann`) est une carte du monde en projection équivalente de Behrmann. Sur cette projection, le plan de projection est sécant à l'ellipsoïde au niveau des parallèles 30°N et 30°S. Sur ces parallèles les cercles présentent donc une forme proche de la réalité puis se déforment à mesure qu'on s'en éloigne. Même si ils sont très déformés vers les hautes latitudes, leurs superficies sont conservées. .. figure:: figures/fig_Tissot_Behrmann.png :width: 40em :align: center :alt: Projection de Behrmann et polygones de Tissot. :figclass: align-center :name: tissot-behrmann Polygones de Tissot sur projection équivalente de Behrmann. Il existe d'autres types de projections comme les projections **équidistantes** et les projections **aphylactiques**. Néanmoins ces projections sont peu utilisées en géomatique, et surtout utilisées pour des représentations cartographiques. La projection équidistante est centrée sur un point précis et conserve ensuite les longueurs entre le point central et tout cercle l'englobant. C'est une représentation qui peut être utilisée dans les cartographies des séismes. le point central est l'épicentre du séisme, puis nous avons des cercles concentriques autour de ce point indiquant les distances à l'épicentre. Dans un autre contexte, cette projection est utilisée pour cartographier les portées des missiles depuis une base de lancement. Les projections aphylactiques, quand à elles, ne conservent ni les formes ni les surfaces mais minimisent les deux erreurs. Quelque soit la nature de la projection, une carte peut se projeter selon des plans de projection de différents types. .. _projections-plans: Le type de plans de projection ******************************* Le plan de projection est une surface plane sur laquelle est projetée l'ellipsoïde. Ce plan de projection peut être de trois grands types différents. - **Plan cylindrique** : ce plan est comme *enroulé* sous forme de cylindre puis positionné autour de l'ellipsoïde. Le cylindre peut être positionné avec des angles différents autour de l'ellipsoïde et peut simplement le *toucher* sur un seul cercle, ou bien le *couper* sur deux cercles. Dans le premier cas le plan est dit **tangent** et dans le second cas **sécant**. Le cylindre ne fera que *toucher* si son diamètre est égale à la circonférence terrestre et *coupera* si son diamètre est inférieur à la circonférence terrestre. Plus de précisions seront données dans les paragraphes qui suivent. Par exemple, les projections *UTM* sont des projections cylindriques sécantes. - **Plan conique** : ce plan est enroulé sous forme de *chapeau chinois* puis posé sur l'ellipsoïde. Il peut être posé n'importe où, sous n'importe quel angle. Il peut simplement *toucher* l'ellipsoïde (plan tangent) sur un seul cercle ou *couper* l'ellipsoïde selon deux cercles (plan sécant). Si les deux cercles sécants sont des parallèles, nous parlons de **parallèles automécoïques**, nouveau mot qui permet de briller en société. Par exemple, la projection *Lambert93* est une projection utilisant un plan conique sécant selon deux parallèles automécoïques. - **Plan azimutal** : ce plan est laissé *à plat* et l'ellipsoïde est soit *posé dessus* soit *coupé* par ce plan. Si l'ellipsoïde est simplement *posé*, le contact entre l'ellipsoïde et le plan ne sera qu'un point mais si l'ellipsoïde est *coupé* par le plan, alors le contact entre les deux sera un cercle. La figure suivante (:numref:`types-plans`) schématise les trois grands types de plans de projection pour des positions fixées. Sur cette figure les plans sont tous *tangents* mais il est possible d'avoir les mêmes projections avec des plans *sécants*. .. figure:: figures/fig_types_plans_projection.png :width: 40em :align: center :alt: Types de plans de projection. :figclass: align-center :name: types-plans Les trois grands types de plans de projection. Comme évoqué, les plans sont de types différents mais présentent également des positions différentes par rapport à l'ellipsoïde. Nous parlons alors de la position des plans de projection par rapport à l'ellipsoïde. .. _projections-positions: La position du plan de projection ********************************** Que le plan soit cylindrique, conique ou azimutal, sa position vis-à-vis de l'ellipsoïde peut varier. D'une manière générale le plan est **tangent** si il ne fait que *toucher* l'ellipsoïde et **sécant** si il le coupe. En plus de cette propriété, nous pouvons définir le plan par rapport à sa position vis-à-vis de l'Équateur. Si le plan *touche* ou *coupe* l'ellipsoïde parallèlement à l'Équateur, le plan est **normal**. Si le plan *touche* ou *coupe* l'ellipsoïde perpendiculairement à l'Équateur, le plan est **transverse**. Pour toute position ni parallèle ni perpendiculaire, le plan est **oblique**. Ces trois positions peuvent s'appliquer pour les trois types de plans, cylindrique, conique et azimutal. La figure suivante (:numref:`positions-plans`) présente les trois positions possibles pour une projection cylindrique. .. figure:: figures/fig_positions_plans.png :width: 40em :align: center :alt: Positions de plans de projection. :figclass: align-center :name: positions-plans Positions du plan de projection pour une projection cylindrique. Le principe est le même pour une projection conique ou azimutale. Au final, ces différentes caractéristiques se croisent pour former une projection proprement dite. Par exemple, la projection *Lambert93* utilisée en France métropolitaine est une projection *équivalente conique sécante*. Le plan de projection est donc conique est coupe l'ellipsoïde selon deux parallèles automécoïques. Les projections *UTM* sont, quant à elles, des projections *conformes cylindriques tangentes transverses* (comme sur le schéma du milieu sur la figure précédente). La figure suivante (:numref:`monde-proj`) présente deux vues du monde selon une projection normale et une projection transverse. Dans le premier cas les altérations sont minimales le long de l'Équateur et dans le second elles sont minimales le long du méridien de Greenwich. .. figure:: figures/fig_proj_normale_transverse.png :width: 40em :align: center :alt: Monde en projection normale ou transverse. :figclass: align-center :name: monde-proj Le monde vu en projection normale centrée sur l'Équateur sur une carte soviétique des années 1920 (A) et vu en projection transverse le long du méridien de Greenwich (B). Par la suite seront présentées plus en détails quelques grands types de projections utilisées à l'échelle de la France, de l'Europe et du monde notamment par les données globales ou les images de télédétection. .. _codes-EPSG: Les codes EPSG **************** Les systèmes de coordonnées de références ont des noms officiels mais également des noms plus courts utilisés de façon plus simple. Par exemple, la zone UTM pour la longitude de Paris est la 31. Le nom complet du SCR de cette zone est *WGS 84 / UTM zone 31N* mais nous parlons couramment de *Zone UTM 31 N*. Mais il existe également une zone UTM 31 basée sur l'ancien ellipsoïde WGS1972. Afin de s'y retrouver plus facilement, un ancien consortium nommé *European Petroleum Survey Group* a proposé une codification des SCR en leur attribuant un code unique, appelée **code EPSG**. Ainsi, chaque SCR possède un code EPSG par lequel il est pratique de se référer dans les logiciels. .. _projections-courantes: Quelques projections courantes ------------------------------- Nous verrons dans cette partie quelques projections et systèmes de coordonnées de référence (SCR) couramment employés en France, en Europe et pour les données globales. Le tableau suivant récapitule les principaux systèmes de coordonnées de référence utilisés dans le contexte français. .. list-table:: Codes EPSG de quelques SCR :widths: 35 25 25 25 :header-rows: 1 * - Système de coordonnées - Code EPSG - Type de projection - Domaine de validité * - RGF93/Lambert93 - 2154 - Conforme - France métropolitaine * - LAEA/ETRS89 - 3035 - Équivalente - Europe * - WGS84 / UTM 30N - 32630 - Conforme - 6°W and 0°W * - WGS84 / UTM 31N - 32631 - Conforme - 0°E et 6°E * - WGS84 / UTM 32N - 32632 - Conforme - 6°E et 12°E * - WGS84 - 4326 - Non projeté - Monde * - WGS84 / Pseudo Mercator - 3857 - Conforme - Monde * - WGS 84 / Equal Earth Greenwich - 8857 - Équivalente - Monde Le *SCR WGS84 (4326)* est un système non projeté à l'échelle globale dont les unités sont des degrés. Par contre, le *SCR WGS84 / Pseudo Mercator (3857)* est une projection cylindrique normale conforme, également utilisable à l'échelle globale, et ses unités sont donc des mètres. Cette projection est notamment utilisée par les services de webmapping mondiaux comme GoogleMaps et OpenStreetMap. .. _Lambert-93: Pour la France ***************** En France métropolitaine, depuis 2000, la projection officielle est le **Lambert 93**, basé sur le **RGF93** (*Réseau Géodésique Français 93*). Le nom complet du SCR est donc **RGF93/Lambert93**. Le code EPSG de ce SCR est le **2154**. Cette projection est **conique conforme sécante**. Le plan de projection est donc de forme conique et est sécant à la surface de la Terre. De plus, étant *conforme*, cette projection conserve les formes mais pas les surfaces. Ainsi, la France aura sa *vraie* forme mais pas sa *vraie* surface. Cependant, la déformation est minime car elle est de l'ordre de `3m/km aux extrémités du domaine d'application`_. Concrètement, la déformation est de 2m/km à Dunkerque et de 3m/km à Bonifacio. Comme `indiqué par l'IGN`_, le Lambert93 présente les paramètres du tableau suivant. .. list-table:: Principaux paramètres de la projection Lambert 93 :widths: 30 25 :header-rows: 1 * - Paramètre - Valeur * - Ellipsoïde de référence - GRS80 * - Type de projection - Conique Conforme directe sécante * - Zone d'application - 41° - 51° Nord * - Méridien central - 3° Est * - Latitude origine - 46°30'Nord * - Décalage est (False easting) - 700 000 m * - Décalage nord (False northing) - 6 600 000 m * - Premier parallèle automécoïque - 44°N * - Second parallèle automécoïque - 49°N Le *méridien central* est le méridien sur lequel le cône de projection est centré. Le méridien 3° Est correspond à peu près au méridien qui coupe la France en deux parties égales dans le sens de la hauteur. La *latitude origine* est à peu près la latitude qui coupe la France en deux parties égales dans le sens de la largeur. Comme le Lambert93 est une projection, nous nous retrouvons sur un plan en deux dimensions qui dispose d'une origine arbitraire. Cette origine est positionnée très loin de la France, quelque part dans le Golfe de Guinée, comme présenté sur la figure suivante (:numref:`L93-repere`). .. figure:: figures/fig_repere_Lambert93.png :width: 25em :align: center :alt: Le repère de la projection Lambert93. :figclass: align-center :name: L93-repere L'origine du repère utilisé dans la projection Lambert 93. Le point d'origine (0,0) du repère de la projection Lambert 93 se situe à 700 000 m à l'ouest du méridien central de la projection (le méridien 3° E) et à 6 600 000 m au sud de la latitude origine 46°30'N. Autrement dit, le point d'intersection entre le méridien central et la latitude origine se trouve à 700 000 m à l'ouest du point d'origine (*False easting*) et à 6 600 000 m au sud (*False northing*) du point d'origine. Ce fort décalage a été choisi pour clairement différencier le Lambert 93 d'autres projections françaises utilisées par le passé. Au final, si nous recentrons le repère sur la France, nous obtenons la figure suivante (:numref:`L93-france`). .. figure:: figures/fig_France_Lambert93.png :width: 40em :align: center :alt: La France dans la projection Lambert93. :figclass: align-center :name: L93-france La France représentée dans la projection Lambert 93. Ainsi, tout point du territoire présente des coordonnées comprises entre 100 000 m et 1 300 000 m en latitude et 6 000 000 m et 7 000 000 en longitude. Nous retrouvons bien les 1000 km (1 000 000 m) d'envergure nord-sud et est-ouest du territoire français. Par exemple, le `bâtiment de la Halle aux Farines de l'Université Paris Cité`_ dans le 13ème arrondissement de Paris a pour coordonnées en Lambert 93 X : 654 612.688 m et Y : 6 858 969.046 m. Si à partir de ces coordonnées nous nous déplaçons de 100 m plein est, nous nous retrouvons au point de coordonnées X : 654 712.688 m et Y : 6 858 969.046 m. De même si nous nous déplaçons de 100 m plein nord, nous nous retrouvons au point de coordonnées X : 654 612.688 m et Y : 6 859 069.046 m. Enfin, les parallèles *automécoïques* sont les parallèles sur lesquels le cône de projection est sécant à la surface terrestre (:numref:`par-auto`). Les erreurs de surface sont minimales le long de ces deux parallèles. .. figure:: figures/fig_paralleles_automecoiques.png :width: 20em :align: center :alt: Exemple de parallèles automécoïques. :figclass: align-center :name: par-auto Représentation en 2 dimensions du concept de parallèles automécoïques. .. tip:: QGIS propose un petit outil de visualisation pour vérifier l'étendue de validité d'un SCR. Cet outil s'affiche dès que vous souhaitez régler un SCR. Par exemple, même si vous ne savez pas à priori quelle portion du globe est couverte par la projection *Pulkovo 1942 / Gauss-Kruger zone 27*, vous pouvez voir dans l'outil dédié, que cette projection est valable pour le Kamtchatka. .. figure:: figures/fig_proj_kamtchatka.png .. _ETRS89: Pour l'Europe ************** L'Union Européenne a défini une projection officielle à utiliser pour les données à l'échelle de l'Europe. Il s'agît de la projection **ETRS89-LAEA** pour *European Terrestrial Reference System 1989 - Lambert Azimuthal Equal Area*. Il s'agit donc d'une projection azimutale **équivalente**. Sur cette projection, les contours sont déformés mais les surfaces sont respectées. Le caractère équivalent de cette projection a été choisi pour pouvoir cartographier des phénomènes en lien avec les surfaces. C'est notamment le cas de tout ce qui est densité : densité de population, rendement agricole, densité de cheptel... Les principaux paramètres de cette projection sont résumés dans le tableau suivant. Le code EPSG de cette projection est le **3035**. .. list-table:: Principaux paramètres de la projection ETRS89 :widths: 30 25 :header-rows: 1 * - Paramètre - Valeur * - Ellipsoïde de référence - GRS80 * - Type de projection - Azimutal * - Zone d'application - Europe * - Méridien central - 10° Est * - Latitude origine - 52° Nord * - Décalage est (False easting) - 4 321 000 m * - Décalage nord (False northing) - 3 210 000 m Sur cette projection, les déformations sont minimales à l'intersection du méridien 10 °E et du parallèle 52°N, puis s'accroissent à mesure que nous nous éloignons de ce point dans toutes les directions. L'Europe prend alors un air arrondi (:numref:`europe-etrs89`). .. figure:: figures/fig_Europe_3035.png :width: 38em :align: center :alt: L'Europe dans la projection ETRS89. :figclass: align-center :name: europe-etrs89 L'Europe dans la projection ETRS89. .. _UTM-zones: Pour le monde par zones UTM ***************************** À l'échelle globale, *i.e.* à l'échelle du monde, nous pouvons utiliser les projections en **Zones UTM** pour *Universal Transverse Mercator*. Il s'agît d'une projection *cylindrique transverse conforme*. Le plan de projection est donc un cylindre *posé* perpendiculairement à l'Équateur et tangent à la Terre selon un méridien central (:numref:`utm-principes`). Nous parlons de *zones UTM* car il s'agit en fait d'un corpus de 60 projections suivant le même principe. Le globe terrestre, mesurant par définition 360° de circonférence, a été divisé en 60 zones de 6° de largeur chacune, numérotées de *01* à *60*. Chaque zone peut être perçue comme une sorte de *quartier* de Terre, comme les différents quartiers d'une orange. Pour chaque zone, un méridien central coupant la zone en deux parties égales est défini et le cylindre de projection est positionné sur ce méridien. Puis, pour la zone suivante, le cylindre est pivoté de 6° et tombe ainsi sur le méridien central de la zone suivante, et ainsi de suite 60 fois. Comme il s'agît d'un type de projection *conforme*, les formes sont respectées mais pas les surfaces. Les images **Landsat** et **Sentinel-2** sont notamment fournies dans ces projections. .. figure:: figures/fig_utm_principes.jpg :width: 38em :align: center :alt: Le principe des zones UTM. :figclass: align-center :name: utm-principes Les 60 zones UTM (à gauche) et le cylindre de projection transverse utilisé pour chaque zone (à droite). Le méridien en rouge est le méridien le long duquel le cylindre est tangent à la surface terrestre. Les 60 zones sont numérotées d'ouest en est avec la zone *01* positionnée sur la pointe extrême orientale de l'Eurasie. Ces zones sont ensuite chacune divisées en deux : une zone au nord de l'Équateur et une zone au sud de l'Équateur (:numref:`utm-decoup`). Nous parlons ainsi de zones UTM Nord et Sud. Par exemple, Paris se trouve sur la zone UTM 31 N et `Hobart`_ (Tasmanie, Australie) sur la zone UTM 55 Sud. .. figure:: figures/fig_utm_zones.jpg :width: 40em :align: center :alt: Les 60 zones UTM Nord et Sud. :figclass: align-center :name: utm-decoup Le découpage en 60 zones UTM de 6° d'est en ouest et subdivisées en nord et sud. La France métropolitaine tombe sur trois zones UTM : 30, 31 et 32 Nord. La Bretagne est sous la zone 30 N, l'Îlde-de-France sous la zone 31 N et la Corse sous la 32 N (:numref:`utm-france`). Les codes EPSG de ces trois projections sont respectivement **32630**, **32631** et **32632**. .. figure:: figures/fig_utm_zone_France.png :width: 27em :align: center :alt: Les 3 zones UTM de la France. :figclass: align-center :name: utm-france La France sous les zones UTM 30N, 31N et 32N. .. tip:: Il est facile de retrouver les codes ESPG des zones UTM. Les 120 zones (60 Nord et 60 Sud) ont toutes un code EPSG à 5 chiffres commençant par *32*. Le troisième chiffre est un *6* pour les zones Nord et un *7* pour les zones Sud. Enfin, les deux derniers chiffres correspondent au numéro de la zone, de *01* à *60*. Ainsi, la zone UTM 31 Nord a pour code *32631* et la zone UTM 22 Sud a pour code *32722*. Le tableau suivant présente les principaux paramètres de la projection UTM 32 Nord. Ces paramètres seront tout à fait identiques pour toutes les zones UTM de l'hémisphère nord à l'exception du méridien central qui changera de 6° en 6° et de la zone d'application. .. list-table:: Principaux paramètres de la projection UTM 32 Nord :widths: 30 25 :header-rows: 1 * - Paramètre - Valeur * - Ellipsoïde de référence - WGS84 * - Type de projection - Cylindrique * - Zone d'application - De 6°E à 12°E * - Méridien central - 9° Est * - Latitude origine - 0° Nord * - Décalage est (False easting) - 500 000 m * - Décalage nord (False northing) - 0 m Nous remarquons la présence d'un *false easting* qui décale l'origine du repère de projection à l'ouest du méridien central, ce qui permet de n'avoir que des coordonnées en *X* positives. Par contre le *false northing* est égal à O car l'Équateur reste l'origine de la projection. Comme nous sommes sur l'hémisphère nord, en nous éloignant vers le nord dans cette projection nous restons de toute façon en coordonénes *Y* positives. Par contre, comme montré dans le tableau suivant présentant les paramètres de la projection UTM 32 Sud, nous remarquons un *false northing* défini à 10 000 000 m. En effet, sans ce *false easting* nous aurions à faire à des coordonnées *Y* négatives au fur et à mesure que nous nous éloignons de l'Équateur vers le sud. .. list-table:: Principaux paramètres de la projection UTM 32 Sud :widths: 30 25 :header-rows: 1 * - Paramètre - Valeur * - Ellipsoïde de référence - WGS84 * - Type de projection - Cylindrique * - Zone d'application - De 6°E à 12°E * - Méridien central - 9° Est * - Latitude origine - 0° Nord * - Décalage est (False easting) - 500 000 m * - Décalage nord (False northing) - 10 000 000 m .. warning:: Les images Landsat de l'hémisphère sud ne sont pas fournies en projection zones UTM Sud, contrairement à ce qu'il serait rigoureux de faire. Elles sont `toutes fournies en zones UTM Nord`_, ce qui engendre des coordonnées *Y* négatives. Dans la plupart des cas cela ne pose aucun problème pour l'affichage et les traitements mais peut être à l'origine de surprises si les images sont croisées à des couches en zones UTM Sud. .. _WGS84-monde: Pour le monde *************** Lorsqu'il est nécessaire de représenter des données à l'échelle globale dans leur entièreté, en mode *carte du monde*, le système de coordonnées de référence **WGS84** est couramment utilisé, dont le code EPSG est **4326**. Il ne s'agît **pas** d'une projection, mais d'un système géographique **non projeté**. Dans ce système, les coordonnées ne sont pas en mètres (comme pour les projections vues plus haut) mais sont en degrés de longitude et latitude. En gros, dans ce système nous sommes à la surface de l'ellipsoïde WGS84. Les longitudes et latitudes peuvent être exprimées en degrés, minutes et secondes, mais également en degrés décimaux. Il est facile de passer de l'un à l'autre et des `convertisseurs existent`_. Par définition, dans 1 degré d'angle il y a 60 minutes d'angle et dans 1 minute d'angle il y a 60 secondes d'angle. Ainsi, dans 1 degré il y a 3600 secondes d'angle (:math:`60 \times 60`). Nous pouvons donc facilement convertir la coordonnée suivante en degrés, minutes, secondes : 42°33'56'' par le calcul qui suit : .. math:: 42 + 33 \times \frac{1}{60} + 56 \times \frac{1}{3600} Ce qui nous donne une coordonnée en degré décimal de 42.565555. La plupart des données à l'échelle globale produites par différentes institutions sont fournies dans ce système de coordonnées. .. warning:: Si nous calculons des superficies en restant dans ce SCR, nous obtenons non pas des mètres carrés mais des degrés carrés, ce qui n'est pas évident à se représenter... Si nous souhaitons travailler en mode *projeté* à l'échelle globale, il existe une projection conforme, qui ne respecte donc pas les surfaces (:numref:`monde-wgs84`), utilisée notamment par GoogleMaps et OpenStreetMap. Il s’agit de la projection **WGS84 / Pseudo Mercator**, dont le code EPSG est **3857**. C'est une projection *cylindrique conforme* et *tangente* à l'Équateur. Le point d'origine est le point d'intersection entre l'Équateur et le méridien de Greenwich sans *false eastin* ni *false northing*. Les points se situant dans le quart nord-est de cette projection auront des coordonnées *X* et *Y* positives, mais partout ailleurs, au moins une des deux coordonnées sera négative. Comme sur une projection de Mercartor classique, les régions polaires sont très exagérées, mais leurs contours respectent la réalité. Par contre, si une projection *équivalente* est nécessaire, à l'échelle globale, il est possible d'utiliser la projection **WGS 84 / Equal Earth Greenwich** de code EPSG **8857**. Les continents seront alors de plus en déformés à mesure que nous nous éloignons de l'Équateur, mais les surfaces seront préservées. Il peut être judicieux d'utiliser cette projection pour représenter des phénomènes se rapportant aux surfaces comme des densités de population ou des rendements agricoles. .. figure:: figures/fig_monde_wgs84.png :width: 45em :align: center :alt: Le monde représenté dans trois SCR. :figclass: align-center :name: monde-wgs84 Une représentation des pays du monde dans trois SCR différents : WGS84 non projeté (A), la projection conforme WGS84 / Pseudo Mercator utilisée par GoogleMaps ou OpenStreetMap (B) et la projection équivalente WGS84 / Equal Earth Greenwich (C) à privilégier pour les représentations cartographiques faisant intervenir des surfaces. .. _scr-sig: Implications dans les SIG et en télédétection ------------------------------------------------ .. _calc_sup: Sur le calcul des surfaces *************************** Ces aspects théoriques ont des implications concrètes et directes sur les calculs que nous pouvons faire sur des couches utilisées dans un SIG. L'implication la plus forte a lieu sur le calcul des superficies pour une couche vecteur de type *polygones*. Pour un même polygone, les résultats seront différents selon que la superficie est calculée sur l'ellipsoïde (donc prenant en compte la courbure de la Terre), ou sur la projection. De plus, selon la nature *conforme* ou *équivalente* de la projection le calcul donnera des résultats différents. Le tableau suivant résume la superficie de la France calculée dans différents SCR. .. list-table:: La superficie de la France dans différents SCR :widths: 30 20 25 :header-rows: 1 * - Référentiel - Nature - Superficie (km\ :sup:`2`) * - Sur l'ellipsoïde - - 548 702,564 * - Lambert 93 - Conforme - 548 411,411 * - WGS84/UTM 31N - Conforme - 548 894,166 * - ETRS89 - Équivalente - 548 702,521 * - WGS84 / Pseudo Mercator - Conforme - 1 163 480,151 * - WGS84 / Equal Earth Greenwich - Équivalente - 548 702,581 * - WGS84 - Non projeté - 64 degrés carrés La superficie la plus proche de la réalité est la première du tableau, celle calculée en mètres carrés directement sur l'ellipsoïde, donc en prenant en compte la courbure de la Terre. De façon logique, les superficies calculées sur des projections *équivalentes* conservant les surfaces se rapprochent grandement de cette valeur. Par exemple, les superficies calculées sur les projections *ETRS89* et *WGS84 / Equal Earth Greenwich* (toutes deux équivalentes) donnent des superficies identiques au kilomètre carré près à la superficie calculée sur l'ellipsoïde. Par contre, sous le kilomètre carré, il y a tout de même un effet distorsion inhérent à la projection utilisée. Dans le même esprit, les superficies calculées sur des projections *conformes* ne conservant pas les surfaces donnent des résultats qui s'éloignent de la réalité. Néanmoins, même si la projection Lambert 93 est conforme, l'écart est faible car le Lambert 93 est vraiment fait pour la France. Par contre, dans la projection conforme *Pseudo Mercator*, la superficie de la France est extrêmement exagérée. Dans cette projection, la France mesure près du double de la réalité. Enfin, nous pouvons calculer la superficie de la France en degrés carrés, si nous souhaitons utiliser les unités géographiques (possible mais peu pertinent...) .. _fonctions-sup: Quelle fonction choisir pour calculer des surfaces ? ******************************************************** Pour calculer des superficies dans les logiciels de SIG ou les librairies géospatiales des langages de programmation, nous faisons appel à différentes fonctions. Ces fonctions ont des comportements différents qui peuvent être modifiés en changeant leurs paramètres. Il est nécessaire de savoir ce que fait chaque fonction afin de pouvoir exploiter ses résultats de superficies. Le tableau suivant répertorie le comportement de quelques fonctions couramment utilisées dans différents outils. .. list-table:: Méthode de calcul de superficies selon différents outils :widths: 30 25 25 25 :header-rows: 1 * - Outil - Fonction - Paramètre - Méthode * - QGIS - $area - - Selon l'ellipsoïde * - QGIS - area(@geometry) - - Selon la projection * - R - terra::expanse - transform=TRUE (1) - Selon l'ellipsoïde * - R - terra::expanse - transform=FALSE - Selon la projection (1) Le paramètre *transform* de *terra::expanse* est mis à *TRUE* par défaut. Au final, souvent par défaut dans QGIS nous utilisons la fonction ``$area`` et dans *R* avec *terra* nous utilisons ``expanse``. Dans les deux cas, les superficies sont calculées le long de l'ellipsoïde et ne prennent donc pas en compte la projection. Ainsi, même si nous disposons d'une couche en projection conforme, Mercator par exemple, sur laquelle les entités sont largement déformées en taille, nous obtiendrons des superficies correctes en utilisant ces fonctions par défaut. .. _implications-rasters: Implications sur les rasters ****************************** Même si nous y pensons moins immédiatement, ces questions de projections ont des implications sur les traitements rasters. En effet, tout comme les données vecteurs, les rasters sont géoréférencés dans un système de coordonnées de référence, qui peut être non projeté (comme le *WGS84*) ou projeté dans une projection conforme ou équivalente. Par exemple, le modèle numérique de terrain global *SRTM* est fourni non projeté en *WGS84* et présente donc une résolution spatiale exprimée en degrés. De leurs côtés, les images satellites Landsat ou Sentinel-2 sont fournies projetées dans le système UTM. Or les projections UTM sont des projections conformes, qui ne conservent donc pas les surfaces. Par exemple, les capteurs à bord des satellites Landsat ont une résolution spatiale de 30 m. Ils sont donc capables de distinguer au mieux des portions de surface terrestre de 30 mètres sur 30 mètres. Les images sont fournies à la résolution du capteur, donc sous formes de rasters présentant des pixels de 30 m de côté (:numref:`info-landsat`). .. figure:: figures/fig_info_raster_Landsat_R.png :width: 32em :align: center :alt: Informations spatiales d'un raster Landsat. :figclass: align-center :name: info-landsat Informations spatiales d'un raster de bande spectrale Landsat. Le SCR est *WGS 84 / UTM zone 31N* avec une *résolution* de 30 m sur 30 m. Cependant, comme la projection UTM est *conforme*, elle ne conserve pas les surfaces. Ainsi, les pixels ne représentent pas en réalité 900 m\ :sup:`2` de territoire. Sur la figure suivante (:numref:`mos-landsat`), nous avons une mosaïque de quatre images Landsat (des *paths and row* 197-35 aux 194-35) sur le nord de l'Algérie. Ces quatre rasters s'étalent sur la totalité de la zone UTM 31 Nord. Comme toutes les projections UTM, celle-ci est une projection conforme transverse donc centrée sur un méridien central, qui est le méridien 3°E (en rouge sur la figure :numref:`mos-landsat`). .. figure:: figures/fig_Landsat_UTM31N_SCR.png :width: 45em :align: center :alt: Mosaïque Landsat zone UTM 31N. :figclass: align-center :name: mos-landsat Mosaïque de quatre tuiles Landsat recouvrant la zone UTM 31 N qui a pour méridien central le 3°E (en rouge). Dans cette projection, chaque pixel de cette mosaïque mesure 30 m sur 30 m. Cependant, il est possible de calculer la superficie *réelle* de chaque pixel directement sur l'ellipsoïde. Cette manipulation peut se faire, par exemple, avec la fonction ``cellsize`` de la librairie *terra* de R. En sortie, nous obtenons un raster où la valeur de la superficie de chaque pixel, en mètres carrés, est directement associée à chaque pixel (:numref:`mos-cellsize`). .. figure:: figures/fig_Landsat_UTM31N_SCR_cellsize.png :width: 45em :align: center :alt: Superficie réelle de chaque pixel. :figclass: align-center :name: mos-cellsize Superficie réelle (le long de l'ellipsoïde) de chaque pixel de la mosaïque Landsat. Sur ce raster, les pixels directement sous le méridien 3°E, utilisé comme méridien central de la projection UTM 31N, ont logiquement une superficie très proche de 900 m\ :sup:`2` mais tout de même légèrement supérieure (900.7203 m\ :sup:`2`). Plus nous nous éloignons de ce méridien central, à l'est et à l'ouest, plus la superficie réelle des pixels s'éloigne des 900 m\ :sup:`2`. En bordure de la zone, les pixels ne mesurent plus que 898 m\ :sup:`2`. La différence est faible mais il est intéressant de regarder l'effet sur des surfaces plus grandes qu'un pixel. Par exemple, sur la figure suivante (:numref:`mos-seuil`), nous avons extrait tous les pixels dont la valeur est supérieure à 0.38. Ces pixels sont représentés en rouge. .. figure:: figures/fig_Landsat_UTM31N_SCR_seuil.png :width: 45em :align: center :alt: Pixels supérieurs à 0.38. :figclass: align-center :name: mos-seuil Extraction des pixels dont la valeur est supérieure à 0.38 (en rouge). Nous trouvons 23 386 234 pixels dont la valeur est supérieure à 0.38. Dans la projection du raster, chaque pixel mesure 900 m\ :sup:`2`. L'ensemble des pixels extraits recouvrent donc 21 047.61 km\ :sup:`2` selon la projection. Par contre, ces mêmes pixels, le long de l'ellipsoïde, recouvrent 21 053.28 km\ :sup:`2`, soit une différence de 5.66 km\ :sup:`2`. Cette différence est certes faible mais tout de même éloignée de zéro. Si nous souhaitons que les résolutions des pixels correspondent réellement à leurs superficies, il faut projeter nos rasters dans des projections équivalentes. Par exemple, si nous projetons notre mosaïque dans la projection *WGS84 / Equal Earth Greenwich (EPSG 8857)*, la résolution des pixels est de 29.35782 m de côté. Sur ce raster reprojeté, nous trouvons 23 916 148 pixels ayant une valeur supérieure à 0.38, représentant ainsi une superficie de 20 612,88 km\ :sup:`2` sur la projection. Or, si nous calculons la superficie de ces pixels supérieurs à 0.38 sur l'ellipsoïde nous trouvons également 20 612,88 km\ :sup:`2`. La projection équivalente conserve donc bien la superficie de chaque pixel. .. tip:: Si nous souhaitons calculer une superficie sur un raster en multipliant un nombre de pixels par la résolution d'un pixel, il est plus juste de le faire si notre raster est dans une projection équivalente. .. _georef: Géoréférencement ---------------------- Une des façons de créer des données est de digitaliser manuellement des informations à partir d'une carte ou de tout autre document géoréférencé comme une photographie aérienne. Dans cette section nous allons voir comment passer d'une carte papier scannée à une carte géoréférencée exploitable dans un logiciel SIG. Lorsque nous disposons d'en document scanné, nous disposons d'un fichier raster (une image composée de pixels) mais ne présentant aucune information spatiale. Chaque pixel est caractérisé par une couleur mais pas par un couple de coordonnées. Le processus de *géoréférencement* permet d'attribuer à chaque pixel du document un couple de coordonnées X et Y. A la fin du processus nous obtenons le même raster qu'au départ mais exploitable dans un logiciel de SIG. Le principe est de repérer visuellement sur la carte des points dont les coordonnées sont connues. Ces points seront appelés *points d'ancrage* ou *points d'amer*. Ensuite, manuellement nous indiquons au logiciel pour chacun de ces points les coordonnées correspondantes. Une fois un certain nombre de points ainsi renseignés, le logiciel se chargera de faire une interpolation entre ces points pour attribuer à tout pixel de l'image un couple de coordonnées. En fonction de la qualité du document, de la précision de l'opérateur, et de la précision des coordonnées initiales, le document final sera forcément imparfait et plus ou moins déformé. Mais il sera exploitable. La question centrale est de bien choisir les points d'ancrage. Dans le cas d'un scan d'une carte topographique, il est facile de s'appuyer sur le quadrillage en longitudes et latitudes que ce type de cartes présentent généralement. Il faudra simplement être attentif au système de coordonnées présentées sur la carte. Dans le cas de cartes anciennes ou de photographies aériennes, un tel quadrillage n'est pas indiqué. Tout l'art consistera dans ce cas à repérer des points reconnaissables dont on pourra trouver les coordonnées par ailleurs. Par exemple, si nous reconnaissons une intersection de routes ou un bâtiment emblématique comme un pont ou un clocher, nous pouvons nous en servir comme points d'ancrage. Il suffira de repérer sa position sur une source annexe comme Open Street Map. .. note:: Le processus de géoréférencement sera beaucoup utilisé par certains et jamais par d'autres. Ce processus demande précision, concentration et méticulosité. Géoréférencer un document peut être long et fastidueux. .. _georef-qgis: Géoréférencer un document dans QGIS ************************************* QGIS possède un menu de géoréférencement simple à prendre en mains. Dans cet exemple, nous allons géoréférencer une carte topographique scannée. Il s'agît d'une carte soviétique de la fin des années 1980 relevée pour la région de `Kaliningrad`_. La carte scannée est disponible sur `cette page`_. Le menu de géoréférencement se trouve dans le menu :menuselection:`Raster --> Géoréférencer...` ou le menu :menuselection:`Couche --> Géoréférencer...`, selon la version de QGIS. La fenêtre suivante apparaît (:numref:`georef-vide-qgis`). .. figure:: figures/fen_georef_vide.png :width: 35em :align: center :alt: alternate text :figclass: align-center :name: georef-vide-qgis Menu de géoréférencement. Nous ajoutons le raster à géoréférencer en cliquant sur l'icône *Ouvrir un raster...* |icone-charger-raster|. La carte apparaît alors dans le panneau central. Nous pouvons naviguer et zoomer / dézoomer à l'aide des icônes de navigation |icone-georef-navigation|. Nous pouvons maintenant passer à la définition des points d'ancrage. .. tip:: Combien de points d'ancrage faut-il définir pour obtenir un résultat correct ? La question n'a pas de réponse absolue. Dans le cas d'un document *moderne* comme une carte topographique, un minimum de quatre points (comme les quatre angles par exemple) pourra faire l'affaire. Dans le cas d'une carte ancienne, il en faudra beaucoup plus pour palier le mieux possible aux déformations inhérentes aux processus. Une dizaine sera bien. La difficulté consiste à trouver des points d'ancrage dont on est sûrs. Dans notre cas, nous nous appuierons sur la quadrillage des coordonnées. Premièrement nous relevons le système de coordonnées employée sur cette carte. En zoomant sur un coin, nous nous apercevons qu'il s'agit d'un système global en degrés, minutes et secondes. Même les soviétiques se référaient au méridien de Greenwich, nous sommes donc dans le système global WGS84 (EPSG 4326). Lorsque nous zoomons sur le coin supérieur gauche, nous voyons les coordonnées suivantes (:numref:`georef-coin`). .. figure:: figures/fen_georef_coin_coord.png :width: 20em :align: center :alt: alternate text :figclass: align-center :name: georef-coin Coordonnées du coin supérieur gauche. Le coin supérieur gauche de la carte a ainsi pour coordonnées 55°20'00'' Nord et 20°00'00'' Est. Nous repérons de la même façon les coordonnées des quatre coins (:numref:`georef-points`). .. figure:: figures/fig_georef_points.png :width: 30em :align: center :alt: alternate text :figclass: align-center :name: georef-points Les quatre points d'ancrage. Pour récapituler nous obtenons le tableau suivant. +-------+------------+------------+ | Point | Latitude | Longitude | +=======+============+============+ | P1 | 55°20'00'' | 20°00'00'' | +-------+------------+------------+ | P2 | 55°20'00'' | 21°00'00'' | +-------+------------+------------+ | P3 | 54°40'00'' | 21°00'00'' | +-------+------------+------------+ | P4 | 54°40'00'' | 20°00'00'' | +-------+------------+------------+ Maintenant que nous avons relevé les coordonnées de nos points d'ancrage, nous allons les sélectionner sur notre scan et renseigner leurs coordonnées. Nous commençons par zoomer sur le premier point, nous retombons sur la figure vue précédemment. Nous allons ajouter un point à cette intersection en cliquant sur l'icône *Ajouter un point* |icone-add-point|. Le curseur devient une croix et nous cliquons, le plus précisément possible, sur l'intersection des deux lignes du quadrillage. Une fenêtre apparaît dans laquelle il nous est demandé de renseigner les coordonnées du point sélectionné (:numref:`georef-coord`). .. warning:: Il faut d'abord renseigner la longitude (champ X) puis la latitude (champ Y). Comme indiqué dans le texte explicatif de la fenêtre, les coordonnées doivent être entrées comme suit : degrés minutes secondes (dd mm ss.ss). .. figure:: figures/fen_georef_coord.png :width: 30em :align: center :alt: alternate text :figclass: align-center :name: georef-coord Les coordonnées du premier point d'ancrage. .. tip:: Dans le cas de la définition d'un point d'ancrage par comparaison avec une donnée existante géoréférencée, il est possible d'utiliser l'outil ``Depuis le canevas de la carte`` disponible sur la fenêtre précédente. Il suffirait alors de cliquer sur le point d'ancrage sur le document géoréférencé de comparaison préalablement chargé dans QGIS. Une fois ce premier point défini et enregistré en cliquant sur :guilabel:`OK`, nous répétons l'opération pour les trois points suivants. Une fois la démarche effectuée pour tous les points, le tableau récapitulatif suivant s'affiche en bas de la fenêtre (:numref:`georef-tableau-points`). .. figure:: figures/fig_georef_tableau_points.png :width: 40em :align: center :alt: alternate text :figclass: align-center :name: georef-tableau-points Tableau récapitulatif des points d'ancrage. Ce tableau présente plusieurs informations : * ID : un simple identifiant unique pour chaque point. * Source X : la coordonnée X du point mais dans le référentiel image dont l'origine est le point supérieur gauche. Par exemple, ici le point 0 se trouve à 205.065 pixels à droite du bord gauche de l'image. * Source Y : la coordonnée Y du point dans le référentiel image. Par défaut, les coordonnées Y sont négatives. * Destination X : la coordonnée géographique X du point telle qu'indiquée par l'utilisateur. Elle est exprimée en degrés décimaux. * Destination Y : la coordonnée géographique X du point telle qu'indiquée par l'utilisateur, exprimée également en degrés décimaux. Les trois autres colonnes seront remplies par la suite. Nous allons maintenant passer au processus de géoréférencement proprement dit. Pour cela nous allons dans le menu :menuselection:`Paramètres --> Paramètres de transformation...`. La fenêtre suivante apparaît (:numref:`georef-param-transfo`). .. figure:: figures/fen_georef_param_transfo.png :width: 20em :align: center :alt: alternate text :figclass: align-center :name: georef-param-transfo Paramétrage du géoréférencement. Dans le champ ``Type de transformation`` nous sélectionnons la façon dont l'interpolation sera effectuée. Ici par exemple nous choisissons ``Thin plate spline``. Dans le champ ``Méthode de rééchantillonnage``, nous pouvons choisir une méthode "lissante" comme ``Cubic spline``. Dans le champ ``SCR``, nous spécifions le système de coordonnées des coordonnées que nous avons entrées. Il s'agît ici du WGS 84 (EPSG 4326). Dans le champ ``Raster de sortie``, nous spécifions le chemin et le nom du raster géoréférencé qui sera produit. Si nous le souhaitons, nous pouvons générer un rapport qui fera un état des lieux sur la qualité du géoréférencement. Nous pouvons également cocher la case ``Charger dans QGIS lorsque terminé`` afin d'afficher automatiquement le raster produit. Une fois ces paramètres renseignés et après avoir cliqué sur :guilabel:`OK`, nous revenons à la fenêtre principale du module de géoréférencement. Pour lancer le processus, nous cliquons sur l'icône *Débuter le Géoréférencement* |icone-georef-lancer|. .. tip:: Il est possible de générer la commande GDAL sous jacente à ce processus de géoréférencement dans le menu :menuselection:`Fichier --> Générer un script GDAL`. Il s'agira alors d'une commande au format texte qu'un utilisateur averti pourra réutiliser et modifier. Une fois le processus achevé, le raster géoréférence apparaît dans QGIS. Il a été déformé lors de la manipulation, ce qui est normal. Pour s'assurer du bon déroulement du processus, il est possible d'afficher en fond les images Google Earth par exemple et de mettre en transparence le raster produit (:numref:`georef-resultat`). .. figure:: figures/fig_georef_resultat.png :width: 35em :align: center :alt: alternate text :figclass: align-center :name: georef-resultat Superposition de la carte géoréférencée et du fond Google Earth. Visuellement le résultat semble tout à fait satisfaisant. Maintenant, nous pouvons nous intéresser à une question scientifique du type *Le trait de côte de cette région a-t-il évolué depuis la fin des années 1980 ?* ou *Le réseau routier s'est-il densifié ici où là ?* ... Nous pouvons quantifier la qualité du processus en regardant les résultats quantitatifs se trouvant dans le rapport PDF que nous avons sorti. Nous y retrouvons le tableau présenté plus haut dans lequel les colonnes de résidus ont été mises à jour. Dans l'idéal, les résidus devraient être égaux à 0. Dans notre cas ils sont très faible, nous pouvons être satisfaits des résultats. .. warning:: Nous avons vu qu'il est nécessaire de définir un type de transformation et une méthode de rééchantillonnage. Nous pouvons nous poser la question légitime de quelles méthodes choisir. Il n'y a pas de réponse absolue, le choix effectué ici donne empiriquement de bons résultats, mais c'est à l'utilisateur de tester différentes méthodes si les résultats ne sont pas satisfaisants. .. _pratique-reprojection: Changement de Système de Coordonnées de Référence -------------------------------------------------- Par *changement de système de coordonnées* nous entendons le fait de changer le système de coordonnées de références (SCR) associé à un fichier. Nous le faisons passer d'un système de coordonnées *A* vers un système de coordonnées *B*. Ce procédé est souvent appelé *reprojection*. C'est une manipulation de base car il est généralement conseillé de ne travailler que dans un seul système de coordonnées au sein d'un projet géomatique. Cette manipulation n'est cependant pas anodine car elle va forcément déformer la donnée initiale. Il est évidemment nécessaire de choisir un système pertinent vis-à-vis de sa zone d'étude mais aussi selon les données mobilisées. Par exemple, si vous travaillez sur une région de France métropolitaine et que vous manipulez des données vecteurs "officielle", comme un découpage administratif fourni par l'IGN par exemple, travailler en Lambert 93 serait pertinent. Par contre, si vous mobilisez également des données Landsat dans votre projet, vous aurez ces images satellites en UTM. Dans ce cas, même si d'un point de vue géographique c'est moins précis, il sera préférable de reprojeter vos données vecteurs Lambert 93 vers l'UTM. En effet, reprojeter des images Landsat en Lambert 93 va s'avérer fastidieux, gourmand en ressources et va déformer les pixels de ces images. Ainsi il sera judicieux de faire passer toutes les données du projet en UTM. Concrètement, lors de ce processus, dans la plupart des cas, une nouvelle couche sera créée. En effet, il est rare de changer directement le système de coordonnées d'une couche donnée. Il s'agît généralement d'une copie de la couche originale avec le nouveau système de coordonnées. La déformation engendrée par la reprojection, dans le cas des données rasters, implique une certaine interpolation dans le processus. Nous verrons qu'il est possible, ou pas, de régler la qualité de cette interpolation. Dans cet exemple, nous allons changer le système de coordonnées d'un MNT (SRTM) dont le SCR est *EPSG:32632 - WGS 84 / UTM zone 32N - Projeté* vers le SCR officiel français, le Lambert 93 (*EPSG:2154 - RGF93 / Lambert-93 - Projeté*). .. _reprojection-raster: Reprojection d'un raster ******************************* Nous allons voir ici différents outils qui permettent de reprojeter un raster. .. _reprojection-raster-qgis: Reprojection d'un raster dans QGIS ++++++++++++++++++++++++++++++++++++++ Version de QGIS : 3.16.1 Dans QGIS, il existe deux méthodes différentes pour effectuer ce changement de SCR de couche raster. La première méthode fait appel à un module propre à QGIS, qui nécessite très peu de paramétrages. La seconde méthode utilise également un module de base de QGIS mais qui interface en fait une fonctionnalité de `GDAL`_. Cette seconde méthode présente plus d'options et permet de régler plus finement la qualité de la reprojection. **Avec QGIS seul** Une fois notre raster chargé dans QGIS, nous faisons un clic droit sur cette couche dans le panneau des couches. Nous cliquons ensuite sur :menuselection:`Exporter --> Enregistrer sous...` . La fenêtre suivante apparaît (:numref:`repro-raster-qgis`). .. figure:: figures/fen_repro_raster_qgis.png :width: 35em :align: center :alt: alternate text :figclass: align-center :name: repro-raster-qgis Changer le SCR d'un raster avec QGIS seul. Dans le champ ``Nom de fichier``, nous spécifions le chemin et le nom sous lequel nous souhaitons enregistrer le nouveau raster reprojeté. Dans le champ ``SCR``, nous indiquons le SCR que nous souhaitons en résultat. Il est possible de cliquer sur l'icône |icone_choix_SCR| pour ouvrir une nouvelle fenêtre qui permet de sélectionner le système désiré (:numref:`choix-scr2`). .. figure:: figures/fen_choix_scr.png :width: 35em :align: center :alt: alternate text :figclass: align-center :name: choix-scr2 Choix d'un système de coordonnées de références. Dans cette fenêtre de sélection de SCR, le plus simple est de chercher le système de coordonnées souhaité en cherchant son code EPSG dans le champ ``Filtre``. Le code EPSG correspondant au Lambert 93 est *2154*. Nous entrons donc ce code. Dans le panneau ``Systèmes de Coordonnées de Références Prédéfinis``, les systèmes de coordonnées contenant ce code apparaissent. Il suffit de cliquer sur celui qui nous intéresse. Notons que les SCR que nous utilisons le plus apparaissent dans le panneau ``Systèmes de Coordonnées de Références récemment utilisés``. En bas à droite de la fenêtre, QGIS nous propose un aperçu du territoire sur lequel ce SCR est valable. Il s'agît bien de la France métropolitaine dans notre cas. Il suffit ensuite de cliquer sur :guilabel:`OK`. Dans la fenêtre précédente, le champ ``SCR`` se met bien à jour. Nous cliquons sur :guilabel:`OK` et la nouvelle couche reprojetée apparaît automatiquement dans QGIS. Si nous regardons ses propriétés, nous constatons que son système de coordonnées est bien le Lambert 93. Nous pouvons constater que cette méthode ne permet pas de choisir et de régler l'interpolation sous-jacente. Il semblerait qu'elle utilise par défaut une méthode de type *plus proche voisin*. Ce type d'interpolation peut avoir des répercussions sur la qualité du raster résultat. **Avec GDAL interfacé par QGIS** Il est également possible d'utiliser GDAL interfacé par QGIS pour changer le SCR d'un raster. Cette méthode nous permet de régler l'interpolation sous-jacente, ce qui peut s'avérer précieux dans certains cas. Ce module se trouve dans le menu :menuselection:`Raster --> Projections --> Projection (warp)...`. La fenêtre suivante s'affiche (:numref:`repro_raster_gdal`). .. figure:: figures/fen_repro_raster_gdal.png :width: 35em :align: center :alt: alternate text :figclass: align-center :name: repro_raster_gdal Changer le SCR d'un raster via GDAL interfacé par QGIS. Dans le champ ``Couche source``, nous spécifions la couche dont nous voulons changer le SCR. À priori, QGIS trouve tout seul dans quel SCR se trouve la couche à reprojeter, par conséquent le champ ``SCR d'origine`` peut rester vide. Par contre, il est nécessaire de préciser le SCR souhaité. De la même façon que précédemment, il est possible de sélectionner le système Lambert 93 - 2154. Un champ important se trouve juste au-dessous : ``Méthode de ré-échantillonage à utiliser``. Par défaut, la méthode choisie est *Plus Proche Voisin*, ce qui revient à celle utilisée par défaut dans le module QGIS dédié vu précédemment. Ici, nous choisissons *Bilinéaire*. Cette méthode permet de "lisser" la repojection. Enfin, dans le champ ``Reprojeté``, nous spécifions un chemin et un nom pour le raster reprojeté résultat. Les autres champs peuvent garder leurs valeurs par défaut. Notons que le panneau du bas ``Console GDAL/OGR`` contient la commande GDAL qui sera utilisée au final. Cette commande prend en argument les différents paramètres que nous avons ajustés dans les différents champs. Après avoir cliqué sur :guilabel:`Exécuter`, le nouveau raster s'affiche dans sa nouvelle projection. .. warning:: **Finalement, quelle méthode dois-je utiliser ?** Au final, le choix ne doit pas réellement se faire entre utiliser le module interne de QGIS ou le module GDAL. La vraie question est "*quelle méthode de ré-échantillonage dois-je utiliser ?*". Si une méthode de type *Plus proche voisin est adaptée*, vous utilisez le module de QGIS ou bien le module GDAL avec l'option "*Plus proche voisin*". Si vous devez utiliser une méthode qui "lisse" le résultat, alors vous devez nécessairement utiliser le module GDAL avec une option de type *Bilinéaire*, *Cubique* ou *Cubic spline*. Ces trois méthodes donneront des résultats très semblables. La question sous-jacente devient donc *quand dois-je utiliser une méthode Plus proche voisin ou une méthode lissante ?* D'une manière générale, lorsque que vous disposez d'un raster discret, comme un raster d'usage du sol par exemple, il faut utiliser une méthode *Plus proche voisin*. Si nous utilisons une méthode lissante dans ce cas, nous risquons de créer de nouvelles valeurs absurdes qui ne correspondent à aucune classe. Imaginons un raster d'usage du sol en trois classes comme suit : 1 : eau, 2 : bâti et 3 : végétation. Après une reprojection lissante nous risquons de nous retrouver avec quelques pixels de valeurs 2.3, 2.7 ... entre les aplats de classes 2 et 3. Les méthodes lissantes sont à réserver aux rasters continus, comme les MNT ou les images satellites par exemple. Pour ce type de rasters, il est également possible d'utiliser les méthodes de type *Plus proche voisin*, mais certains artefacts peuvent apparaître, notamment sur les MNT. .. _reprojection-raster-qgis-lot: Reprojection par lot de rasters dans QGIS ++++++++++++++++++++++++++++++++++++++++++++ Dans certains cas il est intéressant de reprojeter plusieurs rasters par lot. C'est par exemple le cas lorsque nous manipulons des images satellites multi-spectrales où nous avons un raster par bande spectrale. Au lieu de reprojeter bande spectrale par bande spectrale, QGIS propose un moyen simple pour reprojeter tous les rasters des bandes spectrales en une fois. Pour cela, le plus simple est de charger dans QGIS les rasters à reprojeter (mais ce n'est pas obligatoire). Ici, nous allons reprojeter les trois rasters correspondant aux trois premières bandes spectrales d'une scène Landsat 8. Une fois les rasters chargés, nous ouvrons le menu :menuselection:`Raster --> Projections --> Projection (warp)...`. Nous retombons sur la fenêtre de la figure (:numref:`repro_raster_gdal`). Au lieu de remplir les champs à ce niveau, nous cliquons sur le menu :guilabel:`Exécuter comme processus de lot...` La fenêtre suivante apparaît (:numref:`repro-raster-lot-init-qgis`). .. figure:: figures/fen_repro_raster_lot_init_qgis.png :width: 55em :align: center :alt: alternate text :figclass: align-center :name: repro-raster-lot-init-qgis Reprojeter des rasters par lot, réglages. Dans cette fenêtre, nous inclurons autant de lignes que de rasters à reprojeter. Par défaut, une ligne est affichée. Dans cette ligne nous allons régler la reprojection de notre premier raster. Dans la colonne *Couche source* nous sélectionnons dans le menu déroulant le premier raster à reprojeter. Notons que nous pouvons également pointer vers un raster non ouvert dans QGIS grâce à l'icône |icone_browse|. Dans la colonne *SCR cible* nous spécifions le nouveau SCR que nous souhaitons. Ici nous sélectionnons le Lambert 93 (EPSG 2154). Dans la colonne *Méthode de ré-échantillonnage* nous sélectionnons la méthode souhaiteé, *Bilinéaire* par exemple. Enfin, dans le champ *Reprojeté* nous spécifions le chemin et le nom du raster reprojeté qui sera créé. Il est maintenant nécessaire d'ajouter autant de lignes que de rasters à reprojeter. Pour cela soit nous cliquons sur |icone_plus| pour ajouter des lignes une à une, soit nous cliquons sur la case *Auto-remplissage...* qui se trouve immédiatement sous le titre de la colonne *Couche source*. Dans le menu déroulant qui apparaît nous pouvons *Sélectionner des fichiers* depuis un répertoire, *Ajouter tous les fichiers d'un répertoire*, ou *Sélectionner à partir des couches chargées*, ce que nous choisissons ici. Un menu nous permet alors de sélectionner les rasters chargés que nous souhaitons reprojeter. La fenêtre précédente se met à jour (:numref:`repro-raster-lot-qgis`). .. figure:: figures/fen_repro_raster_lot_qgis.png :width: 55em :align: center :alt: alternate text :figclass: align-center :name: repro-raster-lot-qgis Reprojeter des rasters par lot dans QGIS. Il suffit alors de régler les colonnes pour les différents rasters et de remplir la colonne *Reprojeté* pour chaque raster. Il est possible de faire des copier-coller (un par un) du chemin de la première ligne vers les autres lignes et de changer juste le nom du fichier. Ces clics et ces copier-coller peuvent s'avérer un peu fastidieux mais toujours moins que de traiter les rasters un par un. En cliquant sur :guilabel:`Exécuter`, les rasters se reprojettent par lot. .. _reprojection-raster-scp: Preprojection d'un raster avec le module SCP +++++++++++++++++++++++++++++++++++++++++++++++++++ Version de QGIS : 3.18.2 Version de SCP : 7.8.16 Le module supplémentaire à QGIS nommé *Semi-Automaic Classification Plugin* et connu sous l'acronyme de *SCP* a été conçu pour faire des classifications d'images satellites. Mais ce plugin propose également de nombreux outils connexes, notamment un pour la reprojection de couches rasters. Dans cet exemple, nous allons utiliser SCP pour reprojeter des bandes spectrales Landsat 8 prises au-dessus du Val d'Oise depuis le SCR *EPSG:32631 - WGS 84 / UTM zone 31N - Projeté* vers le SCR *Lambert 93* (EPSG 2154). Pour utiliser les fonctionnalités de reprojection de SCP, il est nécessaire d'avoir au préalable défini un *Jeu de bandes* SCP (voir la section dédiée à la classification supervisée avec SCP). Ainsi, par défaut SCP fait des reprojections par lot. Le menu dédié se trouve dans :menuselection:`SCP --> Pré-traitement --> Reproject raster bands`. Le menu suivant apparaît (:numref:`repro_scp`). .. figure:: figures/fen_scp_reproject.png :width: 40em :align: center :alt: alternate text :figclass: align-center :name: repro_scp Changer le SCR de couches rasters avec SCP. Dans le champ ``Sélectionner un jeu de bandes``, nous indiquons le jeu de bandes à reprojeter. Dans le champ ``Use EPSG code``, nous indiquons le SCR vers lequel nous souhaitons convertir nos rasters, ici le Lambert 93 (EPSG 2154). Enfin, dans le champ ``Resampling method`` nous pouvons choisir une méthode de ré-échantillonnage. Il semblerait que SCP utilise ici GDAL, mais nous ne savons pas avec quels paramètres. Pour finir, nous pouvons indiquer un ``Préfixe de sortie`` puis cliquer sur :guilabel:`Lancer`. Il sera alors nécessaire de spécifier un chemin où sauver les rasters reprojetés. Pour le choix de la méthode de ré-échantillonnage, reportez vous à la discussion plus haut sur cette page. Quoi qu'il en soit, ce choix a des répercussions sur le raster reprojété (:numref:`repro_scp_compa`).. .. figure:: figures/fig_scp_repro_compa.png :width: 50em :align: center :alt: alternate text :figclass: align-center :name: repro_scp_compa Effet de la méthode de ré-échantillonnage choisie sur un zoom de la bande 5 d'une image Landsat 8. Les styles sont les mêmes sur les six vignettes. La méthode *Nearest neighbour* semble créer le moins d'artefacts. Les méthodes *Average*, *Maximum* et *Median* apportent un flou et le *First quartile* fait "baver" les pixels. En plus de la reprojection, ce menu propose deux outils très intéressants : **ré-échantillonnage spatial** et **changement d'encodage**. Comme nous voyons sur le menu (:numref:`repro_scp`), il est possible de changer les résolutions en X (``X resolution``) et en Y (``Y resolution``). Nous pouvons même nous servir de cette option sans changer le SCR. La figure (:numref:`scp_resamp`) présente un ré-échantillonnage spatial à 200 mètres en X et Y d'une bande Landsat. .. figure:: figures/fig_scp_resamp.png :width: 40em :align: center :alt: alternate text :figclass: align-center :name: scp_resamp Ré-échantillonnage à 200 mètres d'une bande spectrale Landsat 8. L'option de changement de **l'encodage** est également très intéressante. Cette option peut également s'utiliser sans reprojection. Avec cette fonctionnalité, il est possible de, par exemple, changer un raster codé en *UInt16 - nombre entier non signé de seize bits* vers un raster codé en *Float32 - nombre à virgule flottante de 32 bits*. Ces conversions peuvent s'avérer nécessaires dans certains traitements, attention toutefois à l'éventuelle altération de la donnée que ça peut engendrer. .. note:: Au final, les fonctionnalités de reprojection de SCP sont similaires à ce que propose GDAL mais sont intéressantes lorsque nous travaillons sur un jeu de bandes SCP. Par contre, les options de ré-échantillonnage spatial et de transformation de l'encodage sont très intéressantes même si nous ne nous servons d'aucune autre fonctionnalité de SCP. .. _repro-raster-GDAL: Reprojection d'un raster avec GDAL ++++++++++++++++++++++++++++++++++++++ Version de GDAL : 3.0.4 Une fois que GDAL (:ref:`logiciels-GDAL`) est correctement configuré dans le *PATH* du système, il est possible de l'utiliser en ligne de commande pour reprojeter un raster. Le menu à utiliser est *gdalwarp*. La commande présentée ici est la commune minimale, suffisante dans la plupart des cas, mais elle peut être étoffée avec des options comme présentées sur la `documentation de gdalwarp`_. Dans l'exemple, nous allons reprojeter une bande spectrale Landsat de l'UTM vers le Lambert 93 dont le code EPSG est 2154. Il est nécessaire de se placer dans le répertoire contenant le raster à reprojeter, ou bien de bien renseigner les chemins en relatif ou en absolu. .. code-block:: sh gdalwarp -t_srs EPSG:2154 -r bilinear -of GTiff LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF B5_gdal_L93.tif avec : * gdalwarp : la commande GDAL pour la reprojection raster * -t_srs : le mot clé pour définir le SCR désiré en sortie (*t* comme *target*) * -r : le mot clé pour choisir la méthode de rééchantiollonnage * -of : le mot clé pour le format du raster en sortie * LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF : le raster à reprojeter * B5_gdal_L93.tif : le raster issu de la reprojection .. _reprojection-raster-R: Reprojection d'un raster avec R ++++++++++++++++++++++++++++++++++++++ Version de R : 4.1.2 Lorsqu'un raster est chargé dans R avec l'un ou l'autre des packages dédiés (:ref:`import-raster-R`), il est possible de reprojeter facilement ce raster. Ici, nous verrons comment reprojeter un raster de type *SpatRaster* chargé avec *terra*. Dans l'exemple ci-arpès, nous allons changer de SCR du MNT du `bassin de la Roya`_ issu du SRTM et nommé *srtm_roya_L93.tif*. Nous le reprojèterons du *Lambert 93* vers le système *UTM zone 31N* dont le code EPSG est 32631. Le tout est résumé dans les lignes suivantes. .. code-block:: R # chargement du raster à reprojeter mnt_roya <- terra::rast('srtm_roya_L93.tif') # reprojection du raster dans le système 32631 mnt_roya_reproj <- terra::project(mnt_roya, 'EPSG:32631') Si nous affichons les propriétés du raster (en entrant simplement son nom dans la console R par exemple), nous constatons bien que le SCR a été modifié. .. _reprojection-raster-Python: Reprojection d'un raster avec Python ++++++++++++++++++++++++++++++++++++++ Version de Python : 3.8.10 Il est tout à fait possible de changer le SCR d'un raster importé dans un script Python. **Avec rasterio** Nous présenterons ici la manipulation avec la librairie *rasterio*. Il faut bien avouer que ça fait pas mal de lignes de code pour une petite manipulation... Si la reprojection est effectuée pour mettre en concordance un raster avec un vecteur, il peut être plus simple de reprojeter le vecteur plutôt que le raster. .. code-block:: Python # chargement de la librairie de gestion des rasters rasterio import rasterio # chargement des modules de rasterio spécifiques à la reprojection from rasterio.warp import calculate_default_transform, reproject, Resampling # choix du SCR de destination en code EPSG dst_crs = 'EPSG:2154' # chargement du raster à reprojeter et on le stocke dans la variable src # on récupère ses paramètres "transform" et "width" # on calcule la hauteur que prendra le raster reprojeté # on stocke les paramètres du raster dans la variable kwargs # on met à jour ces paramètres pour le raster de destination with rasterio.open('./Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF') as src: transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs = src.meta.copy() kwargs.update({ 'crs': dst_crs, 'transform': transform, 'width': width, 'height': height }) # on créé le raster de destination en lui attribuant les paramètres mis à jour précédemment with rasterio.open('./Landsat_13/B5_L93.tif', 'w', **kwargs) as dst: for i in range(1, src.count + 1): reproject( source=rasterio.band(src, i), destination=rasterio.band(dst, i), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=dst_crs, resampling=Resampling.nearest) Le script présenté ci-dessous est également utilisable pour reprojeter un raster multi-bandes. .. tip:: Il peut être plus simple d'utiliser une simple commande GDAL que nous pouvons exécuter dans un script Python. **Avec une commande GDAL** GDAL (:ref:`logiciels-GDAL`) propose une commande très simple pour reprojeter un fichier raster via l'outil *gdalwarp* (:ref:`repro-raster-GDAL`). L'idée est de construire une chaîne de caractères qui sera la commande GDAL puis d'exécuter cette commande en tant que commande externe (:ref:`cmd_gdal_python`). Ça suppose bien sûr que l'exécutable de GDAL soit bien reconnu dans le *PATH* du système. La commande présentée ici est la commune minimale, suffisante dans la plupart des cas, mais elle peut être étoffée avec des options comme présentées sur la `documentation de gdalwarp`_. .. code-block:: Python # librairie pour exécuter des commandes système import os # le chemin vers le fichier à reprojeter src = './Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF' # le nom du fichier reprojeté et son chemin de stockage dst = './Landsat_13/B5_gdal_python_L93.tif' # le SCR de reprojection (en EPSG) scr_repro = '2154' # la méthode de rééchantillonnage ech = 'bilinear' # on construit la commande gdalwarp par concaténation cmd = 'gdalwarp -t_srs EPSG:' + scr_repro + ' -r ' + ech + ' -of GTiff ' + src + ' ' + dst os.system(cmd) Cette solution est moins "pythonesque" mais plus rapide à coder et à lire. Cette commande fonctionne aussi pour un raster multi-bandes. .. _reprojection-vecteur: Reprojection d'une couche vecteur ************************************ Nous allons voir ici différents outils permettant de changer le SCR d'une couche vecteur. .. _reprojection-vecteur-qgis: Reprojection d'un vecteur dans QGIS ++++++++++++++++++++++++++++++++++++++ Version de QGIS : 3.16.1 Il est très simple de changer le SCR d'une couche vecteur dans QGIS. Que cette couche soit de type ponctuelle, linéaire ou surfacique la manipulation est la même. Comme pour le cas des rasters, la couche vecteur résultante sera déformée, il est donc nécessaire de choisir un SCR pertinent. Dans cet exemple, nous allons changer le SCR du réseau hydrographique de la Vésubie ("*vesubie_hydro.gpkg*") du Lambert 93 (EPSG 2154) vers le WGS84 UTM Zone 32 Nord (EPSG 32632). Après avoir chargé la couche dans QGIS, nous faisons un clic droit sur cette couche dans le panneau des couches. Puis nous allons dans :menuselection:`Exporter --> Sauvegarder les entités sous...`. La fenêtre suivante s'affiche (:numref:`repro_vecteur`). .. figure:: figures/fen_repro_vecteur_qgis.png :width: 30em :align: center :alt: alternate text :figclass: align-center :name: repro_vecteur Changer le SCR d'une couche vecteur avec QGIS. Dans le champ :guilabel:`Nom de fichier` nous spécifions un chemin d'export et un nom pour la couche reprojetée, ici *vesubie_hydro_UTM32.gpkg*. Dans le champ ``SCR`` nous choisissons le système de coordonnées de la couche résultante, ici l'UTM 32 Nord. Il est possible de sélectionner ce système de coordonnées en cliquant sur l'icône |icone_choix_SCR|. Dans la fenêtre qui apparaît alors il est possible de sélectionner le SCR souhaité en filtrant les résultats par le code EPSG 32632. Ensuite, en cliquant sur :guilabel:`OK` le processus se lance et la nouvelle couche reprojetée apparaît automatiquement. En allant voir ses propriétés il est possible de vérifier que le changement de SCR a bien été pris en compte. .. _reprojection-vecteur-gdal: Reprojection d'un vecteur avec GDAL/OGR ++++++++++++++++++++++++++++++++++++++++ Version de GDAL : 3.0.4 Il est possible de reprojeter une couche vecteur avec GDAL/OGR. Ici, nous allons reprojeter une couche des communes d'Île-de-France (*communes_IDF.gpkg*) en Lambert 93 vers du WGS84 UTM Zone 31 Nord (EPSG 32631). Cette manipulation se fait avec une commande ``ogr2ogr``, comme montré ci-après. .. code-block:: sh ogr2ogr -t_srs epsg:32631 communes_IDF_32631.gpkg communes_IDF.gpkg Avec les éléments suivants : * *ogr2ogr* : le penchant vecteur de GDAL * *-t_srs epsg:32631* : le mot clé pour définir le SCR de destination * *communes_IDF_32631.gpkg communes_IDF.gpkg* : bien mettre en premier la couche de destination et en second la couche à reprojeter .. _reprojection-vecteur-R: Reprojection d'un vecteur dans R ++++++++++++++++++++++++++++++++++ Reprojeter un vecteur dans R dépend si vous travaillez avec la librairie *sf* ou *terra*. Nous présentons ici les deux cas de figure. .. _reprojection-vecteur-terra: Reprojection d'un vecteur avec terra ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Version de R : 4.1.2 Version de terra : 1.7.29 Si nous avons chargé une couche vecteur avec *terra* nous avons un objet *SpatVector*. Il est facile de changer son SCR en lui spécifiant un nouveau SCR via son code EPSG. .. code-block:: R # chargement de la couche vecteur dep <- terra::vect("departements_france_L93.gpkg") # reprojection en UTM Zone 31N (EPSG 32631) dep_reporj <- terra::project(dep, 'EPSG:32631') .. _reprojection-vecteur-sf: Reprojection d'un vecteur avec sf ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Version de R : 4.1.2 Version de sf : 1.0-12 Une fois un vecteur chargé dans R via le package *sf* par exemple, il est possible de le reprojeter vers un nouveau système de coordonnées de référence. Il suffit d'utiliser la commande *st_transform()* qui prend en argument la couche à reprojeter ainsi que le SCR souhaité. Ici, nous allons reprojeter une couche depuis le Lamber 93 en UTM Zone 31 N (EPSG 32631). La commande est la suivante. .. code-block:: R # chargement de la couche vecteur dep <- sf::st_read("departements_france_L93.gpkg") # reprojection en UTM Zone 31N dep_reproj <- sf::st_transform(dep, crs = 32631) La couche résultat est bien maintenant en UTM Zone 31 N. .. |icone_choix_SCR| image:: figures/icone_choix_SCR.png :width: 25 px .. |icone-georef-lancer| image:: figures/icone_georef_lancer.png :width: 15 px .. |icone-add-point| image:: figures/icone_georef_add_point.png :width: 25 px .. |icone-charger-raster| image:: figures/icone_charger_raster_georef.png :width: 25 px .. |icone-georef-navigation| image:: figures/icone_georef_navigation.png :width: 7 em .. |icone_browse| image:: figures/icone_browse.png :width: 20 px .. |icone_plus| image:: figures/icone_export_3D_add_frame.png :width: 15 px .. _documentation de gdalwarp: https://gdal.org/programs/gdalwarp.html .. _Kaliningrad: https://www.openstreetmap.org/relation/1674442#map=12/54.7056/20.4716 .. _cette page: https://maps.vlasenko.net/smtm200/n-34-09.jpg .. _GDAL: https://gdal.org/ .. _European Petroleum Survey Group: https://en.wikipedia.org/wiki/EPSG_Geodetic_Parameter_Dataset .. _Universal Transverse Mercator: https://fr.wikipedia.org/wiki/Transverse_universelle_de_Mercator .. _Corine Land Cover: https://land.copernicus.eu/pan-european/corine-land-cover .. _bassin de la Roya: https://fr.wikipedia.org/wiki/Roya .. _comme la mission Grace: https://earth.esa.int/eogateway/missions/grace .. _-106 m de sa hauteur moyenne: https://www.scientificamerican.com/article/giant-gravity-hole-in-the-ocean-may-be-the-ghost-of-an-ancient-sea1/ .. _International Meridian Conference de 1884: https://www.history.navy.mil/our-collections/photography/numerical-list-of-images/nhhc-series/nh-series/NH-96000/NH-96688.html .. _Place de l'Étoile Rouge à Cotonou: https://fr.wikipedia.org/wiki/Place_de_l%27%C3%89toile_rouge .. _Union Soviétique: https://i.pinimg.com/originals/ac/6f/62/ac6f625116a3b7385723897142b968db.jpg .. _sextant: https://www.maxisciences.com/sciences/sextant-decouvrez-en-plus-sur-cette-invention-qui-a-revolutionne-lhistoire-de-la-navigation_art50732.html .. _chronomètre de marine: https://www.hautehorlogerie.org/fr/watches-and-culture/connaissances-horlogeres/encyclopedie/chronometre-de-marine .. _baromètre: https://www.montre-cardio-gps.fr/altimetre-barometrique-ou-gps-les-differences/ .. _la Chine: https://www.youtube.com/watch?v=LxMtasUnYY8 .. _3m/km aux extrémités du domaine d'application: http://www.geomag.fr/sites/default/files/68_91.pdf .. _indiqué par l'IGN: https://geodesie.ign.fr/contenu/fichiers/Lambert93_ConiquesConformes.pdf .. _bâtiment de la Halle aux Farines de l'Université Paris Cité: https://www.geoportail.gouv.fr/carte?c=2.381454224423593,48.82939927639973&z=19&l0=ORTHOIMAGERY.ORTHOPHOTOS::GEOPORTAIL:OGC:WMTS(1)&permalink=yes .. _Hobart: https://www.openstreetmap.org/search?query=hobart+australie#map=7/-41.968/147.711 .. _toutes fournies en zones UTM Nord: https://www.usgs.gov/faqs/why-do-landsat-scenes-southern-hemisphere-display-negative-utm-values .. _convertisseurs existent: https://www.rapidtables.org/fr/convert/number/degrees-minutes-seconds-to-degrees.html .. _très bon livre de Dava Sobel: https://www.babelio.com/livres/Sobel-Longitude/2797s