# Argan Tree (Argania spinosa (L.) Skeels) Mapping Based on Multisensor Fusion of Satellite Imagery in Essaouira Province, Morocco

Authors: Aicha Moumni; Tarik Belghazi; Brahim Maksoudi; Abderrahman Lahrouni
Journal: Journal of Sensors
 (2021)

---

## Abstract

Tree species identification and their geospatial distribution mapping are crucial for forest monitoring and management. The satellite-based remote sensing time series of Sentinel missions (Sentinel-1 and Sentinel-2) are a perfect tool to map the type, location, and extent of forest cover over large areas at local or global scale. This study is focused on the geospatial mapping of the endemic argan tree (Argania spinosa (L.) Skeels) and the identification of two other tree species (sandarac gum and olive trees) using optical and synthetic aperture radar (SAR) time series. The objective of the present work is to detect the actual state of forest species trees, more specifically the argan tree, in order to be able to study and analyze forest changes (degradation) and make new strategies to protect this endemic tree. The study was conducted over an area located in Essaouira province, Morocco. The support vector machine (SVM) algorithm was used for the classification of the two types of data. We first classified the optical data for tree species identification and mapping. Second, the SAR time series were used to identify the argan tree and distinguish it from other species. Finally, the two types of satellite images were combined to improve and compare the results of classification with those obtained from single-source data. The overall accuracy (OA) of optical classification reached 86.9% with a kappa coefficient of 0.84 and declined strongly to 37.22% (kappa of 0.29) for SAR classification. The fusion of multisensor data (optical and SAR images) reached an OA of 86.51%. A postclassification was performed to improve the results. The classified images were smoothed, and therefore, the quantitative and qualitative results showed an improvement, in particular for optical classification with a highest OA of 89.78% (kappa coefficient of 0.88). The study confirmed the potential of the multitemporal optical data for accurate forest cover mapping and endemic species identification.


---

## Body

## 1. Introduction

The monitoring of forest cover plays a crucial role in biodiversity, feedstock, and water cycle, etc. Therefore, tree species discrimination is necessary and fundamental for this process. Satellite imagery together with machine learning techniques have become an irreplaceable tool for tree species mapping. Passive and/or active sensors provide valuable geospatial information to identify tree types. Recently, machine learning algorithms are techniques that have been successfully used based on remote sensing data for tree type classification [1]. Current satellite sensors, such as Sentinel missions, facilitated the tree species mapping at local as well as national scale.In literature, the most remote sensing data used for forest type mapping is optical imagery. For such data, the forest type classification is based on reflected spectral features that are acquired by optical sensors. However, the different forest type classes (tree crops) can be characterized by similar spectral signatures. Thus, tree species identification becomes difficult and complicated using only spectral features [2]. Weather is one of the factors that caused such a confusion due to the dependence of optical data on sunlight. SAR imagery is a supplement of optical data due to its capacity to acquire images in all-day and in all-weather conditions through penetrating rain and clouds. In 2004, Touzi et al. [3] used the C band of SAR data to discriminate forest tree species and figured out that the information of VV, HH, and VH polarizations identified better the forest trees without leaves. Another study [4] showed that information extracted from SAR imagery can discriminate easily between forest and nonforest types.In addition, the texture features [5, 6] extracted from gray-level cooccurrence matrix (GLCM) have been widely investigated in forest type classification using SAR imagery [6–9]. However, the results of land cover/land use (LC/LU) classification using SAR data are not consistent [9–11]. Therefore, many studies indicated that combining SAR features with optical data can improve the forest type classification results [12, 13]. Machine learning algorithms are powerful methods for forest and crop classification. SVM classification algorithm has proven its effectiveness for classifying forest tree species using satellite imagery [14–16].The endemic Moroccan tree called argan (Argania spinosa) covered 950,000 ha in 2010 [17]. It is part of the semiarid Mediterranean domain in southwest Morocco in transition towards the Sahara [18]. The zone was declared a UNESCO MAB (Man and the Biosphere Reserve) in 1998. Argan trees have provided multiple ecosystem products and services including the provision of fruits from which argan oil is produced. Argan is a slow-growing tree with shrubby architecture which has a lifespan of around 200 years [19]. Argan trees have traditionally provided multiple ecosystem products and services including the provision of fruits from which argan oil is produced; leaves and young shoots eaten by sheep, goats, and camels; and wood for fiber and fuel [20].During the last century, more than half of the argan forest of Morocco disappeared, mainly in the plains, and its mean density fell from 100 to less than 30 stumps/ha to meet the growing fuel needs of major Moroccan cities [21]. Recent work distinguishes between lowland and mountain argan plantations, where ecological and economic conditions differ significantly [22]. The former seems very threatened, as it has great agricultural potential and absorbs most of the emigration from mountain areas. Socioeconomic changes and the farms’ modernisation since the 1980s have led to a cumulative 2.6% reduction in forest cover over the last 17 years [23]. In mountain argan groves, where the user population retains a traditional pattern of use, several trends can be discerned. Therefore, these several factors caused degradation, desertification, and problems in the natural regeneration of the argan population [24, 25]. To overcome these problems, specific management strategies are needed, in such a way that the state of forest resources should be assessed with detailed and accurate geospatial distribution of different forest tree species and, more specifically, argan trees.The estimate of forest cover and its mapping relies mainly on fieldwork and ground surveys, which require cost and time. However, high-resolution satellite data such as Landsat and Sentinel images are available and provide huge amounts of data regularly on large areas in short time and remotely. Satellite image-based land cover land use mapping, generally, still have some challenges despite the improvement of spatial and temporal resolution [26–29]. Particularly, forest tree species discrimination can face the presence of similar classes in the study area and/or the presence of several features within the same pixel, contributing to spectral confusion between different forest over types [30, 31]. In this contest, we attempted to identify the argan tree from other forest species based on the classification of time series of SAR data and its derived products (GLCM texture features: correlation variance and mean), the classification NDVI time series derived from optical images, and finally, the classification of combined multisensor time series data using an SVM classifier.The achievement of the objective of the present work is illustrated in four steps, namely, (i) to model and discriminate the phenological evolution of the Moroccan argan tree from other tree species using S2-derived NDVI time series; (ii) to map the spatial extent of forest cover and identify its tree species; (iii) to evaluate the tree species mapping using optical data, SAR data, and the fusion of the two types of satellite imagery; and (iv) to evaluate and assess the potential of SAR data to improve the tree species classification performance and overcome the confusion constraints.


## 2. Background

Morocco is one of the few countries in North Africa to have a diverse range of endemism and biodiversity [32]. The only species in its genus, Argania spinosa (L.) Skeels, belongs to the Sapotaceae, a tropical and subtropical tree family. It is endemic only to Morocco, where it grows in arid and semiarid climates with 150 to 400 mm of rainfall per year. It is a slow-growing spiny tree that can reach a maximum height of around 10 m (Figure 1). The leaves are small (20 mm), with a spatulated shape and an entire margin. Despite having evergreen foliage, leaves may be shed entirely or partially in response to summer stress.Figure 1
Argan tree in Essaouira region (07/09/2019).In 1998, M’Hirit et al. classified it as a facultative summer semideciduous tree. Flowering takes place mostly in the spring, with a lot of variation between individuals, and there is a secondary peak in the autumn. As a result, fruit production continues over a long period of time, with fruits appearing in various stages of ripening from April to September [20]. The fruits are stalkless drupes that look like olives and have a hard endocarp that protects an oil-rich kernel. In different cameras, one endocarp can contain up to four kernels, although the normal number is one.It is a multipurpose forest tree (forest, fruit, and fodder) of great importance to the country, both biologically, phylogenetically, and ecologically, as well as economically and socially [33]. It has significant medicinal and therapeutic properties due to the oil extracted from its fruits [34, 35]. Furthermore, it is highly sought after in cosmetics as a skin and hair revitalizing agent. This ecosystem is based on a balance of resources and human exploitation, and it also plays an important role in the fight against desertification and erosion.Regardless of the various roles played by this species, the threat of argan degradation is a major concern for both the general public and scientists [36]. Artificial regeneration is the main way to ensure the survival of the argan tree in particular and forest cover in general [37, 38]. Figure 2 illustrates some reforested areas around Smimou in Essaouira province.Figure 2
Reforested areas in Essaouira (07/09/2019).With the absence of accurate and detailed mapping of the argan forest in the region, as well as an assessment of recent changes, we use the potential of satellite images acquired from different sensors to detect the actual state of forest species trees and, more specifically, the argan tree. Forest cover mapping is a practical tool for decision making in forest management. Moreover, it is a means for the definition of the area’s objectives and management operations, to establish a framework of intervention for the manager.


## 3. Materials and Methods

### 3.1. Experimental Site

The test site is an area of approximately 2000 ha across Smimou village in Essaouira of Morocco (Figure3). This rural district is considerably large and varied, with the mountain range covering a large part of it. The most famous mountain and the highest peak in Essaouira province (912 m), is Jbel Amsittene [39]. The Smimou zone is characterized by an arid to semiarid climate. During the year, the temperature generally varies from 10°C to 22°C and is rarely below 8°C or above 26°C [40, 41]. The mean temperature is around 20°C. The rains are irregular with an annual average of 300 mm.Figure 3
Study area location.Indeed, the climate is also diverse in terms of both temperature and precipitation. This is due to the oceanic (Atlantic) setting on one side and the height of the mountains on the other. Precipitation varies by 48.6 mm between the driest and the wettest month.


