INTRODUCTION

Pour mener à bien ce rapport, nous avons tout d’abord arrêté le choix de notre territoire. Après une première sélection, de la Réserve naturelle nationale de l’estuaire de la Seine, se révélant trop pauvre en termes de géodiversité, notre choix s’est porté sur le territoire de l’Opération Grand Site “Falaises d’Étretat – Côte d’Albâtre”. Le Grand Site, situé dans le département de la Seine Maritime en Normandie, s’étend sur 13 communes, de Fécamp à Saint-Jouin-Bruneval, avec Etretat en son centre. Ce territoire possède une mosaïque de paysages, des villages habités au plateau agricole et verdoyant, des hautes falaises de craies blanches entrecoupées de valleuses typiques de ce littoral, qui offrent parfois des accès aux plages de galets. Ce tout indissociable forme l’esprit des lieux de ce Grand Site. Ce «monument falaise» peut culminer à près de 100 m d’altitude. Le site est composé de 3 sites classés et de plusieurs propriétés du Conservatoire du littoral, intégrées aux espaces naturels sensibles. Celui-ci s’est ainsi avéré offrir un degré suffisant de géodiversité, tout en pouvant être analysé par la télédétection grace à une seule image satellite, de par la largeur de fauchée des satellites utilisés.

Nous allons maintenant présenter notre méthode, ainsi que nos résultats obtenus, et l’analyse dans le temps, que nous avons pu y porter.

Le présent rapport est rédigé avec le logiciel Quarto

PREMIERE PARTIE : TRAITEMENT DE L’OCCUPATION DU SOL PAR TELEDETECTION

La télédetecion est l’ensemble des techniques mises en œuvre à partir de satellites, qui ont pour objet d’étudier la surface de la Terre, en utilisant les propriétés des rayonnements électromagnétiques (REM), émis ou réfléchis par les différents objets observés. Le principe fondamental ici, est qu’il est possible d’établir une correspondance entre le REM mesuré et la nature de l’objet observé. Le concept clé en est la signature (l’identité) spectrale.

A) Méthodes de sélection et de traitement de l’image employées et premier niveau d’analyse

1. Outils utilisés

Pour mener ce travail à bien, nous avons utilisé la version 3.34.12 du logiciel QGIS, support à ce travail, couplé à la version 2.12.99 de l’extension d’OrfeoToolbox Processing provider (OTB), qui permet de conduire le traitement. OTB a été utilisé en lieu et place de SCP, après une impasse technique découverte pendant les cours, alors que la classification non supervisée du territoire pour l’image la plus récente avait été réalisée. Les sites utilisés pour trouver et télécharger les images satellites ont été l’EO Browser, de Copernicus, pour 2024, et CREODIS pour 1995 et 2001.

2. La recherche de données satellites disponibles et pertinentes

Il s’agira d’effectuer une analyse de trois images satellite, récupérées sur un site internet en libre accès, doté de fonctionnalités suffisantes pour le bien de notre recherche (voir “Outils utilisés”, Infra). Ces images devront être capturées chacune par le satellite ayant la meilleure qualité d’image possible aux dates choisies et les images devront présenter une absence de couvert nuageux. Si c’est possible, trouver une seule image pour couvrir tout le territoire d’étude et ainsi éviter des fusions de dalles sur QGIS, qui s’avèreraient lourdes à traiter. Ces images devront être réparties entre les années 1990 et 2020, à des périodes saisonnales les plus rapprochées possibles (pour une apparence de sol la plus approchante possible (notamment pour les champs cultivés et les espaces boisés)). En l’occurrence, les dates des images trouvées pour respecter ces critères sont celles du 19 septembre 2024 (capturée par le satellite Sentinel-2), du 25 août 2001 (par Landsat 7) et du 1er août 1995 (par Landsat 5).

Une fois nos images choisies en fonction de ces différents critères, il faudra les télécharger dans un format travaillable (géoTIF et non JPEG, comme le propose EOS landviewer, par exemple), les dézipper, en charger tous les fichiers correspondant à une bande spectrale sur QGIS, puis les fusionner pour créer un stack raster. Pour celui-ci, il convient de choisir un jeu de 3 bandes, dans « symbologie », pour obtenir le meilleur rendu possible des différentes composantes terrestres et maritimes du territoire.

