APPORT DES IMAGES LANDSAT À L'ÉTUDE DE L'ÉVOLUTION DE L'OCCUPATION DU SOL
DANS LA PLAINE DE SAÏSS AU MAROC, POUR LA PÉRIODE 1987-2018.

EL GAROUANI Abdelkader et AHARIK Kamal

Université Sidi Mohamed Ben Abdallah,
Faculté des Sciences et Techniques de Fès,
Route d'Imouzzer, B.P. 2202, FES, MAROC,
Tel : +21 25 35 60 80 14 - Fax : +21 25 35 60 82 14

e-mail : el_garouani@yahoo.fr ; https://orcid.org/0000-0001-5379-6469

Résumé

Cet article concerne la plaine de Saïss au Maroc et porte sur l'évolution de l'occupation et de l'utilisation des sols pour la période allant de 1987 à 2018. Cette plaine s'avère très importante au niveau économique pour le pays. La méthodologie adoptée comporte successivement le calcul d'indices spectraux à partir d'images Landsat (NDVI : Normalized Difference Vegetation Index, NDWI : Normalized Difference Water Index et NDBI : Normalized Difference Built-up Index), puis l'utilisation de l'algorithme de vraisemblance afin de réaliser quatre classifications thématiques pour les années 1987, 2003, 2014 et 2018. La précision globale de ces classifications est déterminée à partir de la matrice de confusion, et varie entre 83 et 87% ; le coefficient kappa est, pour les quatre années, supérieur à 0,80. Entre 1987 et 2018, les surfaces correspondant aux terres irriguées, aux oliviers et au milieu urbain, ont progressé respectivement de 123%, 136% et 115%. À l'inverse, les forêts, les parcours et les terres arables ont vu leur surface diminuer respectivement de 10%, 6% et 29%.

**Mots-clés** : occupation du sol, classification, images Landsat, plaine de Saïss, Maroc.

Abstract

The objective of this article is to study the evolution of land cover / use for the period (1987 - 2018). The study area concerns the Saïss plain which is a very important economic region of the country. The adopted methodology begins first with the calculation of indices (NDVI: Normalized Difference Vegetation Index, NDWI: Normalized Difference Water Index and NDBI: Normalized Difference Built-up Index) then the classification of 4 Landsat images from 1987, 2003, 2014 and 2018 with the maximum likelihood algorithm. The precision of image segmentation is determined by the kappa coefficient using the confusion matrix (the precision varies between 0.83 and 0.87). Areas of irrigated land, olive trees and urban areas increased respectively by 123%, 136% and 115% between 1987 and 2018. On the other hand, for areas of forest, rangelands and arable land, they decreased respectively by 10%, 6% and-29%.

**Keywords** : land use, classification, Landsat images, Saïss plain, Morocco.
# 1. Introduction

La plaine de Saïss est caractérisée par un développement agricole important. Ce développement nécessite des ressources croissantes en eaux de surface et souterraine. Vue l'aspect très irrégulier des précipitations dans la région, les ressources souterraines sont de plus en plus sollicitées. Les activités anthropiques croissantes (agricoles, industrielles, touristiques et domestiques) ont affecté l'état de ces ressources souterraines pour l'ensemble de la région. La piézométrie du complexe aquifère de Saïss est en diminution permanente. Cette diminution influence la quantité et la qualité des réserves en eau disponibles (Essahlaoui et al., 2001).

La recharge annuelle de l'aquifère présente le paramètre du bilan hydrique le plus sensible aux effets de la variabilité climatique. De plus, l'infiltration des eaux de surface qui alimentent la nappe phréatique, dépend de plusieurs éléments tels que, d'une part, les divers flux hydriques (précipitations, l'évapotranspiration, etc.) et d'autre part les caractéristiques bio-géophysiques de la surface (occupation des sols, pédologie, altitude, pente, etc. (Bateelan et De Smedt 2007 ; Richard, 1999 ; Taylor et al., 2013).

Actuellement, un projet de transfert d'eau (700 Millions de m³/an) depuis le barrage Mdez qui se trouve à 60 km à l'est de la plaine de Saïss est en cours de réalisation. Le but de ce projet est de fournir une alternative à la surexploitation actuelle de l'aquifère et des ressources en eau de surface (ABHS, 2016). La disponibilité d'une carte d'occupation du sol de la plaine à jour est donc d'une grande importance pour l'optimisation de l'usage et de la gestion des ressources en eau.

Les images de télédétection sont extrêmement utiles pour comprendre l'évolution spatio-temporelle de l'occupation du sol. La classification est l'une des activités crucales en matière de télédétection. Les procédures de classification tentent de réaliser l'interprétation d'une image numérique en s'appuyant sur la puissance de calcul informatique. La plus usuelle de ces procédures est la classification supervisée.

Classifier une image en télédétection consiste à rassembler tous les pixels d'une image en un nombre réduit de classes correspondant aux principaux millieux du paysage observé, pour aboutir à une cartographie thématique approchée (Lassale 2015 ; Le Hégarat, 2005). On distingue ainsi deux types de classification : supervisée et non supervisée.

