TP7 — Opérations Raster

Au cours de ce TP, tu vas te familiariser avec certains outils de géotraitement pour les données au format raster (=image). Tu travailleras avec des données satellites au format raster (provenant de la NASA), une couche thématique raster qui identifie les classes d’occupations du sol (créée par nos soins avec les données satellites de la NASA) et enfin un MNT (modèle numérique de terrain) que tu extrapoleras (grâce à la méthode de l’IDW “Inverse Distance Weighting”) à partir de données ponctuelles disponibles sur Swisstopo.

Ce TP n’aurait pas été possible sans les ressources listées ci-dessous:

• TP7 du cours “Géomatique et SIG” du Privat-docent Dr. Marj Tonini
• (Image de couverture) “Sardinia in satellite imagery” by M. Jurzyk

1a) Télécharge à cet hyperlien un projet ArcGIS qui contient déjà les couches nécessaires ainsi qu’une mise en page finale de la carte que tu pourras utiliser pour l’habillage.

1b) Ouvre le projet que tu viens de télécharger. Il contient une géodatabase regroupant les trois couches raster que nous allons utiliser :

• MapTicino1990_8classes.tif : occupation du sol en 1990 au Tessin avec une classification en 8 classes (“Forêt”, “Prés”, “Eau”, “Neige”, “Sols nus”, “Aires urbaines”, “Nuages”, “Ombres”) et calculées à partir des données de Landsat4.
• Landsat4_1990_194028_bx.tif : les 7 bandes spectrales individuelles collectées par le satellite Landsat 4 de la NASA (b1: Bleu ; b2:Vert ; b3: Rouge ; b4: Proche infrarouge ; b5: Infrarouge moyen -1 ; b6: Thermique ; b7: Infrarouge moyen -2).
• MNT_Ticino.tif : un modèle numérique de terrain du Canton Tessin. Nota bene : ce n’est pas ce MNT qu’il faudra utiliser pour la carte finale ! Cette couche nous servira uniquement comme base de calcul.

1c) Vérifie dans les propriétés de la carte que le système de référence est bien CH1903+ LV95. Tu peux maintenant commencer le TP.

1d) Prends le temps d’explorer les métadonnées des couches (on parle de métadonnées, mais il s’agit, dans ArcGis Pro, des propriétés de la couche) et réponds aux premières questions sur le quiz Moodle. Choisis bien les métadonnées de la bande 4 pour répondre aux questions.

Pour la suite du quiz Moodle, nous allons travailler sur les bandes spectrales qui apparaissent comme 7 couches séparées. Afin d’effectuer les opérations d’analyse d’image, il faudra les fusionner dans un fichier multibandes unique. Pour ce faire, il existe un outil nommé Composite Bands.

1e) Une fois le menu de l’outil ouvert, ajoute les fichiers qui correspondent à chaque bande en ordre croissant (de b1 à b7).

1f) Sauvegarde la nouvelle couche dans la géodatabase et affiche le résultat.

Solution – Outil composite band

2a) Pour commencer, il faut uniquement visualiser la couche multibandes que tu viens de créer, en décochant toutes les autres couches depuis la fenêtre Contents. Cela évitera qu’ArcGis ne doive rafraîchir toutes les couches à chaque fois que tu bouges sur la carte.

2b) Ouvre le menu de symbologie de la couche multibandes et essaye les différentes combinaisons de bandes (“Band”) dans les canaux (“Channel”) correspondant aux couleurs rouge, vert et bleu.

2c) Après t’être familiarisé.e avec le fonctionnement des bandes, essaye de reproduire les compositions du tableau ci-dessous, qui indique la composition recommandée pour les images Landsat.

Des opérations peuvent être effectuées sur les différentes bandes de façon à mettre en évidence certains éléments. Dans le menu “Analysis” en haut de la page, tu trouveras l’outil “Raster Functions” qui te permet de faire différentes analyse basées sur les images multibandes.