Le but de ce travail, sera de définir les principaux modes d’occupation du sol (MOS) pour ce territoire, discriminés par une réflectance différente des rayons solaires qui sont renvoyés au capteur du satellite, et dont la signature spectrale correspond à une longueur d’onde ; puis d’en étudier l’évolution à travers le temps.

Pour cela, nous effectuerons ce que l’on appelle un “traitement de haut niveau” en télédétection, qui consiste en la conduite d’une classification non supervisée, puis supervisée, de chaque image satellite choisie, jusqu’à parvenir à un résultat, en termes de précision, établi à l’aide d’un outil appelé “matrice de confusion”.

3. L’image du 19 septembre 2024 : sélection et présentation intégrale de son traitement

Tout d’abord, créons une image correspondant au visible (couleurs naturelles) en utilisant les couleurs correspondant au rouge, vert et bleu (RVB), pour Sentinel-2, soit les bandes 4, 3 et 2. Pour chaque date, c’est l’image satellite avec la combinaison des bandes dans le visible, qui a été choisie.

Figure 1 : Caractéristiques des bandes spectrales de Sentinel-2

Figure 2 : Stack raster de l’image satellite du 19 septembre 2024 (bandes RVB), Satellite Sentinel-2, résolution 10m.

Maintenant, créons un autre jeu de bande, pour faire une composition colorée avec de fausses couleurs, qui donne du contraste à l’image et fasse encore davantage ressortir les différents modes d’occupation du sol. Pour cela, combinons la bande 8 (proche infra-rouge, ou NIR), la bande 4 (le rouge), et la bande 3 (le vert). Le NIR et le rouge sont utilisés pour calculer l’indice de présence de végétation NDVI, et le NIR et le vert le sont pour l’indice de présence d’eau NDWI. Cette combinaison est donc très pertinente pour créer et analyser une composition colorée.

Nous obtenons l’image suivante :

Figure 3 : Composition colorée de l’image satellite du 19 septembre 2024 (bandes 8, 4, 3), Satellite Sentinel-2, résolution 10m.

Nous voyons, qu’en effet, la végétation ressort vivement (en rouge), ainsi que les aires urbaines et l’eau en bleu sombre.

3.1 Calcul d’indices d’occupation du sol

Nous allons maintenant calculer trois indices : l’indice de végétation, l’indice d’eau, et l’indice du bâti, que nous avons identifiés comme pertinents au vu de notre territoire d’étude et qui sont des indices classiquement utilisés en télédétection. Pour ce faire, nous allons utiliser l’extension OrfeoToolbox Processing provider (OTB), sous QGIS, et sa fonction « radiometrics », qui calcule ces indices et produit les cartes correspondantes.

3.1.1 L’indice de présence de végétation (ou NDVI)

Figure 4 : Indice de présence de végétation sur le territoire (NDVI).

Comme écrit précédemment, la végétation ressort vivement en rouge, ce qui permet de bien l’identifier. Pour un beau rendu dans QGIS, aller dans « Symbologie », cliquer sur « Palette/Valeurs uniques », palette de couleur « Turbo » (la plus contrastée en couleurs), puis l’inverser (pour obtenir une échelle colorée continue, allant du bleu le plus foncé au rouge le plus vif), et enfin la classer, pour qu’elle se répartisse en fonction de la signature spectrale des bandes de longueur d’ondes choisies.

3.1.2. L’indice de présence d’eau (NDWI)

Figure 5 : Indice de présence d’eau sur le territoire (NDWI)

3.1.3. L’indice de présence de bâti (BuiltUpISU)

Figure 6 : Indice de présence de bâti sur le territoire (BuiltUpISU)

Le rendu des deux derniers indices diffère très peu. C’est OTB qui a effectué les calculs et a produit les cartes.

3.2. Classification non-supervisée (K-means)

Avant de se lancer dans une classification supervisée de sa zone d’études, il est souvent utile de commencer par une classification non supervisée. Cette approche, un peu “en aveugle”, permet de se faire une idée sur la façon dont nous pourrons différencier les usages du sol. Concrètement, l’utilisateur stipule simplement un nombre de classes désirées et laisse la main à l’algorithme pour « faire au mieux », le but étant de minimiser la variance intra-classe et de maximiser la variance inter-classes (plus de cohésion en interne versus plus de contraste avec l’externe).