### 3.2. Reference Data

The field data are used to extract profiles and train the classifiers. Moreover, those samples are exploited to validate the classification results. In collaboration with the team of the Forestry Research Center (CRF), samples were collected during two field campaigns in 2019, and completed by samples extracted using the archive of very high spatial resolution images from Google Earth. The surveyed plots are accompanied by descriptions informing about the land cover type, and they are randomly split into two groups: calibration samples and validation samples. The land cover typology was grouped into eight main classes that were selected in terms of abundance over the area. This typology contains argan, sandarac gum, fallow, bare soil, olive, argan + sandarac gum, noncrop, and agricultural lands. Table1 illustrates the number of training and validation samples for each land cover type.Table 1
Detailed description of reference data used for tree species classification.
CalibrationValidationTotalSample numbersPixel numbersSample numbersPixel numbersArgan tree56255103134159Olive tree7920138136117Sandarac gum tree233434424367Argan + sandarac gum176625526672Agricultural land132722413837Fallow92392210931Bare soil286572610154Non-crop91892813337Total234340574Reference date was used also for the extraction of the temporal NDVI profiles of land cover types, in order to determine the different phenological states of the argan tree and other tree species


### 3.3. Satellite Data and Data Processing

#### 3.3.1. Optical Imagery

The most convenient satellite data is determined by a number of factors, including (i) spatial resolution, (ii) revisit time, (iii) image availability and cost over the study area, and (v) study area characteristics [42]. S2 is primarily aimed at freely offering high spatiotemporal resolution optical imagery for operational land cover and change detection mapping [43]. These data can be downloaded on the Theia Land Service website [44], which produces and distributes level 2A products, corrected atmospherically using the MAJA/MACCS software developed thanks to the coordination between Cnes/Cesbio and DLR [45–48].For the present study, S2 images have been downloaded from December 2018 to December 2019 over the Smimou experimental site. The 35 downloaded images, covering the whole year of 2019 and well distributed over the agricultural seasons, have been used to derive the Normalized Difference Vegetation Index (NDVI) time series.This index is the most commonly used vegetation index in remote sensing. NDVI [49] is a normalized index that generates an image showing the degree of vegetation cover (biomass). The absorption of chlorophyll pigment in the red channel and the high reflectivity of plant matter in the near infrared (NIR) channel are used to create this index, which is based on the comparison in two channels of a multispectral raster data set [50, 51]. NDVI is calculated based on the normalized difference of the NIR and R spectral bands and illustrated as follows:
(1)NDVI=NIR−RNIR+R.This index generates values between -1 and 1, mainly representing the vegetation cover, where negative values are mainly generated by clouds, water, and snow and values close to zero are mainly generated by rock and bare soil. The very low NDVI values (0.1 and below) correspond to waste rock, sand, or snow surfaces. Intermediate values (0.2 to 0.3) represent areas of shrubs and grasslands, while high values (0.6 to 0.8) indicate temperate or tropical rainforests. The transformation of the images into NDVI layers was carried out for the whole series of 35 images, in order to extract the phenological evolution of the tree species and the other land cover classes selected to be identified and mapped in this study.


#### 3.3.2. SAR Imagery

Sentinel-1 (S1) is a two-satellite constellation that orbits the Earth in 6 days. The twin Sentinel-1A/B constellations use synthetic aperture radar (SAR) imaging to acquire imagery regardless of weather conditions.S1 images were acquired in the Ground Range Detected (GRD) format and interferometric wide (IW) swath mode with dual polarization (VV and VH) through the PEPS CNES websitehttps://peps.cnes.fr/. The S1 products offered by the PEPS platform are calibrated and orthorectified on the S2 grid in order to facilitate the joint use of the two missions [52–57]. Two SAR images were downloaded each month to coincide with key dates in terms of crop vegetative cycles. This produced a set of 22 SAR images (Figure 4).Figure 4
Distribution of radar and optical image acquisition dates.S1 images were preprocessed using the Sentinel Application Platform (SNAP) and ENVI software as follows: (i) the removal of the speckle noise using the enhanced Lee filter, (ii) the conversion of SAR bands' pixel values to dB values to facilitate the interpretation and the analysis of data, and (iii) the extraction of texture features using a gray-level cooccurrence matrix (GLCM).(1) Speckle Noise Filtering and Conversion to dB. The main issue in SAR imaging classification and interpretation is the speckle noise effect, which is caused by coherent interference of the wave reflected from several basic scatterers [58]. Speckle noise removal filtering is a technique for improving image quality by decreasing the speckle noise effect.The enhanced Lee filter has been chosen in this study to minimize the speckle effect on the images without loss of information. According to the literature, it is one of the most used filters to reduce speckle noise in SAR remote sensing data [59, 60].We have applied the enhanced Lee filter [61] implemented in ENVI software to the S1 images acquired in VV and VH polarizations, with a window size of 5×5 [62, 63]. The filtered images were then converted into decibels (dB) using the following:
(2)X=10×log10x,where x is the value of each pixel, X is the value converted to dB.(2) Texture Feature Extraction. SAR imagery contains texture information [5, 64]. Image classification will benefit greatly from SAR texture information features. New texture images can be reconstructed using the spatial information in the SAR image. Texture features provide details from neighboring pixels, which is useful in discriminating the various crop types found in agricultural fields [65].A gray-level cooccurrence matrix (GLCM) is considered one of the most popular techniques used for texture analysis. It is calculated based on a statistical method, also known as the gray-level spatial dependence matrix. GLCM is capable of assessing texture features by taking into account the spatial relationship between pixels and their surroundings [66, 67].From many proposed statistical texture methods, the gray-level cooccurrence matrix (GLCM) is one of the most powerful for land cover/land use (LC/LU) mapping [68, 69].The texture analysis was applied with VV and VH polarizations using the Sentinel toolbox called SNAP (SeNtinel Application Platform) based on GLCM statistics. GLCM is a probability measure of two gray levels separated by a given distance occurring in the same orientation. According to the literature [70–73], the most common and popular GLCM texture features used for tree type discrimination are as follows: correlation, mean, and variance. These three texture measurements were calculated and used for the present study.Thecorrelation function defines the similarity between the elements of columns or rows in GLCM and is a measure of the image’s gray-tone linear dependencies. The GLCM correlation feature can be calculated as follows:
(3)∑i,j=0N−1iPi,j−μxμyσxσy.The GLCMmean value is more than just the average of all original pixel values; each pixel value is weighted based on how frequently it occurs in combination with a specific neighbor pixel value. The expression of the mean is illustrated as follows:
(4)∑i,j=0N−1iPi,j.Thevariance texture parameter targets the SAR image’s partial characteristic. It can be computed as follows:
(5)∑i,j=0N−1iPi,ji−μ2, where
(i)
Pi,j is a normalized gray-tone spatial dependency matrix such as ∑i,j=0K−1iPi,j=1(ii)
i and j represent, respectively, the rows and columns for the mean, variance, and correlation measurements(iii)
μ is the mean for the variance texture measurement(iv)
N is the number of distinct gray levels in the quantized image(v)
μx,μy and σx,σy are the means and standard deviations of Px and Py respectively, for the correlation texture measure [74]


### 3.4. Methodology

The methodological approach followed in this work (Figure5) consists of 5 main steps:
(1)
Data acquisition(2)
Data preprocessing(3)
Extraction and analysis of profiles and classification(4)
Accuracy assessment and postclassification(5)
Comparison, analysis, and map productionFigure 5
Methodological approach followed for the present study.


#### 3.4.1. SVM Classification

SVM is one of the most appealing algorithms for remotely sensed data classification because of the benefits of generalization even with limited training samples, which is common in remote sensing data processing [26, 75]. Furthermore, because SVM is a supervised nonparametric statistical learning method, it does not require a training set that strictly follows the standard independent and identical distribution. The benefits of SVM stem from two factors: transforming the original space training set into a very high dimensional new space and discovering a large linear margin boundary in the new space [76]. SVM is a classifier based on the structural risk minimization theory that is aimed at reducing generalization error by maximizing margins on training data [77]. SVM classifiers can use a variety of kernel functions. However, only four of them have been widely used to classify satellite data: polynomial, linear, radial basis function (RBF), and sigmoid kernels [77]. In the present study, The SVM-RBF kernel was selected and used for all image classifications as one of the most effective machine learning algorithms for classifying forest tree species using satellite imagery [14–16].


#### 3.4.2. Classification Scenarios