Une des fonctions qui peuvent être appliquées à un raster multibandes est le calcul du NDVI. Cet index de végétation normalisé (“Normalized Difference Vegetation Index”) met en évidence la couverture végétale sur le territoire. Il est très utilisé en agriculture et sylviculture pour connaître l’état de santé des plantes, comme le montre l’image ci-dessous.

Source à consulter pour plus d’informations : https://eos.com/blog/ndvi-faq-all-you-need-to-know-about-ndvi/

Sur ArcGIS le calcul de cet indicateur est automatique, et effectué grâce à la formule ci-dessous :

Où R désigne la réflectance spectrale dans la bande rouge (la bande 3 dans le TP) et PIR (ou NIR, Near InfraRed en anglais) indique la réflectance spectrale dans la bande du Proche-Infrarouge (la bande 4).

2d) Applique la fonction NDVI (en noir et blanc) et dans les propriétés de l’opération, indique comme bande visible la bande 3 (correspondant au rouge) et comme bande infrarouge la bande 4 (correspondant au proche infrarouge). Assure-toi de sélectionner l’option “Scientific Output” qui permet de maintenir les valeurs du NDVI scientifiquement correctes, comme décrit par la formule ci-dessus (entre -1 et 1, sinon les valeurs sont comprises entre 0 et 200, plus d’infos sur le site d’ArcGis Pro).

Remarque : Il est important de noter qu’après avoir sélectionné le rendu en noir et blanc, l’image apparaît en couleurs. Ceci est normal et est fait pour faciliter la visualisation du phénomène sans toutefois compromettre les valeurs calculées à partir du NDVI. Ce n’est, en revanche, pas le cas avec l’option “colorized NDVI” qui, elle, n’offre pas les valeurs scientifiques du NDVI (Elle affiche les valeurs de 0-255, plus d’infos sur le site d’ArcGis Pro).

2e) Explore désormais les autres rendus possibles dans le menu des fonctions raster, et réponds aux questions de la deuxième page du quiz Moodle.

Solution NDVI

2f) Sauvegarde bien la nouvelle couche du NDVI dans ta géodatabase ! Tu devras utiliser ce rendu à la fin du travail.

Le raster “MapTicino1990_8classes.tif” est une cartographie de l’occupation du sol au Tessin en 1990. Les données sont recueillies par le satellite Landsat 4 de la NASA, et ensuite analysées de sorte à créer 8 classes d’occupation du sol : “Forêt”, “Prés”, “Eau”, “Neige”, “Sols nus”, “Aires urbaines”, “Nuages”, “Ombres”.

La fonction qui permet de faire cette différenciation en se basant uniquement sur les données satellitaires (les 7 bandes spectrales vues plus haut) est par contre très sensible et engendre des erreurs.

L’erreur principale qu’on retrouve est due à la similitude entre les longueurs d’ondes émises par les sols nus et les sols anthropiques (les zones urbaines). Ainsi, on retrouve des pixels classés comme “aires urbaines” en haute montagne à la place d’une classification comme “sols nus”.

Pour corriger ces erreurs, il existe dans ArcGIS un outil permettant d’effectuer des calculs et des modifications d’une image raster : l’outil “Raster Calculator” (dans “Spatial Analyst Tools”).

3a) Cherche l’outil “Raster Calculator” dans le paquet d’outil “Spatial Analyst Tools”.

⚠️ Fais attention à ne pas le confondre avec un autre outil du même nom mais contenu dans un autre paquet d’outils.

3b) Ouvre l’outil en question et essaye de trouver l’expression qui te permet de modifier les valeurs erronées dont on a parlé avant (lis d’abord la suite avant de t’y attaquer). En pratique, il faudra créer une requête permettant de modifier l’affectation des pixels classés comme “zone urbaine” et situés au-dessus de 1’400 mètres d’altitude en les classants comme “sols nus”.

