Auteur : Paul Passy

Licence :

Les Systèmes de Coordonnées de Références

Dans cette section nous allons découvrir quelques points pratiques sur les systèmes de coordonnées de références (SCR), aussi appelés abusivement projections. La théorie se sera présentée rapidement. De plus amples explications théoriques sont disponibles dans vos cours magistraux ou sur des sites très bien faits comme celui-ci. Nous verrons ensuite quelques manipulations pratiques faisant intervenir les SCR sur des données vecteur et raster.

Table des matières

Un (tout petit) peu de théorie 

La planète Terre est connue pour être une sphère. Dans la réalité c’est un géoïde. Un géoïde est une sphère légèrement bosselée, comme si elle s’était pris des coups de marteaux. En plus d’être bosselé, ce géoïde est légèrement aplati au niveau des pôles. Un rayon allant du centre de la Terre à un pôle est à peu près 20 km plus court qu’un rayon allant du centre la Terre à l’Équateur.

En géomatique, les données que nous utilisons doivent être géoréférencées. C’est à dire qu’à tout point de notre donnée doit être associé un couple de coordonnées X et Y permettant de localiser ce point à la surface du globe. Le « souci » c’est que nous travaillons sur des surfaces planes, soit du papier soit des écrans d’ordinateur. Il est donc nécessaire de passer de la surface géodésique à une surface plate.

Ce processus se fait en deux temps. En premier lieu, le géoïde est approximé par un ellipsoïde. Un ellipsoïde est à l’ellipse ce que la sphère est au cercle. Cet ellipsoïde est un volume ellispsoïdal qui se rapproche au mieux du géoïde. Il existe une multitude d’ellipsoïdes possibles. À chaque ellipsoïde est associé un système géodésique, aussi appelé datum, qui définit les paramètres de l’ellipsoïde comme son aplatissement, sa forme et son orientation.

En second lieu, il est nécessaire de passer de cette surface ellipsoïdale à une surface plane. Ce passage se fait par une projection. Là aussi, il existe une infinité de possibilités de projections différentes. Il est mathématiquement impossible de parfaitement projeté un ellipsoïde sur un plan. Toute projection se rapprochera donc au mieux de la réalité mais conservera plus ou moins une part d’erreur.

Certaines projections conservent les angles, c’est-à-dire les formes, mais distordent les surfaces. Ces projections sont qualifiées de conformes. D’autres, au contraire, conservent les surfaces mais pas les angles. Elles ont tendance à déformer les objets. Ces projections sont qualifiées d”équivalentes (Fig. 139).

alternate text — Fig. 139 Les pays du monde en projection conforme WGS 84 (4326) (A) et en projection équivalente WGS 84 / Equal Earth Greenwich (8857) (B).

Les SCR en pratique 

Pour s’y retrouver dans tous ces systèmes de coordonnées, un code unique a été associée à chacune d’entre elles. Ce code est connu sous l’acronyme EPGS (European Petroleum Survey Group). Les codes EPSG des SCR sont utilisés par la grande majorité des logiciels de géomatique ainsi que par les bibliothèques géospatiales des langages de programmation comme R ou Python.

Dans cette partie pratique nous verrons quelques SCR couramment utilisés et leurs principales caractéristiques.

Les SCR globaux 

Le SCR le plus utilisé pour travailler à l’échelle globale, c’est-à-dire à l’échelle du globe est le WGS 84. Le WGS 84 repose sur l’ellipsoïde IAG GRS80. Ce SCR n’est pas une projection, car son unité est le degré. Dans ce SCR, chaque point est repéré par un couple de coordonnées correspondant à la longitude et à la latitude exprimées en degrés, minutes et secondes ou en degrés décimaux. Par conséquent, il faut être prudent lorsque nous calculons des surfaces sur des données exprimées dans ce SCR car le résultat sera exprimé en degrés carrés, ce qui est compliqué à interpréter.

C’est dans ce SCR que nous trouverons la plupart des données diffusées à l’échelle du globe. C’est notamment le cas des données diffusées par le site Natural Earth.

Le code EPSG associé à ce SCR est le 4326.

Il existe d’autres SCR globaux, notamment ceux utilisés par les grands diffuseurs de données en ligne comme Google Maps ou OpenStreetMap (EPSG 3857). Ce SCR est également un SCR global de type Mercator mais projeté. C’est-à-dire que son unité est le mètre.

Astuce

Moins utilisés en géomatique, il existe néanmoins des SCR globaux à projections équivalentes qui respectent les surfaces. Dans QGIS, le SCR WGS 84 / Equal Earth Greenwich (code EPSG 8857) en est un exemple.

Les SCR semi globaux 

Dans cette catégorie nous retrouvons les SCR de type UTM (Universal Transverse Mercator). Les zones UTM peuvent être vues comme des « quartiers » de globe de 6 degrés de longitude de large. Il y a donc 60 zones (6 * 60° = 360°) UTM. De plus ces zones UTM sont divisées en hémisphère nord et hémisphère sud. Ainsi, chaque zone est divisée en Nord et Sud. Il y a donc au final 120 (60 * 2) zones UTM à l’échelle du globe (Fig. 140).

À l’échelle de la France, nous trouvons trois zones UTM différentes. Pour la partie occidentale du pays, nous trouvons la zone UTM 30 Nord (EPSG 32630), pour la partie centrale et orientale nous trouvons la zone UTM 31 Nord (EPSG 32631) et pour la partie la plus à l’est et la Corse nous trouvons la zone UTM 32 Nord (EPSG 32632).