ENVI Classification Scenarios is a professional software for classifying satellite imagery using both unsupervised and supervised algorithms. For this work, we have limited ourselves to using supervised methods, in particular the support vector machines (SVM). The techniques generally provide better results in the production of land cover maps [78–80]. As a first step, we started by classifying the time series of NDVI, VV, VH, and VH/VV. We then combined the data from the two sensors to explore different possible scenarios (combinations). Table 2 describes these scenarios in the order followed during the image classification step.Table 2
Scenarios created for the classification step.
ScenarioBandsSatellite sensorType1NDVIS2Optical2VVS1SAR3VHS1SAR4VH/VVS1SAR5VV, VHS1SAR6VV, VH/VVS1SAR7VH, VH/VVS1SAR8VV, VH, VH/VVS1SAR9TextureS1SAR10Texture, VH, VH/VVS1SAR11Texture, VH, VVS1SAR12Texture, VV, VH/VVS1SAR13NDVI, VVS1 and S2Optical and SAR14NDVI, VHS1 and S2Optical and SAR15NDVI, VH/VVS1 and S2Optical and SAR16NDVI, VV, VHS1 and S2Optical and SAR17NDVI, VV, VH/VVS1 and S2Optical and SAR18NDVI, VH, VH/VVS1 and S2Optical and SAR19NDVI, VV, VH, VH/VVS1 and S2Optical and SAR20NDVI, textureS1 and S2Optical and SAR21NDVI, texture, VH, VH/VVS1 and S2Optical and SAR22NDVI, texture, VVS1 and S2Optical and SAR


#### 3.4.3. Performance Assessment

Evaluation of classification performance was done by comparing images and by calculating the accuracy of the land use classification using statistical tools such as the confusion matrix and the kappa index (K) introduced by Cohen [81], which is used to validate and estimate the precision in the vegetation classification scheme [80]. Its formula is as follows:
(6)Kappa=P0−PcPp−Pc, where
(i)
P0: is equal to the actual percentage obtained from classification of land covers; it is equal to the quotient of the sum of the diagonal digits of the matrix with the total number of observations(ii)
Pc is the estimate of the probability of obtaining a correct classification. For a correct classification, the value of Pc is generally lower than P0(iii)
Pp: is the percentage obtained when the classification is perfectConfusion matrices have been calculated in order to evaluate the results, these matrices reveal not only the general errors made at the level of each class during interpretation but also errors due to confusion between land use classes. The mathematical expression used to calculate the OA of the classification is given by the average of the percentages of correctly classified pixels:(7)OA=∑i=1rXiiN×100, where
(i)
Xii is the number of diagonal pixels (i.e., correctly classified)(ii)
N is the total number of pixels(iii)
R is the total number of classes


#### 3.4.4. Postclassification: Smoothing

After the classification, the obtained images showed a lack of sharpness in the definition of the classified plots; the images always seem “noisy.” We then apply a smoothing filter, called majority/minority analysis [82], to the images using the ENVI software. This analysis filters the image by replacing the value of the center pixel in a window with a size of n×n (where n has to be defined) with the majority value in that window. In order to avoid “oversmoothing,” we have chosen a window with a size of 3×3.


## 3.1. Experimental Site

The test site is an area of approximately 2000 ha across Smimou village in Essaouira of Morocco (Figure3). This rural district is considerably large and varied, with the mountain range covering a large part of it. The most famous mountain and the highest peak in Essaouira province (912 m), is Jbel Amsittene [39]. The Smimou zone is characterized by an arid to semiarid climate. During the year, the temperature generally varies from 10°C to 22°C and is rarely below 8°C or above 26°C [40, 41]. The mean temperature is around 20°C. The rains are irregular with an annual average of 300 mm.Figure 3
Study area location.Indeed, the climate is also diverse in terms of both temperature and precipitation. This is due to the oceanic (Atlantic) setting on one side and the height of the mountains on the other. Precipitation varies by 48.6 mm between the driest and the wettest month.


## 3.2. Reference Data

The field data are used to extract profiles and train the classifiers. Moreover, those samples are exploited to validate the classification results. In collaboration with the team of the Forestry Research Center (CRF), samples were collected during two field campaigns in 2019, and completed by samples extracted using the archive of very high spatial resolution images from Google Earth. The surveyed plots are accompanied by descriptions informing about the land cover type, and they are randomly split into two groups: calibration samples and validation samples. The land cover typology was grouped into eight main classes that were selected in terms of abundance over the area. This typology contains argan, sandarac gum, fallow, bare soil, olive, argan + sandarac gum, noncrop, and agricultural lands. Table1 illustrates the number of training and validation samples for each land cover type.Table 1
Detailed description of reference data used for tree species classification.
CalibrationValidationTotalSample numbersPixel numbersSample numbersPixel numbersArgan tree56255103134159Olive tree7920138136117Sandarac gum tree233434424367Argan + sandarac gum176625526672Agricultural land132722413837Fallow92392210931Bare soil286572610154Non-crop91892813337Total234340574Reference date was used also for the extraction of the temporal NDVI profiles of land cover types, in order to determine the different phenological states of the argan tree and other tree species


## 3.3. Satellite Data and Data Processing

### 3.3.1. Optical Imagery

The most convenient satellite data is determined by a number of factors, including (i) spatial resolution, (ii) revisit time, (iii) image availability and cost over the study area, and (v) study area characteristics [42]. S2 is primarily aimed at freely offering high spatiotemporal resolution optical imagery for operational land cover and change detection mapping [43]. These data can be downloaded on the Theia Land Service website [44], which produces and distributes level 2A products, corrected atmospherically using the MAJA/MACCS software developed thanks to the coordination between Cnes/Cesbio and DLR [45–48].For the present study, S2 images have been downloaded from December 2018 to December 2019 over the Smimou experimental site. The 35 downloaded images, covering the whole year of 2019 and well distributed over the agricultural seasons, have been used to derive the Normalized Difference Vegetation Index (NDVI) time series.This index is the most commonly used vegetation index in remote sensing. NDVI [49] is a normalized index that generates an image showing the degree of vegetation cover (biomass). The absorption of chlorophyll pigment in the red channel and the high reflectivity of plant matter in the near infrared (NIR) channel are used to create this index, which is based on the comparison in two channels of a multispectral raster data set [50, 51]. NDVI is calculated based on the normalized difference of the NIR and R spectral bands and illustrated as follows:
(1)NDVI=NIR−RNIR+R.This index generates values between -1 and 1, mainly representing the vegetation cover, where negative values are mainly generated by clouds, water, and snow and values close to zero are mainly generated by rock and bare soil. The very low NDVI values (0.1 and below) correspond to waste rock, sand, or snow surfaces. Intermediate values (0.2 to 0.3) represent areas of shrubs and grasslands, while high values (0.6 to 0.8) indicate temperate or tropical rainforests. The transformation of the images into NDVI layers was carried out for the whole series of 35 images, in order to extract the phenological evolution of the tree species and the other land cover classes selected to be identified and mapped in this study.


### 3.3.2. SAR Imagery

Sentinel-1 (S1) is a two-satellite constellation that orbits the Earth in 6 days. The twin Sentinel-1A/B constellations use synthetic aperture radar (SAR) imaging to acquire imagery regardless of weather conditions.S1 images were acquired in the Ground Range Detected (GRD) format and interferometric wide (IW) swath mode with dual polarization (VV and VH) through the PEPS CNES websitehttps://peps.cnes.fr/. The S1 products offered by the PEPS platform are calibrated and orthorectified on the S2 grid in order to facilitate the joint use of the two missions [52–57]. Two SAR images were downloaded each month to coincide with key dates in terms of crop vegetative cycles. This produced a set of 22 SAR images (Figure 4).Figure 4
Distribution of radar and optical image acquisition dates.S1 images were preprocessed using the Sentinel Application Platform (SNAP) and ENVI software as follows: (i) the removal of the speckle noise using the enhanced Lee filter, (ii) the conversion of SAR bands' pixel values to dB values to facilitate the interpretation and the analysis of data, and (iii) the extraction of texture features using a gray-level cooccurrence matrix (GLCM).(1) Speckle Noise Filtering and Conversion to dB. The main issue in SAR imaging classification and interpretation is the speckle noise effect, which is caused by coherent interference of the wave reflected from several basic scatterers [58]. Speckle noise removal filtering is a technique for improving image quality by decreasing the speckle noise effect.The enhanced Lee filter has been chosen in this study to minimize the speckle effect on the images without loss of information. According to the literature, it is one of the most used filters to reduce speckle noise in SAR remote sensing data [59, 60].We have applied the enhanced Lee filter [61] implemented in ENVI software to the S1 images acquired in VV and VH polarizations, with a window size of 5×5 [62, 63]. The filtered images were then converted into decibels (dB) using the following:
(2)X=10×log10x,where x is the value of each pixel, X is the value converted to dB.(2) Texture Feature Extraction. SAR imagery contains texture information [5, 64]. Image classification will benefit greatly from SAR texture information features. New texture images can be reconstructed using the spatial information in the SAR image. Texture features provide details from neighboring pixels, which is useful in discriminating the various crop types found in agricultural fields [65].A gray-level cooccurrence matrix (GLCM) is considered one of the most popular techniques used for texture analysis. It is calculated based on a statistical method, also known as the gray-level spatial dependence matrix. GLCM is capable of assessing texture features by taking into account the spatial relationship between pixels and their surroundings [66, 67].From many proposed statistical texture methods, the gray-level cooccurrence matrix (GLCM) is one of the most powerful for land cover/land use (LC/LU) mapping [68, 69].The texture analysis was applied with VV and VH polarizations using the Sentinel toolbox called SNAP (SeNtinel Application Platform) based on GLCM statistics. GLCM is a probability measure of two gray levels separated by a given distance occurring in the same orientation. According to the literature [70–73], the most common and popular GLCM texture features used for tree type discrimination are as follows: correlation, mean, and variance. These three texture measurements were calculated and used for the present study.Thecorrelation function defines the similarity between the elements of columns or rows in GLCM and is a measure of the image’s gray-tone linear dependencies. The GLCM correlation feature can be calculated as follows:
(3)∑i,j=0N−1iPi,j−μxμyσxσy.The GLCMmean value is more than just the average of all original pixel values; each pixel value is weighted based on how frequently it occurs in combination with a specific neighbor pixel value. The expression of the mean is illustrated as follows:
(4)∑i,j=0N−1iPi,j.Thevariance texture parameter targets the SAR image’s partial characteristic. It can be computed as follows:
(5)∑i,j=0N−1iPi,ji−μ2, where
(i)
Pi,j is a normalized gray-tone spatial dependency matrix such as ∑i,j=0K−1iPi,j=1(ii)
i and j represent, respectively, the rows and columns for the mean, variance, and correlation measurements(iii)
μ is the mean for the variance texture measurement(iv)
N is the number of distinct gray levels in the quantized image(v)
μx,μy and σx,σy are the means and standard deviations of Px and Py respectively, for the correlation texture measure [74]