L’interprétation des classes, sorties automatiquement, peut se faire de façon empirique en comparant la classification avec une composition colorée ou une image à plus haute résolution spatiale. Mais l’interprétation peut également se faire à l’aide des signatures spectrales moyennes calculées pour chaque classe. Nous avons ainsi choisi une image satellite (raster multi-bandes) et créé 10 classes, dans lesquelles se sont réparties les signatures de bandes spectrales des objets au niveau du sol, classes auxquelles nous nous sommes efforcés d’attribuer la couleur la plus proche possible du rendu réel.

Figure 7 : Résultat de la classification non-supervisée (K-means).

Nous pouvons maintenant afficher la signature spectrale des objets des classes retenues in fine, pour la suite et la fin de notre traitement d’image dans le cadre de la classification supervisée, à savoir les champs, l’eau, les zones boisées et les zones bâties. Celles-ci correspondent, en effet; à la très grande majorité des modes d’occupation du sol du territoire de l’étude.

Figure 8 : Signature spectrale des quatre classes choisies.

L’eau et le bâti apparaissent respectivement dans la bande bleue et dans la bande verte, alors que la végétation (bois et champs) s’inscrit dans les bandes rouges et proche infra-rouge (NIR). Nous observons des chevauchements entre les signatures de bandes spectrales, ce qui indique que le niveau de réflectance des classes “zone batie” et “zone boisée”, ainsi que les classes “zone boisée” et “champs” peuvent se confondre et envoyer des résultats erronés.

Par conséquent, il va convenir, dans le cadre de la classification supervisée que nous allons maintenant présenter, de bien délimiter nos polygones d’entraînement et de validation sur des parties contrastées de territoire d’étude, pour obtenir le résultat le plus précis possible.

3.3. Classification supervisée

3.3.1. Le pré-traitement

La procédure avec OTB se passe en deux temps, d’abord une phase d’apprentissage basée sur des régions d’intérêts puis une phase d’application du modèle. Ces deux temps s’effectuent avec deux modules différents : trainImagesClassifier et ImageClassifier. La première étape est de digitaliser dans Qgis, sur la zone d’emprise, une couche d’apprentissage de type polygones, ainsi qu’une couche de validation de type polygones également. Ces polygones doivent être tous géographiquement distincts. Nous faisons également en sorte qu’il n’y ait pas un grand déséquilibre de taille de polygones (et donc de pixels en leur sein) entre les différentes classes. En l’occurrence, nous avons créé quatre macro-classes, avec 5 polygones à chaque fois, à savoir les champs, l’eau, les zone boisées et les zones bâties. Ces prétraitements présentent les avantages de pouvoir équilibrer le nombre de pixels par classe, grace à la smallest class strategy, ainsi que d’extraire dans une table les signatures radiométriques de chaque échantillon d’apprentissage.

Figure 9 : Exemples de polygones d’entraînement digitalisés sur le territoire, avec le résultat de l’échantillonnage d’entraînement

Chaque rond coloré apparaissant sur un polygone représente un pixel. L’ultime étape de ce prétraitement est d’associer, pour les pixels d’échantillons sélectionnés, les valeurs sous-jacentes du raster à classifier. Cela se fait à l’aide du module sampleExtraction. Cette étape est intéressante, car elle va nous permettre de tracer les signatures radiométriques de chaque polygone d’apprentissage et les signatures radiométriques moyennes de chaque classe.

L’étape suivante est de « faire apprendre » à un modèle à reconnaître les classes identifiées à partir des bandes spectrales, c’est-à-dire du raster multi-bandes, à classifier. Pour chaque classe, l’algorithme va extraire une sorte de signature radiométrique moyenne et va apprendre à la reconnaître. Pour cet apprentissage nous utilisons le module TrainImagesClassifier de OTB, que nous avons conduit selon deux algorithmes différents (Random Forest et SVM), pour pouvoir comparer les résultats de nos matrices de confusion, qui seront obtenues à la phase de post-traitement de cette classification supervisée.

Enfin, La phase d’application se fait à l’aide du module ImageClassifier. Nous obtenons un raster au format .tif présentant notre classification en 4 classes, selon les algorithmes Random Forest et Support-vector machine (SVM). Nous attribuons à chaque classe une couleur lui correspondant, afin d’obtenir l’image la plus réaliste possible (champs = vert ; mer = bleu ; zones bâties = gris ; zones boisées = vert très foncé).