⚠️ Utilise les variables ainsi que les opérateurs offerts par l’outil. Tu pourrais écrire toi-même le tout, mais il arrive que les symboles ne soient pas identiques, ce qui pourrait engendrer une erreur. En outre, si tu fais une erreur minime dans l’écriture des couches, l’outil produira un autre message d’erreur.

Voici quelques exemples de requêtes qui peuvent être menées :

• (“MapTicino1990_8classes” == 4) : sélection de tous les pixels de la couche “MapTicino1990_8classes” ayant une valeur égale à 4 (c.à.d. qui appartiennent à la classe “neige”).
• (“MNT25_Ticino” < 600) : sélection de tous les pixels situés à une altitude inférieure à 600 mètres dans le MNT “MNT25_Ticino”.

Les deux requêtes ci-dessus donneront comme résultat un nombre, puisqu’elles indiqueront combien de pixels répondent à la requête en question. Elles donnent aussi la position de ces pixels.

Il existe aussi des fonctions conditionnelles, qui s’appliquent uniquement dans le cas où une condition préalablement fixée est respectée.

• Con ( condition, valeur-si-vrai, valeur-si-faux) : cette requête appliquera la valeur “valeur-si-vrai” dans le cas où la condition initiale est respectée ; et appliquera la valeur “valeur-si-faux” dans le cas où la condition initiale n’est pas respectée. Il fonctionne de la même manière que le “=SI()” dans Excel.
• Con ( (“MNT25_Ticino” > 600) , 600 , “MNT25_Ticino”) : cette requête définit un plafond d’altitude à 600 mètres. La condition est “si le MNT25_Ticino présente une valeur supérieure à 600 mètres d’altitude, assigne la valeur de 600m” ; “si le MNT25_Ticino présente une valeur inférieure à 600 mètres d’altitude, laisse la valeur d’origine”.

Les requêtes ci-dessus ont une seule condition, mais on peut très bien indiquer plusieurs conditions à remplir avec le symbole “&”.

Con ( ( “MNT25_Ticino” > 600) & (“MapTicino1990_8classes” == 7) , 4 , “MapTicino1990_8classes” ) : cette requête modifie la classification des pixels 7 (“ombres”) en pixels 4 (“neige”) si l’altitude est supérieure à 600 mètres. Donc : “si le MNT indique une valeur supérieure à 600 m, et que le pixel appartient à la classes 7 “Ombres”, alors assigne la valeur du pixel à 4 (“neige”). Si une des condition n’est pas remplie, laisse la valeur d’origine”.

⚠️ C’est une requête de ce type qui te permettra d’effectuer le reclassement des pixels ayant la valeur 6 “Aires urbaines” et situés au-dessus de 1’400 mètres d’altitude en pixels avec valeur 5 “Sols nus”.

3c) Effectue la requête, puis, une fois la requête effectuée, sauvegarde le résultat dans une nouvelle couche que tu nommeras avec la mention “ReClass” (ex. ReClass_MapTicino1990_8classes).

Il se peut qu’il ne te sorte que 2 classes en sortie au lieu des 8: dans ce cas, relance l’outil (en changeant le nom dans “Output raster”) jusqu’à obtenir le total des classes. Garde cette couche dans ta géodatabase car elle te sera redemandée plus tard pour le rendu final.

Solution reclassement

Si tu ne visualises pas la vidéo ci-dessous, active le mode plein écran si tu es sur MacOS, ou tourne ton appareil si tu es sur iOS ou iPadOS. La vidéo est toujours visible depuis Windows (sur les machines virtuelles).

Étant donné la qualité de la vidéo, on t’offre un zoom sur la condition qui a été utilisée. Essaye toutefois de l’écrire toi-même avant de regarder la solution.

Solution Requête

Bravo, tu as presque fini ! Dans cette dernière partie, tu créeras un MNT à partir de données ponctuelles disponibles sur le site de l’Office fédéral de topographie (Swisstopo). Il s’agit du modèle numérique de base du terrain de la Suisse, utilisé pour la production du MNT avec une maille de 25 m. Pour ce faire, on se basera sur une méthode déterministe classique, à savoir l’interpolation par moyenne mobile pondérée ou Inverse Distance Weighting (IDW).