Ces SCR sont des projections avec une unité métrique. Ces SCR sont très couramment utilisés car ils présentent une précision homogène sur l’ensemble du globe. Il s’agît des SCR choisis par les fournisseurs d’images satellites. Lorsque nous téléchargeons une image Landsat, Sentinel ou SPOT, cette image sera projetée dans la zone UTM locale. Par exemple, si nous téléchargeons une image Landsat ou Sentinel sur la Corse, son SCR sera en UTM zone 32 N (32632).

Travailler en UTM zone est souvent un bon choix par défaut lorsque nous travaillons sur une région dont nous ignorons le SCR officiel. Par exemple, si nous travaillons sur le Bangladesh, nous ne connaissons pas le SCR officiel de ce pays à priori. Nous pourrons alors choisir de travailler en UTM zone 45 Nord.

Note

Cette série de SCR reposant sur une projection de type Mercator, les déformations engendrées dans les régions polaires sont donc importantes. Si vous êtes amenés à travailler sur ces régions, il faudra sans doute penser à un SCR plus approprié.

Les SCR continentaux 

Nous ne parlerons que d’un seul SCR continental, à savoir celui utilisé lorsque nous travaillons avec des données à l’échelle de l”Europe. Le SCR dédié est le ETRS89-extended / LAEA Europe de code EPSG 3035 (Fig. 141). Il s’agît du SCR officiel choisi par les instances de l’Union Européenne pour la diffusion de données géographiques.

C’est par exemple dans ce SCR ci que nous trouvons les données du programmes européen d’observation de la Terre Copernicus. Les données d’usage du sol européennes Corine Land Cover sont notamment diffusées sous ce SCR.

Les SCR nationaux 

Chaque pays possède son (ou ses) SCR officiels. En France, il s’agît du Lampert 93, dont le code EPSG est 2154. Le Lambert 93 est une projection conique conforme, reposant, comme le WGS 84, sur l’ellipsoïde IAG GRS80. Comme il s’agît d’une projection conforme, les surfaces ne sont pas conservées. Les déformations sont les plus grandes aux extrémités nord et sud de la France, mais restent totalement satisfaisantes pour la plupart des usages (Fig. 142).

Les données nationales sont, normalement, toujours fournies dans ce référentiel. C’est également le SCR à privilégier lorsque nous travaillons sur un territoire français.

Astuce

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 bien pour le Kamtchatka.

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.

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 Raster ‣ Géoréférencer… ou le menu Couche ‣ Géoréférencer…, selon la version de QGIS. La fenêtre suivante apparaît (Fig. 143).

Nous ajoutons le raster à géoréférencer en cliquant sur l’icône Ouvrir un raster… . La carte apparaît alors dans le panneau central. Nous pouvons naviguer et zoomer / dézoomer à l’aide des icônes de navigation . Nous pouvons maintenant passer à la définition des points d’ancrage.

Astuce

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 (Fig. 144).

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 (Fig. 145).

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 . 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é (Fig. 146).

Avertissement

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

Astuce

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 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 (Fig. 147).

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 Paramètres ‣ Paramètres de transformation…. La fenêtre suivante apparaît (Fig. 148).

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

Astuce

Il est possible de générer la commande GDAL sous jacente à ce processus de géoréférencement dans le menu 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 (Fig. 149).

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.

Avertissement

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.

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 d’un raster 

Nous allons voir ici différents outils qui permettent de reprojeter un raster.

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 Exporter ‣ Enregistrer sous… . La fenêtre suivante apparaît (Fig. 150).

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 pour ouvrir une nouvelle fenêtre qui permet de sélectionner le système désiré (Fig. 151).

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

Dans la fenêtre précédente, le champ SCR se met bien à jour. Nous cliquons sur 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 Raster ‣ Projections ‣ Projection (warp)…. La fenêtre suivante s’affiche (Fig. 152).

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 Exécuter, le nouveau raster s’affiche dans sa nouvelle projection.

Avertissement

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 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 Raster ‣ Projections ‣ Projection (warp)…. Nous retombons sur la fenêtre de la figure (Fig. 152). Au lieu de remplir les champs à ce niveau, nous cliquons sur le menu Exécuter comme processus de lot… La fenêtre suivante apparaît (Fig. 153).

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 .

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 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 (Fig. 154).

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 Exécuter, les rasters se reprojettent par lot.

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 SCP ‣ Pré-traitement ‣ Reproject raster bands. Le menu suivant apparaît (Fig. 155).

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 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é (Fig. 156)..

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 (Fig. 155), 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 (Fig. 157) présente un ré-échantillonnage spatial à 200 mètres en X et Y d’une bande Landsat.

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.

Reprojection d’un raster avec GDAL 

Version de GDAL : 3.0.4

Une fois que GDAL (GDAL/OGR) 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.

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 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 (Importer un raster dans 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.

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

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

Astuce

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 (GDAL/OGR) propose une commande très simple pour reprojeter un fichier raster via l’outil gdalwarp (Reprojection d’un raster avec 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 (Lancer GDAL avec 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.

# 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 d’une couche vecteur 

Nous allons voir ici différents outils permettant de changer le SCR d’une couche vecteur.

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 Exporter ‣ Sauvegarder les entités sous…. La fenêtre suivante s’affiche (Fig. 158).

Dans le champ 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 .

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

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

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

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