Figure 10 : classification en 4 classes, avec l’algorithme “Random Forest”

Figure 11 : classification en 4 classes, avec l’algorithme “SVM”

La différence notable, entre ces deux cartes, est que l’algorithme SVM fait apparaître davantage de zones en vert franc, voire foncé, que Random Forest. Nous verrons, à la lecture des matrices de confusion, si cela a une incidence sur les résultats.

3.3.2. Le post-traitement

Concrètement, une validation de classification résulte de la comparaison entre le raster de classification et des polygones de contrôle digitalisés manuellement par l’utilisateur. OTB propose un module dédié au calcul d’une matrice de confusion suite à une classification supervisée. En entrée, il est logiquement nécessaire de disposer d’un raster de classification d’occupation du soL et de polygones de validation. OTB va confronter ces deux données pour calculer une matrice de confusion qui renseignera sur la qualité de la classification. Nous allons présenter, en les commentant, les deux matrices de confusion obtenues (soit deux par date étudiée pour l’image du 19 septembre 2024), chacune avec un des deux algorithmes retenus pour les déterminer (“Random Forest” et “SVM”).

3.3.2.1. Avec Random Forest

Figure 12 : Matrice de confusion avec Random Forest

3.3.2.2. Avec SVM

Figure 13 : matrice de confusion avec SVM.

Au final, les différents indicateurs, qu’ils soient quantitatifs ou qualitatifs, sont unaninements “presque parfaits”, avec un indice Kappa de 0,9776 pour la matrice réalisée à l’aide de Random Forest, et 0,9994 pour celle produite avec SVM, qui considère encore davantage de “vrais positifs” que Random Forest (ce qui est normal, puisque ce dernier algorithme est spécialisé sur les questions de végétation et les valide de manière plus pointue). Le résultat est vraiment très bon, s’élevant entre 92,91% et 100% de bonnes classifications, et ce quel que soit l’indicateur utilisé. D’ailleurs, si nous effectuons la moyenne des indicateurs Kappa des deux algorithmes, il s’élève à 0,9885, soit très au-dessus des 60 à 70% de bonnes prédictions fixés pour ce dossier. Intéressons-nous tout de même aux quelques erreurs de classification, identifiées par nos matrices de confusion. Si l’on considére la première, faite avec Random Forest, qui fait apparaître le plus d’erreur (même si tout reste très relatif à ce niveau de résultats), sur les 1638 pixels que nous avons classés comme étant situés dans des champs, 53 se sont trouvés être en zone bâtie et 63 en zone boisée. Cela peut aisément s’expliquer, au regard des cartes présentées plus haut, car, dans le premier cas, de la végétation est présente en ville et a dû se retrouver comprise dans un polygone que nous avons dessiné en milieu urbain ; alors que dans le second cas, des pixels classés comme faisant partie d’un champ, ont été validés comme étant situés dans une zone boisée. Ici, un ou plusieurs polygone a pu déborder sur une zone boisée, à moins que certaines parties d’un, ou de plusieurs champs, aient renvoyé une réflectance ayant une signature spectrale s’apparentant à une zone boisée. Nous avons effectivement vu, plus tôt dans cette étude, que les signatures spectrales se chevauchaient entre ces classes.

4. L’image du 25 août 2001

SENTINEL-2 n’existant pas encore dans les années 2000, nous avons dû rechercher un autre satellite, donnant des images correctes, et étant passé au dessus de notre zone d’étude dans les années 2000. Notre choix s’est porté sur Landsat-7, dont les images sont constituées par 7 bandes spectrales, ainsi que d’une bande panchromatique, qui recouvre le spectre du vert jusqu’à l’infrarouge, en passant par le rouge.

Figure 14 : Bandes constituant l’image de Landsat 7, avec longueur d’onde et résolution.