Les phénomènes spatio-continus sont définis en tout point de l’espace géographique (ex. l’altitude et la température) mais sont généralement étudiés à travers des données ponctuelles. Entre les points d’échantillonnage, les valeurs de ces phénomènes ne sont pas mesurées. L’objectif des méthodes d’interpolation consiste à prédire ces valeurs inconnues sur la base de l’autocorrélation spatiale :

« Deux objets proches ont plus de chance [d’interagir] que deux objets éloignés » (première loi de la géographie de Waldo Tobler).

Pour qu’une modélisation soit satisfaisante, il est primordial qu’elle soit basée sur une analyse exploratoire des données et sur une analyse des erreurs (quelle que soit la méthode d’interpolation choisie).

Le but de ce TP est d’interpoler assez rapidement (au détriment de la qualité du résultat) une surface raster à partir de données ponctuelles de l’altitude.

4a) Pour télécharger les données, rends-toi sur le site de swisstopo et appuie sur « DHM25 – Modèle de base ESRI Shapefile » qui se trouve sous l’onglet Géodonnées et applications > Modèles d’altitude > MNT25.

4b) Ensuite, importe le fichier shape « dhm25_p » dans ton projet (sans oublier de dézipper le fichier, deux dossiers en seront extraits).

Pour prédire l’altitude en tout point de l’espace géographique Suisse on utilisera l’outil IDW.

4c) Cherche cet outil dans la barre Find Tools ou trouve-le en cliquant Tools > Geostatistical Analyst Tools > Interpolation > IDW.

4d) Grâce à cet outil, estime le MNT à partir des données Swisstopo.

⚠️ Fais bien attention à choisir le bon champ à interpoler, i.e., l’élévation Shape.Z et non le code identifiant la commune OBJECTID dans “Z value field”. Il se peut que l’outil affiche une erreur (un thème récurrent avec ce logiciel), dans ce cas : sauvegarde le projet, quitte ArcGis Pro et relance le projet.

Solution

4e) Finalement, génère une représentation 3D de la surface du MNT grâce à l’outil Hillshade (Ombrage) dans Raster Function.

╔═╗
║═╬╦╦═╦═╦═╗ ║
╠═║║║╬║╩╣╔╝ ║
╚═╩═╣╔╩═╩╝. . ═

6a) Tu peux d’ores et déjà soumettre les trois rasters au format .pdf sur Moodle: Rendus_Cartes_TP7. N’oublie pas le format du rendu : nom_prenom_cartes_TP7.pdf (tu peux consulter la grille d’évaluation sur Moodle).

6b) Crée un paquetage de projet à partager en suivant les instructions à ce lien.

N’oublie pas d’inclure un résumé (summary) et des balises ou mots-clés (tags) contenant les métadonnées identifiées précédemment. Tu peux cliquer la case “Share outside of organization” et le bouton “Analyze” pour optimiser ton paquetage pour un partage public.

Il se peut que l’historique (History Items) de ton projet crée une erreur à l’analyse ou pendant le paquetage, tu peux désélectionner cette option si c’est le cas.

Solution

6c) Une fois le projet paqueté, vérifie que le fichier .ppkx que tu viens de créer te permet bien de réouvrir le projet avec la géodatabase complète (par exemple en double-cliquant dessus).

Solution

6d) Copie le fichier .ppkx de la machine virtuelle sur ton OneDrive et crée un lien de partage.

6e) Tu peux maintenant te rendre à nouveau sur Moodle pour soumettre à ce lien : Rendus_projet_TP7. N’oublie pas le format du rendu : nom_prenom_TP7.ppkx

Félicitations pour avoir terminé le TP et à la semaine prochaine pour l’analyse de distance !