## 3.3.1. Optical Imagery

The most convenient satellite data is determined by a number of factors, including (i) spatial resolution, (ii) revisit time, (iii) image availability and cost over the study area, and (v) study area characteristics [42]. S2 is primarily aimed at freely offering high spatiotemporal resolution optical imagery for operational land cover and change detection mapping [43]. These data can be downloaded on the Theia Land Service website [44], which produces and distributes level 2A products, corrected atmospherically using the MAJA/MACCS software developed thanks to the coordination between Cnes/Cesbio and DLR [45–48].For the present study, S2 images have been downloaded from December 2018 to December 2019 over the Smimou experimental site. The 35 downloaded images, covering the whole year of 2019 and well distributed over the agricultural seasons, have been used to derive the Normalized Difference Vegetation Index (NDVI) time series.This index is the most commonly used vegetation index in remote sensing. NDVI [49] is a normalized index that generates an image showing the degree of vegetation cover (biomass). The absorption of chlorophyll pigment in the red channel and the high reflectivity of plant matter in the near infrared (NIR) channel are used to create this index, which is based on the comparison in two channels of a multispectral raster data set [50, 51]. NDVI is calculated based on the normalized difference of the NIR and R spectral bands and illustrated as follows:
(1)NDVI=NIR−RNIR+R.This index generates values between -1 and 1, mainly representing the vegetation cover, where negative values are mainly generated by clouds, water, and snow and values close to zero are mainly generated by rock and bare soil. The very low NDVI values (0.1 and below) correspond to waste rock, sand, or snow surfaces. Intermediate values (0.2 to 0.3) represent areas of shrubs and grasslands, while high values (0.6 to 0.8) indicate temperate or tropical rainforests. The transformation of the images into NDVI layers was carried out for the whole series of 35 images, in order to extract the phenological evolution of the tree species and the other land cover classes selected to be identified and mapped in this study.


## 3.3.2. SAR Imagery

Sentinel-1 (S1) is a two-satellite constellation that orbits the Earth in 6 days. The twin Sentinel-1A/B constellations use synthetic aperture radar (SAR) imaging to acquire imagery regardless of weather conditions.S1 images were acquired in the Ground Range Detected (GRD) format and interferometric wide (IW) swath mode with dual polarization (VV and VH) through the PEPS CNES websitehttps://peps.cnes.fr/. The S1 products offered by the PEPS platform are calibrated and orthorectified on the S2 grid in order to facilitate the joint use of the two missions [52–57]. Two SAR images were downloaded each month to coincide with key dates in terms of crop vegetative cycles. This produced a set of 22 SAR images (Figure 4).Figure 4
Distribution of radar and optical image acquisition dates.S1 images were preprocessed using the Sentinel Application Platform (SNAP) and ENVI software as follows: (i) the removal of the speckle noise using the enhanced Lee filter, (ii) the conversion of SAR bands' pixel values to dB values to facilitate the interpretation and the analysis of data, and (iii) the extraction of texture features using a gray-level cooccurrence matrix (GLCM).(1) Speckle Noise Filtering and Conversion to dB. The main issue in SAR imaging classification and interpretation is the speckle noise effect, which is caused by coherent interference of the wave reflected from several basic scatterers [58]. Speckle noise removal filtering is a technique for improving image quality by decreasing the speckle noise effect.The enhanced Lee filter has been chosen in this study to minimize the speckle effect on the images without loss of information. According to the literature, it is one of the most used filters to reduce speckle noise in SAR remote sensing data [59, 60].We have applied the enhanced Lee filter [61] implemented in ENVI software to the S1 images acquired in VV and VH polarizations, with a window size of 5×5 [62, 63]. The filtered images were then converted into decibels (dB) using the following:
(2)X=10×log10x,where x is the value of each pixel, X is the value converted to dB.(2) Texture Feature Extraction. SAR imagery contains texture information [5, 64]. Image classification will benefit greatly from SAR texture information features. New texture images can be reconstructed using the spatial information in the SAR image. Texture features provide details from neighboring pixels, which is useful in discriminating the various crop types found in agricultural fields [65].A gray-level cooccurrence matrix (GLCM) is considered one of the most popular techniques used for texture analysis. It is calculated based on a statistical method, also known as the gray-level spatial dependence matrix. GLCM is capable of assessing texture features by taking into account the spatial relationship between pixels and their surroundings [66, 67].From many proposed statistical texture methods, the gray-level cooccurrence matrix (GLCM) is one of the most powerful for land cover/land use (LC/LU) mapping [68, 69].The texture analysis was applied with VV and VH polarizations using the Sentinel toolbox called SNAP (SeNtinel Application Platform) based on GLCM statistics. GLCM is a probability measure of two gray levels separated by a given distance occurring in the same orientation. According to the literature [70–73], the most common and popular GLCM texture features used for tree type discrimination are as follows: correlation, mean, and variance. These three texture measurements were calculated and used for the present study.Thecorrelation function defines the similarity between the elements of columns or rows in GLCM and is a measure of the image’s gray-tone linear dependencies. The GLCM correlation feature can be calculated as follows:
(3)∑i,j=0N−1iPi,j−μxμyσxσy.The GLCMmean value is more than just the average of all original pixel values; each pixel value is weighted based on how frequently it occurs in combination with a specific neighbor pixel value. The expression of the mean is illustrated as follows:
(4)∑i,j=0N−1iPi,j.Thevariance texture parameter targets the SAR image’s partial characteristic. It can be computed as follows:
(5)∑i,j=0N−1iPi,ji−μ2, where
(i)
Pi,j is a normalized gray-tone spatial dependency matrix such as ∑i,j=0K−1iPi,j=1(ii)
i and j represent, respectively, the rows and columns for the mean, variance, and correlation measurements(iii)
μ is the mean for the variance texture measurement(iv)
N is the number of distinct gray levels in the quantized image(v)
μx,μy and σx,σy are the means and standard deviations of Px and Py respectively, for the correlation texture measure [74]


## 3.4. Methodology

The methodological approach followed in this work (Figure5) consists of 5 main steps:
(1)
Data acquisition(2)
Data preprocessing(3)
Extraction and analysis of profiles and classification(4)
Accuracy assessment and postclassification(5)
Comparison, analysis, and map productionFigure 5
Methodological approach followed for the present study.


### 3.4.1. SVM Classification

SVM is one of the most appealing algorithms for remotely sensed data classification because of the benefits of generalization even with limited training samples, which is common in remote sensing data processing [26, 75]. Furthermore, because SVM is a supervised nonparametric statistical learning method, it does not require a training set that strictly follows the standard independent and identical distribution. The benefits of SVM stem from two factors: transforming the original space training set into a very high dimensional new space and discovering a large linear margin boundary in the new space [76]. SVM is a classifier based on the structural risk minimization theory that is aimed at reducing generalization error by maximizing margins on training data [77]. SVM classifiers can use a variety of kernel functions. However, only four of them have been widely used to classify satellite data: polynomial, linear, radial basis function (RBF), and sigmoid kernels [77]. In the present study, The SVM-RBF kernel was selected and used for all image classifications as one of the most effective machine learning algorithms for classifying forest tree species using satellite imagery [14–16].


### 3.4.2. Classification Scenarios

ENVI Classification Scenarios is a professional software for classifying satellite imagery using both unsupervised and supervised algorithms. For this work, we have limited ourselves to using supervised methods, in particular the support vector machines (SVM). The techniques generally provide better results in the production of land cover maps [78–80]. As a first step, we started by classifying the time series of NDVI, VV, VH, and VH/VV. We then combined the data from the two sensors to explore different possible scenarios (combinations). Table 2 describes these scenarios in the order followed during the image classification step.Table 2
Scenarios created for the classification step.
ScenarioBandsSatellite sensorType1NDVIS2Optical2VVS1SAR3VHS1SAR4VH/VVS1SAR5VV, VHS1SAR6VV, VH/VVS1SAR7VH, VH/VVS1SAR8VV, VH, VH/VVS1SAR9TextureS1SAR10Texture, VH, VH/VVS1SAR11Texture, VH, VVS1SAR12Texture, VV, VH/VVS1SAR13NDVI, VVS1 and S2Optical and SAR14NDVI, VHS1 and S2Optical and SAR15NDVI, VH/VVS1 and S2Optical and SAR16NDVI, VV, VHS1 and S2Optical and SAR17NDVI, VV, VH/VVS1 and S2Optical and SAR18NDVI, VH, VH/VVS1 and S2Optical and SAR19NDVI, VV, VH, VH/VVS1 and S2Optical and SAR20NDVI, textureS1 and S2Optical and SAR21NDVI, texture, VH, VH/VVS1 and S2Optical and SAR22NDVI, texture, VVS1 and S2Optical and SAR