Nous aurions souhaité trouver une image au milieu ou à la fin des années 2000, afin de bénéficier d’un bon espacement entre l’année 1995 et l’année 2024, et ainsi mieux mesurer la potentielle évolution des modes d’occupation des sols sur cette durée. Malheureusement, en image estivale (entre le 15 juin et le 30 septembre) et ce sans nébulosité au-dessus du territoire, nous n’avons pas trouvé plus récent, dans les années 2000, que cette date. Nous avions indentifié une image datant du 14 juillet 2006 ; cependant, même avec uniquement 7% de taux de nébulosité, plusieurs nuages se retrouvaient au-dessus du territoire ce qui nous a posé le double problème du nuage et de son ombre portée. Après avoir cherché une méthode pour supprimer ces nuages et masquer leur ombre, à l’aide des extensions GRASS et Cloudmasking, utilisables sur Qgis, nous avons décidé de nous mettre à la recherche d’une image vierge de tout nuage. Cela a été assez fastidieux, ayant recherché une image lors d’un été, en partant de 2010, puis rétrochronologiquement, pour arriver à trouver celle-ci, datant de 2001. L’image présentée est un stack raster combinant les bandes 3, 2 et 1, afin d’être, comme pour l’image de 2024, dans le domaine du visible, avec des couleurs réelles. Il en sera de même pour notre troisième et dernière image, ce qui aidera à les comparer.

Figure 15 : Image satellite prise par le satellite LANDSAT 7, le 25 août 2001.

Sans reprendre tout le processus du traitement d’image sur OTB et Qgis, quel résultat obtenons-nous, dans notre classification supervisée, avec nos deux matrices de confusion réalisées à l’aide des algorithmes “Random Forest” et “SVM” ?

4.1. Matrice de confusion avec Random Forest

Figure 16 : Matrice de confusion avec Random Forest

4.2. Matrice de confusion avec SVM

Figure 17 : Matrice de confusion avec SVM

L’indice de précision Kappa, était, en 2024, de 97% pour la matrice de confusion Random Forest et de 99% pour la matrice de confusion SVM. Avec notre image satellite de 2001, il n’est plus que de 83% pour la première, et 88% pour la seconde. Ceci représente tout de même un très bon résultat, puisqu’étant toujours considéré, et ce dans les deux cas, d’“accord presque parfait”. Nous retrouvons à peu près le même écart dans l’indice Kappa, entre les deux matrices, qu’en 2024, et dans le même sens. La confusion dans la classifications des champs est la même, pour la même raison (Random Forest est plus pointu dans sa validation de zones végétalisées, certains pixels devenant des zones bâties ou boisées au moment de leur validation. En revanche, la nouveauté réside dans la grande baisse de la qualité producteur pour la zone bâtie, qui chute à 66% dans la matrice Random Forest, alors qu’elle était de 95% en 2024. Cela semble avoir deux facteurs explicatifs, à savoir une urbanisation qui s’est accrue en 23 ans (notamment les zones pavillonaires en entrée/sorties de ville), ainsi qu’une moindre qualité de l’image utilisée, du fait d’un satellite moins fin et peut-être de notre ignorance, du fait d’un manque de temps d’autoformation, à corriger les images Landsat pour obtenir la réflectance “Top of atmosphere” (TOA), qui rendrait, sans doute, la carte plus nette, lisible et plus aisée à travailler.

5. L’image du 1er août 1995

En 1995, Landsat 7 n’avait pas encore été lancé. Le plus proche satellite alors existant était Landsat 5. Nous en voyons effectivement la proximité dans le tableau ci-dessous. Il manque la bande panchromatique numéro 8 ; à part cela, les autres bandes, leur longueur d’onde et leur résolution sont identiques.

Figure 18 : Bandes constituant l’image de Landsat 5, avec longueur d’onde et résolution.

Figure 19 : Image satellite prise par le satellite LANDSAT 5, le 19 septembre 2024.

5.1. Matrice de confusion avec Random Forest

Figure 20 : Matrice de confusion avec Random Forest

5.2. Matrice de confusion avec SVM

Figure 21 : Matrice de confusion avec SVM

Ici, pour la première fois, la matrice Random Forest présente des résultats (légèrement) supérieurs à celle faite à l’appui de SVM. Ce qu’il ressort de ces deux matrices est la poursuite de la dégradation de la qualité utilisateur pour les champs (un peu plus de la moitié de classement validé) ; alors, qu’à l’inverse, la classification des zones boisées regagne les 10 points de précision sur perdus en 2001 dans la matrice Random Forest.

Le résultat pour les deux matrices reste tout à fait acceptable, pour le traitement d’une image datant d’il y 29 ans, prise avec un satellite équivalent deux générations de plus que le Sentinel-2.

B) Analyse globale des résultats : qualité des trois images satellites, hypothèses sur l’évolution des modes d’occupation du sol et prise de recul sur les résultats issus des matrices de confusion

