SIG et coordonnées géographiques === *"Si on était capable de voir à 40 000 km, on se verrait de dos."* *~Philippe Geluck* ### :warning: À lire avant de commencer :warning: 1. Copiez (ne coupez pas!) tout le répertoire partagé nommé `IGAST_QGIS_AVANCE` sur votre machine. Ce sera votre base de travail. 2. Faites les exercices seuls, on ne mémorise rien en regardant seulement les autres faire. 3. Posez des questions dès que vous doutez d'avoir bien compris quelque chose! 4. Ce TP part du principe que vous avez déjà manipulé QGIS et que vous savez utiliser un terminal. Si vous n'êtes pas à l'aise, deux réflexes (dans l'ordre d'utilisation): -- https://stackoverflow.com/search -- https://qgis.org/fr/docs/index.html 5. Certaines réponses peuvent être trouvées directement sur le Web. Toutefois le but de cette partie n'est pas tant de donner les bonnes réponses que de s'exercer à jongler avec les coordonnées géographiques, de vous inviter à réfléchir en mettant *les mains dans le cambouis* sur les implications concrètes d'un modèle de la Terre. Répondre "bêtement" aux questions, c'est passer à coté de l'exercice. 6. Posez des questions dès que vous doutez d'avoir bien compris quelque chose! # I] Coordonnées géographiques Dans cette première partie vous manipulerez des coordonnées géographiques basées sur un modèle de la Terre sphérique ou ellipsoïdal de paramètres suivants : > **Modèle sphérique** Rayon $R = 6 371km$ (rayon moyen de la Terre) > **Modèle ellipsoïdal** (SCR : WGS84) Demi grand axe $a= 6378.137km$ Aplatissement $f= \frac{1}{298.257223563}$ Circonférence équatoriale $2\pi a=40075.017km$ Première exentricité $e=0.08181919132$ Deuxième exentricité $e'=0.082094438$ Vous aurez aussi besoin d'une calculatrice, par exemple https://www.desmos.com/scientific ## Ex.1 - Échauffement planétaire. 1. Un *grand cercle* est un cercle tracé à la surface d'une sphère et qui a le même diamètre qu'elle. -- Dans le modèle sphérique, les méridiens sont-ils des grands cercles? Et dans le modèle ellipsoïdal? -- Notre Terre ellipsoïdale compte un seul grand cercle. Quel est-il? ``` Sur une sphère, sont des grand-cercles tous ceux qui résultent de la section de la sphère par un plan passant par le centre de cette sphère. On peut aussi dire que ce sont tous les cercles de périmètre maximal tracés à la surface de la phère. Les méridiens sont donc de grand cercles, comme l'équateur ainsi que tous les cercles obliques dont le centre coincide avec celui de la sphère. Dans le modèle ellipsoïdal, seul l'équateur correspond à cette définition. ``` 2. Voici un tableau comprenant des coordonnées géographiques exprimées en degrés décimaux et sexagésimaux. -- décrivez la méthode de passage des degrés décimaux aux degrés sexagésimaux, et vice-versa. -- complétez les valeurs manquantes. | Degrés décimaux | Degrés sexagésimaux | | -------- | -------- | | | 30° | | 15.41° | | | | 27°12'58" | | 64.1816884° | | | | 48°16'24.55" | ``` | Degrés décimaux | Degrés sexagésimaux | | -------- | -------- | | 30° | 30° | | 15.41° | 15°24'36" | | 27+0.2+0.016111 =27.21611° | 27°12'58" | | 64.1816884° | 64°10'54.08" | | 48.2734° | 48°16'24.55" | ``` 3. Si la Terre était réduite à un ellipsoïde de demi grand-axe $a=1m$, quelle serait la distance entre un pôle et la position de ce même pôle si la Terre était une sphère de rayon $R=a$ ?. À titre de comparaison, quelle serait alors l'altitude du mont Everest (8448m)? Pensez-vous que l'on peut considérer la Terre comme sphérique pour réaliser des mesures grossières sur SIG? ``` Pour rappel le facteur d'aplatissement de l'ellipsoïde est donné par f = (a-b)/a où a est son demi grand-axe et b son demi petit-axe. La distance est donnée par l=a-b = a*f. Réduisant la Terre à un ellispoïde de demi grand axe a=1m, on trouve l=3.35mm. L'Everest (8448m) se trouve ainsi réduit à 1.38mm. L'aplatissement est si faible qu'il serait imperceptible à l'oeil nu si l'on regarde la Terre depuis l'espace. Toutefois, il n'est quand même pas négligeable et affecte les mesures de précision. Le modèle sphérique est néammoins largement utilisé dans les SIG pour réaliser des calculs rapides et dont la précision est satisfaisante pour des applications ne nécessitant pas une précision très grande et à condition de travailler à grande échelle (l'approximation sphérique affecte les calculs d'autant plus que les distance sont grandes). ``` ## Ex. 2 - Longueur de degré de latitude $\Delta^1_{\phi}$. En 1795 naissait le mètre, unité de longueur voulue universelle destinée en France à remplacer toutes les anciennes mesures de longueur jusqu'ici en usage, notamment la *toise*. Il fut décidé que sa définition s'appuiera sur un objet naturel, commun à tous les Hommes, connu de tous les peuples et qui permette de re-mesurer le mètre en tout point de la planète. Cet objet naturel, bien sûr, c'est la Terre. Le mètre fut alors défini arbitrairement comme le dix-millionième ($10^{-7}$) du quart du méridien terrestre. Il fut alors nécessaire de mesurer la longueur d'une portion de méridien, ce qui fut exécuté en pleine guerre par les citoyens astronomes Delambre et Méchain entre Dunkerque et Barcelone par triangulation. Après près de 3 ans de mesures, ils trouvèrent finalement que les deux villes étaient séparées de 9°40'25.40", ou 551584.72 toises et calculèrent qu'un quart de méridien mesurait 5130740 toises. On obtenait donc que 1 toise = 1.9490625m. On s'attendrait du coup à ce que la longueur d'un quart de méridien soit de dix millions de mètres. Pas de chance, les satellites modernes qui ont mesuré entièrement ce quart de méridien avec très grande précision nous apprennent qu'il est en fait long de 10002290m, soit 2290 mètres d'écart! Le mètre initial était donc trop court de 0.229 millimètres! De nos jour, le mètre est défini d'après une grandeur (cf. [Wikipedia](https://fr.wikipedia.org/wiki/M%C3%A8tre)). 1. L'affirmation *la longueur d'un degré de méridien est plus grande près des pôle que près de l'équateur* est elle vraie? Pouvez-vous trouver une explication intituive à cela? ``` La différence est causée par la variation de courbure des ellipses méridiennes. Lorsque l'on s'approche des pôles, la courbure de l'arc méridien est plus faible : la longueur de l'arc intercepté par 2 rayons espacés de 1° augmente d'autant plus qu'on approche de la zone la plus "plate" de l'ellipse, c'est à dire les pôles. Voici un exemple pour en donner l'intuition : Imaginez que vous marchez le long d'un méridien avec en main un bâton toujours pointé vers le centre de la Terre (la Terre est supposée lisse). Pour déterminer la longueur en km d'un 1° d'arc méridien, il vous suffit de choisir un point de départ puis de marcher le long du méridien passant par ce point jusqu'à ce que votre bâton ait pivoté de 1°. Naturellement, plus le méridien est "applati", plus il faudra marcher longtemps pour le bâton ait suffisamment pivoté. ``` 2. Calculer la longueur d'une portion de méridien quelconque sur un ellipsoïde est [assez complexe](https://en.wikipedia.org/wiki/Meridian_arc), mais les choses sont plus simples si on mesure seulement 1° d'arc. On a alors : $\Delta^1_{\phi} = \frac{\pi a (1-e^2)}{180°(1-e^2\sin^2\phi)^{\frac{3}{2}}}$. Où $a$ est le demi-grand axe de l'ellipsoïde, $e$ son excentricité, $\phi$ la latitude à laquelle est mesuré le degré (en degrés décimaux). Calculez la longueur d'un degré d'arc méridien à l'équateur $\phi=0°$, à Paris $\phi\approx 48.85°$ et au pôle nord $\phi=90°$. Est-ce cohérent avec votre réponse à la question précédente? ``` Delta^1_0 = 110.57km Delta^1_48.85 = 111.21km Delta^1_90 = 111.69km On constate bien que la longueur d'un degré de latitude augmente en se rapprochant des pôles. ``` 3. La différence de latitude entre Paris et Amiens est de 1°, idem entre Pretoria et Johannesburg en Afrique du Sud. Ces villes sont deux à deux situées approximativement sur le même méridien. On vous dit un jour que la distance Paris-Amiens est la même que Pretoria-Johannesburg. Que pensez-vous de cette affirmation? ``` C'est bien sûr faux, la distance entre Paris et Amiens est nécessairement plus grande que celle entre Pretoria et Johannesburg, cf. questions précédentes. ``` 4. Qu'en serait-il si notre planète était une sphère parfaite? ``` Dans ce cas l'affirmation serait vrai puisque dans ce cas un méridien est un cercle. La longueur d'un arc de 1° sir ce cercle est constante. ``` 5. L'erreur dans la détermination initiale du mètre provenait du fait que Delambre et Méchain utilisaient un fil à plomb pour déterminer la verticale locale et calculer la latitude de leur position lors des opérations de triangulation. -- À votre avis, pourquoi cela a-t-il causé des erreurs? Pour vous aider à répondre, rappelez vous la définition du goïde, et pensez que les erreurs principales étaient dans les mesures de Méchain effectuées dans les Pyrénées. ``` Le fil à plomb que Méchain utilisait était dévié par la masse des Pyrénées : il ne pointaint donc plus vers le centre des masses de la Terre. Comme ce fil à plomb était utilisé pour déterminer la latitude du lieu où étaient effectuées les mesures, elles s'en sont trouvées faussées. ``` ## Ex. 3 - Longueur d'un degré de longitude $\Delta^1_{\lambda}$. On appelle *longueur d'un degré de longitude* la longueur (en mètres, km, etc.) d'une portion de parallèle entre deux méridiens séparés d'une distance angulaire de 1°. 1. L'affirmation suivante est-elle vraie : *"Un degré de longitude est plus long à l'équateur qu'à la latitude de Paris"*? Expliquez pourquoi. 2. Considérons un instant la Terre comme parfaitement sphérique ($R=6371km$). Calculez la longueur d'un degré d'arc de longitude à 0° de latitude, 40° et 65° avec la formule ci-dessous. Reportez vos réponses dans le tableau ci-dessous. $\Delta^1_{\lambda} = \frac{\pi R \cos \phi}{180°}$. 3. Faites de même avec un modèle de Terre ellipsoïdique en utilisant la formule suivante: $\Delta^1_{\lambda} = \frac{\pi a \cos \phi}{180°\sqrt(1-e^2\sin^2\phi)}$. | Latitude $\phi$| km/° - Sphère| km/° - Ellipsoïde|Emplacement approximatif| | -------- | -------- |-------- |-------- | | 0 | ||Quito| | 40 | ||Naples| | 65 | ||Reykjavik| | 75 | ||Station Concordia, Antarctique| 4. Si l'on dessine les contours d'une zone de 1° de coté en latitude et longitude, quelle forme obient-on... -- ...si la zone se trouve proche du pôle Sud? -- ...du pôle Nord? -- ...de l'équateur? 5. La différence de longitude entre le centre de Paris et l'ENSG est d'environ 14'24". En supposant que ces deux lieux sont situés sur le même parallèle, quelle distance les sépare? ## Ex. 4 - Une affaire de précision. 1. Complétez le tableau ci-dessous avec la longueur en mètres de chaque portion de parallèle pour une latitude constante de 48.84°. | Longitude $\lambda$ | Longueur en m. | | -------- | -------- | | 1° | | | 1' | | | 1" | | | 0.1" | | | 0.01" | | | Longitude $\lambda$ | Longueur en m. | | -------- | -------- | | 1° | | | 0.1° | | | 0.01° | | | 0.001° | | | $10^{-6}$° | | 2. Sur Google Maps, l'ENSG est localisée aux coordonnées géographiques $(\phi=48.841023°,\lambda=2.5873236°)$. -- Les coordonnées sont données avec une précision arithmétique de $10^{-6}$°. À quelle longueur sur le terrain correspond $10^{-6}$°? Utilisez les résultats de la question précédente pour répondre. -- Sachant que la BDTOPO, base de référence de l'IGN, a une exactitude planimétrique de l'ordre du mètre, pensez-vous vous que le degré de précision de Google Maps reflète la réelle précision des données du site? 3. [Wikipedia](https://en.wikipedia.org/wiki/Wikipedia:WikiProject_Geographical_coordinates)(https://en.wikipedia.org/wiki/Wikipedia:WikiProject_Geographical_coordinates) fournit un guide de saisie de coordonnées géographiques pour ses utilisateurs intitulé "WikiProject Geographical coordinates". Ce document est précieux car les conseils qu'il donne sont valables dans le cas général. Il donne des tables de précision pour différentes valeurs d'angles de latitude et différentes tailles d'entités géographiques. Lisez la partie intitulée *Precision guidelines*, puis répondez aux questions suivantes : - Combien de décimales sont nécessaires pour positionner suffisamment précisément un objet de 1000m de coté à 45° de latitude? - Cela a-t-il un sens de localiser un polygone représentant Paris (approx. 7km entre les limites Est et Ouest de la ville) avec une longitude en degrés décimaux ayant 5 chiffres après la virgule? Pourquoi? - Quelle serait la précision la plus adaptée dans ce cas là? 4. Vous avez traversé la rue et trouvé un travail de géomaticien à Paris. Un jour, quelqu'un vous envoie un fichier contenant des points GPS localisant des cabanes à oiseaux dont les coordonnées en degrés décimaux sont arrondies à 3 chiffres après la virgule. Ces données sont-être utilisables? Un tel objet mesurant une dizaine de centimètres, quelle serait la précision adaptée d'après Wikipedia? --- # II] GDAL, QGIS et les coordonnées géographiques ## Ex.5 : Reprojections avec GDAL et OGR. Dans le répertoire *Exercice 5* se trouvent deux jeux de données : un raster au format GeoTiff nommé `natural_earth_1.tif` et un vecteur au format Shapefile nommé `countries`. Si vous travaillez sous Windows, vous aurez besoin de la console *OSGeo4WShell*. Vous la trouverez en principe dans le `Menu Démarrer->QGIS`. Sous les distributions MacOS ou Linux, GDAL et OGR sont directement accessibles depuis un terminal. Pour les questions suivantes, utilisez la documentation des commandes GDAL et les exemples qu'elle donne pour construire vos requêtes. 1. À l'aide de l'outil `gdalinfo` (https://www.gdal.org/gdalinfo.html), trouvez le SCR `natural_earth_1.tif` et son code EPSG. Ce SCR est-il associé à une projection cartographique? 2. Utilisez `ogrinfo` (https://www.gdal.org/ogrinfo.html) avec l'option `al` pour vérifier que le fichier `countries.shp` utilise le même SCR. 3. Chargez `natural_earth_1.tif` dans QGIS et vérifiez que le SCR du projet est celui du raster. Celui-ci utilise des coordonnées géographiques, et devrait donc vous présenter une représentation en 3D sur ellipsoïde. Mais la vue de QGIS est en 2D, donc nécessairement en projection. Observez les distorsions de la carte, et identifiez quelle projection utilise QGIS pour ce SCR à l'aide de [cette page Wikipedia ](https://fr.wikipedia.org/wiki/Liste_de_projections_cartographiques). Est-ce une bonne idée de l'utiliser pour cartographier la France? 3. A l'aide de l'outil `gdalwarp` (https://www.gdal.org/gdalwarp.html), reprojetez ce raster dans le SCR Lambert 93 (EPSG:2154). Ouvrez le résultat dans QGIS. Quel est l'effet ce cette projection sur les formes des pays? Quel est à votre avis le pays qui se trouve le moins déformé? Utiliseriez-vous cette projection pour cartographier le Brésil? Est-ce, de manière générale, une projection pertinente pour un planisphère? 4. Ouvrez le fichier `countries.shp` dans QGIS 5. A l'aide de l'outil `ogr2ogr` (https://gdal.org/programs/ogr2ogr.html), projetez le jeu de données vectoriel `countries` dans le SCR Lambert 93. Attention, l'ordre des fichiers d'entrée et de sortie sont inversés par rapport à gdalwarp! `ogr2ogr` devrait produire une erreur. Quelle est-elle et, à votre avis, pourquoi s'est-elle produite? Pour vérifier votre hypothèse, ouvrez le fichier vecteur reprojeté dans QGIS : quel pays manque-t-il? ## Ex.6 - Préparation du projet QGIS 1. Créez un nouveau projet dans QGIS. 2. Allez dans le menu `Préférences -> Options -> SCR`. Assignez le SCR par défaut "EPSG:3857 WGS 84 /Pseudo Mercator" et cochez l'option "Demander le SCR". Cliquez sur "OK" pour valider. 3. Ce SCR est aussi nommé "Web Mercator" car il est très largement utilisé pour la carographie sur le Web. Lisez le paragraphe *Properties* de sa page Wikipedia : https://en.wikipedia.org/wiki/Web_Mercator_projection. Que signifie qu'elle n'est pas *conforme*? Comparez à la projection utilisée par QGIS pour le SCR EPSG:4326. Celle-ci est-elle conforme? 4. Allez dans `Extensions->Installer des extensions`, recherchez "QuickMapServices" et installez le si ce n'est déjà fait. 5. Activer l'extension "Saisie de coordonnées". 6. Dans le menu `Internet->QuickMapServices`, charger le flux WMS de OpenStreetMap Standard. 7. Sauvegardez votre projet. Les prochains exercices sont à réaliser avec ce projet. N'oubliez pas de sauvegarder régulièrement, d'autant plus si vous travaillez sous Windows où QGIS n'est pas toujours très stable. ## Ex.7 - Mesures de distances Mesurer la distance entre deux points sur Terre est l'opération la plus basique que l'on peut demander à un SIG. Nous allons voir pourtant qu'une fois encore, la forme de la Terre introduit quelques subtilités... Passez le projet en projection `WGS 84 / Mercator (EPSG:3395)`. Les projections de Mercator sont conformes : les angles sont conservés localement. Une conséquence de cette propriété est qu'une ligne tracée sur une projection Mercator coupe tous les méridiens à angle constant. ![](https://i.imgur.com/hrV7m5Y.png) Sous QGIS, créez une nouvelle couche temporaire en mémoire de type `Polyligne/Courbe` et de SCR `EPSG:3395`, puis tracez une polyligne dont l'une extrémité est la Tour Eiffel à Paris et l'autre la statue de la liberté à New York. Sélectionnez la couche dans le panneau des couches, et selectionnez la polyligne créée avec l'outil de sélection d'entité ![](https://i.imgur.com/rbHwRpj.png). Ouvrez ensuite la console Python (`Extension->console Python`) et entrez y le code suivant : ```python layer = qgis.utils.iface.activeLayer() # récupère la couche actuellement sélectionnée features = layer.selectedFeatures() # récupère les éléments selectionnés dans la couche layer geom = features[0].geometry() # récupère la géométrie du premier élément sélectionné print(geom.length()/1000.0) # calcule la longueur en km de cette géométrie # Note : cette ligne fonctionne car la géométrie a été créée dans un SCR projeté dont l'unité est le mètre. Si vou s aviez créé la couche avec le SCR WGS84 EPSG:4326, geom.length() renverrait une longueur en degrés. ``` 1. Notez la valeur que vous obtenez. Utilisez maintenant l'outil de mesure de QGIS ![](https://i.imgur.com/uoUt12V.png) pour mesurer la polyligne *à la main*. Assurez vous que dans la fenêtre de mesure le radio-bouton "Ellipsoïdale" est coché. Pour vous faciliter la vie vous pouvez activer l'accrochage du curseur aux données dans le menu `Projet->Options d'accrochage`. 2. Notez la nouvelle valeur obtenue. Pourquoi les deux valeurs sont-elles différentes? Cochez le radio-bouton "Cartésien", et comparez aux valeurs déjà obtenues. Comment expliquez-vous le résultat? 3. Pensez-vous que la polyligne tracée décrit le chemin le plus court entre ses deux extrémités? Si vous pensez que non, à quoi ressemblerait la route la plus courte? 4. Trouvez sur la page https://fr.wikipedia.org/wiki/Loxodromie le nom du type de route tracée par la polyligne. Déduisez-en le nom de la route mesurée par l'outil de mesure de QGIS. 5. Est-ce que tracer des polylignes pour mesurer des distances sur de longues distances vous semblent une bonne idée? 6. Est-ce qu'il est possible sous QGIS de mesurer des distances indépendament de la projection cartographique utilisée? 7. À l'aide de la fenêtre de propriétés du projet, trouvez quel modèle d'ellispoïde est utilisé par QGIS pour effectuer des mesures de distance. Pourquoi est-ce cet ellipsoïde qui est sélectionné par défaut? 8. On dit que l'ellipsoïde utilisé par les SCR français légaux est une densification locale de l'ellipsoïde standard WGS84, et qu'on peut utiliser WGS84 sans affecter les calculs. Trouvez la page du CRS français légal Lambert 93 sur le site https://epsg.io, et le nom de l'ellipsoïde qu'utilise ce SCR. Changez l'ellipsoïde utilisé par QGIS par celle-ci, et refaite la mesure de distance précédente. Pouvez-vous confirmer l'absence d'effet sur le calcul? 9. Certes, QGIS permet d'effectuer des mesures de distances sur ellipsoïde quelque soit la projection choisie, mais si les projections ne permettaient pas de mesurer des distances sur des cartes planes, cela réduirait considérablement leur utilité! En fait, chaque projection est adaptée à une région donnée, une échelle donnée et un but précis. Ainsi, la projection Lambert 93 a été spécialement conçue pour permettre d'effectuer des mesures planes en France métropolitaine, en minimisant l'erreur résultante. Cette erreur est appelée *altération linéaire*. -- Cherchez sur le Web l'amplitude de l'altération linéaire de Lambert93. Trouvez sur la figure ci-dessous l'altération linéaire à Paris. -- Passez le projet en Lambert 93. Mesurez la longueur d'un bâtiment de l'ENSG sur QGIS, en repère cartsésien et sur ellipsoïde. Est-ce cohérent avec l'indication de la figure? Même question pour un baiment à Dunkerque. ![](https://i.imgur.com/GKjOpgU.png) Figure : altération linéaire comparée entre projections Lambert. 10. Sachant tout cela, est-il acceptable d'utiliser une mesure plane en Lambert 93 pour estimer la longueur en km d'un trajet dans Paris? Et pour calculer la longueur de fibre optique nécessaire pour relier directement l'ENSG à l'IGN? **Question bonus 1** : Pensez-vous que le même problème se pose pour les surfaces? S'il vous reste du temps en fin de TP, vérifiez cela en effectuant la même procédure pour cette fois mesurer une aire (`utilisez la fonction geom.area() en python`). **Question bonus 2** : Il existe une classe de projection, dites gnomoniques (en fait des azimutales), qui ont la propriété interessante de transformer les grands cercles en droites. Dans ce type de projection, on voit la Terre comme si on la regardait depuis son centre, en imaginant qu'elle est creuse (cf. la figure ci-dessous). ![](https://i.imgur.com/6N2zgQQ.png) Figure: projection gnonomique, le plan de projection est ici tangent au pôle Nord. Image issue de https://gisgeography.com/azimuthal-projection-orthographic-stereographic-gnomonic/. Aucun des SCR de QGIS ne propose de projection gnomonique. C'est l'occasion de voir qu'il est possible de créer nous-même des projections personnalisées. Allez dans le menu `Préférences->Projections personnalisées`, puis ajoutez un nouveau SCR dont la chaîne PROJ4 est la suivante : ```+proj=gnom +lat_0=90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs``` -- En lisant cette chaîne, pouvez-vous dire où est situé le plan tangent? Attention, dans ce type de projection les surfaces situées à plus de 90° de latitude du centre de la projection ont une aire infinie. QGIS risque de planter si vous affichez un raster. Vous pouvez le remplacer par la couche vecteur "countries.shp" qui évitera ce problème. Voici mainetnant une façon de vérifier visuellement la propriété des projections gnonomiques sur les grands cercles : -- commencez par changer la projection du projet pour la projection que vous venez de créer. Profitez-en pour vérifier que le rendu ressemble bien à la figure ci-dessus (à ceci-près qu'on visualise toute la Terre. -- créez une couche en mémoire tempraire de type polyligne et ajouter une polyligne allant de Paris à New-York. Il faut ensuite ajouter des points dans la polyligne. Pour cela, utilisez l'outil "densifier par le nombre de sommets" dans le menu vecteur et ajouter 2000 points. -- Ne conservez que la polyligne densifiée, puis revenez à une projection Mercator classique (EPSG:3395). Pouvez-vous confirmer que la polyligne décrit bien une route orthodromique? ## Ex.8 - Sauvetage de fichier. Les fichiers GeoTiff peuvent directement stocker leurs métadonnées de géoréférencement à l'interieur de leur structure, ou les stocker dans des fichiers externes : - un fichier .prj contenant la chaîne proj4 (syntaxe WKT) de leur SCR; - un fichier .wld (on trouve aussi l'extension .tfw) contenant les paramètres [GeoTransform](https://www.gdal.org/gdal_datamodel.html). GeoTransform est le nom de la transformation affine qui permet au SIG d'associer les coordonnées géographiques ou projetées à chaque pixel de l'image. Bien que le fichier .wld seul suffise à un SIG pour calculer les coordonnées d'un raster, encore faut-il qu'il sache dans quel SCR sont exprimées ces coordonnées..ce que contient le fichier .prj. Les deux sont donc essentiels lorsque le raster ne stocke pas lui même ses métadonnées de géoréférencement. Vous trouverez un fichier raster géoréférencé `champs_sur_marne.tif`, stocké dans le répertoire *Exercice 8*. Il s'accompagne d'un fichier .wld, mais, pas de chance, le fichier .prj est manquant! 1. Allez dans `Préférences->Options->SCR` et cochez `Aucune projection`. 2. Ouvrez le raster `champs_sur_marne.tif` dans QGIS, ainsi que la couche WMS OpenStreetMap. Que se passe-t-il et pourquoi? 3. Ouvrez le fichier .wld associé au raster avec un éditeur de texte. A l'aide de la documentation du format *World File* (https://fr.wikipedia.org/wiki/World_file), calculez la transformation les coordonnée du point $(0,0)$ de l'image GeoTiff. Activez la projection Lambert 93 dans QGIS, et regardez à l'aide de la souris où 'tombent' ces coordonnées. Que pouvez-vous en conclure sur le SCR inital de `champs_sur_marne.tif`? Vérifiez votre hypohtèse en copiant les coordonnées calculées sur le site http://geofree.fr/gf/projguess.asp. 4. Modifiez de nouveau les préférences pour que QGIS demande le SCR des nouvelles couches. 5. Il se trouve que ce raster est un extrait de la BDORTHO de l'IGN et est donc projeté en Lambert 93 (EPSG : 2154), ce que vous devriez d'ailleurs avoir déterminé dans la question 3. Reste donc à corriger ce raster. Pour cela, deux stratégies sont possibles : ### Stratégie 1 Par chance, il y a dans le répertoire *Exercice 12* le raster `saint_maur.tif`, qui n'a rien à voir avec Champs sur Marne mais qui est projeté en Lambert 93 avec des métadonnées internes. Une façon de faire est alors d'extraire les métadonnées de projection de ce raster et de les copier vers le raster `champs_sur_marne.tif`. 5.1. Ouvrez le raster `saint_maur.tif` dans QGIS, vérifiez qu'il est bien projeté en Lambert 93. 5.2. Pour extraire les métadonnées de projections, allez dans le menu `Raster->Projections->Extraire la projection `, sélectionnez le fichier `saint_maur.tif` et cochez **Créer le .prj**. 5.3. Rendez vous dans le dossier du raster, vous devriez vous deux nouveaux fichiers : `saint_maur.wld` et `saint_maur.prj`. Laissons le premier pour le moment, et ouvrez le fichier `.prj` avec un éditeur de texte. Vous voilà devant la représentation [Well Known Text](https://en.wikipedia.org/wiki/Well-known_text) (WKT) de la chaîne proj4 du SCR Lambert 93. 5.4. Quel est le nom complet du SCR? Quel est le nom de l'ellipsoïde qui lui est associé? Quel est le type de cette projection? Est-ce que la projection associée à ce SCR conserve les formes des entités géographiques? 5.5. Copiez maintenant le fichier `saint_maur.prj` et placez sa copie dans le même dossier que le raster `champs_sur_marne.tif`. Nommez cette copie `champs_sur_marne.prj`.QGIS associe les fichier en se basant sur leurs noms, ils doivent être exactement similaires, à l'exception de l'extension. 5.6. Ré-ouvrez le fichier `champs_sur_marne.tif` dans QGIS. Que constatez-vous? ### Stratégie 2 Conserver les métadonnées des GeoTiff dans des fichiers annexes a l'avantage de la lisibilité. Cependant, il n'est pas toujours pratique de traîner 3 fichiers pour un seul raster. Nous allons voir comment, plus simplement, on peut assigner une projection à un raster. 5.7. Cette fois, allez dans le menu `Raster->Projections->Assigner une projection`. Choisissez `champs_sur_marne.tif` comme fichier source, et `RGF93 / Lambert-93` comme SCR. Notez au passage que QGIS utilise GDAL et vous met la commande exacte en dessous. Décochez "Charger dans le canvas une fois terminé" et cliquez sur OK. 5.8. QGIS a renseigné les métadonnées de géoréférencement à l'intérieur du raster; vous n'avez plus besoin du fichier .wld qui l'accompagne. Supprimez-le, réouvrez le raster dans QGIS et vérifiez que tout c'est bien passé : il devrait être maintenant parfaitement superposé au fond de carte. ## Ex.9 - Manipulation des SCR avec QGIS Dans le répertoire *Exercice 9*, vous trouverez deux jeux de données raster au format GeoTiff et un jeu de données vecteur au format Shapefile. 1. Cochez "Aucune projection" dans les propriétés du projet, volet SCR. 3. Chargez les deux GeoTiff et le Shapefile. Pourquoi sont-ils localisés à des endroits différents? Notez le SCR de chacun à partir des propriétés de chaque couche. 4. En vue d'une publication Web de ces données, vous aimeriez commencer par les harmoniser en les projetant toutes en Pseudo Mercator, la projection la plus utilisée sur le Web. Pour commencer, nous allons transformer les données vecteur. Cette opération est très simple : faites un clic droit sur la couche vectorielle et choisissez `Enregistrer sous`. Enregistrez le nouveau fichier sous un nouveau nom au format ESRI Shapefile et choisissez le bon SCR (EPSG:3857). 5. Vérifiez que le nouveau vecteur est bien superposé à la couche WMS de OSM, elle aussi en projection Pseudo-Mercator. 6. Pour les rasters, il y a de nouveau deux stratégies possibles : ### Stratégie 1 **Pour le raster nord.tif.** 6.1. Comme pour les données vecteur, faites un clic droit sur la couche 'nord' et choisissez `Enregistrer sous`. Le menu est un peu différent. Choisissez "Image" comme mode de sortie, le format GTiff et choisissez de nouveau le bon SCR. 6.2. Vérifiez que le raster ainsi exporté est bien en projection Pseudo-Mercator. ### Stratégie 2 **Pour le raster sud.tif.** 6.3. Cette fois, allez dans le menu `Raster->Projections->Projection`. Le fichier source doit être sud.tif, qui est en projection Lambert Zone I (EPSG:27571). Notez que dans cette outil les options sont bien plus nombreuses qu'avec `Enregistrer sous`; mais leur utilisation sort du cadre de ce TP. Choisissez toujours le même SCR cible. 6.4. En principe, les deux rasters et les données vecteurs reprojetées devraient tous être situés au même endroit sans la projection à la volée. 6.5. On peut en profiter pour fusionner les deux rasters en un seul. Pour ceci, allez dans le menu `Raster->Divers->Fusionner`, choisissez les deux fichiers rasters, enregistrez le résultat dans un nouveau fichier. Vous avez maintenant un raster général de Paris, projeté en Pseudo-Mercator et donc directement superposable à des flux WMS venant de Google ou OSM.L'IGN par contre diffuse ses flux en Lambert 93. --- # Troisième partie : Géoréférencement Pour cette troisième partie, les ressources à utiliser se trouvent dans le répertoire *Exercice 10*. Les trois exercices qui suivent portent sur 3 situations de géoréférencement que vous risquez de rencontrer régulièrement, chacun permettant une stratégie de géoréférencement différente. ## Ex.10 - Les coordonnées des bords du raster sont connues. Le répertoire *Exercice 10* contient une photographie aérienne nommée `14.jpg`. Si vous l'ouvrez dans QGIS, vous pourrez constater qu'elle n'est pas géoréférencée : le coin supérieur gauche de l'image est aux coordonnées $(0,0)$. On vous fournit les informations suivantes concernant l'emplacement géographique de cette image : ``` En projection Lambert zone CC 45 (EPSG:3945), les coordonnées des 4 bords sont : Upper Left ( 1414428.599, 4190808.712) Lower Left ( 1414428.599, 4179607.583) Upper Right ( 1430271.987, 4190808.712) Lower Right ( 1430271.987, 4179607.583) ``` Avec ces informations, il est possible d'utiliser l'outil de la suite GDAL `gdal_translate` pour assigner des coordonnées exprimées dans un SCR aux coins d'une image. 1. Allez sur la page de documentation de l'outil `gdal_translate` [https://www.gdal.org/gdal_translate.html](https://www.gdal.org/gdal_translate.html). Sachant que "bounds" signifie "emprise/cadre", quel paramètre va vous être particulièrement utile ici? 2. Que signifient les termes `ulx`, `uly`, `lrx` et `lry`? 3. A l'aide du paramètre `a_srs` et de celui trouvé dans la question précédente, assignez les coordonées Lambert CC45 des coins de `14.jpg` et enregistrez le résultat dans un fichier `14.tif`. 4. Est-ce que l'opération que vous venez de faire est réellement un géoréférencement au sens vu dans le cours? 5. A l'aide de [ce document ](https://www.ideobfc.fr/upload/gedit/1/Documents/FicheT3Les_projections_coniques_conformes_9_zones.pdf), expliquez quel est l'interêt des projections Lambert CC. Quelle est la projection Lambert CC à utiliser pour Paris? ## Ex.11 - Le raster n'a aucun repère géodésique utilisable. Le plus souvent -et presque systématiquement pour les cartes anciennes-, vous n'aurez aucune information sur les coordonnées géographiques ou projetées d'une image à géoréférencer, ni de repère géodésique exploitable dans la carte elle-même. Dans ce cas, il faut passer par une étape de géoréférencement indirecte utilisant des paires de *points de contrôle* définis manuellement. Chaque paire est constituée de deux points : le premier localisé dans léimage (et donc exprimé en pixels) et le second dans un SCR géographique ou projeté cible, dans lequel la carte sera transformée. Cette tâche est particulièrement fastidieuse, mais heureusement QGIS fournit un outil graphique facilitant grandement la tâche. Dans cet exercice, nous allons géoréfencer une carte historique de Londres réalisée par le Docteur John Snow en 1854 et [ayant révolutionné la médecine](https://fr.wikipedia.org/wiki/%C3%89pid%C3%A9mie_de_chol%C3%A9ra_de_Broad_Street). **Conservez bien le résultat final : vous en aurez besoin pour un prochain TP.** 1. Si ce n'est pas déjà fait, chargez le flux WMS OpenStreetMap depuis le plugin QuickMapService. 2. Choisissez comme SCR du projet : OSGB 1936 British National Grid (EPSG : 27700). 3. Placez vous à Londre, au niveau de [Golden Square](https://goo.gl/maps/HXKHudVza6y) (aidez vous de Google Maps pour le localiser). 1. Ouvrez le géoréférenceur de QGIS, via le menu `Raster->Géoréférencer->Géoréférencer`. S'il n'est pas disponible, vous devrez d'abord l'activer dans les extensions de QGIS. 2. Dans la fenêtre du géoréférenceur, chargez l'image `snow_map.jpg` à l'aide du bouton ![](https://i.imgur.com/Wb9RIUG.png) (ou `Fichier->Ouvrir Raster`). QGIS va vous demander le SCR de cette carte. Dans le cas de la carte de Snow, on ne le connait pas. Par défaut, choisissez le même SCR que celui visé, soit OSGB 1936 British National Grid. 3. Commencez à placer des paires de points de contrôle avec l'outil ![](https://i.imgur.com/5PjEQL8.png) sur des éléments qui semblent communs sur la carte OpenStreetMap et la carte de Snow : coin d'immeuble, fond d'impasse, etc. La méthode générale est la suivante : -- Naviguez dans la carte de Snow et sur OpenStreetMap en parallèle avec les outils zoom (molette de la souris) et pan (![](https://i.imgur.com/1CbSw5H.png)) et trouvez le croisement des rues *New Brulington Street* et *Savile row* (en bas à gauche de la carte). -- Placez un premier point dans la carte de Snow, sur l'immeuble formant l'angle sud de cet embranchement : ![](https://i.imgur.com/gI6rGlw.png) -- Cliquez devrait ouvrir un popup qui vous demandera les coordonnées correspondantes dans le SCR cible. Choisissez l'option `Depuis le canvas de la carte` qui vous permet de selectionner un emplacement sur la vue générale de QGIS. La fenêtre de géoréférencement est alors momentanément cachée. Cliquez alors sur le même élément dans OSM, ce qui fera réapparaitre le géorférenceur et le popup. Cliquez sur "OK" : votre première paire de points de contrôle est placée! 7. Placez encore au moins **4 points de contrôle** dans le quart "bas-gauche" de la carte de Snow. Note 1 : Essayez de placer tous les points en étant à la même échelle sur OSM comme sur la carte de Snow. Sur OSM, une échelle de 1:2000e fait l'affaire (utilisez la molette de la souris pour naviguer entre plusieurs niveaux de zoom). Note 2 : Vous pouvez annuler le placement en cours d'un point sur le canvas QGIS avec un clic droit. Note 3 : Utilisez les boutons ![](https://i.imgur.com/O8QTW1Q.png) pour supprimer ou déplacer un point de contrôle, sur l'image à géoréférencer ou sur le canvas QGIS. 8. Ouvrez le menu ![](https://i.imgur.com/63KCOSa.png) : celui-ci vous permet de paramétrer différents éléments du géoréférencement, notamment (1) le type de transformation géométrique qui sera utilisé, (2) le SCR dans lequel géoréférencer l'image, (3) la méthode de rééchantillonage des pixels lors de la transformation et enfin (4) le nom du fichier géoréférencé. Décochez l'option "Utiliser 0 pour la transparence si nécessaire". 9. Enregistrez les points de contrôle avec le bouton ![](https://i.imgur.com/G1oLR3i.png) à coté de l'image `snow_map.jpg`. Ainsi, vous pourrez réouvrir le géoréférenceur plus tart et améliorer le géoréférencement en ajoutant de nouveaux points. 10. Choisissez une transformation de Helmert, choisissez "OSGB 1936 British National Grid (EPSG : 27700)" comme SCR cible et choisissez le nom du fichier GeoTiff géoréférencé qui sera créé. 11. Lancez le géoréférencement avec ![](https://i.imgur.com/S7HCktt.png). Observez le résultat en jouant sur la transparence du raster géoréférencé (`clic droit sur le nom du raster dans la liste des couches -> Propriétés -> Transparence`). Observez plus précisément la correspondance entre le fond OSM et la carte géoréférencée dans sa partie Sud-Ouest, puis Nord-Est. Trouvez-vous le résultat satisfaisant? A votre avis, pourquoi obtenez vous ce résultat? 12. Reprenez le géoréférencement et ajoutez encore au moins encore 8 paires de points, réparties cette fois sur toute la carte. Enregistrez ces points dans un fichier nommé `snow_map.png.points`. 13. Exécutez le géoréférencement en choisissant une transformation **polynomiale d'ordre 3**. Vérifiez le que le résultat est meilleur que précédemment. 14. Dans le géoréférenceur, quelle est la valeur de la RMSE pour cette transformation? Qu'est-ce qu'elle signifie? 16. Executez un nouveau géoréférencement avec une transformation **Thin Plate Spline (TPS)**. Mais, cette fois, vous allez exécuter le géoréférencement en ligne de commande avec GDAL pour stocker le résultat dans un fichier nommé `map_snow_tps_osgb1936.tif`. Enregistrez préalablement les point de controle dans le fichier `map_snow.png.points`. Cliquez sur le bouton ![](https://i.imgur.com/5jKA0tO.png) puis copiez le script généré. -- combien de commandes sont exécutées par ce script? -- exécutez les deux commandes dans un terminal, et ouvrez le résultat dans QGIS. -- le géoréférencement est un type de reprojection particulier dans lequel la projection de départ est inconnue. `gdalwarp` est la commande qui effectue cette reprojection. Quels sont les paramètres de la commande à utiliser pour reprojeter un raster depuis le SCR ESPG:2154 vers le SCR EPSG:3395 (utilisez la documentation de `gdalwarp` https://gdal.org/programs/gdalwarp.html)? Vérifiez avec l'outil QGIS `Raster->Projections->projection (warp)`. Vous pouvez constater au passage qu'un certain nombre d'outils de QGIS ne sont qu'une surcouche graphique par dessus GDAL! ## Ex.12 - Géoréférencement vecteur. Parfois, on géréoférence des cartes ou des images pour pouvoir en extraire le contenu par *vectorisation*, c'est à dire dessinant des objets géographiques vectoriels par dessus l'image pour constituer des bases de données géographiques. C'est encore la base de la production des bases de données vectorielles de l'IGN, à ceci près que les images aériennes support ne sont pas géoréfencées à la main... Lorsque l'objectif est de vectoriser un raster, on a deux possibilités : 1) Géoréférencer le raster, puis vectoriser son contenu dans le système de projection du raster géoréférencé. Ainsi, les coordonnées des données vectorielles seront directement dans ce système. 2) Vectoriser le raster contenu dans l'espace image original, non géoréférencé. Géoréférencer ensuite le raster, puis réutiliser les points de contrôle créés pour appliquer le même géoréférencement aux données vectorielles. Les deux méthodes ont leurs avantages et leurs inconvénients. Le principal avantage de la première est la simplicité, mais en contrepartie il est très difficile de retoucher le géoréférencement du raster une fois les données vectorielles créées. La seconde pose problème si ce que l'on veut vectoriser est réparti sur plusieurs images (ou bien il faut les fusionner préalablement, ce qui n'est pas toujours possible). En contrepartie, elle permet de modifier le géoréférencement des données même après vectorisation. Le but de cet exercice est d'expérimenter la seconde méthode pour vectoriser la carte du docteur Snow. Vous aurez l'occasion d'expérimenter la première méthode, largement plus commune, dans d'autres TP. 1. Dans QGIS, chargez la carte de Snow non géoréférencée `snow_map.png`. Regardez les coordonnées de la carte indiquées par QGIS en bas de la fenêtre principale. Etes-vous en projection? Pourquoi les coordonnées du coin haut-gauche de l'image sont $(0,0)$? 2. Créez ensuite une nouvelle couche vectorielle de type point, puis saisissez les positions des 13 pompes. **Placez chaque point au centre du cercle représentant la pompe!** Sauvegardez la couche dans un fichier shapefile nommé `pompes.shp` ![](https://i.imgur.com/d3W48FO.png) 3. La commande [`ogr2ogr`](https://gdal.org/programs/ogr2ogr.html) permet d'effectuer de nombreuses opérations sur des données vectorielles, notamment de les reprojeter et de les géoréférencer avec des points de contrôles. Préparez dans un fichier texte la commande complète permettant de géoréférencer le fichier Shapefile `pompes.shp` à l'aide des points de contrôle déjà créés pour le raster `snow_map.png`. La transformation à utiliser est `Thin Plane Spline` et le resultat doit être stocké dans un second fichier Shapefile nommé `pompes_tps_osgb1936.shp`. Note: vous aurez besoin des paramètres `t_srs` et `-gcp` et `-tps`. `-gcp` permet de passer une paire de point de contrôle à la commande sous la forme : `-gcp ungeoref_x ungeoref_y georef_x georef_y` Il vous faut utiliser les points de contrôle stockés dans le fichier de points de contrôles `snow_map.png.points` dont chaque ligne contient une paire de points de contrôle 2D sous la forme `georef_x,georef_y,ungeoref_x,ungeoref_x,enable,dX,dY,residual`. **:warning:** L'ordre des points de contrôles demandé par le paramètre `-gcp` est inversé par rapport au fichier! **:warning:** Ainsi, si le fichier `snow_map.png.points` contient 2 paires de points de contrôle, par exemple : ``` mapX,mapY,pixelX,pixelY,enable,dX,dY,residual 804.23,-2114.26,996.33,-1265.64,1,0,0,0 5011.52,-1838.92,2778.93,-3279.98,1,0.00,0.00,0.00 ``` Alors vous devrez passer ces deux paires à la commande `ogr2ogr` ainsi : ```ogr2ogr -gcp -2114.26,996.33026874116058025 804.23,-2114.26 -gcp 2778.93,-3279.98 5011.52,-1838.92 [suite des paramètres...] pompes_tps_osgb1936.shp pompes.shp ``` Chargez le fichier `pompes_tps_osgb1936.shp` dans QGIS et vérifiez qu'il est correctement géoréférencé. :::warning **:warning: IMPORTANT :warning:** **Conservez les fichiers `snow_map.png`, `snow_map.png.points`, `map_snow_tps_osgb1936.tif`, `pompes.*` et `pompes_tps_osgb1936.*`, vous en aurez besoin pour un prochain TP! ::: ## Ex.13 - Le raster possède un carroyage. Un carroyage est une grille dont les coordonnée sont données dans la projection de la carte. Pour les cartes IGN en projection Lambert, ce carroyage est aligné avec le méridien origine de la projection : toutes les lignes verticales du carroyage sont parallèles à ce méridien. Vous trouverez dans le répertoire *Exercice 10* une image non géoréférencée nommée `16.jpg`. Il s'agit d'un extrait d'une carte historique du *Petit Nanterre* datant de 1968. 1. Ouvrez-la avec un éditeur d'images (Paint ou autre) ou même QGIS et cherchez le carroyage sur la carte. Identifiez la façon dont est symbolisé sur la carte les points où se croisent les lignes verticales et horizontales du carroyage. 2. Sachant qu'avant 1972, l'IGN utilisait les SCR 'Lambert Nord', 'Lambert Centre', 'Lambert Sud' et 'Lambert Corse' (EPSG:27561, EPSG:27562, EPSG:27563 et EPSG:27564), aidez vous des pages du site epsg.io se rapportant à chaque système pour déterminer quel est le SCR de cette carte. Observez notamment les emprises cartographiques de chaque SCR. 3. Dans cet exercice, vous allez utiliser l'outil de géoréférencement de QGIS selon un mode quelque peu particulier. En effet, vous connaissez déjà les coordonnées Lambert des points du carroyage puisqu'elles sont indiquées sur la carte ancienne. Ouvrez le géoréférenceur de QGIS. 4. Dans la fenêtre du géoréférenceur, chargez l'image `16.jpg` à l'aide du bouton ![](https://i.imgur.com/Wb9RIUG.png) (ou `Fichier->Ouvrir Raster`). 5. Identifiez 4 croix symbolisant les intersections des lignes du carroyage à proximité des 4 coins de la carte. Notez les coordonnées de chacun. 6. Placez maintenant un point de contrôle au point central de chacune de ces croix. Lorsque le popup s'ouvre, renseignez les coordonnées Lambert du point à la place de choisir "Depuis le canvas de la carte", puis cliquez sur "OK". 7. Une fois les 4 paires de points définies, choisissez une transformation globale rigide adaptée pour seulement 4 paires de points (menu ![](https://i.imgur.com/63KCOSa.png)), et choisissez `NTF (Paris) / Lambert Nord France` comme SCR cible. Choisissez enfin dans quel fichier enregistrer le résultat (et cochez "Charger dans QGIS losque terminé"). 8. Lancez la transformation (![](https://i.imgur.com/S7HCktt.png)) et admirez le résultat! 9. A votre avis, pourquoi une transformation rigide suffit pour ce type de géoréférencement? 10. Observez les résidus des points de contrôle et la RMSE du géoréférencement. Pourquoi ne sont-ils pas nuls? 11. Enregistrez les points de contrôle à coté du fichier JPG source. --- # Exercices bonus ## Comparaison des distances en projection, grand-cercle et sur ellipsoïde. Vous avez vu jusqu'ici que QGIS permet de mesurer des distances, des aires et même des angles soit en projection, soit sur l'ellipsoïde. La première solution est approximative, mais très rapide à calculer, tandis que la seconde est particulièrement couteuse en calculs. En réalité, tous les outils du type Géoportail, Google Maps et autres ne calculent pas les distances sur ellipsoïde mais sur sphère pour accélérer les calculs. On parle de distance grand-cercle, puisque le plus court chemin entre deux points sur une sphère est une portion de grand cercle. Le site https://www.movable-type.co.uk/scripts/latlong.html donne la possibilité de calculer ce type de distance en utilisant la formule la plus répandue nommée *Haversine*. Sous QGIS, placez 2 points et mesurez la distance sur ellipsoïde entre eux. Utiliser le calculateur de movable-type.co.uk pour effectuer le calcul sur sphère. Quel est l'écart entre les deux mesures? On peut aller plus loin dans l'analyse en choisissant un point à la surface du globe, puis en calculant l'écart entre les 2 distances entre ce point et une grille de points couvrant la planète. Utilisez QGIS et PostGIS pour réaliser cette analyse. ## Explorons le répertoire EPSG. Les sites Web [https://epsg.io](https://epsg.io) et [http://spatialreference.org](http://spatialreference.org) contiennent l'intégralité du répertoire EPSG. epsg.io fournit une description détaillée et conviviale de chaque SRC. Les deux sites permettent de récupérer la description de n'importe quel SCR sous de nombreux formats. 1. Allez sur le site [https://epsg.io](https://epsg.io) et recherchez le Lambert 93. A l'aide des informations de la page, répondez aux questions suivantes : - quel est son identifiant EPSG? - pour quel type de cartographie est-il le plus adéquat? - quelle est sa précision? - Cela a-t-il du sens d'utiliser Lambert 93 pour une carte de Grèce? 2. Sur le même site, cherchez la projection "WGS 84/pseudo mercator". Sur la page du SCR, lisez la remarque. Que pensez-vous de l'utilisation de ce SCR pour des mesures de précision? A votre avis, pourquoi les créateurs de ce SRC ont fait le choix d'utiliser une sphère plutôt qu'un ellipsoïde? ## Lecture d'une chaîne Proj4. Voici la chaîne proj4 du SCR Pseudo-Mercator : ``` +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs ``` A l'aide de la documentation de Proj4 (https://proj4.org/usage/projections.html), répondez aux questions suivantes: - Quel est le type de projection utilisée pour ce SCR? Cherchez sur Wikipedia quel effet elle a sur les angles et les surfaces des objets cartographiés. - Qu'est-ce qui vous indique que le datum de ce SCR est bel et bien une sphère? Voici maintenant la chaîne Proj4 du SCR Lambert 93 : ``` +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ``` - Que représente le paramètre `+x_0`, quel conséquence a-t-il sur les coordonnées projetées en Lambert 93? Cliquez sur la carte accompagnant la description de Lambert 93 sur le site epsg.io, et regardez où tombe la coordonnée $x=0$. A votre avis, pourquoi cette valeur de 700000m a été choisie par les géodésiens français?