### 3.4.3. Performance Assessment

Evaluation of classification performance was done by comparing images and by calculating the accuracy of the land use classification using statistical tools such as the confusion matrix and the kappa index (K) introduced by Cohen [81], which is used to validate and estimate the precision in the vegetation classification scheme [80]. Its formula is as follows:
(6)Kappa=P0−PcPp−Pc, where
(i)
P0: is equal to the actual percentage obtained from classification of land covers; it is equal to the quotient of the sum of the diagonal digits of the matrix with the total number of observations(ii)
Pc is the estimate of the probability of obtaining a correct classification. For a correct classification, the value of Pc is generally lower than P0(iii)
Pp: is the percentage obtained when the classification is perfectConfusion matrices have been calculated in order to evaluate the results, these matrices reveal not only the general errors made at the level of each class during interpretation but also errors due to confusion between land use classes. The mathematical expression used to calculate the OA of the classification is given by the average of the percentages of correctly classified pixels:(7)OA=∑i=1rXiiN×100, where
(i)
Xii is the number of diagonal pixels (i.e., correctly classified)(ii)
N is the total number of pixels(iii)
R is the total number of classes


### 3.4.4. Postclassification: Smoothing

After the classification, the obtained images showed a lack of sharpness in the definition of the classified plots; the images always seem “noisy.” We then apply a smoothing filter, called majority/minority analysis [82], to the images using the ENVI software. This analysis filters the image by replacing the value of the center pixel in a window with a size of n×n (where n has to be defined) with the majority value in that window. In order to avoid “oversmoothing,” we have chosen a window with a size of 3×3.


## 3.4.1. SVM Classification

SVM is one of the most appealing algorithms for remotely sensed data classification because of the benefits of generalization even with limited training samples, which is common in remote sensing data processing [26, 75]. Furthermore, because SVM is a supervised nonparametric statistical learning method, it does not require a training set that strictly follows the standard independent and identical distribution. The benefits of SVM stem from two factors: transforming the original space training set into a very high dimensional new space and discovering a large linear margin boundary in the new space [76]. SVM is a classifier based on the structural risk minimization theory that is aimed at reducing generalization error by maximizing margins on training data [77]. SVM classifiers can use a variety of kernel functions. However, only four of them have been widely used to classify satellite data: polynomial, linear, radial basis function (RBF), and sigmoid kernels [77]. In the present study, The SVM-RBF kernel was selected and used for all image classifications as one of the most effective machine learning algorithms for classifying forest tree species using satellite imagery [14–16].


## 3.4.2. Classification Scenarios

ENVI Classification Scenarios is a professional software for classifying satellite imagery using both unsupervised and supervised algorithms. For this work, we have limited ourselves to using supervised methods, in particular the support vector machines (SVM). The techniques generally provide better results in the production of land cover maps [78–80]. As a first step, we started by classifying the time series of NDVI, VV, VH, and VH/VV. We then combined the data from the two sensors to explore different possible scenarios (combinations). Table 2 describes these scenarios in the order followed during the image classification step.Table 2
Scenarios created for the classification step.
ScenarioBandsSatellite sensorType1NDVIS2Optical2VVS1SAR3VHS1SAR4VH/VVS1SAR5VV, VHS1SAR6VV, VH/VVS1SAR7VH, VH/VVS1SAR8VV, VH, VH/VVS1SAR9TextureS1SAR10Texture, VH, VH/VVS1SAR11Texture, VH, VVS1SAR12Texture, VV, VH/VVS1SAR13NDVI, VVS1 and S2Optical and SAR14NDVI, VHS1 and S2Optical and SAR15NDVI, VH/VVS1 and S2Optical and SAR16NDVI, VV, VHS1 and S2Optical and SAR17NDVI, VV, VH/VVS1 and S2Optical and SAR18NDVI, VH, VH/VVS1 and S2Optical and SAR19NDVI, VV, VH, VH/VVS1 and S2Optical and SAR20NDVI, textureS1 and S2Optical and SAR21NDVI, texture, VH, VH/VVS1 and S2Optical and SAR22NDVI, texture, VVS1 and S2Optical and SAR


## 3.4.3. Performance Assessment

Evaluation of classification performance was done by comparing images and by calculating the accuracy of the land use classification using statistical tools such as the confusion matrix and the kappa index (K) introduced by Cohen [81], which is used to validate and estimate the precision in the vegetation classification scheme [80]. Its formula is as follows:
(6)Kappa=P0−PcPp−Pc, where
(i)
P0: is equal to the actual percentage obtained from classification of land covers; it is equal to the quotient of the sum of the diagonal digits of the matrix with the total number of observations(ii)
Pc is the estimate of the probability of obtaining a correct classification. For a correct classification, the value of Pc is generally lower than P0(iii)
Pp: is the percentage obtained when the classification is perfectConfusion matrices have been calculated in order to evaluate the results, these matrices reveal not only the general errors made at the level of each class during interpretation but also errors due to confusion between land use classes. The mathematical expression used to calculate the OA of the classification is given by the average of the percentages of correctly classified pixels:(7)OA=∑i=1rXiiN×100, where
(i)
Xii is the number of diagonal pixels (i.e., correctly classified)(ii)
N is the total number of pixels(iii)
R is the total number of classes


## 3.4.4. Postclassification: Smoothing

After the classification, the obtained images showed a lack of sharpness in the definition of the classified plots; the images always seem “noisy.” We then apply a smoothing filter, called majority/minority analysis [82], to the images using the ENVI software. This analysis filters the image by replacing the value of the center pixel in a window with a size of n×n (where n has to be defined) with the majority value in that window. In order to avoid “oversmoothing,” we have chosen a window with a size of 3×3.


## 4. Results and Discussion

### 4.1. Analysis of Temporal Profiles

#### 4.1.1. NDVI Profiles

For each of the selected classes, the NDVI, VV, VH, and VH/VV profiles were calculated, in order to study the possible confusions and/or separability to be expected between these different classes, and to be able to interpret and analyze the classification results. The temporal evolution of the NDVI profiles allows the modeling of land cover type dynamics, and particularly the phenological evolution of the different vegetative classes, which makes them relatively interpretable. On the other hand, the SAR profiles are more complex for the vegetative classes, reflecting the evolution of the proportion of the backscattered signal, mainly influenced by the surface roughness of the canopy as well as its water content. Figure6 shows, for each of the selected classes, the temporal evolution of NDVI, and that of the C-band backscattered signals (in VV, VH, and VH/VV polarizations).Figure 6
NDVI, VV, VH, and VH/VV profiles of the different selected classes present in the Smimou zone.NDVI is a measurement of the normalized difference between the energy emitted in infrared and red by surface objects, when applied to vegetation communities. This index establishes a value indicating the importance or dominance of vegetation in a given area. For annual crops, as an example, the start of the season (i.e., the start date of the vegetation growth phase) is the point at which the rate of increase in NDVI values is higher than previous successive observations during the vegetation growth period. While the end of the season is defined as the point in the maturity period when a significant drop in NDVI is observed.These changes can be seen in Figure6 for fallow and agricultural land, where the NDVI spectral profiles of these crops illustrate their seasonal aspects. The amplitude of the graphs becomes important during periods of vegetation development. Agricultural lands in this region grow vegetables that develop in January, but they represent only small plots with produce intended for local consumption. Fallow land resembles cereals that are cultivated in early December, reach their maximum development in late March, and are harvested in mid-May and early June at the latest, which is due to the inputs of winter rains.Argan, sandarac gum, argan + sandarac gum, and olive form a group of tree species present in the study area; however, we have separated them. The NDVI value of these classes is generally all year round greater than 0.2 (Figure6), given that they are perennial crops. Decreases in the NDVI values of these trees in summer are due to water stress resulting from lack of water and increased temperatures during this period, and also sometimes due to loss of leaves through pruning. The temporal profile of the mixture (argan + sandarac gum) is approximately the same as that of the olive trees.The classes of bare soil and noncrop are relatively easy to detect due to the absence of vegetation translated by NDVI values of no more than 0.20 (Table3). These two classes have been combined into one class called “bare soil.”Table 3
NDVI ranges and amplitudes of all classes during 2019.
ClassMaxMinAmplitudeArgan0.480.170.31Sandarac gum0.810.430.38Olive0.530.220.31Fallow0.370.110.26Argan + sandarac gum0.480.280.20Agricultural land0.620.120.50Bare soil0.200.100.10


#### 4.1.2. VV, VH, and VH/VV Profiles