Comparons maintenant les images satellites prises sur les trois années de l’étude, en les disposant au sein d’un carton.

Figure 22 : Images satellites prises pendant les étés 1995 (Landsat 5), 2001 (Landsat 7) et 2024 (Sentinel-2).

La comparaison entre ces trois images est visuellement éloquente. La différence de netteté et de précision est frappante entre les deux images plus anciennes et la plus récente, avec également des variations entre les deux plus anciennes. Pourtant, la saisonnalité, qui peut influer sur l’image, est globalement la même.

Les images plus anciennes ont été capturées par des satellites de génération antérieure, dotés de moins de diversité de lecture de bandes spectrales, d’une résolution spatiale passant de 10 mètres (pour le satellite SENTINEL-2), à 30 mètres, pour les satellites Landsat 5 et 7, et enfin d’une absence de correcteur atmosphérique intégré, contrairement à SENTINEL-2. La différence se voit, mais ce n’est peut-être pas là le seul facteur explicatif.

Nous n’avons pas eu le temps d’effectuer une recherche climatologique, pour savoir si les conditions étaient très différentes entre les trois étés. Toujours est-il que l’occupation du sol semble assez semblable, mais si elle est difficile à lire empiriquement, étant donné la mosaique de paysage représentée sur le territoire de l’Observatoire du Grand Site Falaises d’Etretat - Côte d’Albâtre.

Enfin, nous posons l’hypothèse de l’évolution de l’anthropisation de ce milieu naturel, en 29 ans, ce qui expliquerait les différences observées entre les années. En termes d’urbanisation, tout d’abord, avec une artificialisation accrue, par la construction de grandes zones pavillonnaires autour des villes ; d’autant que ce territoire, situé autour d’Etretat et avec Fécamp en son Nord, est très attractif. Cela peut également concerner les pratiques culturales, avec un champ qui était cultivé et qui ne l’est plus, ou qui ne fait plus pousser le même végétal. Dans les deux cas, la réflectance renvoyée évolue et peut indiquer un autre objet, ce qui influe, souvent négativement, sur les résultats d’une classification supervisée.

Même si nos résultats sont très corrects et oscillent de 77 à 99% de précision, selon l’indice Kappa, nous avons souhaité tenter de porter une analyse sur les erreurs par omission ou commission qui ont été trouvées. L’anthropisation par l’Homme de son milieu en est une.

Par ailleurs, comme nous l’avons vu lors de la classification non-supervisée de l’image satellite de 2024, nous observons des chevauchements entre les signatures de bandes spectrales, ce qui indique que le niveau de réflectance des classes concernées (à savoir “zone bâtie” et “zone boisée”, ainsi que les classes “zones boisées” et “champs”), peuvent se confondre et envoyer le signal d’un objet erroné.

Par ailleurs, une résolution spatiale trois fois moindre en 1995 et 2001 implique moins de pixels à exploiter. Un pixel doit être occupé au moins à ¾ par un mode d’occupation du sol (MOS). Il faut des pixels petits, pour être adaptés aux objets terrestres simples. En effet, plus la surface du pixel est petite, plus la valeur radiométrique correspondante est fiable, localisée et donc représentative du MOS considéré. En ayant su ceci au début de notre travail, nous aurions dessiné des polygones d’entraînement et de validation plus grands, afin d’obtenir un nombre de pixel qui serait demeuré satisfaisant en 1995 et 2001.

Pour conclure, ce territoire est une mosaique de paysages, en évolution, et ce sans compter sur les aléas climatiques pèsent, et vont peser de plus en plus à l’avenir, que ce soit l’érosion des côtes, l’élévation du niveau de la mer avec les risques de submersion, les inondations continentales et l’érosion des sols agricoles. Cela va impacter les activités humaines et le peuplement d’un territoire, toujours aussi attractif, mais qu’il faudra peut-être davantage protéger, en l’aménageant et le règlementant. C’est le sens de l’existence de notre territoire d’études, l’Observation du Grand Site Falaises d’Etretat - Côte d’Albâtre, qui prône depuis plusieurs années son inscription en tant que Grand Site de France.

Après cette première partie consacrée à la télédétection, nous allons maintenant nous intéresser à la géocomputation d’une indice de géodiversité, grace au langage R.

