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 vas extrapoler (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”).
• 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. NB. !! ce n’est pas ce MNT qu’il faudra utiliser pour la carte finale !! Cette couche nous servira uniquement comme base de calcul.

1c) Prends le temps d’explorer les métadonnées des couches 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.

Les bandes spectrales 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 Band”.

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

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

2a) Pour commencer, il faut uniquement visualiser la couche multibandes que tu viens de créer, en décochant toutes les autres couches depuis le bandeau “Contents”.

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 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, indiques 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 de l’NDVI scientifiquement correctes, comme on les a vues précédemment.

Remarque : Il est juste que malgré avoir sélectionné le rendu en noir et blanc, l’image apparaît en couleurs. Ceci est fait pour faciliter la visualisation du phénomène sans toutefois compromettre sa scientificité, contrairement à l’option “colorized NDVI” qui n’offre pas les valeurs scientifiques.

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.

2f) Sauvegarde bien la nouvelle couche du NDVI dans ta géodatabase ! Tu devra ajouter 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 passant uniquement sur les données satellitaires (les 7 bandes spectrales qu’on a vu avant) est par contre très sensible et engendre des erreurs.

L’erreur principale qu’on retrouve est du à 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”.

Il existe dans ArgGIS un outil permettant d’effectuer des calculs et des modifications d’une image raster : l’outil “Raster Calculator” (in “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. 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 soit 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éléction de tous les pixels de la couche “MapTicino1990_8classes” ayant une valeur égale à 4 (c.à.d. qui appartient à 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 ou la condition initiale est respectée ; et appliquera la valeur “valeur-si-faux” dans le cas ou la condition initiale n’est pas respectée.
• Con ( (“MNT25_Ticino” > 600) , 600 , “MNT_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, applique une 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 si-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 change 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) Une fois la requête effectuée, sauvegarde le résultat dans une nouvelle couche que tu vas appeler avec la mention “ReClass” (ex. ReClass_MapTicino1990_8classes). Garde cette couche dans ta géodatabase car elle te sera redemandé plus tard pour le rendu final.

Bravo, tu as presque fini ! 🥳 Dans cette dernière partie, tu vas créer 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 va se baser 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 (quelque soit la méthode d’interpolation choisie).

Le but de ce TP est d’interpoler très 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 (tu devras dézipper deux dossier).

Pour prédire l’altitude en tout point de l’espace géographique Suisse on va utiliser 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”.

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

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

Il ne te reste plus qu’à rendre tes résultats sur Moodle. Pour ce faire, nous avons préparé trois layouts : NDVI, occup_sol, et IDW, avec lesquels tu peux exporter les rasters respectifs au format .pdf.

4f) Pour chaque layout, fais correspondre l’information présente sur la carte à son titre sans oublier le fond de carte. 

Conseil : Pour modifier la couche affichée dans une mises en page (Layout), il te suffit d’activer la carte et sélectionner, désélectionner ou ajouter (Map > Add Data) les bonnes couches dans la fenêtre contenu (Contents).

4g) Pour la carte concernant l’occupation du sol, insère une légende qui illustre les 8 classes.

4h) Finalement, pour chaque carte, indique l’auteur.e en utilisant ton nom, ton prénom, ainsi que les sources des données.

4i) N’oublie pas de soumettre les trois rasters au format .pdf sur Moodle: Rendus_TP7.

Félicitations pour avoir terminé le TP et à la semaine prochaine pour (ré)apprendre les bases de Python utiles aux SIG 🐍