The VV, VH, and VH/VV profiles presented in Figure6 illustrate the temporal evolution in intensity (radar backscatter coefficient σ expressed in dB of the VV, VH, and VH/VV polarizations) of the 7 cover types selected: agricultural land, olive, sandarac gum, fallow, argan, argan + sandarac gum, and bare soil. First of all, we can clearly see that for all the classes studied (Figure 6), the values of intensity in parallel polarization VV (values between -3 dB and -10 dB) are above those in cross polarization VH (values between -9 dB and -17 dB), while the values of intensity of the VH/VV ratio (in dB) are positive and above those of the two previous polarizations (varying between 4 dB and 7 dB), and this for all the classes. Moreover, even if the graphs look similar, the ranges of the signal values are different. We can also notice the oscillatory behavior of the argan and sandarac gum classes and their mixture, the reasons for which are still unknown.In contrast to NDVI profiles, radar profiles seem to be more complicated to interpret. The two main physical variables influencing the backscattering process are water content and roughness, and since this latter varies from one class to another, in addition to the fact that the region is dry, it can be said that the similarity between the backscattered signals of the different classes is due to the low water content of the observed surfaces. This hypothesis is supported by several bibliographical works [83–85] aimed at mapping the rice fields, and which have shown that the temporal change in rice field backscatter is very high compared to other classes, given the amounts of water present on flooded fields.


### 4.2. Classification Accuracies

Table4 presented the results corresponding to the scenarios of uncombined Sentinel products. The classification results of the products from each of the two sensors S1 and S2 confirm the superiority of the optical data and the SVM algorithm in terms of LC/LU classification and performance.Table 4
Classification results of uncombined sentinel products.
ScenariosSVMOA (%)KNDVI86.900.84VV27.140.15VH26.400.15VH/VV28.950.17VH, VV33.850.24VH, VH/VV32.150.22VV, VH/VV31.860.21VH, VV, VH/VV32.680.22Texture36.650.28Texture, VH, VH/VV37.220.28Texture, VH, VV35.910.27Texture, VV, VH/VV36.580.27The OA resulting from the classification of the NDVI time series is 86.90% (K=0.84), while the best values of OA that have been obtained by classifying the SAR time series and more precisely scenarios (VV and VH) and (VV, VH, and VH/VV) are only 33.85% (K=0.24) and 32.68% (K=0.22), respectively. We note that the addition of the VH/VV to the VH or VV band increases the accuracy by about 5% to 6%, while the inclusion of textural features improved the classification to 37.22% (K=0.29). Taking into consideration that SAR sensors are not primarily intended for land cover mapping, the classification of SAR products in this area presents much lower results than those found by optical data (NDVI).The classification results of combined optical and SAR products are grouped in Table5.Table 5
Classification results of combined S1 and S2 images.
ScenariosSVMOAKNDVI, VH82.460.79NDVI, VV86.510.84NDVI, VH/VV84.980.82NDVI, VH, VV83.030.80NDVI, VH, VH/VV81.970.78NDVI, VV, VH/VV84.270.80NDVI,VH, VV, VH/VV81.470.78NDVI, texture58.620.52NDVI, texture, VH, VH/VV67.240.62NDVI, texture, VV68.020.63The combination of S1 and S2 products, compared to S2 classification, does not improve the classification accuracy. On the contrary, in the study area, the integration of the two types of satellite data gave practically the same results as the classification of optical data.In order to visually compare the quality of the classified images, we present in Figure7 the images obtained for the scenarios reached good OA.Figure 7
Maps resulting from the best classifications.Indeed, we observed differences between the classifications of S1 and S2 products. The quality of the classified SAR images is lower than that of the classified NDVI image in terms of sharpness even if the speckle effect has been reduced after applying the filter. This can be explained by the huge confusion of the argan class with the olive tree, which appears to be dominant on the maps resulting from the classification of SAR products. The combined products, as expected, had quite sharp images, which is due to the complementary contribution of the optical data. The presence of isolated pixels on homogeneous plots over the entire area contributes to the “noisy” appearance of the classified SAR products. Therefore, we have tried to improve the quality of these images and the accuracy of the results by applying a postprocessing which consists of smoothing the classified products.Visually, the quality of image classification has improved as shown in Figure8. The smoothed product looks more like the classified image from the NDVI time series. This process has been applied to all scenarios (optical and SAR classified images) and has given satisfactory results. Table 6 shows the OA and K coefficient for these scenarios before and after the smoothing process has been applied.Figure 8
Smoothing of scenario classifications (NDVI), (VH and VV), (VV, VH, and VV/VH), and (VH, VH/VV, and texture).Table 6
Improvement of OA andK after application of smoothing.
ScenariosBefore smoothingAfter smoothingOA (%)KOA (%)KNDVI86.900.8489.780.88VV, VH33.850.2435.940.26NDVI, VV86.510.8488.140.86NDVI, VH/VV84.980.8287.360.85NDVI, VH, VV77.890.7481.380.78VH, VV, VH/VV32.680.2234.450.24NDVI, VV, VH/VV80.990.7883.220.80Texture, VH, VH/VV37.220.2838 .040.29Table7 shows the producer accuracy (PA) and user accuracy (UA) of classifications after smoothing.Table 7
UA and PA of the best results of the smoothed image.
VV, VHVV, VH, VH/VVNDVI, VH/VVNDVI, VVNDVI, VV, VH/VVNDVI, VV, VHNDVITexture, VH, VH/VVClassPA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)Argan32.1643.1629.8042.7089.8072.9386.6776.4783.4970.1364.7367.8990.2082.1446.2730.18Bare soil12.3322.2512.0219.8593.0092.1695.2892.4785.5388.4385.6987.2394.8292.994.2622.58Noncrop52.9148.0850.7944.8686.2488.5989.4296.0267.6374.8792.6292.0089.9593.9262.9634.00Olive49.7516.2654.7316.6447.2648.9771.1452.5778.1950.4064.0938.9355.2259.6863.1823.30Fallow46.0340.8944.7741.9690.3898.6395.8298.2882.2987.8063.0892.4892.47100.0056.9047.06Agricultural areas32.7232.1326.8431.2099.6391.2598.5392.1094.6387.3198.5887.7699.6391.2532.3528.48Sandarac gum54.5238.6451.6037.5899.4297.9999.1397.4298.9595.1998.4298.42100.0096.8952.4848.26Argan + sandarac gum39.8864.2338.2262.1680.9789.4873.7291.9074.1896.7979.8590.1184.7489.1941.6963.16We ended with the elaboration of three tree species maps (Figure9), the first one derived from NDVI, the second one derived from the best scenario of SAR data (VH, VH/VV, and texture), and the third one derived from the best scenario of combined data (multisensor data) (NDVI, VV). Finally, as the goal of the present work, we developed the spatial distribution maps of the argan tree and other tree species present in the Smimou zone.Figure 9
Maps of the spatial distribution of the argan tree using (a) optical data (OA=89.78%), (b) SAR data (OA=38.04%), and (c) combined optical and SAR data (OA=88.14%).


### 4.3. Discussion

In general, the NDVI profiles model the phenological behavior of the classes, and compared to the VV, VH, and VH/VV profiles, they show a high separability of the classes which facilitates their discrimination when applying the classification algorithms. By quantifying the results using the confusion matrices, the separability was translated by the difference shown by the accuracy statistics (OA andK indices). The highest classification accuracy was recorded for the optical data (OA=89.78% and K=0.88) compared to fused multisensor data (OA=88.14% and K=0.86), while the classification of the best case scenario of SAR data (VH, VH/VV, and texture features) gave an OA of 38.04% and a K of 0.29. The PA values of a single SAR classification before the integration of texture features, shown in Table 7, illustrate that the most misclassified cover types are bare soil, agricultural areas, and argan trees. These confusions can be explained by the fact that the argan trees can be found on cultivated land, in pure forests at different densities, with or without undergrowth, and in mixed forests. However, the integration of texture features (scenario combining texture, VH, and VH/VV) increased the discrimination of the argan tree (PA=46.27% instead of PA=29.80% and PA=32.16% for the scenarios VV, VH, VH/VV, and VV, VH, respectively). This confirms the contribution of texture features in the improvement of the tree species classification accuracy using SAR data [6, 8, 68, 71].On the other hand, concerning the UA values, the olive tree is the most confused tree type (UA value between 16% and 23%). Olive trees in the study area are generally planted around the small villages for local consumption with different densities and with or without annual crops. Therefore, the percentage of UA for bare soil and agricultural areas presents also low values.For optical classification, the PA and UA values (NDVI time series derived from S2 data) were greater than 70% for all land covers except for olive trees. The misclassification of this tree type can be explained by the fact that there are some few zones in the study area where argan is mixed with olive trees.For the fused multisensor data (NDVI, VV which is the best case scenario of combined SAR and optical data), the results of classifications (OA=88.14% and K=0.86) were lower than those of optical classification (89.78% and K=0.88), contrary to what was expected [1, 6, 85]. As shown in Table 7, the PA and UA values for the argan tree (PA=90.20% and UA=82.14%), using NDVI classification, are better than those obtained with combined multisensor data (PA=86.67% and UA=78.47%).However, Sirro et al. [13] used optical and SAR data for forest and land cover classification and revealed that the classification accuracy of optical data was better compared to fused multisensor data. This study aligns with the results that we found in the present study.S1 and S2 were combined in several studies [86–88] for crop type classification and confirmed that the fusion of the two types of data outperformed the classification based on optical sensor images.


## 4.1. Analysis of Temporal Profiles

### 4.1.1. NDVI Profiles