SECONDE PARTIE : GEOCOMPUTATION D’UN INDICE DE GEODIVERSITE GRACE AU LANGAGE R

Définition du concept de géodiversité

Il résume la variété géologique, géomorphologique, pédologique et hydrologique (sous-indices) d’un territoire. L’idée est d’extraire automatiquement (par SIG) la géorichesse dans des mailles régulières.

Méthode employée

Afin de mener notre étude à bien sur ce territoire, nous en calculerons tout d’abord l’indice de diversité lithologique (shannon couches lithologiques), puis l’indice de diversité hydrologique (Shannon ordre de Strahler (l’ordre de Strahler varie en fonction de la présence de cours d’eau et d’affluents)), et enfin nous couplerons les deux, pour définir l’indice de géodiversité du territoire.

Outils utilisés

R studio, avec la version 4.4.2 de R, pour, grace au langage R, établir l’indice de géodiversité de notre territoire,

1. Préparation de l’analyse

1.1. Installation des packages, qui seront utiles, et leur activation

install.packages(“sf”, dependencies = TRUE) install.packages(“mapsf”, dependencies = TRUE) install.packages(“tmap”, dependencies = TRUE) install.packages(“viridis”, dependencies = TRUE) install.packages(“rmarkdown”) install.packages(“tinytex”)

library(tidyr) library(dplyr) library(sf) library(terra) library(tmap) library(vegan) library(viridis) library(rmarkdown) library(“tinytex”)

1.2. Recherche et téléchargement des données nécessaires à l’étude

Il s’agira d’utiliser la base de données “BD-TOPO” de la Seine-Maritime, pour sa couche “hydrographie”, ainsi que la BD-CHARM 50 (couche FGEOL) pour la lithologie et enfin télécharger un fichier en format “Shapefile”, pour délimiter l’emprise spatiale de notre aire d’étude.

1.3. Création de la grille et définition de l’emprise

cellsize <- 1000 # En mètres si EPSG:2154 emprise <- st_read(“data/FGEOL/OGS_Falaises.gpkg”, quiet = TRUE) View(emprise)

grid <- st_make_grid(x = emprise, cellsize = cellsize) %>% st_as_sf() %>% st_intersection(emprise) %>% mutate(ID = row_number())

plot(st_geometry(grid))

2. Calcul de l’indice de diversité lithologique

Import de la donnée

litho <-st_read(“data/FGEOL/GEO050K_HARM_076_S_FGEOL_2154.shp”, quiet = TRUE)

Intersection entre les mailles et la lithologie

litho_grid <- st_intersection(litho, grid) %>% mutate(surface_intersection = as.numeric(st_area(.))) # Surface après intersection

Agrégation par maille et lithologie

part_surface_litho <- litho_grid %>% st_drop_geometry() %>% group_by(ID, CODE) %>% summarise(surface_litho = sum(surface_intersection), .groups = “drop”)

Calcul des pourcentages de chaque lithologie par maille

part_litho <- part_surface_litho %>% group_by(ID) %>% mutate(surface_totale = sum(surface_litho), part_litho = (surface_litho / surface_totale) * 100) %>% ungroup() %>% select(ID, CODE, part_litho) %>% pivot_wider(names_from = CODE, values_from = part_litho, values_fill = 0)

Calcul de l’indice de Shannon

C’est l’indice principal utilisé en géodiversité.

part_litho$shannon_litho <- diversity(part_litho %>% select(-ID), index = “shannon”)

Jointure des résultats avec la grille

grid <- grid %>% left_join(part_litho, by = “ID”) %>% mutate(across(where(is.numeric), ~ replace_na(.x, 0)))

Visualisation cartographique

tmap_mode(“view”) tm_shape(grid) + tm_basemap(“OpenTopoMap”) + tm_polygons(fill = “shannon_litho”, fill.scale = tm_scale_intervals(style = “quantile”, n = 6, values = “viridis”), fill.legend = tm_legend(title = paste(“Diversité lithologique”)))

3. Calcul de l’indice de diversité hydrologique

3.1. Diversité des ordres de Strahler