La classification supervisée (par l'algorithme du maximum de vraisemblance, par exemple) nécessite un échantillonnage manuel effectué par l'utilisateur, grâce auquel il différencie les classes thématiques dans une région donnée (Bolstad et Lillesand 1991 ; El Garouani et al., 2010 ; Martin et al., 1998). Le résultat obtenu cartographie les classes d'occupation du sol retenues lors de l'échantillonnage. Cette opération facilite en la

simplifiant la compréhension des phénomènes naturels qui s'avèrent être en évolution continue dans l'espace et dans le temps (Blaschke 2010 ; Campbell et Wynne, 2011 ; Lillesand et al., 2014 ; Richards, 1999 ; Rodes, 2016).

Les méthodes de classification supervisée sont ainsi utilisées en télédétection pour la production de cartes d'occupation des sols et pour la reconnaissance d'objets sur les images satellitaires (Chandran et Christy, 2015 ; Ilenco et al., 2017). Au fil des années, plusieurs algorithmes de classification supervisées ont été développés. Citons, par exemple, les réseaux de neurones, les approches bayésiennes, les arbres de décision. (Matthieu, 2017 ; Lu et Weng, 2007).

L'objectif principal de nos travaux est de produire des cartes d'occupation des sols pour plusieurs dates (1987, 2003, 2014 et 2018) afin d'analyser les variations spatio-temporelles de cette occupation sur une période de 31 ans. Nous considérons que ce type d'analyse constitue un outil essentiel pour la planification du développement et pour la gestion du territoire concerné.

# 2. Matériel et Méthode

## 2.1 Zone d'étude

La plaine de Saïss est localisée dans le nord central du Maroc et s'étend sur une superficie de 2100 km² (Figure1). Notre zone test couvre une surface rectangulaire d'une longueur de 95 km et d'une largeur de 30 km. Cette superficie représente environ 6% de la superficie totale du bassin du Sebou. Le nombre d'habitants est d'environ 2 408 584 selon les résultats de recensement général de la population et de l'habitat (HCP, 2014).

La structure géologique des formations aquifères de la plaine a été étudiée grâce à des sondages électriques (Amraoui et al., 2003 ; Bargach et al., 2004 ; Chabli et al., 2005 ; et Essahlaoui et El Ouali, 2003).
<figure><img src="image_1.png" /><figcaption>**Figure 1** : Situation géographique de la région d'étude (le fond de carte correspond à une image en composition colorée (RVB = Bandes 6, 5, 4) de Landsat 8 du 29 Juin 2018. La couleur verte correspond à la végétation)</figcaption></figure>

Une activité sismique a affecté la région après le dépôt de formations quaternaires. Le résultat de cette activité est la division de la plaine de Saïss en deux parties (la plaine de Fès et le plateau de Meknès) et la création de fissurations au niveau du toit de la nappe profonde. Ces fissures ont créé des points d'échange entre la nappe profonde karstique et la nappe phréatique alluviale (Amraoui 2005 ; Chamayou et al., 1975). La plaine comporte 51621 ha (valeur en 2003) de terres irriguées, dont 93% sont alimentés par des pompages privés à partir de la nappe phréatique de Saïss. Les grandes surfaces irriguées sont concentrées dans la région d'El Hajeb (68%), Meknès (20%), Fès (9%), et sont irriguées à partir des eaux de surface et des eaux souterraines. Les exploitations agricoles adoptent des modes d'irrigation économiques, comme l'aspersion et le goutte à goutte. Malgré les efforts effectués pour généraliser ces nouvelles techniques d'irrigation, le mode gravitaire demeure le plus utilisé pour 57% des surfaces irriguées.

## 2.2 Données LANDSAT

Nous avons utilisé les images Landsat pour la cartographie de l'occupation des sols et du couvert végétal. L'analyse concerne quatre scènes pour les années 1987, 2003, 2014 et 2018. Ces images sont acquises à partir des sites GLCF (Global Land Cover Facility) et USGS (United State Geological Survey). Nous avons retenu des images acquises pendant la saison d'été (tableau 1). Ces images se caractérisent par une couverture nuageuse quasi nulle (le taux de la couverture nuageuse est indiqué dans le fichier de métadonnées de chaque image).

<table><caption>**Tableau 1** : Images Landsat utilisées</caption><thead><tr><td>Satellites</td><td>Dates d'acquisition</td><td>Bandes spectrales utilisées</td></tr></thead><tbody><tr><td>Landsat 5</td><td>07 Mai 1987</td><td>1, 2, 3, 4, 5</td></tr><tr><td>Landsat 7</td><td>27 Mai 2003</td><td>1, 2, 3, 4, 5</td></tr><tr><td>Landsat 8</td><td>02 Juin 2014</td><td>2, 3, 4, 5, 6, SWIR 1 et 2</td></tr><tr><td>Landsat 8</td><td>29 Juin 2018</td><td>2, 3, 4, 5, 6, SWIR 1 et 2</td></tr></tbody></table>
Les données des satellites Landsat-5, 7 et 8 utilisées correspondent à des images multispectrales de 30 m de résolution spatiale. Les images fournies ont subi les corrections radiométriques et géométriques partielles (sans corrections atmosphériques) grâce à l'usage d'un modèle numérique de terrain (MNT). Ce prétraitement correspond au niveau L1TP (USGS, 2019). Ensuite, un mosaïquage (figure 2) rendu nécessaire par le fait que la région d'étude est couverte par deux *Path-Row* différents (201-36 et 201-37). Enfin, l'extraction du secteur d'étude et la projection des images dans le système de coordonnées conforme conique de Lambert ont été effectuées.

<figure><img src="image_2.png" /><figcaption>**Figure 2** : Mosaïquage des scènes couvrant la région d'étude.</figcaption></figure>

Les deux scènes se trouvent sur la même trajectoire descendante du satellite Landsat et sont acquises le même jour, avec un décalage d'environ 25 secondes. L'assemblage des bandes spectrales identiques des deux scènes forme une seule image qui couvre toute la région d'étude, image pour laquelle les niveaux des bandes spectrales ont été corrigées afin d'éviter les discontinuités entre images.

## 2.3 Données d'apprentissage et de validation

Pour entrainer et valider les classifications, des photographies aériennes (réalisées en 1987), des images *Google Earth* (réalisées en 2003, 2014 et 2018) et des visites sur le terrain (effectuées en 2018) ont été utilisées. Des missions de reconnaissance ont été effectuées en 2018 pour l'identification des classes d'occupation du sol et des sites d'entrainement. Un GPS a alors été utilisé pour se localiser. L'expertise de terrain a essentiellement consisté à interpréter et à géoréférencer des parcelles de toutes les classes d'occupation du sol identifiées pour la classification. Huit classes d'occupation du sol ont été retenues (surface d'eau, urbanisme, parcours, terres irriguées, vergers, oliviers, forêts et terres arables) et cinq sites ont été visités pour chaque classe. Les sites visités sont des parcelles de quelques pixels de surface. Le choix

et le nombre des classes d'occupation du sol retenues sont imposés, d'une part par le contexte géographique sud-méditerranéen caractérisé par un climat semi-aride, et d'autre part par la résolution spatiale des images Landsat (30 m), qui ne permet pas d'identifier plus de classes (Figure 3). Les pixels de 900 m² peuvent donc présenter des hétérogénéités d'occupation du sol, et par conséquent, conduire à des confusions entre classes.

L'identification et la délimitation des parcelles d'apprentissage (sous forme de polygones contenant environ une dizaine de pixels) sur les compositions colorées des images, sont assistées par des photographies aériennes d'une échelle de 1/10 000 réalisées en 1987, par des images *Google Earth* pour les années 2003, 2014, et par des campagnes de terrain et des images *Google Earth* pour l'année 2018. Entre cinq et 15 sites d'apprentissage pour chaque classe ont été identifiés.

Après avoir numérisé les sites d'apprentissage, un fichier de classes spectrales correspondant à ces différents sites a été créé. Pour faciliter le choix des classes de réflectance, nous avons effectué une analyse graphique de la séparabilité des classes thématiques selon leurs signatures spectrales dans les différentes bandes (Figure 4). Cette analyse nous a permis d'identifier les classes présentant des confusions et nous a conduit à choisir un échantillonnage plus pertinent.

Après la réalisation des classifications, une campagne de terrain (vérification des sites aléatoires) pour la validation a été effectuée. Environ 100 points d'échantillonnage aléatoires ont été utilisés pour chaque image. La collecte de ces échantillons a été basée sur les photos aériennes, les images *Google Earth* et la reconnaissance de terrain effectuée (lors de cette campagne de validation).

## 2.4. Méthodologie

Afin de générer des cartes d'occupation des sols correspondant aux années 1987, 2003, 2014 et 2018, nous avons lancé le processus de classification présenté précédemment sur les images des années correspondantes. Ce processus a utilisé les trois indices NDVI, NDWI et NDBI, qui ont été calculés après traitement de correction des scènes brutes pour les transformer en réflectance ; cela dans le but de caractériser la surface indépendamment des conditions atmosphériques et d'ensoleillement. Cette opération de correction a été menée grâce au logiciel ArcGIS 10.6, développé par ESRI (Environmental Systems Research Institute).

- Indice NDVI (Normalized Difference Vegetation Index) : l'indice de végétation par différence normalisée, initialement proposé par Rouse et al. (1973), étroitement lié à la biomasse verte et à la surface foliaire, constitue l'un des trois indices spectraux les plus largement utilisés
pour la surveillance de l'agriculture. Des études antérieures ont montré que le potentiel de développement, le stress hydrique et le potentiel de production des plantes sont exprimés dans la réflectance spectrale des couvertures des cultures et que les conditions de croissance des cultures peuvent ainsi être quantifiées à l'aide de NDVI (Tucker, 1979 ; Jackson et al., 1986 ; Weigand et Richardson, 1990). Les valeurs de NDVI varient entre -1 et +1. Les valeurs négatives correspondent aux surfaces d'eau. Pour les terrains nus, l'indice est proche de 0 et pour les surfaces de végétation, il peut atteindre +1. Le NDVI permet d'identifier les surfaces construites, les terrains agricoles et la densité du couvert végétal.

- Indice NDWI (Normalized Difference Water Index) : le NDWI est un indicateur sensible au changement de la teneur en eau des feuilles (Gao, 1996). Le NDWI est calculé en utilisant les réflectances du proche infrarouge et de l'infrarouge moyen Les valeurs de NDWI varient entre -1 et +1. Les valeurs situées entre -1 et 0 correspondent aux surfaces brillantes sans végétation ni surface d'eau et celles qui se trouvent entre 0.3 et 1 correspondent aux surfaces d'eau (McFeeters, 2013). Cet indice est utilisé pour identifier les surfaces d'eau continues (McFeeters, 1996 ; Xu, 2006).

- Indice NDBI (Normalized Difference Built-up Index) : le NDBI, qui est basé sur les caractéristiques spectrales des zones bâties et de couvert végétal, a été largement utilisé pour extraire les zones bâties des images de télédétection. Le NDBI a été utilisé pour automatiser le processus de cartographie des zones bâties en intégrant la méthode d'Otsu (Tamim, 2015). Le NDBI a été calculé en utilisant la bande proche infrarouge (NIR) et la bande infrarouge à ondes courtes 1 (SWIR1) (Zha et al., 2003) : le résultat génère des valeurs entre -1 et +1. Cet indice est utilisé pour séparer les zones construites. Il permet de différencier les zones de constructions des terrains nus.

- Ces trois indices spectraux sont calculés grâce aux formules suivantes :

$$NDVI = \frac{NIR-Red}{NIR+Red}$$

$$NDWI = \frac{Green-NIR}{Green+NIR}$$

$$NDBI = \frac{SWIR-NIR}{SWIR+NIR}$$

Avec : Green : bande verte, NIR : bande proche infrarouge, Red : bande rouge et SWIR : bande moyen infrarouge.

Note : bien que les corrections effectuées pour calculer la réflectance reposent sur des hypothèses simplificatrices et

portent uniquement sur les effets supposés invariables sur chaque image, cette transformation permet de comparer dans de meilleures conditions les signatures spectrales des cibles selon leur position et leur évolution temporelle (Bru, 2015). Les informations nécessaires pour effectuer ces corrections sont disponibles dans des fichiers de métadonnées fournis avec les images.
<figure><img src="image_3.png" /><figcaption>**Figure 3** : Zoom sur les classes d'occupation des sols identifiées à partir de l'image Landsat8 de 2018 (composition colorée : RVB = Bandes 6, 5, 4) et leurs correspondances thématiques sur le terrain (photos)</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 4** : Test de la séparabilité des classes des sites d'apprentissage par nuage de points pour les bandes spectrales bleue et proche infra-rouge</figcaption></figure>

Le processus de la classification a commencé par l'analyse visuelle des compositions colorées, réalisées en attribuant à chaque couleur primaire (rouge, vert et bleu) trois bandes spectrales d'une image. Ceci nous a permis de faciliter la photo-interprétation et de mettre en évidence des phénomènes environnementaux. En pratique, les combinaisons suivantes ont été utilisées : composition des bandes 5-4-3 (composition colorée infrarouge ; la végétation apparait rouge) et composition des bandes 4-3-2 (composition en couleurs naturelles). Après la différenciation visuelle des classes, les sites d'entrainement ont été créés en comparant les compositions colorées et les trois indices spectraux normalisés calculés. L'étape suivante a consisté à produire la classification en utilisant l'algorithme du maximum de vraisemblance (MLC) (Bolstad et Lillesand, 1991). Cet algorithme repose sur deux principes :

* Les pixels de chaque classe dans l'espace multidimensionnel sont répartis selon une loi normale.

* Le théorème de Bayes guide la prise de décision.

L'outil tient compte à la fois des variances et des covariances des signatures de classe lors de l'attribution de chaque pixel à l'une des classes correspondant aux parcelles d'apprentissages (Shivakumar et Rajashekaradhyab, 2018). Il suppose que la distribution d'un exemple de classe est caractérisée par une loi normale, et que cette classe peut être identifiée par le vecteur moyen et la matrice de covariance. Compte tenu de ces deux caractéristiques pour chaque famille de pixels, la probabilité statistique est calculée pour chaque classe,

ce qui permet de connaître la probabilité d'appartenance des pixels aux diverses classes.

Lors de la première visite de terrain, les photographies aériennes et les images *Google Earth* permettent l'identification des sites pour effectuer les classifications. La résolution spatiale des images fournies par *Google Earth* dépend de l'endroit géographique. Cette résolution n'est donc pas homogène sur toute la terre. Sur notre site d'étude, elle est de 65 cm.

Sur l'image en composition colorée, des zones d'apprentissage sont définies. Ce sont des régions de l'image que l'on juge représentatives pour chaque classe. Le choix des régions d'apprentissage mérite beaucoup d'attention parce que pour une même classe, il peut y avoir une certaine variabilité spectrale, due par exemple au changement de la densité des arbres dans une forêt ou à la variation de la profondeur de l'eau au niveau d'un barrage. Dans ce cas, il faut, soit créer des sous-classes, ce qui augmente la complexité des calculs, soit définir au moins un site d'apprentissage par type de paysage.

Dans notre cas, pour les classes thématiques plus ou moins homogènes du point de vue spectral (surface d'eau et forêt), cinq sites d'apprentissage ont été identifiés. Pour les autres classes plus hétérogènes, quinze sites d'apprentissage ont été utilisés.

La Figure 5 présente de manière synthétique le processus de classification.

<figure><img src="image_5.png" /><figcaption>**Figure 5** : Processus de classification supervisée à partir des trois indices NDVI, NDWI, NDBI</figcaption></figure>
En raison de la résolution décamétrique de l'image (30 mètres) et du processus de traitement (chaque pixel est traité individuellement), certaines classifications de pixels par maximum de vraisemblance sont incorrectes, ce qui génère un bruit aléatoire. Certaines petites régions le sont également. Pour limiter ce phénomène et améliorer le résultat, nous avons procédé à une étape de post-classification. Un filtre dit « majoritaire » a été appliqué. Le filtre majoritaire corrige le classement obtenu en prenant en compte la distribution des classes dans leurs voisinages contigus et en attachant le pixel au groupe de plus grand effectif. L'utilisation de Quatre voisins permet de conserver les angles des zones anguleuses alors que l'utilisation de Huit voisins permet de lisser ces mêmes angles. Dans notre cas, le filtre majoritaire (3x3) est appliqué en utilisant les huit pixels les plus proches. Il requiert qu'au moins la moitié des valeurs (quatre des huit pixels) soient identiques pour y affecter la classe correspondante. (ESRI, 2020).

La carte résultant de l'application de cet outil de nettoyage des données semble plus esthétique (Figure 6), car lissée. Ce processus de filtrage qui supprime les pixels isolés (ou bruit) de la sortie de classification présente malgré tout un inconvénient. Il peut en effet avoir pour conséquence une perte d'information géographique, d'où la nécessité de bien contrôler visuellement et avec beaucoup de précaution l'emploi de ce type de filtre.

<figure><img src="image_6.png" /><figcaption>**Figure 6 :** Traitements de post-classification (A avant filtrage, B après filtrage)</figcaption></figure>

L'évaluation de la qualité de la classification est réalisée par le calcul de la précision globale et du coefficient Kappa de Cohen. Lorsque le coefficient de Kappa est compris entre 0.61 et 1, la classification peut être jugée comme bonne (Cheng et al., 2018; Mubako et al., 2018; Vázquez-Jiménez et al., 2018).

Le référentiel de Landis et Koch (1977) propose une interprétation du coefficient Kappa en fonction de sa valeur (Tableau 2).

<table><caption>**Tableau 2.** Référentiel pour interpréter la valeur de Kappa.
Source : Landis et Koch (1977).</caption><thead><tr><td>Valeur de Kappa</td><td>Interprétation</td></tr></thead><tbody><tr><td>1.00 - 0.81</td><td>Excellente</td></tr><tr><td>0.80 - 0.61</td><td>Bonne</td></tr><tr><td>0.60 - 0.41</td><td>Faible</td></tr><tr><td>0.20 - 0.00</td><td>Négligeable</td></tr><tr><td>&lt; 0.00</td><td>Mauvaise</td></tr></tbody></table>

La précision globale et le coefficient de Kappa sont obtenus à partir de la matrice de confusion qui compare les classes obtenues par classification automatique avec des points de référence dont les caractéristiques sont connues (Girard et Girard, 2010).

La matrice de confusion est un tableau à double entrée. Chaque ligne représente le nombre d'occurrences d'une classe de la donnée de référence (réelle), tandis que chaque colonne représente le nombre d'occurrences d'une classe prédite par le système de classification (Pelletier, 2017).

Les échantillons sur la diagonale représentent le nombre d'échantillons correctement prédit par l'algorithme de classification. De cette matrice de confusion, il est possible de calculer un ensemble de métriques caractérisant la performance de l'algorithme de classification utilisé. Les échantillons test utilisés doivent référencer la réalité sur le terrain.

La précision globale de la classification est égale au nombre de pixels bien classés (pixels de la diagonale de la matrice de confusion) par rapport au nombre total de pixels sondés. Afin de s'assurer de la représentativité de cet indice et d'affiner ainsi cette information, la précision de chaque classe prise séparément peut également être calculée (erreurs ommission et commission).

$$ \text{Précision globale} = \left( \frac{\text{Nombre total d'échantillons corrects}}{\text{Nombre total d'échantillons}} \right) * 100 $$

Le coefficient de Kappa (K) traduit la réduction de l'erreur commise lors de l'utilisation de la classification par rapport à l'erreur qui se produirait dans le cas d'une attribution aléatoire des classes. Il indique comment les données à classer s'accordent aux données de référence. Il constitue une mesure fiable dans l'évaluation des classifications thématiques, par l'examen de tous les éléments de la matrice de confusion, et prend en compte à la fois les erreurs d'omissions et de commissions (Rosenfiel et Filztpatrick-lins, 1986). Il se calcule de la manière suivante :
$$K = \frac{Tp * (Tp \text{ bien classés}) - \text{Somme}(Tl * Tc)}{Tp^2 - \text{Somme}(Tl * Tc)}$$

Avec :

Tp : Total pixels

Tl : Total lignes

Tc : Total colonnes

Des cartes de sites (réalité terrain) répartis de manière aléatoire sur toute la zone d'étude pour chaque date sont réalisées. Ces cartes vont servir de vérité terrain pour le calcul des matrices de confusion.

### 3. Résultats et discussions

L'analyse des images, pour les années 1987, 2003, 2014 et 2018, a permis d'établir une cartographie de l'évolution de l'occupation des sols de la plaine de Saïss (Figures 7 à 10).

Avant d'entamer l'analyse et l'interprétation de ces cartes, il est judicieux d'évaluer les résultats de la classification.

Cette validation s'est effectuée par la comparaison visuelle des résultats de la classification avec des données externes (Photographies aériennes, images Google Earth et la réalité du terrain) et par l'intermédiaire de l'analyse de la matrice de confusion (Cheng et al., 2017 ; El Garouani et al., 2017 ; Khatami et al., 2017).

En pratique, une carte de sites répartis de manière aléatoire sur toute la zone d'étude est élaborée. Environ 100 échantillons pour les huit classes sont établis en se basant sur les images de Google Earth. Cette carte, rastérisée par la suite, va servir de vérité terrain. Une opération d'intersection de cette couche avec celle issue de la classification a été réalisée pour obtenir une couche dans laquelle les points sont affectés à leur classe d'occupation du sol. A l'aide d'une requête SQL sous ArcGIS, les données de la matrice de confusion ont été générées. À titre d'exemple, la matrice de confusion pour l'image de 2018 est présentée par le tableau 3.

<caption>**Tableau 3 :** Matrice de confusion pour la classification de l'image de 29-06-2018.</caption>

Précision globale = 84% Coefficient Kappa = 0,80

<table><thead><tr><th>Classification /Référence</th><th>Surface d'eau</th><th>Urbanisme</th><th>Parcours</th><th>Terres irriguées</th><th>Vergers</th><th>Oliviers</th><th>Forêt</th><th>Terres arables</th><th>Total</th><th>Errore Omission</th></tr></thead><tbody><tr><td>Surface d'eau</td><td>20</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>21</td><td>0,05</td></tr><tr><td>Urbanisme</td><td>0</td><td>20</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>23</td><td>0,13</td></tr><tr><td>Parcours</td><td>0</td><td>0</td><td>14</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>17</td><td>0,18</td></tr><tr><td>Terres irriguées</td><td>0</td><td>0</td><td>0</td><td>27</td><td>3</td><td>1</td><td>0</td><td>1</td><td>32</td><td>0,16</td></tr><tr><td>Vergers</td><td>1</td><td>1</td><td>0</td><td>2</td><td>26</td><td>2</td><td>2</td><td>0</td><td>34</td><td>0,24</td></tr><tr><td>Oliviers</td><td>0</td><td>0</td><td>1</td><td>2</td><td>2</td><td>27</td><td>1</td><td>2</td><td>35</td><td>0,23</td></tr><tr><td>Forêt</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>2</td><td>25</td><td>1</td><td>29</td><td>0,14</td></tr><tr><td>Terres arables</td><td>0</td><td>2</td><td>1</td><td>1</td><td>2</td><td>1</td><td>0</td><td>41</td><td>48</td><td>0,18</td></tr><tr><td>Total</td><td>21</td><td>23</td><td>16</td><td>33</td><td>35</td><td>34</td><td>29</td><td>48</td><td>239</td><td></td></tr><tr><td>Errore Commission</td><td>0,05</td><td>0,13</td><td>0,13</td><td>0,18</td><td>0,26</td><td>0,21</td><td>0,14</td><td>0,20</td><td></td><td>84</td></tr></tbody></table>

Le pourcentage de la précision globale et le coefficient Kappa pour les quatre classifications dépassent respectivement 83% et 0.80. Les classifications sont donc jugées fiables et exploitables pour nos domaines d'intérêt, qui sont la planification du développement et la gestion du territoire. Le tableau 4 présente les précisions globales et

les coefficients Kappa obtenus par les différentes matrices de confusion pour les quatre années d'étude.

En analysant les résultats des matrices de confusion classe par classe, nous remarquons que les plus fortes imprécisions ont lieu entre la classe des vergers, la classe des cultures irriguées et la classe des oliviers. Cela
s'explique par le fait que l'attribution de certaines parcelles à l'une des classes n'est pas évidente lors du choix des parcelles d'entrainement, tant par photo-interprétation que sur le plan spectral.

<table><caption>**Tableau 4 :** Précisions globales et coefficients Kappa des matrices de confusion pour les quatre années d'étude</caption><thead><tr><td>Années</td><td>Précision globale %</td><td>Coefficient de Kappa</td></tr></thead><tbody><tr><td>1987</td><td>83,15</td><td>0,80</td></tr><tr><td>2003</td><td>85,17</td><td>0,82</td></tr><tr><td>2014</td><td>84,20</td><td>0,81</td></tr><tr><td>2018</td><td>85,32</td><td>0,81</td></tr></tbody></table>

Les principaux changements dans l'occupation du sol concernent la transformation des superficies occupées par les terres arables en surfaces irriguées, en surfaces d'oliviers et d'arbres fruitiers ainsi qu'en zones urbaines et périurbaines. Les différentes conversions sont détaillées dans le tableau 5. Ce tableau présente la superficie de chacune des classes pour les quatre années d'étude et la variation de l'occupation du sol pendant cette période de 31 ans.

<table><caption>**Tableau 5 :** Évolution des occupations du sol entre 1987 et 2018 dans la plaine de Saïss.</caption><thead><tr><th rowspan="2">Occupation du sol</th><th colspan="4">Surfaces (km²)</th><th rowspan="2">Variation (%) (1987-2018)</th></tr><tr><th>1987</th><th>2003</th><th>2014</th><th>2018</th></tr></thead><tbody><tr><td>Urbanisme</td><td>62</td><td>105</td><td>123</td><td>133</td><td>115</td></tr><tr><td>Vergers</td><td>287</td><td>420</td><td>514</td><td>577</td><td>101</td></tr><tr><td>Oliviers</td><td>77</td><td>145</td><td>168</td><td>182</td><td>136</td></tr><tr><td>Surface d'eau</td><td>4</td><td>5.5</td><td>5.3</td><td>5</td><td>20</td></tr><tr><td>Terres irrigués</td><td>35</td><td>64</td><td>72</td><td>78</td><td>123</td></tr><tr><td>Forêts</td><td>21</td><td>20</td><td>19</td><td>19</td><td>-10</td></tr><tr><td>Terres arables</td><td>1430</td><td>1165</td><td>1024</td><td>958</td><td>-29</td></tr><tr><td>Parcours</td><td>118</td><td>115</td><td>112</td><td>111</td><td>-6</td></tr></tbody></table>

L'analyse des changements de l'occupation du sol a permis de conclure que la classe d'urbanisme a progressé de 62 km² en 1987, à 133 km² en 2018. Selon le recensement de la population publié par l'HCP en 2014, la population a évolué de 1 986 000 habitants en 2004 à 2 408 580 habitants en 2014, ce qui explique l'augmentation des zones urbaines observée entre 1987 et 2018.

Depuis 2008, le pays a adopté une stratégie de développement dont l'objectif est de faire de l'agriculture l'un des moteurs du développement économique national. Cette stratégie, appelée Plan Maroc Vert (FAO, 2016), a permis le développement de la productivité des petites et moyennes exploitations familiales. Les résultats obtenus se manifestent par le développement de la surface des cultures irriguées et des arbres fruitiers (y compris l'olivier). Les ressources en eau utilisées pour l'irrigation sont essentiellement d'origine souterraine (pompage au niveau de la nappe phréatique). L'une des conséquences visibles de l'extension de l'arboriculture est que les forêts et les terres arables ont diminué respectivement de 10% et de 29% entre 1987 et 2018. La superficie de l'eau libre augmente aussi et ce changement dépend non seulement des précipitations, mais également de la réalisation de deux barrages de protection autour de la ville de Fès, pour prévenir les inondations (Barrage El Gaada et Barrage Mehraz). Les activités d'irrigation effectuées par la construction de plusieurs bassins de stockage des eaux de pompage ont aussi contribué à l'augmentation des surfaces d'eau libre. Les parcours occupent une faible superficie, mais ils montrent une baisse significative (d'environ 6 %), due à la mise en culture et/ou à l'urbanisation.
<figure><img src="image_7.png" /><figcaption>Figure 7 : Catre d'occupation du sol du 25 mai 1987</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8 : Catre d'occupation du sol du 27 mai 2003</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 9** : Catre d'occupation du sol du 2 juin 2014</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 10** : Carte d'occupation du sol du 29 juin 2018</figcaption></figure>
Pour mettre en évidence les changements spatio-
temporels de l'occupation du sol observés entre 1987 et
2018 sur le nord-est de la zone d'étude, nous proposons

(Figure 11) un agrandissement de cette zone (repérable
par le cadre rouge situé au NE de la figure 7).

<figure><img src="image_11.png" /><figcaption>**Figure 11 :** Agrandissement de la zone nord-est pour mettre en évidence la variation de l'occupation du sol</figcaption></figure>

4. Conclusion

Les résultats obtenus par la classification des images
Landsat en utilisant l'algorithme du maximum de
vraisemblance nous apparaissent très satisfaisants
(valeurs élevées obtenues pour les précisions de
classification et des coefficients Kappa). Ces résultats
confirment que la plaine de Saïss a connu des
changements notables en termes d'occupation du sol sur
cette période de 31 ans (entre 1987 et 2018). Ces
changements concernent un étalement urbain important,
une extension des surfaces implantées en oliviers et en
arbres fruitiers (vergers) au détriment des terres arables et
des parcours, qui sont en régression. La zone de forêt a

subi un faible changement (-10%) en partant d'une surface
initiale (21 km²) réduite par rapport à la surface totale de la
plaine (2100 km²).

Ces cartes d'occupation du sol constituent une base pour
mener des travaux qui visent à la gestion de l'eau et des
terres à l'échelle de la plaine de Saïss. Les résultats
obtenus permettront par la suite d'évaluer les composants
du bilan hydrologique dans le système aquifère de Saïss
(potentiel de recharge). Ils seront aussi utilisés dans
l'évaluation des besoins en eau potable et en eau pour
l'agriculture. Notre programme de recherche est d'étudier
d'éventuelles corrélations entre l'évolution de l'occupation
du sol et l'alimentation des réserves hydriques souterraines.

Références

ABHS (Agence du Bassin Hydraulique du Sebou), 2016. Projet de Sauvegarde de la plaine du Saïss : Etude et Plan d'Action Environnemental et Social. Rapport, OPID 48414, 80p.

Amraoui F., 2005. L'aquifère karstique du Moyen Atlas tabulaire : problème de la turbidité affectant les sources du complexe Ribâa-Bittit (Maroc). Sciences et techniques de l'environnement, Mémoire hors-série, 5-8.

Amraoui F., Razack M., Bouchaou L., 2003. Turbidity dynamics in karstic systems. Example of Ribaa and Bittit springs in the Middle Atlas (Morocco). Hydrological Sciences Journal, **48** (6), 971-984

Bargach, K., Ruano P., Chabli A., Galindo-Zaldívar J., Chalouan A., Jabaloy A., Akil M., Ahmamou M., Sanzde Galdeano C, Benmakhlouf M., 2004. *Recent Tectonic Deformations and Stresses in the Frontal Part of the Rif Cordillera and the Saïss Basin (Fes and Rabat Regions, Morocco).* Pure and Applied Geophysics, **161**, 521-540.

Bateelan O., De Smeedt F., 2007. GIS-based recharge estimation by coupling surface–subsurface water balances. Journal of hydrology, **337**, 337-355.

Blaschke T., 2010. Object based image analysis for remote sensing. Journal of Photogrammetry and Remote Sensing, **65**(1), 2-16.

Bolstad P., Lillesand T.-M., 1991. Rapid maximum likelihood classification. Photogrammetric Engineering and Remote Sensing, **57**(1), 67-74.

Bru D., 2015. Corrections atmosphériques pour capteurs à très haute résolution spatiale en zone littorale. Thèse de doctorat, Université de Bordeaux, 203 p.

Campbell J.-B., Wynne R.-H., 2011. *Introduction to Remote Sensing*. Guilford Press, New York, 667 pp.

Chabli A., Akil M., A. Chalouan A., Ruano P., Bargach k., 2005. Déformations Néotectoniques dans Les Dépôts Plio- Quaternaires De La Région De Casablanca – Mohammedia (Meseta Côtière, Maroc). Revista Sociedad Geológica de España, **18**, 169-178.

Chandran A. J., Christy J., 2015. A survey of cloud detection techniques for satellite images. *International Research Journal of Engineering and Technology*, **2**(9), 2485-2490.

Chamayou J., Combe M., Genetier B., Leclerc C., 1975. Plaines et bassins du Maroc atlantique Bassin de Meknès-Fès et couloir de Fès-Taza. In *Ressources en Eau du Maroc, Tome 2, Editions du Service Géologique du Maroc*, Rabat, 41-71.

Cheng G., Han J., Lu X., 2017. Remote sensing image scene classification: Benchmark and state of the art. Proceedings of the IEEE, **105**(10), 1865-1883.

Cheng J., Liang S., Liu H., Nie A., Liu Q., 2018. A disaggregation approach for estimating high spatial resolution broadband emissivity for bare soils from Landsat surface reflectance. International Journal of Digital Earth, **11**(7), 691-702.

El Garouani A., Mulla D. J., El Garouani S., Knight J., 2017. *Analysis of urban growth and sprawl from remote sensing data: Case of Fez, Morocco*. International Journal of Sustainable Built Environment, **6**(1), 160-169.

El Garouani A., Tribak A., Abahrour M., 2010. Assessment the effects of land use-cover changes on regional soil loss susceptibility using RUSLE model and remote sensing data. Red Books, International Association of Hydrological Sciences Publication, **340**, 343-349

Essahlaoui A., Sahbi H., Bahi L., EL-Yamine N., 2001. Preliminary Survey of the Structure and Hydrogeology of the Western Sais Basin, Morocco, Using Electrical Resistivity. Journal of African Earth Sciences, **32**, 777-789.

Essahlaoui A., El Ouali A., 2003. Determination of the Geological Structure of the Southern Part of the Saïss Plain (Meknes-Fes Basin, Morocco) Using Geo-Electric Techniques. Bulletin of Engineering Geology and the Environment, **62**, 155-166.

ESRI (Environmental Systems Research Institute) 2020. Environnements d'analyse et Spatial Analyste : Filtre Majoritaire. https://www.esrifrance.fr

FAO, 2016. Plan Maroc Vert : les grands principes et avancées de la stratégie agricole marocaine. Bulletin de Synthèse Souveraineté Alimentaire, **20**, 1-8.

Gao B.-C., 1996. NDWI - A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment, **58**, 257-266.

Girard M.C. et Girard C.M., 2010. Traitement des données de télédétection : Environnement et ressources naturelles. DUNOD, 2ème édition, 553p.

HCP (Haut-Commissariat au Plan), 2014. Recensement général de la population et de l'habitat du Maroc en 2014. https://www.hcp.ma/

Ienco D., Gaetano R., Dupaquier C., Maurel P., 2017. Land Cover Classification via Multitemporal Spatial Data by Recurrent Neural Networks. IEEE Geoscience and Remote Sensing Letters, **14** (10), 1685-1689. DOI: 10.1109/LGRS.2017.2728698

Jackson R. D., Pinter P. J., Reginato R. J., 1986. Detection and evaluation of plant stress for crop management decisions. IEEE Transactions on Geoscience and Remote Sensing, **1**, 99-106.

Khatami R., Mountrakis G., Stehman S.V., 2017. Mapping per-pixel predicted accuracy of classified remote sensing images. Remote Sensing of Environment, **191**, 156-167.
Landis J. R., Koch G. G., 1977. *The measurement of observer agreement for categorical data*. Biometrics, 1,159–174

Lassale P., 2015. Etude du passage à l'échelle des algorithmes de segmentation et de classification en télédétection pour le traitement de volumes massifs de données. Thèse PhD, Université de Toulouse

Le Hégarat M. S., 2005. Classification d'images de télédétection pour l'estimation et le suivi de paramètres géophysiques. Rapport d'Habilitation à Diriger des Recherches, l'Université de Versailles Saint Quentin, CETP IPSL, 111p.

Lillesand T., Kiefer R-W., Chipman J., 2014. *Remote Sensing and Image Interpretation*. Wiley Publisher, 804 pp.

Lu D., Weng Q., 2007. *A survey of image classification methods and techniques for improving classification performance*. International Journal of Remote Sensing, 28 (5), 823-870.

Martin M. E., Newman S. D., Aber J. D., Congalton R. G., 1998. *Determining forest species composition using high spectral resolution remote sensing data*. Remote Sensing of Environment, 65(3), 249-254.

Matthieu M. G., 2017. Techniques d'analyse de contenu appliquées à l'imagerie spatiale. Thèse PhD, Université de Toulouse, 241p.

McFeeters S. K., 2013. *Using the Normalized Difference Water Index (NDWI) within a Geographic Information System to Detect Swimming Pools for Mosquito Abatement: A Practical Approach*. Remote Sensing, 5(7), 3544-3561. DOI:10.3390/rs5073544

McFeeters S. K., 1996. *The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features*. International Journal of Remote Sensing, 17(7), 1425-1432.

Mubako S., Belhaj O., Heyman J., Hargrove W., Reyes C., 2018. *Monitoring of Land Use/Land-Cover Changes in the Arid Transboundary Middle Rio Grande Basin Using Remote Sensing*. Remote Sensing, 10(12), 1-17

Pelletier Ch., 2017. Cartographie de l'occupation des sols à partir de séries temporelles d'images satellitaires à hautes résolutions : Identification et traitement des données mal étiquetées. Thèse de Doctorat de l'Université de Toulouse, 289p.

Richards J.-A., 1999. *Remote Sensing Digital Image Analysis: An Introduction*. Springer-Verlag, Berlin, Germany, 365p.

Rodes I., 2016. Estimation de l'occupation des sols à grande échelle pour l'exploitation d'images d'observation de la Terre à hautes résolutions spatiale, spectrale et temporelle. Thèse PhD, Toulouse

Rouse J.W., Haas R.H., Schell J.A., Deering D.W., 1973. *Monitoring vegetation systems in the great plains with*

ERTS. Third Earth Resources Technology Satellite-1 Symposium Washington D.C., 309-317.

Rosenfiel G.H., Filzpatrick-lins K., 1986. *A Coefficient of Agreement as a Measure of Thematic Classification Accuracy*. Photogrammetric Engineering and Remote Sensing, 52, (2), 223-227.

Shivakumar B. Ra., Rajashekaradhyab S. V., 2018. *Investigation on Land Cover Mapping Capability of Maximum Likelihood Classifier: A Case Study on North Canara, India*. Procedia Computer Science, 143, 579-586

Tamim A., 2015. Segmentation et classification des images satellitaires : application à la détection des zones d'upwelling côtier marocain et mise en place d'un logiciel de suivi spatiotemporel. Thèse de Doctorat, Université Mohammed V, Faculté des Sciences, Rabat, 93 p.

Taylor R.G., Scanlon B., Döll P., Rodell M., Van Beek R., Wada Y., Konikow L., 2013. *Ground water and climate change*. Nature Climate Change, 3(4), 322-329.

Tucker C.J., 1979. *Red and photographic infrared linear combinations for monitoring vegetation*. Remote Sensing of the Environment, 8, 127-150.

USGS 2019. Landsat collection 1 level 1 product definition. EROS Sioux Falls, South Dakota, USA, 32 p.

Vázquez-Jiménez R., Ramos-Bernal R.N., Romero-Calcerrada R., Arrogante-Funes P., Tizapa, S.S., Novillo C.J., 2018. *Thresholding Algorithm Optimization for Change Detection to Satellite Imagery*. In Colorimetry and Image Processing, INTECH, 163-182.

Weigand C.L., Richardson A.J., 1990. *Use of spectral vegetation indices to infer leaf area, evapotranspiration and yield*. Journal of Agronomy, 82, 623-629.

Xu H., 2006. *Modification of Normalised Difference Water Index (NDWI) to enhance open water features in remotely sensed imagery*. International Journal of Remote Sensing, 27(14), 3025-3033.

Zha Y., Gao J., Ni S., 2003. *Use of normalized difference built-up index in automatically mapping urban areas from satellite imagery*. International Journal of Remote Sensing, 24, 583-594.
**Titre en anglais:** CONTRIBUTION OF LANDSAT IMAGES TO THE STUDY OF LAND-USE DEVELOPMENT OVER THE PERIOD 1987-2018, IN THE SAÏSS PLAIN IN MOROCCO

**Légendes en anglais**

Figure 1: Geographical location of the study area (The base map corresponds to a composite image (RGB = Bands 6, 5, 4) of Landsat8 of June 29, 2018. The green color corresponds to the vegetation)

Figure 2: Scenes mosaic covering the study area

Figure 3: Test of classes separability (learning sites) by point cloud for the blue and near infrared spectral bands

Figure 4: Zoom in on land cover classes identified from the Landsat8 image of 2018 (composite image: RGB = Bands 6, 5, 4) and their thematic correspondences in the field (photos)

Figure 5: Supervised classification process

Figure 6: Post-classification processes (before filtering in the left, after filtering in the right)

Figure 7: Land cover map of May 25, 1987

Figure 8: Land cover map of May 27, 2003

Figure 9: Land cover map of June 02, 2014

Figure 10: Land cover map of June 29, 2018

Figure 11: Zoom in of the NE area to highlight the variation in land cover

**Tables:**

Table 1: Used Landsat images

Table 2. Standard for interpreting the Kappa value. Source: Landis and Koch (1977)

Table 3: Confusion matrix for a classified image (06-29-2018)

Table 4: Global accuracy and Kappa coefficients of confusion matrices for four years of study

Table 5: Land use evolution between 1987 and 2018 in the Saïss plain.