For each of the selected classes, the NDVI, VV, VH, and VH/VV profiles were calculated, in order to study the possible confusions and/or separability to be expected between these different classes, and to be able to interpret and analyze the classification results. The temporal evolution of the NDVI profiles allows the modeling of land cover type dynamics, and particularly the phenological evolution of the different vegetative classes, which makes them relatively interpretable. On the other hand, the SAR profiles are more complex for the vegetative classes, reflecting the evolution of the proportion of the backscattered signal, mainly influenced by the surface roughness of the canopy as well as its water content. Figure6 shows, for each of the selected classes, the temporal evolution of NDVI, and that of the C-band backscattered signals (in VV, VH, and VH/VV polarizations).Figure 6
NDVI, VV, VH, and VH/VV profiles of the different selected classes present in the Smimou zone.NDVI is a measurement of the normalized difference between the energy emitted in infrared and red by surface objects, when applied to vegetation communities. This index establishes a value indicating the importance or dominance of vegetation in a given area. For annual crops, as an example, the start of the season (i.e., the start date of the vegetation growth phase) is the point at which the rate of increase in NDVI values is higher than previous successive observations during the vegetation growth period. While the end of the season is defined as the point in the maturity period when a significant drop in NDVI is observed.These changes can be seen in Figure6 for fallow and agricultural land, where the NDVI spectral profiles of these crops illustrate their seasonal aspects. The amplitude of the graphs becomes important during periods of vegetation development. Agricultural lands in this region grow vegetables that develop in January, but they represent only small plots with produce intended for local consumption. Fallow land resembles cereals that are cultivated in early December, reach their maximum development in late March, and are harvested in mid-May and early June at the latest, which is due to the inputs of winter rains.Argan, sandarac gum, argan + sandarac gum, and olive form a group of tree species present in the study area; however, we have separated them. The NDVI value of these classes is generally all year round greater than 0.2 (Figure6), given that they are perennial crops. Decreases in the NDVI values of these trees in summer are due to water stress resulting from lack of water and increased temperatures during this period, and also sometimes due to loss of leaves through pruning. The temporal profile of the mixture (argan + sandarac gum) is approximately the same as that of the olive trees.The classes of bare soil and noncrop are relatively easy to detect due to the absence of vegetation translated by NDVI values of no more than 0.20 (Table3). These two classes have been combined into one class called “bare soil.”Table 3
NDVI ranges and amplitudes of all classes during 2019.
ClassMaxMinAmplitudeArgan0.480.170.31Sandarac gum0.810.430.38Olive0.530.220.31Fallow0.370.110.26Argan + sandarac gum0.480.280.20Agricultural land0.620.120.50Bare soil0.200.100.10


### 4.1.2. VV, VH, and VH/VV Profiles

The VV, VH, and VH/VV profiles presented in Figure6 illustrate the temporal evolution in intensity (radar backscatter coefficient σ expressed in dB of the VV, VH, and VH/VV polarizations) of the 7 cover types selected: agricultural land, olive, sandarac gum, fallow, argan, argan + sandarac gum, and bare soil. First of all, we can clearly see that for all the classes studied (Figure 6), the values of intensity in parallel polarization VV (values between -3 dB and -10 dB) are above those in cross polarization VH (values between -9 dB and -17 dB), while the values of intensity of the VH/VV ratio (in dB) are positive and above those of the two previous polarizations (varying between 4 dB and 7 dB), and this for all the classes. Moreover, even if the graphs look similar, the ranges of the signal values are different. We can also notice the oscillatory behavior of the argan and sandarac gum classes and their mixture, the reasons for which are still unknown.In contrast to NDVI profiles, radar profiles seem to be more complicated to interpret. The two main physical variables influencing the backscattering process are water content and roughness, and since this latter varies from one class to another, in addition to the fact that the region is dry, it can be said that the similarity between the backscattered signals of the different classes is due to the low water content of the observed surfaces. This hypothesis is supported by several bibliographical works [83–85] aimed at mapping the rice fields, and which have shown that the temporal change in rice field backscatter is very high compared to other classes, given the amounts of water present on flooded fields.


## 4.1.1. NDVI Profiles

For each of the selected classes, the NDVI, VV, VH, and VH/VV profiles were calculated, in order to study the possible confusions and/or separability to be expected between these different classes, and to be able to interpret and analyze the classification results. The temporal evolution of the NDVI profiles allows the modeling of land cover type dynamics, and particularly the phenological evolution of the different vegetative classes, which makes them relatively interpretable. On the other hand, the SAR profiles are more complex for the vegetative classes, reflecting the evolution of the proportion of the backscattered signal, mainly influenced by the surface roughness of the canopy as well as its water content. Figure6 shows, for each of the selected classes, the temporal evolution of NDVI, and that of the C-band backscattered signals (in VV, VH, and VH/VV polarizations).Figure 6
NDVI, VV, VH, and VH/VV profiles of the different selected classes present in the Smimou zone.NDVI is a measurement of the normalized difference between the energy emitted in infrared and red by surface objects, when applied to vegetation communities. This index establishes a value indicating the importance or dominance of vegetation in a given area. For annual crops, as an example, the start of the season (i.e., the start date of the vegetation growth phase) is the point at which the rate of increase in NDVI values is higher than previous successive observations during the vegetation growth period. While the end of the season is defined as the point in the maturity period when a significant drop in NDVI is observed.These changes can be seen in Figure6 for fallow and agricultural land, where the NDVI spectral profiles of these crops illustrate their seasonal aspects. The amplitude of the graphs becomes important during periods of vegetation development. Agricultural lands in this region grow vegetables that develop in January, but they represent only small plots with produce intended for local consumption. Fallow land resembles cereals that are cultivated in early December, reach their maximum development in late March, and are harvested in mid-May and early June at the latest, which is due to the inputs of winter rains.Argan, sandarac gum, argan + sandarac gum, and olive form a group of tree species present in the study area; however, we have separated them. The NDVI value of these classes is generally all year round greater than 0.2 (Figure6), given that they are perennial crops. Decreases in the NDVI values of these trees in summer are due to water stress resulting from lack of water and increased temperatures during this period, and also sometimes due to loss of leaves through pruning. The temporal profile of the mixture (argan + sandarac gum) is approximately the same as that of the olive trees.The classes of bare soil and noncrop are relatively easy to detect due to the absence of vegetation translated by NDVI values of no more than 0.20 (Table3). These two classes have been combined into one class called “bare soil.”Table 3
NDVI ranges and amplitudes of all classes during 2019.
ClassMaxMinAmplitudeArgan0.480.170.31Sandarac gum0.810.430.38Olive0.530.220.31Fallow0.370.110.26Argan + sandarac gum0.480.280.20Agricultural land0.620.120.50Bare soil0.200.100.10


## 4.1.2. VV, VH, and VH/VV Profiles

The VV, VH, and VH/VV profiles presented in Figure6 illustrate the temporal evolution in intensity (radar backscatter coefficient σ expressed in dB of the VV, VH, and VH/VV polarizations) of the 7 cover types selected: agricultural land, olive, sandarac gum, fallow, argan, argan + sandarac gum, and bare soil. First of all, we can clearly see that for all the classes studied (Figure 6), the values of intensity in parallel polarization VV (values between -3 dB and -10 dB) are above those in cross polarization VH (values between -9 dB and -17 dB), while the values of intensity of the VH/VV ratio (in dB) are positive and above those of the two previous polarizations (varying between 4 dB and 7 dB), and this for all the classes. Moreover, even if the graphs look similar, the ranges of the signal values are different. We can also notice the oscillatory behavior of the argan and sandarac gum classes and their mixture, the reasons for which are still unknown.In contrast to NDVI profiles, radar profiles seem to be more complicated to interpret. The two main physical variables influencing the backscattering process are water content and roughness, and since this latter varies from one class to another, in addition to the fact that the region is dry, it can be said that the similarity between the backscattered signals of the different classes is due to the low water content of the observed surfaces. This hypothesis is supported by several bibliographical works [83–85] aimed at mapping the rice fields, and which have shown that the temporal change in rice field backscatter is very high compared to other classes, given the amounts of water present on flooded fields.


## 4.2. Classification Accuracies