cellsize <- 1000 emprise <- st_read(“data/FGEOL/OGS_Falaises.gpkg”, quiet = TRUE) f <-list.files(path = “data/dalles”, pattern = ‘.*\.asc$’, full.names = TRUE) r <- lapply(f, rast) mnt <- do.call(“merge”,r) mnt <- crop(mnt, emprise, mask = TRUE) writeRaster(mnt, “data/dalles/dallesOGS_Falaises.tif”) mnt <- rast(“data/dalles/dallesOGS_Falaises.tif”) plot(mnt)

strahler <- rast(“data/dalles/Strahler_OGSFalaises.sdat”)

Conversion de la grille en SpatVector pour extraction des valeurs raster

grid_vect <- vect(grid)

Extraction des valeurs du raster de lithologie par maille

strahler_values <- extract(strahler, grid_vect)

Calcul des proportions des classes lithologiques par maille

strahler_counts <- strahler_values %>% group_by(ID, Strahler_OGSFalaises) %>% # litho est la valeur du raster (ID de la lithologie) summarise(n_pixels = n(), .groups = “drop”) %>% # Nombre de pixels de chaque litho mutate(total_pixels = sum(n_pixels)) %>% mutate(part_strahler = (n_pixels / total_pixels) * 100) %>% select(ID, Strahler_OGSFalaises, part_strahler)

Restructuration en format large (pivot_wider)

part_strahler <- strahler_counts %>% pivot_wider(names_from = Strahler_OGSFalaises, values_from = part_strahler, values_fill = 0)

Calcul de l’indice de Shannon

part_strahler$shannon_strahler <- diversity(part_strahler %>% select(-ID), index = “shannon”)

Jointure avec la grille spatiale

grid <- grid %>% left_join(part_strahler, by = “ID”) %>% mutate(across(where(is.numeric), ~ replace_na(.x, 0)))

tm_shape(grid) + tm_basemap(“OpenTopoMap”) + tm_polygons(col = “shannon_strahler”, style = “quantile”, n = 6, palette = “viridis”, title = paste(“Diversité de Strahler”))

3.2. Diversité des plans d’eau

Renommer la variable ID pour éviter la confusion avec l’ID de grid

eau <- st_read(“data/BD_TOPO_76/BDTOPO_3-4_TOUSTHEMES_SHP_LAMB93_D076_2024-12-15/BDTOPO_3-4_TOUSTHEMES_SHP_LAMB93_D076_2024-12-15/BDTOPO/1_DONNEES_LIVRAISON_2024-12-00129/BDT_3-4_SHP_LAMB93_D076-ED2024-12-15/HYDROGRAPHIE/COURS_D_EAU.shp”, quiet = TRUE) %>% rename(ID_EAU = ID)

Intersection entre les mailles et les plans d’eau

eau_grid <- st_intersection(eau, grid) %>% mutate(surface_intersection = as.numeric(st_area(.))) # Surface après intersection

Agrégation par maille et plans d’eau

part_surface_eau <- eau_grid %>% st_drop_geometry() %>% group_by(ID, ID_EAU) %>% summarise(surface_eau = sum(surface_intersection), .groups = “drop”)

Calcul des pourcentages de chaque plan d’eau par maille

part_eau <- part_surface_eau %>% group_by(ID) %>% mutate(surface_totale = sum(surface_eau), part_eau = (surface_eau / surface_totale) * 100) %>% ungroup() %>% select(ID, ID_EAU, part_eau) %>% pivot_wider(names_from = ID_EAU, values_from = part_eau, values_fill = 0)

Calcul de l’indice de Shannon

part_eau$shannon_eau <- diversity(part_eau %>% select(-ID), index = “shannon”)

Jointure des résultats avec la grille

grid <- grid %>% left_join(part_eau, by = “ID”) %>% mutate(across(where(is.numeric), ~ replace_na(.x, 0)))

3.3 Diversité hydrologique totale

grid$shannon_hydro <- grid$shannon_eau + grid$shannon_strahler

tm_shape(grid) + tm_basemap(“OpenTopoMap”) + tm_polygons(col = “shannon_hydro”, style = “quantile”, n = 6, palette = “viridis”, title = paste(“Diversité hydrologique”))

4. Géodiversité totale

grid$geodiversite <- grid$shannon_litho + grid$shannon_hydro

tm_shape(grid) + tm_basemap(“OpenTopoMap”) + tm_polygons(col = “geodiversite”, style = “quantile”, n = 6, palette = “viridis”, title = paste(“Géodiversité”))