Table4 presented the results corresponding to the scenarios of uncombined Sentinel products. The classification results of the products from each of the two sensors S1 and S2 confirm the superiority of the optical data and the SVM algorithm in terms of LC/LU classification and performance.Table 4
Classification results of uncombined sentinel products.
ScenariosSVMOA (%)KNDVI86.900.84VV27.140.15VH26.400.15VH/VV28.950.17VH, VV33.850.24VH, VH/VV32.150.22VV, VH/VV31.860.21VH, VV, VH/VV32.680.22Texture36.650.28Texture, VH, VH/VV37.220.28Texture, VH, VV35.910.27Texture, VV, VH/VV36.580.27The OA resulting from the classification of the NDVI time series is 86.90% (K=0.84), while the best values of OA that have been obtained by classifying the SAR time series and more precisely scenarios (VV and VH) and (VV, VH, and VH/VV) are only 33.85% (K=0.24) and 32.68% (K=0.22), respectively. We note that the addition of the VH/VV to the VH or VV band increases the accuracy by about 5% to 6%, while the inclusion of textural features improved the classification to 37.22% (K=0.29). Taking into consideration that SAR sensors are not primarily intended for land cover mapping, the classification of SAR products in this area presents much lower results than those found by optical data (NDVI).The classification results of combined optical and SAR products are grouped in Table5.Table 5
Classification results of combined S1 and S2 images.
ScenariosSVMOAKNDVI, VH82.460.79NDVI, VV86.510.84NDVI, VH/VV84.980.82NDVI, VH, VV83.030.80NDVI, VH, VH/VV81.970.78NDVI, VV, VH/VV84.270.80NDVI,VH, VV, VH/VV81.470.78NDVI, texture58.620.52NDVI, texture, VH, VH/VV67.240.62NDVI, texture, VV68.020.63The combination of S1 and S2 products, compared to S2 classification, does not improve the classification accuracy. On the contrary, in the study area, the integration of the two types of satellite data gave practically the same results as the classification of optical data.In order to visually compare the quality of the classified images, we present in Figure7 the images obtained for the scenarios reached good OA.Figure 7
Maps resulting from the best classifications.Indeed, we observed differences between the classifications of S1 and S2 products. The quality of the classified SAR images is lower than that of the classified NDVI image in terms of sharpness even if the speckle effect has been reduced after applying the filter. This can be explained by the huge confusion of the argan class with the olive tree, which appears to be dominant on the maps resulting from the classification of SAR products. The combined products, as expected, had quite sharp images, which is due to the complementary contribution of the optical data. The presence of isolated pixels on homogeneous plots over the entire area contributes to the “noisy” appearance of the classified SAR products. Therefore, we have tried to improve the quality of these images and the accuracy of the results by applying a postprocessing which consists of smoothing the classified products.Visually, the quality of image classification has improved as shown in Figure8. The smoothed product looks more like the classified image from the NDVI time series. This process has been applied to all scenarios (optical and SAR classified images) and has given satisfactory results. Table 6 shows the OA and K coefficient for these scenarios before and after the smoothing process has been applied.Figure 8
Smoothing of scenario classifications (NDVI), (VH and VV), (VV, VH, and VV/VH), and (VH, VH/VV, and texture).Table 6
Improvement of OA andK after application of smoothing.
ScenariosBefore smoothingAfter smoothingOA (%)KOA (%)KNDVI86.900.8489.780.88VV, VH33.850.2435.940.26NDVI, VV86.510.8488.140.86NDVI, VH/VV84.980.8287.360.85NDVI, VH, VV77.890.7481.380.78VH, VV, VH/VV32.680.2234.450.24NDVI, VV, VH/VV80.990.7883.220.80Texture, VH, VH/VV37.220.2838 .040.29Table7 shows the producer accuracy (PA) and user accuracy (UA) of classifications after smoothing.Table 7
UA and PA of the best results of the smoothed image.
VV, VHVV, VH, VH/VVNDVI, VH/VVNDVI, VVNDVI, VV, VH/VVNDVI, VV, VHNDVITexture, VH, VH/VVClassPA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)PA (%)UA (%)Argan32.1643.1629.8042.7089.8072.9386.6776.4783.4970.1364.7367.8990.2082.1446.2730.18Bare soil12.3322.2512.0219.8593.0092.1695.2892.4785.5388.4385.6987.2394.8292.994.2622.58Noncrop52.9148.0850.7944.8686.2488.5989.4296.0267.6374.8792.6292.0089.9593.9262.9634.00Olive49.7516.2654.7316.6447.2648.9771.1452.5778.1950.4064.0938.9355.2259.6863.1823.30Fallow46.0340.8944.7741.9690.3898.6395.8298.2882.2987.8063.0892.4892.47100.0056.9047.06Agricultural areas32.7232.1326.8431.2099.6391.2598.5392.1094.6387.3198.5887.7699.6391.2532.3528.48Sandarac gum54.5238.6451.6037.5899.4297.9999.1397.4298.9595.1998.4298.42100.0096.8952.4848.26Argan + sandarac gum39.8864.2338.2262.1680.9789.4873.7291.9074.1896.7979.8590.1184.7489.1941.6963.16We ended with the elaboration of three tree species maps (Figure9), the first one derived from NDVI, the second one derived from the best scenario of SAR data (VH, VH/VV, and texture), and the third one derived from the best scenario of combined data (multisensor data) (NDVI, VV). Finally, as the goal of the present work, we developed the spatial distribution maps of the argan tree and other tree species present in the Smimou zone.Figure 9
Maps of the spatial distribution of the argan tree using (a) optical data (OA=89.78%), (b) SAR data (OA=38.04%), and (c) combined optical and SAR data (OA=88.14%).


## 4.3. Discussion

In general, the NDVI profiles model the phenological behavior of the classes, and compared to the VV, VH, and VH/VV profiles, they show a high separability of the classes which facilitates their discrimination when applying the classification algorithms. By quantifying the results using the confusion matrices, the separability was translated by the difference shown by the accuracy statistics (OA andK indices). The highest classification accuracy was recorded for the optical data (OA=89.78% and K=0.88) compared to fused multisensor data (OA=88.14% and K=0.86), while the classification of the best case scenario of SAR data (VH, VH/VV, and texture features) gave an OA of 38.04% and a K of 0.29. The PA values of a single SAR classification before the integration of texture features, shown in Table 7, illustrate that the most misclassified cover types are bare soil, agricultural areas, and argan trees. These confusions can be explained by the fact that the argan trees can be found on cultivated land, in pure forests at different densities, with or without undergrowth, and in mixed forests. However, the integration of texture features (scenario combining texture, VH, and VH/VV) increased the discrimination of the argan tree (PA=46.27% instead of PA=29.80% and PA=32.16% for the scenarios VV, VH, VH/VV, and VV, VH, respectively). This confirms the contribution of texture features in the improvement of the tree species classification accuracy using SAR data [6, 8, 68, 71].On the other hand, concerning the UA values, the olive tree is the most confused tree type (UA value between 16% and 23%). Olive trees in the study area are generally planted around the small villages for local consumption with different densities and with or without annual crops. Therefore, the percentage of UA for bare soil and agricultural areas presents also low values.For optical classification, the PA and UA values (NDVI time series derived from S2 data) were greater than 70% for all land covers except for olive trees. The misclassification of this tree type can be explained by the fact that there are some few zones in the study area where argan is mixed with olive trees.For the fused multisensor data (NDVI, VV which is the best case scenario of combined SAR and optical data), the results of classifications (OA=88.14% and K=0.86) were lower than those of optical classification (89.78% and K=0.88), contrary to what was expected [1, 6, 85]. As shown in Table 7, the PA and UA values for the argan tree (PA=90.20% and UA=82.14%), using NDVI classification, are better than those obtained with combined multisensor data (PA=86.67% and UA=78.47%).However, Sirro et al. [13] used optical and SAR data for forest and land cover classification and revealed that the classification accuracy of optical data was better compared to fused multisensor data. This study aligns with the results that we found in the present study.S1 and S2 were combined in several studies [86–88] for crop type classification and confirmed that the fusion of the two types of data outperformed the classification based on optical sensor images.


## 5. Conclusions

LC/LU mapping is a crucial tool to understand, analyze, and monitor land cover dynamics in order to better manage territories. Moreover, the availability of satellite imagery, which has the advantage of covering surfaces at all scales (local, regional, and continental), offers the possibility of carrying out this mapping as the only way to monitor large areas, and the principle technical tool put at the service of this discipline. The use of remote sensing data has obvious advantages over land cover classification. The accuracy of these data depends on the quality of information extraction from them, and the combination of optical remote sensing, which includes radiometric indices, with SAR remote sensing is generally accepted as a means of improving classification [89]. The argan tree is an endemic and endangered species in southwestern Morocco. It is a large tree with thorny twigs and attenuated leaves with a short petiole. Since 1998, an area of 830,000 hectares between Agadir and Essaouira has had the status of “biosphere reserve” granted by UNESCO to protect the argan tree [90]. In this context, the objective of this study was to map the spatial distribution of the argan tree in the rural community of Smimou using remote sensing data and evaluation of the contribution of the combination of optical and SAR satellite images in improving the identification of the tree species and particularly the argan tree. To do this, we started with data acquisition. 35 Sentinel-2 and 22 Sentinel-1 images were downloaded from THEIA and PEPS websites, respectively, due to the preprocessing already applied to the “ready-to-use” products they offer. In addition, ground truth information was acquired from two field campaigns conducted in 2019, and was supplemented by Google Earth data. The methodology consists of applying the SVM classification of the S2 products, the S1 products, and the combined products of the two sensors, in the following order: the NDVI time series; the time series of VV, VH, texture, and their combinations; and finally, the time series of the combinations of NDVI, VV, VH, and the rest. The quantitative evaluation of the results was carried out by the OA and kappa coefficient obtained by construction of the confusion matrices, while the qualities of the images were visually compared in terms of sharpness. A visual interpretation of the quality of the results showed that the classifications of multiband images containing S1 data remain “noisy” despite the fact that these images were filtered before classification. They contain huge spatial expanses of olive trees confused with argan trees. In order to improve the accuracy further, we used the postclassified product smoothing technique. Indeed, the results improved and the classification of the data derived from S1 reached an accuracy of 38.04% (kappa=0.29) for the best scenario (VH, VH/VV, and texture). Moreover, the synergistic use of the two Sentinel constellation products (1 and 2) led to an increase in terms of accuracy and quality of the map. In perspective, it would be necessary to test the use of S1 and S2 sensors for mapping the spatial distribution of the argan tree over a larger area, i.e., the Essaouira region, and subsequently at the national scale.



---

*Source: 6679914-2021-10-05.xml*