Next Article in Journal
Riparian Plant Evapotranspiration and Consumptive Use for Selected Areas of the Little Colorado River Watershed on the Navajo Nation
Next Article in Special Issue
Dynamic Maize Yield Predictions Using Machine Learning on Multi-Source Data
Previous Article in Journal
Monitoring Seasonal Movement Characteristics of the Landslide Based on Time-Series InSAR Technology: The Cheyiping Landslide Case Study, China
Previous Article in Special Issue
Plant Density Estimation Using UAV Imagery and Deep Learning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Spatio-Temporal Assessment of Olive Orchard Intensification in the Saïss Plain (Morocco) Using k-Means and High-Resolution Satellite Data

1
Bonn International Centre for Conflict Studies, 53121 Bonn, Germany
2
Department of Geography, University of Bergen, 5020 Bergen, Norway
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(1), 50; https://doi.org/10.3390/rs15010050
Submission received: 27 October 2022 / Revised: 7 December 2022 / Accepted: 19 December 2022 / Published: 22 December 2022
(This article belongs to the Special Issue Machine Learning Methods Applied to Optical Satellite Images)

Abstract

:
Olive orchard intensification has transformed an originally drought-resilient tree crop into a competing water user in semi-arid regions. In our study, we used remote sensing to evaluate whether intensive olive plantations have increased between 2010 and 2020, contributing to the current risk of aquifer depletion in the Saïss plain in Morocco. We developed an unsupervised approach based on the principles of hierarchical clustering and used for each year of analysis two images (5 m pixel size) from the PlanetLabs archive. We first calculated area-based accuracy metrics for 2020 with reference data, reaching a user’s accuracy of 0.95 and a producer’s accuracy of 0.89. For 2010, we verified results among different plot size ranges using available 2010 Google Earth Imagery, reaching high accuracy among the 50 largest plots (correct classification rate, CCR, of 0.94 in 2010 and 0.92 in 2020) and lower accuracies among smaller plot sizes. This study allowed us to map super-intensive olive plantations, thereby addressing an important factor in the groundwater economy of many semi-arid regions. Besides the expected increase in plantation size and the emergence of new plantations, our study revealed that some plantations were also given up, despite the political framework encouraging the opposite.

Graphical Abstract

1. Introduction

Olive trees are a drought-tolerant [1] species of high economic importance, traditionally grown in the Mediterranean Region, where 95 percent of global olive cultivation is located [2]. Over the last three decades, super-high density (SHD) planting systems have been taking over the olive sector worldwide, since their smaller varietal sizes and tree-to-tree distances allow for the complete mechanization of pruning and harvest processes, lowering production costs in the long term and increasing yields [3]. However, olive trees also suffer under increasing droughts and higher temperatures [4], and the recommended climate change adaptation would consist of increasing spacing between trees [5]. Besides contrarily reducing tree-to-tree distance, an increasing number of plantations are also adopting (often backup) irrigation measures that are strongly encouraged by governmental subsidies on drip irrigation [6]. Hereby, the limited water resources of many drought-prone areas dedicated to olive cultivation risk being more quickly exploited and in extreme cases depleted, increasing the environmental risk of desertification but also affecting the local population dependent on such water resources [7,8,9].
Over the last ten years, Moroccan olive production has more than doubled, placing the country among the five main producers worldwide [2]. This development has been strongly encouraged by a national strategy called the Green Morocco Plan (GMP) [10], in which important subsidies were allocated to efficient irrigation systems [11] and in the conversion to arboriculture [12]. Over 30 percent of olive cultivation occurs in the Fès and Meknès Region (FMR) [13]. Here, the Saïss plain, a very fertile area of economic importance, is one of the regions with the most intensive use of groundwater in Morocco [8]. It is a source that will be depleted within the coming 25 years [14], as a 100 Mm³ deficit is registered yearly [15]. A spatio-temporal assessment of SHD olive plantations can thus help evaluate the impact of olive orchard intensification on local water resources. For this, using remote sensing (RS) data and methods can be a cost and time-efficient approach for large-scale mapping.
Among RS studies on tree crop mapping, the use of different spatial resolutions of satellite imagery has had varying levels of success depending on the tree spacing and tree crown size that were mapped [16,17]. Sparsely planted orchards are mostly mapped using object-based image analysis (OBIA) on very high spatial resolution (VHR) imagery [18,19]. However, VHR data have high acquisition costs [20] and high computational power requirements [21]. For homogeneous and intensive orchards, some studies achieved promising results using 10 m spatial resolution [17,22].
The adoption of multi-temporal approaches in land use and land cover (LULC) classification studies allows for capturing crop-specific phenology patterns and helps increase classification accuracies. However, the trade-offs between spatial and temporal resolution concerning commercial high spatial resolution data remain. The high acquisition costs of such imagery often limit the feasibility of their application, especially in less-developed countries [23].
Besides spatial and temporal resolution, spectral resolution also plays a crucial role in discriminating crop types. Spectral properties of leaves can provide information about their physical and chemical composition. There are many factors that influence a leaf’s optical properties, such as the structure of the epidermis, waxes, and the cutin [24]. Evergreen leaves, for example, have lower reflectance in the infrared wavelengths than deciduous leaves [25]. Given these characteristics, different vegetation indices (VIs), which reflect the relationship between combinations of spectral bands, have been developed through the years and allow for determining the greenness of vegetation and their water content, and thus, they have different capacities for tree crop discrimination based on their leaves and planting characteristics.
The most widely used VI in crop classification studies is the Normalized Difference Vegetation Index (NDVI) [26,27]. However, it possesses limitations such as soil background brightness [28] and saturation [29,30,31]. The Modified Soil Adjustment Vegetation Index (MSAVI-2) [32] allows one to simultaneously increase the vegetation signal while decreasing the soil-induced variations, thus minimizing the influences of soil on vegetation spectra. Specifically, studies on tree crop classification, such as [33], concluded that NDVI performs worse on the discrimination of tree types than of annual crops. Furthermore, due to the high similarity amongst temporal reflectance profiles of fruit trees, recent work shifted from VI- to full-band-based multi-temporal approaches that consider the entire spectral resolution available in the image data set [34].
In the case of olive trees, their evergreen phenology presents a discriminatory advantage against deciduous trees and other land use classes with high seasonal variability. The lack of seasonal changes in olive trees’ phenology makes the use of multi-temporal data redundant. As [16] concluded, multi-temporal datasets did not improve accuracy in olive detection itself, but rather, they improved the classification accuracy of surrounding classes, reducing misclassification errors. We thus argue that, for olive orchard detection, the use of two images per year allows us to eliminate surrounding ground classes with changing phenology. This would further ease the use of higher spatial resolution with lower temporal availability or associated with high acquisition costs.
In LULC studies, besides both higher data availability and accessibility, increasing computational capacity has allowed for the development and adoption of machine learning algorithms (MLAs). MLAs can be subcategorized into supervised methods, which consist of training a classifying algorithm using reference data [35,36], and unsupervised methods. While supervised approaches are preferred due to their higher robustness and precision [37], when no labelled data is available, as in the present study, unsupervised methods are required.
Unsupervised methods recognize unclassified data without any prior knowledge of LULC types, with an interpreter assigning a class to each cluster of pixels via visual interpretation [38]. They group pixels based on the similarity of their values into clusters or spectral classes [37]. There is the possibility of applying the clustering algorithm in a hierarchical form, either divisive or agglomerative, starting with a large cluster, which is iteratively subdivided, or many clusters that are iteratively merged [35,39]. k-means [40,41] is one of the simplest and most well-known clustering algorithms and is often used in RS studies [42,43].
In this regard, the overall objective of this study was to develop a two-step classification approach using k-means and high-resolution satellite data to enable a spatio-temporal assessment of olive orchard intensification in the Saïss Plain (Morocco). Specifically, we aimed:
(i)
To map the extent of SHD olive plantations in the Saïss plain;
(ii)
To assess the development of this planting pattern from the implementation of the GMP in 2010 until 2020.
Given the limited availability of HR imagery and the lack of reference data over the past years, we developed an unsupervised approach for single-class detection, which uses the advantage of olive trees’ evergreen phenology and the semi-arid climatic conditions in the study area. We hereby aimed to prove that super-intensive olive plantations are a suitable crop for mapping with satellite imagery of low temporal resolution and unsupervised classification techniques. This can be valuable when reference data is scarce and image availability is low or acquisition costs are high.

2. Materials and Methods

2.1. Study Area

Located at 583 m above sea level between the cities of Fès and Meknès in Northern Morocco, the Saïss plain has a semi-continental Mediterranean climate, with moderate winters and hot dry summers [44]. Local rainfall is about 400 mm per year, which mainly occurs during winter. However, the proximity to the Middle Atlas Mountains, where rainfall can reach up to 800 mm [45], provides the Saïss Plain with additional river discharge that also feeds the underlying aquifers, making it a very fertile region [46,47].
As a result of the liberalization and privatization of the Moroccan agricultural sector [48] the Saïss plain’s main characteristic is structural duality caused by the coexistence of smallholder family farms together with an increasing share of large-scale investor-led farming systems [49]. Besides large-scale intensive agriculture, [50] reported that the introduction of olive orchards substituting cereal crops has pushed some farmers in the Saïss area to exploit the rows between the trees with annual crops. This creates a highly heterogeneous agricultural landscape.
Based on regional crop statistics, 65 percent of all tree crops and almost 95 percent of evergreen tree crops cultivated in Morocco are olive (1,098,000 ha), which includes the FMR, where the Saïss plain is located, the leading producer [13,51]. The second most cultivated evergreen tree crops are citrus trees, which represent about 5 percent nationally [13]. However, the main citrus cultivation sites are not in the FMR [51]. This allows us to assume that most evergreen tree crops in the study area are olive.
The Saïss aquifer has an approximate extent of 220,000 ha and provides drinking water to the main urban areas of the region with around 6 million inhabitants and is used for irrigated agriculture [47], which takes place on around 40,000 ha. The extent of the aquifer was used for study area definition, adding a buffer of 1.3 km to reduce contour irregularities, obtaining a total extent of 241,390 ha. The availability of Planet data delimitated the study area within this area to an extent of 133,377 ha (Figure 1).

2.2. Methods Overview

The method we developed consists of applying a hierarchical divisive clustering and masking technique on HR PlanetScope and RapidEye images (one for the dry and one for the rainy season) alternately. First, we preprocessed RapidEye (RE) and PlanetScope (PS) images, generating one rainy and one dry season composite for each year. Then, we applied a k-means clustering and a masking sequence alternatingly on 2020 satellite imagery to extract intensive olive plantations. We collected reference data using a self-programmed NDVI-time series inspector tool and by using bi-temporal photointerpretation in Google Earth Engine (GEE) and Google Earth, respectively. We performed area-based accuracy assessment using reference data in 2020 for 12 different LULC classes, which allowed us to optimize the clustering and masking sequence. We then applied the best performing parameter settings from 2020 tests on 2010 HR imagery to map the olive plantation extent and subsequently evaluate the main changes in intensive olive plantations between 2010 and 2020. Additionally, we sample-wise verified 2010 and 2020 results among different plot size ranges from the detected olive orchards (Figure 2).

2.3. Image Sources and Pre-Processing

Through PlanetLabs’ Education and Research Program [57], we freely acquired and used a limited quota of HR imagery of 3 and 5 m pixel size, dating back to 2009. While open-source satellite alternatives, such as 10 m resolution Sentinel-1 and -2 data, have unlimited access and more spectral bands, they only date back to 2015 and 2017, respectively. Likewise, Landsat, which is the longest-running enterprise for Earth Observation with high temporal availability, was not suitable for our research purpose due to its 30 m pixel size. Thus, we used two seasonally representative images of RapidEye (RE, 5 m pixel size) and PlanetScope (PS, 3 m) for each year of analysis (2010 and 2020, respectively). The images corresponding to the rainy season were from January 2011 and 2020, and those for the dry season were from August 2010 and 2020. There are different formats of imagery available, such as the PlanetScope Analytic Ortho Scene Product (L3B) and the RapidEye Ortho Tile Product (L3A), that are recommended for land cover mapping since they are orthorectified, scaled to Surface Reflectance, and projected to a cartographic projection [56]. Thus, we used these products herein (Table 1):
All images were merged and clipped in QGIS to the defined extent and uploaded into GEE. Further, we resampled PS data to 5 metres and normalized RE data to PS data using the following formula based on [58]:
I n e w = I M i n M a x M i n ( M a x n e w M i n n e w ) + M i n n e w
where for each band of the RE image, I   is the original pixel value; I n e w   is the normalized pixel value; M a x and M i n   are the maximum and minimum pixel values of the RE image, respectively; M a x n e w and M i n n e w are the maximum and minimum pixel values of the corresponding band of PS imagery.
While mapping tree crops usually benefits from using full-band composites instead of merely conducting VI-based classifications [33], the comparatively low spectral resolution of PS and RE images (Table 1) and the nature of the here-developed approach invite opting for a VI-based approach. Based on the limited possibilities and the success of other studies in semi-arid areas [59], we calculated MSAVI-2 and NDVI from the preprocessed imagery to compare their performance in the study area (Table 2).

2.4. Reference Data Collection and Crop Separability

For cluster validation and method evaluation and improvement, we collected reference data for 2020. The aim was to detect whether the method also captured other land use classes as commission errors. Due to the prevailing travel restrictions connected to the COVID SAR-19 pandemic while we conducted this study, we could not perform reference data collection in situ. For this reason, an NDVI time series inspector tool was developed in GEE (https://rebeccanavarrokhoury.users.earthengine.app/view/referencedatacollector, accessed on 21 October 2022) (Figure 3), which allowed for obtaining NDVI time-series information per click on a selected pixel and collecting reference data for approach validation. The satellite data used for this task were 10 m Sentinel-2 (S2) L2A imagery. We filtered S2 time series (TS) available for NDVI-profiling based on cloud cover (5 percent) and applied a cloud mask the imagery using the cloud mask band. Then, we calculated NDVI from the pre-processed S2 imagery and used it for time series plots.
The NDVI time series inspector tool allowed for categorizing annual crops into spring, winter, summer, and double cropping, based on the time of highest NDVI values (Figure 4).
However, the 10 m pixel size of Sentinel-2 imagery was not suitable for NDVI time series calculation on orchards because it was biased by soil background reflectance. For tree crops, we thus used the photointerpretation of Google Earth imagery using the time slider function to visually discriminate deciduous from evergreen trees, which, based on crop statistics [13,51], were assumed to be olive trees (Figure 5).
Collecting training data via photointerpretation revealed that different planting forms based on tree-to-tree distances are present in the study area. The olive class was, thus, subdivided into three classes: traditional, intensive, and super-intensive olive plantations (Table 3, Figure 6). We use the terms super-intensive and SHD in this work interchangeably.
Based on these criteria, we defined the 12 LULC classes and collected the corresponding reference data in the self-developed GEE app with the following listed areas (Table 4).
We plotted monthly NDVI time series and frequency distributions per land-use class computed from the Sentinel-2 data from all sites to visualize land use class separability (Figure 7).
From Figure 7, we can see that deciduous orchards (Figure 7d) had very similar temporal profiles to evergreen olive orchards (Figure 7a,b). The visual examination of multi-temporal VHR Google imagery showed that this was often caused by intercropping practices or natural grass growth during the winter months, which biased the reflectance profile (see Figure 8). Riparian vegetation also showed high NDVI values throughout the year (Figure 7f).
To further investigate crop separability, we plotted values in selected wavelengths for six crop classes (Figure 9). In the visible range (S2 bands blue, green, and red in July 2020), all samples except bare land had similar values, especially in the green range. Riparian vegetation and deciduous orchards had almost identical values in the visible range. Deciduous orchards had the highest NIR values, followed by riparian vegetation and summer crops. Olive orchards had lower NIR values, which agrees with findings from the literature [25] and further proves the high accuracy of the sampling approach based on photointerpretation. SHD olive orchards had higher NIR values than traditional orchards due to mixed pixels in the latter; however, both converged in the SWIR wavelengths. The analysis of the spectral signatures thus confirmed our assumption that infrared values allow for the discrimination of land use classes with similar reflectance values in the visible range.
We later used the collected reference data to evaluate and improve the unsupervised method applied on 2020 imagery, described in the following sections.

2.5. Unsupervised Classification

2.5.1. Cluster and Masking Sequence

The k-means clustering algorithm locates cluster centers evenly distributed in the data space and uses a distance measure in an iterative procedure to associate each data point to the nearest, statistically similar cluster center. Then, cluster centroids are moved to the average of their class values until the smallest possible distance is reached [39].
The here-developed unsupervised approach was based on the principles of hierarchical divisive clustering and only requires the use of two images per year, which are alternatingly clustered and masked as follows (Figure 10).
The developed approach consists of several consecutive steps: (i) First, a two-cluster k-means algorithm was applied to a dry season composite of the tested VI (NDVI or MSAVI-2) for 2020 on Planet Scope imagery and later for 2010 on RapidEye data. This first step allows for the differentiation of non-vegetated from vegetated areas during summer (August) which is a dry season in the study area. The mean VI of each cluster was computed and used to automatically assign the cluster with the highest VI mean to the vegetated cluster, and the lower VI mean to the unvegetated clusters. The non-vegetated areas in the study region include rainfed annual crops, built-up areas, and natural vegetation, such as shrubs. The vegetated cluster includes deciduous orchards, summer annual crops, and (mainly SHD) olive plantations. Most traditionally planted olive orchards with large spaces between trees and rows were included in the non-vegetated area because of predominant soil reflectance in the 5 m pixel of the image composite. The vegetated cluster was then used to mask the winter season composite (January). Then, (ii) a second two-cluster k-means algorithm was applied on the masked winter image, dividing it again into vegetated and non-vegetated areas for the winter season. The non-vegetated areas are now deciduous trees and annual summer crops that were harvested. The remaining vegetated areas in the final cluster are expected be olive orchards of high-density planting patterns (with no or very low soil reflectance between trees).
First attempts revealed high commission errors in this final olive cluster. These were mainly caused by undercropping practices during the winter months on deciduous orchards, which made them have high VI values also on the winter composite. Based on the findings from Section 2.4, we added a (iii) third clustering step using the NIR band of the summer imagery to re-cluster the remaining orchards contained in the second step cluster into deciduous and evergreen tree crops. Here, the cluster with the lower NIR mean was automatically assigned to the olive classes, and the higher was assigned to the remaining commission errors (deciduous orchards and riparian vegetation). This third step aimed to remove the remaining non-olive classes during the peak vegetation period of deciduous and riparian vegetation (see Figure 11).

2.5.2. Area-Based Accuracy Assessment for 2020

The collected reference data (Section 2.4) were used to assess the accuracy of the final SHD olive cluster. In GEE, we calculated the area where the cluster result overlapped different land use classes. Since the approach was a single class approach, only the user’s and producer’s accuracy for the target class (SHD olive plantations) were calculated, from which omission and commission errors relative to the 12 LULC types from the reference data were estimated.
For user’s accuracy ( U i ), the proportion of the area mapped as olive within the reference class olive ( p i i ) and the total area mapped as olive ( p i ) were used:
U i = p i i / p i
For producer’s accuracy ( P j ), we used the area proportion of the reference class olive that was mapped as olive ( p j j ) and the total area of the reference class olive ( p j ), as follows:
P j = p j j / p j
Commission errors could be attributed to the different land use classes, using
C o m m i s s i o n = p i h / p h
where p i h is the mapped olive area within a non-olive reference class, and p h is the total area of the non-olive reference class in question.
Omission errors were calculated using the producer’s accuracy complementary measure:
O m i s s i o n = 1 p j j / p j
For intercomparability, we combined both metrics within the F-score [61]:
F s c o r e = 2 × U i × P i U i + P i

2.5.3. Count-Based Binary Accuracy Assessment for 2010 and 2020

Additionally, for 2010 and 2020, we sample-wise verified results among different plot size ranges with VHR Google Earth Imagery. For this, we verified the 50 largest plots, the 50 smallest plots, and 50 plots around the median area for being or not being tree crop plantations. We calculated the correct classification rate (CCR) for each batch as the sum of the number of correctly classified plots divided by the sample size n ( n   = 50):
C C R =   i = 1 n δ ( y i ,   y ^ i ) n
where y i is the real ground class observed on VHR Google Earth Imagery, y ^ i is the vector of predicted class labels, and δ is an indicator variable such that δ ( y i ,   y ^ i ) = 1   if y i = y ^ i and zero otherwise.

2.6. Post-Processing

To reduce commission errors from small erroneous pixels or road trees included in the cluster, we “sieved” the data in QGIS (version Białowieża) using the GDAL Sieve function [62]. The sieving tool removes smaller polygons below a defined threshold and replaces them with the pixel value of the largest neighbor polygon. This allowed for the removal of a large share of single trees that did not belong to large plantations. At the same time, speckled areas within super-intensive plantations caused by plantation heterogeneity were closed and homogenized.

2.7. Change Detection Analysis

Finally, we used the results from both years as inputs for the change detection analysis. For this, the digital number (DN) value of a pixel from 2010 was subtracted from the DN value of the same pixel in 2020 [63,64], resulting in pixels with new orchards (DN = 1), no change (DN = 0), and orchards given up to other land uses (DN = −1).

3. Results

3.1. Planting Patterns in 2010 and 2020

The accuracy assessment revealed that a large share of intensive olive plantations (~65 percent of the reference data) was included in the final cluster. Main areas of super-intensive and intensive olive plantations were found in the central part of the plain and the region nearby Meknès in 2010. This distribution prevailed also in 2020, with new plantations appearing also in the East near Fès and the Southwest (see Figure 12).
The total area detected in 2010 was 6267 ha and in 2020, 6785 ha. The minor increase in area can be attributed (i) to changes from smaller plot sizes to large-scale plantations, and (ii) based on the sample-based accuracy assessment (see Section 3.3.2), also to greater commission errors in 2010, especially among smaller plot sizes.

3.2. Change Detection Analysis

Algebra-based change detection was performed as described in Section 2.7, and it revealed that, from the 6267 ha in 2010, only 1470 ha remained in 2020 (“No change”), and 4797 ha were given up since then (“Converted to other land uses”). There were 5315 ha of “new” olive plantations in 2020, mainly in the form of large-scale plantations (Figure 13).
VHR Google Earth imagery was also used to sample-wise observe which land use classes had occurred. The land-use change from SHD olive to other classes (here labelled as “losses”) was due to conversions to annual crops or bare land, which, as later fieldwork revealed, was often the previous land preparation for urbanization purposes (Figure 14). Most of the areas converted to olive plantations (here labelled as “new orchards”), were originally mostly annual crops (Figure 15a).

3.3. Accuracy Assessment

3.3.1. Area-Based Accuracy Assessment for 2020: NDVI vs. MSAVI-2

The differences between using NDVI and MSAVI-2 were subtle, with the former reaching higher UA than the latter (see Table 5) due to MSAVI-2′s capacity to reduce soil reflectance, which caused higher commission errors among traditional and intensive olive plantations (see Table 6). Both indices reached an F-Score of 0.78 when including a third step with the NIR band.
The last clustering step that used the NIR band considerably reduced errors of commission from riparian vegetation (by 0.85 and 0.82) and deciduous orchards (by 0.49 and 0.43) for MSAVI-2 and NDVI, respectively. Commission errors from traditional orchards were around 30 percent for MSAVI-2 and 20 percent for NDVI, which were due to the large tree crown size typical among such planting patterns and which were captured by the 5 m pixel size. Additionally, the effect of MSAVI-2 on soil background reflectance increased commission errors. Around 70 percent of intensive olive plantations were mapped by our approach as well, which was due to the larger tree crown sizes and the narrower rows compared to traditional orchards.
Super-intensive olive plantations were less affected by soil background reflectance due to their narrow planting patterns, thus using NDVI or MSAVI-2 would not considerably change the PA. For mapping SHD olive plantations, using NDVI is better than MSAVI-2. Additionally, since 70 percent of the reference intensive olive class were captured by our approach, we merged both classes. Thus, the approach we developed was rather suitable for intensive olive mapping in general. Once merged, an F-Score of 0.89 was reached by MSAVI-2 and 0.87 by NDVI (see Table 7).
To further improve accuracy, we removed erroneous pixels using the GDAL Sieve function in QGIS (version Białowieża), for which we tested different thresholds from 10 to 1000 and found that the best tradeoff was reached using a threshold of 500 pixels. Then, we vectorized the sieved layer and calculated per-plot-size statistics. This allowed us to remove a large share of single trees from traditional olive plantations. At the same time, speckled areas within super-intensive plantations caused by plantation heterogeneity were closed and homogenized (Figure 16), further improving area-based accuracy (see Table 8 and Table 9).
The best results were reached when merging intensive and super-intensive olive classes and using the MSAVI + NIR set. The difference in the NDVI+NIR set was due to commission errors in the land use classes intensive and traditional olive plantations, which were 10 percent higher in the MSAVI than in the NDVI set. We considered it better to have a lower number of traditional olive plantations in the final cluster than to have a higher number of intensive plantations included, and thus we decided that the NDVI set was most representative for our purposes.

3.3.2. Sample-Wise Verification and Comparison of 2010 and 2020

We applied the NDVI-based sequence on 2010 normalized RE imagery and assigned clusters automatically, as we also did for 2020, and verified them using VHR Google Earth imagery available for the year of analysis. After sieving and vectorizing, we sample-wise verified plots among different plot sizes, which included the 50 largest, the 50 smallest, and 50 around the median size. Available Google Earth Imagery for 2010 and 2020 revealed that all large plantations mapped were correct in both years, reaching a high number of true positives and a correct classification rate (CCR) of 0.96 and 0.92 in 2010 and 2020, respectively (Table 10). In the smaller plot size range, results from 2010 performed worse than those for 2020 (reaching a CCR of 0.4 and 0.76, respectively), and most commission errors took place on smaller fields with annual crops. Intensive irrigated annual cropping was among the main classification errors.

4. Discussion

4.1. Methodological Achievements and Shortcomings

In our work, super-intensive and intensive olive plantations were successfully mapped in a highly heterogeneous study area in Morocco. For this, we implemented an unsupervised classification approach based on the principles of hierarchical clustering, which allowed us to detect SHD and intensive olive plantations on 5 m resolution imagery. Hereby, two frequent cost-related problems in remote sensing were addressed: (i) while in most land use analyses the cost- and time-intensive gathering of reference data is necessary, e.g., to train supervised classification algorithms, in our approach no labelled data were required; (ii) the design of the clustering and masking sequence and the phenological characteristics of the target crop allowed us to use a very low temporal resolution, requiring only two images per year of analysis. This can be especially advantageous to reducing acquisition costs when using commercial imagery of higher spatial resolution.
Compared to supervised approaches, which can map intra-class variability, our method could only detect dense and vigorous olive plantations, including mainly super-intensive, but also a large share of intensive, and some traditional olive plantations. Younger plantations with smaller tree crown sizes could not be detected or led to unmapped areas within the same plantation. However, this is a problem that other studies using supervised classification for tree crop mapping with high-resolution data of 3 m pixel size, such as [16], also faced during their analyses. On the other hand, high commission errors were scored on intensive olive plantations due to their similar reflectance values, something that supervised classification may overcome by training the classifying algorithm on different planting intensities as separate land use classes.
Previous work on tree crop mapping, specifically on olive trees, found that using multi-temporal information would not improve olive detection accuracy itself, but rather improve the accuracy of other crop classes, and hereby overall accuracy [16]. Additionally, the difficulty of mapping perennial crops merely based on their phenological profiles was already discussed in [17] and [34], who found that using all reflectance bands helped discriminate different tree crop types. Our study confirmed these findings, since using the low temporal resolution of only two images per year and making use of the NIR band allowed us to reach a remarkably high accuracy in detecting olive plantations. However, limitations of NIR wavelengths were found in the discrimination between some annual crops and olive plantations, which confirms previous findings from the crop separability analysis conducted in Section 2.4 (Figure 9).
Based on previous studies on semi-arid environments, we tested different indices (NDVI vs. MSAVI-2) to see whether they helped improve mapping on sites with higher soil background reflectance. While MSAVI-2 successfully reduced soil background reflectance, to map super-intensive olive plantations, which are less affected by this phenomenon, we found that NDVI performed better, and using MSAVI-2 led to higher commission errors among traditional olive orchards.
Our sample-wise verification of 2010 results using VHR Google Earth Imagery revealed that a higher number of smaller plots of annual crops were misclassified in 2010 than in 2020, and that main commission errors in 2020 results were deciduous tree plantations. This may confirm the expected shift from cereal-based agriculture in 2010 towards arboriculture in 2020, incentivized by the GMP. Hereby, k-means cluster centers may also have shifted, creating such a difference among commission errors, both qualitatively and quantitatively. While applying the GDAL sieving tool removed most of the smaller misclassified spots in both years and improved accuracy, further research may investigate whether other unsupervised clustering methods are more appropriate for our proposed approach.
As opposed to most studies on olive mapping, which were mostly performed on homogeneous environments with a low diversity of surrounding ground classes [19], our study succeeded in extracting olive plantations, mostly of super-intensive and to a great extent also of intensive planting patterns, on a very heterogeneous ground. However, it should be noted that no other evergreen tree crops are known to be cultivated in the study area. Thus, for replicability, this approach needs to be adapted if other evergreen tree crops are also present.

4.2. Implications for Agricultural Management and Policy Monitoring in the Study Region

This study addressed a subject that is currently affecting the olive sector not only in Morocco but also globally and will be of growing concern in times of climate change. Although irrigated olive cultivation can be an important carbon sink and can fight erosion while improving the profitability of the sector, the limited water resources of many drought-prone areas dedicated to olive cultivation may be also exploited and depleted more quickly, presenting not only an environmental risk but also affecting the local population who depend on such water resources.
Mapping super-intensive olive plantations can help locate important actors in the olive oil sector. However, although our study successfully mapped the extent and evolution of an important actor of the groundwater economy in the Saïss plain, we found that super-intensive olive plantations only represent a share of different land uses in the study area and that other actors also need to be addressed to improve water resource management.
Finally, our approach was helpful to analyze policy-induced land use change in the study area. It revealed that the average plantation size increased between 2010 and 2020, reflecting large agribusinesses taking over the sector to the detriment of smaller family farms and medium-size investors. Another remarkable finding was the detection of several larger SHD olive farms that were given up, despite the encouragement of the GMP to adopt intensive forms of arboriculture, especially olive trees. The reasons behind this can be attributed to the short lifecycle of SHD olive plantations. Other contributing factors, such as access to water or better revenues from other crops, should also be considered and could be further investigated.

5. Conclusions

Detecting super-intensive olive plantations may help to quantify the potential risks to local water resource management to address stakeholders and implement climate-smart agriculture technologies when needed.
While the use of remote sensing data and methods for mapping olive and tree crops has long been challenged by high spatial resolution requirements, we found that mapping super-intensive olive plantations comes with advantages compared to other tree crops and planting patterns. Using 5 m resolution allowed us to map all super-intensive olive plantations in the study area, except for areas with younger trees with smaller tree crown sizes.
Furthermore, the olive’s evergreen nature allowed us to use a reduced number of images per year. Especially when the aim was to extract super-intensive olive plantations within a single-class classification approach, as was the case in our study, we concluded that two satellite images were enough to suppress all other surrounding land use classes with higher variability in their phenology. This was made possible by using an unsupervised approach, which also allowed us to overcome the challenge of gathering labelled data, usually necessary for supervised classification methods.
The approach we developed, based on the principles of hierarchical clustering, consisted of applying a two-cluster k-means algorithm alternatingly on one image from the dry season and one from the winter season, extracting the vegetated cluster and creating a mask that was applied on the subsequent image to cluster. Using the NIR band was crucial to refine results and remove commission errors such as deciduous orchards with undercropping practices during the winter season, confirming findings from previous studies on tree crop mapping. In addition, we also concluded that NDVI was better to map SHD olive plantations than MSAVI-2, since they are less affected by soil background reflectance than other orchards of lower tree density.
Finally, this study leaves three open questions that may be addressed by future research: (i) while the performance of the approach we developed was quite promising, sample-wise verification of 2010 results revealed a larger share of commission errors among annual crops, suggesting that k-means may not be the best clustering algorithm for the proposed method, which further research may investigate; (ii) our study revealed that water-intensive SHD olive plantations, among other forms of land use, are a major player in the groundwater economy of the Saïss plain; however, further research is required to investigate which other land use classes are putting the Saïss aquifer at risk of depletion; (iii) lastly, our results showed that, despite the considerable commission errors among the small-sized plots in 2010, there was also a considerable number of larger plantations that appeared to have been given up and converted to other land uses. Future research may investigate the reasons for some areas converting from arboriculture to annual crops, despite the political framework encouraging the opposite.
We further would like to encourage interested researchers to test our approach in other areas with Sentinel-2 data using our GEE Olive Detection App for super-intensive olive plantation mapping: https://rebeccanavarrokhoury.users.earthengine.app/view/olivedetectionapp (accessed on 21 October 2022).

Author Contributions

Conceptualization, R.N. and L.W.; methodology, R.N. and O.D.; software, R.N.; validation, R.N., O.D. and L.W.; formal analysis, R.N.; investigation, R.N.; resources, R.N.; data curation, R.N.; writing—original draft preparation, R.N.; writing—review and editing, O.D.; visualization, R.N.; supervision, O.D. and L.W.; project administration, L.W.; funding acquisition, L.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Federal Ministry for Education and Research, grant number 01LZ1807D.

Data Availability Statement

The Google Earth Engine Code can be found online on https://github.com/rbcnavarro/olive_intensification.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Sofo, A.; Manfreda, S.; Fiorentino, M.; Dichio, B.; Xiloyannis, C. The olive tree: A paradigm for drought tolerance in Mediterranean climates. Hydrol. Earth Syst. Sci. 2008, 12, 293–301. [Google Scholar] [CrossRef] [Green Version]
  2. FAOSTAT. FAOSTAT, 2021. Available online: http://www.fao.org/faostat/en/#data (accessed on 10 September 2021).
  3. Fernández-Escobar, R.; Gil-Ribes, J.A.; Quesada-Moraga, E.; Trapero, A.; Msallem, M. Evolution and sustainability of the olive production systems. Options Mediterr. 2013, 106, 11–42. [Google Scholar]
  4. Brito, C.; Dinis, L.-T.; Moutinho-Pereira, J.; Correia, C.M. Drought Stress Effects and Olive Tree Acclimation under a Changing Climate. Plants 2019, 8, 232. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Boulet, G.; Jarlan, L.; Olioso, A.; Nieto, H. Chapter 2—Evapotranspiration in the Mediterranean region. In Water Resources in the Mediterranean Region; Zribi, M., Brocca, L., Tramblay, Y., Molle, F., Eds.; Elsevier: Amsterdam, The Netherlands, 2020; pp. 23–49. ISBN 978-0-12-818086-0. [Google Scholar]
  6. Rodríguez Sousa, A.A.; Barandica, J.M.; Rescia, A. Ecological and Economic Sustainability in Olive Groves with Different Irrigation Management and Levels of Erosion: A Case Study. Sustainability 2019, 11, 4681. [Google Scholar] [CrossRef] [Green Version]
  7. Martínez-Valderrama, J.; Guirado, E.; Maestre, F. Unraveling Misunderstandings about Desertification: The Paradoxical Case of the Tabernas-Sorbas Basin in Southeast Spain. Land 2020, 9, 269. [Google Scholar] [CrossRef]
  8. Molle, F.; Tanouti, O.; Faysse, N. Morocco. In Irrigation in the Mediterranean; Molle, F., Sanchis-Ibor, C., Avellà-Reus, L., Eds.; Springer International Publishing: New York, NY, USA, 2019; pp. 51–88. ISBN 978-3-030-03696-6 978-3-030-03698-0. [Google Scholar]
  9. Molle, F.; Closas, A. Why is state-centered groundwater governance largely ineffective? A review. WIREs Water 2020, 7, e1395. [Google Scholar] [CrossRef]
  10. ADA. Foundations. Available online: https://www.ada.gov.ma/en/foundations (accessed on 25 November 2021).
  11. ADA. Approches de Mise en Œuvre des Deux Piliers du PMV. Available online: https://www.ada.gov.ma/fr/approches-de-mise-en-oeuvre-des-deux-piliers-du-pmv (accessed on 25 November 2021).
  12. ADA. Main Achievements of the Green Morocco Plan. Available online: https://www.ada.gov.ma/en/main-achievements-green-morocco-plan (accessed on 25 November 2021).
  13. Fellah Trade. Les Chiffres Clés de la Filière Oléiculture—Fellah Trade, 2021. Available online: https://www.fellah-trade.com/fr/filiere-vegetale/chiffres-cles-oleiculture?filiere=filiere_vegetale (accessed on 25 November 2021).
  14. European Bank for Reconstruction; Development EBRD. Saiss and Garet Water Conservation Project, 2020. Available online: https://www.ebrd.com/what-we-do/project-information/board-documents/1395289064192/Saiss_and_Garet_Water_Conservation_Project_(SSBR).pdf?blobnocache=true (accessed on 21 October 2022).
  15. Agence du Bassin Hydraulique du Sebou ABHS. Système Aquifère du Saiss, 2016. ABHS. Available online: http://www.abhsebou.ma/presentation-du-bassin/eaux-souterraines/systeme-aquifere-du-saiss/ (accessed on 4 December 2021).
  16. Lin, C.; Jin, Z.; Mulla, D.; Ghosh, R.; Guan, K.; Kumar, V.; Cai, Y. Toward Large-Scale Mapping of Tree Crops with High-Resolution Satellite Imagery and Deep Learning Algorithms: A Case Study of Olive Orchards in Morocco. Remote Sens. 2021, 13, 1740. [Google Scholar] [CrossRef]
  17. Brinkhoff, J.; Vardanega, J.; Robson, A.J. Land Cover Classification of Nine Perennial Crops Using Sentinel-1 and -2 Data. Remote Sens. 2019, 12, 96. [Google Scholar] [CrossRef] [Green Version]
  18. Khan, A.; Khan, U.; Waleed, M.; Khan, A.; Kamal, T.; Marwat, S.N.K.; Maqsood, M.; Aadil, F. Remote Sensing: An Automated Methodology for Olive Tree Detection and Counting in Satellite Images. IEEE Access 2018, 6, 77816–77828. [Google Scholar] [CrossRef]
  19. Waleed, M.; Um, T.-W.; Khan, A.; Khan, U. Automatic Detection System of Olive Trees Using Improved K-Means Algorithm. Remote Sens. 2020, 12, 760. [Google Scholar] [CrossRef] [Green Version]
  20. Apollo Mapping Price List. Available online: https://apollomapping.com/image_downloads/Apollo_Mapping_Imagery_Price_List.pdf (accessed on 22 December 2021).
  21. Díaz-Varela, R.A.; La Rosa, R.D.; León, L.; Zarco-Tejada, P.J. High-Resolution Airborne UAV Imagery to Assess Olive Tree Crown Parameters Using 3D Photo Reconstruction: Application in Breeding Trials. Remote Sens. 2015, 7, 4213–4232. [Google Scholar] [CrossRef]
  22. Htitiou, A.; Boudhar, A.; Lebrini, Y.; Hadria, R.; Lionboui, H.; Benabdelouahab, T. A comparative analysis of different phenological information retrieved from Sentinel-2 time series images to improve crop classification: A machine learning approach. Geocarto Int. 2020, 37, 1426–1449. [Google Scholar] [CrossRef]
  23. Dubovyk, O.; Menz, G.; Conrad, C.; Thonfeld, F.; Khamzina, A. Object-based identification of vegetation cover decline in irrigated agro-ecosystems in Uzbekistan. Quat. Int. 2013, 311, 163–174. [Google Scholar] [CrossRef]
  24. Ehleringer, J.; Björkman, O.; Mooney, H.A. Leaf Pubescence: Effects on Absorptance and Photosynthesis in a Desert Shrub. Science 1976, 192, 376–377. [Google Scholar] [CrossRef] [PubMed]
  25. Ustin, S.L.; Jacquemoud, S. How the Optical Properties of Leaves Modify the Absorption and Scattering of Energy and Enhance Leaf Functionality. In Remote Sensing of Plant Biodiversity; Cavender-Bares, J., Gamon, J.A., Townsend, P.A., Eds.; Springer International Publishing: New York, NY, USA, 2020; pp. 349–384. ISBN 978-3-030-33157-3. [Google Scholar]
  26. Huang, S.; Tang, L.; Hupy, J.P.; Wang, Y.; Shao, G. A commentary review on the use of normalized difference vegetation index (NDVI) in the era of popular remote sensing. J. For. Res. 2021, 32, 1–6. [Google Scholar] [CrossRef]
  27. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring vegetation systems in the Great Plains with ERTS. NASA Spec. Publ. 1974, 351, 309–317. [Google Scholar]
  28. Elvidge, C.D.; Lyon, R.J. Influence of rock-soil spectral variation on the assessment of green biomass. Remote Sens. Environ. 1985, 17, 265–279. [Google Scholar] [CrossRef]
  29. Carlson, T.N.; Ripley, D.A. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sens. Environ. 1997, 62, 241–252. [Google Scholar] [CrossRef]
  30. Turner, D.P.; Cohen, W.B.; Kennedy, R.E.; Fassnacht, K.S.; Briggs, J.M. Relationships between Leaf Area Index and Landsat TM Spectral Vegetation Indices across Three Temperate Zone Sites. Remote Sens. Environ. 1999, 70, 52–68. [Google Scholar] [CrossRef]
  31. Stroppiana, D.; Boschetti, M.; Brivio, P.A.; Nizzetto, L.; Di Guardo, A. Forest leaf area index in an Alpine valley from medium resolution satellite imagery and in situ data. J. Appl. Remote Sens. 2012, 6, 1–17. [Google Scholar] [CrossRef]
  32. Qi, J.; Kerr, Y.; Chehbouni, A. External factor consideration in vegetation index development. In Proceedings of the 6th International Symposium on Physical Measurements and Signatures in Remote Sensing, Val d’Isere, France, 17–24 January 1994. [Google Scholar]
  33. Peña, M.A.; Brenning, A. Assessing fruit-tree crop classification from Landsat-8 time series for the Maipo Valley, Chile. Remote Sens. Environ. 2015, 171, 234–244. [Google Scholar] [CrossRef]
  34. Peña, M.A.; Liao, R.; Brenning, A. Using spectrotemporal indices to improve the fruit-tree crop classification accuracy. ISPRS J. Photogramm. Remote Sens. 2017, 128, 158–169. [Google Scholar] [CrossRef]
  35. Richards, J.A. Supervised Classification Techniques. In Remote Sensing Digital Image Analysis; Springer: Berlin/Heidelberg, Germany, 2013; pp. 247–318. ISBN 978-3-642-30061-5 978-3-642-30062-2. [Google Scholar]
  36. Wang, S.; Azzari, G.; Lobell, D.B. Crop type mapping without field-level labels: Random forest transfer and unsupervised clustering techniques. Remote Sens. Environ. 2019, 222, 303–317. [Google Scholar] [CrossRef]
  37. Richards, J.A. Clustering and Unsupervised Classification. In Remote Sensing Digital Image Analysis; Springer: Berlin/Heidelberg, Germany, 2013; pp. 319–341. ISBN 978-3-642-30061-5 978-3-642-30062-2. [Google Scholar]
  38. Mather, P.; Tso, B. Classification Methods for Remotely Sensed Data; CRC Press: Boca Raton, FL, USA, 2016; ISBN 978-1-4200-9074-1. [Google Scholar]
  39. Data Mining: Practical Machine Learning Tools and Techniques; Witten, I.H. (Ed.) Elsevier: Amsterdam, The Netherlands, 2017; ISBN 978-0-12-804291-5. [Google Scholar]
  40. Forgy, E.W. Cluster analysis of multivariate data: Efficiency versus interpretability of classifications. Biometrics 1965, 21, 768–769. [Google Scholar]
  41. Macqueen, J. Some Methods for Classification and Analysis of Multivariate Observations. In Multivariate Observations; John Wiley & Sons: Hoboken, NJ, USA, 1967; p. 17. [Google Scholar]
  42. Gulácsi, A.; Kovács, F. Sentinel-1-Imagery-Based High-Resolution Water Cover Detection on Wetlands, Aided by Google Earth Engine. Remote Sens. 2020, 12, 1614. [Google Scholar] [CrossRef]
  43. Li, C.; Wang, J.; Wang, L.; Hu, L.; Gong, P. Comparison of Classification Algorithms and Training Sample Sizes in Urban Land Classification with Landsat Thematic Mapper Imagery. Remote Sens. 2014, 6, 964–983. [Google Scholar] [CrossRef] [Green Version]
  44. Debolini, M.; Valette, E.; François, M.; Chéry, J.-P. Mapping land use competition in the rural–urban fringe and future perspectives on land policies: A case study of Meknès (Morocco). Land Use Policy 2015, 47, 373–381. [Google Scholar] [CrossRef]
  45. Funk, C.; Peterson, P.; Landsfeld, M.; Pedreros, D.; Verdin, J.; Shukla, S.; Husak, G.; Rowland, J.; Harrison, L.; Hoell, A.; et al. The climate hazards infrared precipitation with stations—A new environmental record for monitoring extremes. Sci. Data 2015, 2, 150066. [Google Scholar] [CrossRef] [Green Version]
  46. Dauteuil, O.; Moreau, F.; Qarqori, K. Structural pattern of the Saïss basin and Tabular Middle Atlas in northern Morocco: Hydrological implications. J. Afr. Earth Sci. 2016, 119, 150–159. [Google Scholar] [CrossRef] [Green Version]
  47. Quarouch, H.; Kuper, M.; Abdellaoui, E.H.; Bouarfa, S. Subterranean waters, a source of dignity as well as a social resource: The case of farmers on the Saïss plain of Morocco. Cah. Agric. 2014, 23, 158–165. [Google Scholar] [CrossRef]
  48. Olivié, I.; Pérez, A. The Difficult Escape from Dualism: The Green Morocco Plan at a Crossroads. New Medit. 2018, XVII, 37–50. [Google Scholar] [CrossRef]
  49. Ameur, F.; Kuper, M.; Lejars, C.; Dugué, P. Prosper, survive or exit: Contrasted fortunes of farmers in the groundwater economy in the Saiss plain (Morocco). Agric. Water Manag. 2017, 191, 207–217. [Google Scholar] [CrossRef]
  50. Daoui, K.; Fatemi, Z.E.A. Agroforestry Systems in Morocco: The Case of Olive Tree and Annual Crops Association in Saïs Region. In Science, Policy and Politics of Modern Agricultural System; Behnassi, M., Shahid, S.A., Mintz-Habib, N., Eds.; Springer: Dordrecht, The Netherlands, 2014; pp. 281–289. ISBN 978-94-007-7956-3 978-94-007-7957-0. [Google Scholar]
  51. Fellah Trade. Les Chiffres Clés de la Filière Rosacées Fruitières—Fellah Trade, 2021. Available online: https://www.fellah-trade.com/fr/filiere-vegetale/chiffres-cles-rosacees-fruitieres?filiere=filiere_vegetale (accessed on 19 August 2021).
  52. OCHA. Morocco—Subnational Administrative Boundaries—Humanitarian Data Exchange, 2021. Available online: https://data.humdata.org/dataset/morocco-administrative-boundaries-populated-places (accessed on 3 December 2021).
  53. ArcGIS. ArcGIS—Ressources en Eau du Maroc, 2020. Available online: https://www.arcgis.com/home/webmap/viewer.html?webmap=af42905a1f74421c9d45a9607d6a2815 (accessed on 8 December 2021).
  54. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The Shuttle Radar Topography Mission. Rev. Geophys. 2007, 45, RG2004. [Google Scholar] [CrossRef]
  55. Geofabrik GmbH. Geofabrik Download Server, 2021. Available online: https://download.geofabrik.de/africa/morocco.html (accessed on 3 December 2021).
  56. Planet Labs Inc. Planet Imagery. Product Specifications, 2018. Available online: https://assets.planet.com/docs/Planet_Combined_Imagery_Product_Specs_letter_screen.pdf (accessed on 21 October 2022).
  57. Planet Labs Inc. Education and Research Program, 2021. Available online: https://www.planet.com/markets/education-and-research/ (accessed on 26 November 2021).
  58. Qian, Y.; Zhou, W.; Yu, W.; Han, L.; Li, W.; Zhao, W. Integrating Backdating and Transfer Learning in an Object-Based Framework for High Resolution Image Classification and Change Analysis. Remote Sens. 2020, 12, 4094. [Google Scholar] [CrossRef]
  59. Ge, G.; Shi, Z.; Zhu, Y.; Yang, X.; Hao, Y. Land use/cover classification in an arid desert-oasis mosaic landscape of China using remote sensed imagery: Performance assessment of four machine learning algorithms. Glob. Ecol. Conserv. 2020, 22, e00971. [Google Scholar] [CrossRef]
  60. Google Developers. Landsat8 Harmonic Modeling—Earth Engine Code Editor, 2021. Available online: https://code.earthengine.google.com/07fe82cb6d4ef2f13f85a1b278b98897 (accessed on 6 December 2021).
  61. Xiong, J.; Thenkabail, P.; Tilton, J.; Gumma, M.; Teluguntla, P.; Oliphant, A.; Congalton, R.; Yadav, K.; Gorelick, N. Nominal 30-m Cropland Extent Map of Continental Africa by Integrating Pixel-Based and Object-Based Algorithms Using Sentinel-2 and Landsat-8 Data on Google Earth Engine. Remote Sens. 2017, 9, 1065. [Google Scholar] [CrossRef] [Green Version]
  62. Sieve. Available online: https://docs.qgis.org/2.8/en/docs/user_manual/processing_algs/gdalogr/gdal_analysis/sieve.html (accessed on 8 July 2022).
  63. Afify, H.A. Evaluation of change detection techniques for monitoring land-cover changes: A case study in new Burg El-Arab area. Alex. Eng. J. 2011, 50, 187–195. [Google Scholar] [CrossRef] [Green Version]
  64. Ke, L.; Lin, Y.; Zeng, Z.; Zhang, L.; Meng, L. Adaptive Change Detection with Significance Test. IEEE Access 2018, 6, 27442–27450. [Google Scholar] [CrossRef]
Figure 1. Overview of the study area: (a) administrative boundaries of Morocco and location of the FMR and the Saïss aquifer; (b) elevation and shaded topography (illumination coming from west at 45° altitude) in the Saïss area; (c) false color image ([R, G, B] = [NIR, R, G]) from a PlanetScope mosaic from August 2020, showing the extent of the study area defined by the intersection of available imagery and a 1300 m buffer around the Saïss aquifer. Sources: [52,53] for (a), [54,55] for (b), and [53,55,56] and Copernicus Sentinel (2020) data for (c).
Figure 1. Overview of the study area: (a) administrative boundaries of Morocco and location of the FMR and the Saïss aquifer; (b) elevation and shaded topography (illumination coming from west at 45° altitude) in the Saïss area; (c) false color image ([R, G, B] = [NIR, R, G]) from a PlanetScope mosaic from August 2020, showing the extent of the study area defined by the intersection of available imagery and a 1300 m buffer around the Saïss aquifer. Sources: [52,53] for (a), [54,55] for (b), and [53,55,56] and Copernicus Sentinel (2020) data for (c).
Remotesensing 15 00050 g001
Figure 2. Workflow overview.
Figure 2. Workflow overview.
Remotesensing 15 00050 g002
Figure 3. Reference data collection tool using NDVI time series inspector developed in GEE with S2 NDVI TS and harmonic fit based on GEE’s implementation example [60]. Per click time-series information allowed for categorizing annual crops into ”winter”, ”spring”, ”summer” and “double cropping”. Based on Copernicus Sentinel data (2020).
Figure 3. Reference data collection tool using NDVI time series inspector developed in GEE with S2 NDVI TS and harmonic fit based on GEE’s implementation example [60]. Per click time-series information allowed for categorizing annual crops into ”winter”, ”spring”, ”summer” and “double cropping”. Based on Copernicus Sentinel data (2020).
Remotesensing 15 00050 g003
Figure 4. NDVI time series (blue line and dots) and harmonic fit (red line and dots) of defined annual crop categories for: spring crops (a), double cropping (b), summer crops (c), and winter crops (d).
Figure 4. NDVI time series (blue line and dots) and harmonic fit (red line and dots) of defined annual crop categories for: spring crops (a), double cropping (b), summer crops (c), and winter crops (d).
Remotesensing 15 00050 g004
Figure 5. Bi-temporal photointerpretation of deciduous orchards. (a) Google Earth Image from February and (b) Google Earth Image from October, showing the visible difference in tree crown size and greening.
Figure 5. Bi-temporal photointerpretation of deciduous orchards. (a) Google Earth Image from February and (b) Google Earth Image from October, showing the visible difference in tree crown size and greening.
Remotesensing 15 00050 g005
Figure 6. VHR images showing varying tree-to-tree and row-to-row distances within different orchard management types: (a) traditional, (b) intensive, (c) and super-intensive or SHD. Source: Google Earth Images 2022.
Figure 6. VHR images showing varying tree-to-tree and row-to-row distances within different orchard management types: (a) traditional, (b) intensive, (c) and super-intensive or SHD. Source: Google Earth Images 2022.
Remotesensing 15 00050 g006
Figure 7. Monthly NDVI 1st, 2nd, and 3rd quartiles from collected training data of perennial and non-perennial land use classes. Based on Copernicus Sentinel-2 data (2020).
Figure 7. Monthly NDVI 1st, 2nd, and 3rd quartiles from collected training data of perennial and non-perennial land use classes. Based on Copernicus Sentinel-2 data (2020).
Remotesensing 15 00050 g007aRemotesensing 15 00050 g007b
Figure 8. Google Earth Imagery showing a deciduous orchard with winter soil greening during the winter (a) and with full crown size for comparison during the summer (b).
Figure 8. Google Earth Imagery showing a deciduous orchard with winter soil greening during the winter (a) and with full crown size for comparison during the summer (b).
Remotesensing 15 00050 g008
Figure 9. Spectral signatures of selected crops collected from the image interpretation of NDVI Sentinel-2 classes reveal larger differences in NIR and SWIR wavelengths than in the visible range during July.
Figure 9. Spectral signatures of selected crops collected from the image interpretation of NDVI Sentinel-2 classes reveal larger differences in NIR and SWIR wavelengths than in the visible range during July.
Remotesensing 15 00050 g009
Figure 10. Hierarchical divisive clustering and masking steps used to extract SHD olive plantations using a two-cluster k-means algorithm on two HR seasonally representative processed satellite images.
Figure 10. Hierarchical divisive clustering and masking steps used to extract SHD olive plantations using a two-cluster k-means algorithm on two HR seasonally representative processed satellite images.
Remotesensing 15 00050 g010
Figure 11. The three-step clustering and masking sequence that reduced commission errors with the dry season NIR band in a third masking step.
Figure 11. The three-step clustering and masking sequence that reduced commission errors with the dry season NIR band in a third masking step.
Remotesensing 15 00050 g011
Figure 12. Spatial distribution of the detected olive cluster in the Saïss plain in 2010 (a) and 2020 (b).
Figure 12. Spatial distribution of the detected olive cluster in the Saïss plain in 2010 (a) and 2020 (b).
Remotesensing 15 00050 g012
Figure 13. Land-use change of detected olive plantations between 2010 and 2020.
Figure 13. Land-use change of detected olive plantations between 2010 and 2020.
Remotesensing 15 00050 g013
Figure 14. Agricultural land use changes from intensive olive plantations to land preparation for urbanization purposes: (a) 2010, (b) 2020. Source: Google Earth Images 2020.
Figure 14. Agricultural land use changes from intensive olive plantations to land preparation for urbanization purposes: (a) 2010, (b) 2020. Source: Google Earth Images 2020.
Remotesensing 15 00050 g014
Figure 15. Land use changes from annual crops in 2010 (a) to SHD olive plantations in 2020 (b), and from SHD olive plantations in 2020 (c) to annual crops in 2011 (d). Source: Google Earth Images 2020.
Figure 15. Land use changes from annual crops in 2010 (a) to SHD olive plantations in 2020 (b), and from SHD olive plantations in 2020 (c) to annual crops in 2011 (d). Source: Google Earth Images 2020.
Remotesensing 15 00050 g015
Figure 16. Post-processing the clustered area using the GDAL Sieve tool [62] in QGIS (version Białowieża) with a threshold of 500 px: (a) the raw result from the clustering sequence (NDVI + NIR); (b) the sieved data.
Figure 16. Post-processing the clustered area using the GDAL Sieve tool [62] in QGIS (version Białowieża) with a threshold of 500 px: (a) the raw result from the clustering sequence (NDVI + NIR); (b) the sieved data.
Remotesensing 15 00050 g016
Table 1. PlanetScope and RapidEye properties and data acquired.
Table 1. PlanetScope and RapidEye properties and data acquired.
DescriptionPlanetScope (L3B)RapidEye (L3A)
Covering period2015—today2009–2020
Number of satellites120 approx.5
Image capture capacity150 million km²/day>6 million km²/day
Blue455–515 nm440–510 nm
Green500–590 nm520–590 nm
Red590–670 nm630–685 nm
Red edge--690–730 nm
Near infrared (NIR)780–860 nm760–850 nm
Pixel size3 m5 m
Acquired data5 January 2020
29 August 2020
3 August 2010
15 January 2011
Table 2. Vegetation indices analyzed in this study.
Table 2. Vegetation indices analyzed in this study.
Index AbbreviationNameFormula
NDVI [27]Normalized difference Vegetation index N I R R e d N I R + R e d
MSAVI-2 [32]Modified soil Adjustment vegetation Index 2 ( 2 N I R + 1 ( 2 N I R + 1 ) 2 8 ( N I R R e d ) 2
Table 3. Definition of subcategories for olive plantations based on tree-to-tree and row-to-row distances measured on VHR Google Earth Imagery.
Table 3. Definition of subcategories for olive plantations based on tree-to-tree and row-to-row distances measured on VHR Google Earth Imagery.
Planting PatternTree-to-Tree Distance (m)Row-to-Row Distance (m)
Traditional1010
Intensive3–66
Super-intensive1–24
Table 4. Land use class definitions based on the highest NDVI month (annual crops) and bi-temporal photointerpretation (orchards) and per class areas sampled.
Table 4. Land use class definitions based on the highest NDVI month (annual crops) and bi-temporal photointerpretation (orchards) and per class areas sampled.
Land Use ClassArea of Collected Reference Data (ha)
Evergreen orchards (traditional)103
Evergreen orchards (intensive)293
Evergreen orchards (super-intensive)628
Deciduous orchards (all types)525
Greenhouses87
Winter crops159
Spring crops245
Summer crops/double cropping38
Shrubland766
Other vegetation (riparian mainly)6
Bare land163
Impervious ground58
Table 5. Accuracy assessment and comparison of MSAVI-2 and NDVI, with and without third (NIR) cluster and masking step.
Table 5. Accuracy assessment and comparison of MSAVI-2 and NDVI, with and without third (NIR) cluster and masking step.
MSAVI-2NDVI
2-Steps+NIR2-Steps+NIR
Producer’s accuracy0.950.940.930.92
User’s accuracy0.490.670.520.68
F-Score0.650.780.660.78
Table 6. Commission errors when using MSAVI-2 and NDVI composites, with and without the third (NIR) cluster and masking step.
Table 6. Commission errors when using MSAVI-2 and NDVI composites, with and without the third (NIR) cluster and masking step.
MSAVI-2 CompositesNDVI Composites
2-Steps+NIR2-Steps+NIR
Evergreen orchards (traditional)0.360.290.240.20
Evergreen orchards (intensive)0.810.730.710.66
Deciduous orchards (all types)0.570.080.520.09
Greenhouses0.300.050.240.04
Winter crops0.020.010.010.01
Spring crops0.010.010.010.00
Summer crops / double cropping0.080.040.050.02
Shrubland0.000.000.000.00
Other vegetation (riparian mainly)0.980.130.960.14
Bare land0.000.000.000.00
Impervious ground0.000.000.000.00
Table 7. User’s and producer’s accuracies and F-Score after merging classes.
Table 7. User’s and producer’s accuracies and F-Score after merging classes.
MSAVI-2 + NIRNDVI + NIR
Producer’s accuracy0.870.84
User’s accuracy0.910.91
F-Score0.890.87
Table 8. Commission errors after “sieving” at a 500-pixel threshold using NDVI and MSAVI-2.
Table 8. Commission errors after “sieving” at a 500-pixel threshold using NDVI and MSAVI-2.
Commission Errors
Land Use Class“Sieved” MSAVI-2 + NIR“Sieved” NDVI + NIR
Evergreen orchards (traditional)0.160.06
Evergreen orchards (intensive)0.750.65
Deciduous orchards (all types)0.040.06
Greenhouses0.020.01
Winter crops0.010.01
Spring crops0.000.00
Summer crops / double cropping0.000.00
Shrubland0.000.00
Other vegetation (riparian mainly)0.070.08
Bare land0.000.00
Impervious ground0.000.00
Table 9. Producer’s and user’s accuracy and F-Score of “sieved” cluster, before and after merging intensive and super-intensive olive plantations.
Table 9. Producer’s and user’s accuracy and F-Score of “sieved” cluster, before and after merging intensive and super-intensive olive plantations.
Super-intensive olive plantations only
“Sieved” MSAVI + NIR“Sieved” NDVI + NIR
Producer’s accuracy0.950.93
User’s accuracy0.690.72
F-Score0.800.81
Intensive and super-intensive olive plantations merged
“Sieved” MSAVI + NIR“Sieved” NDVI + NIR
Producer’s accuracy0.890.84
User’s accuracy0.950.95
F-Score0.920.89
Table 10. CCR and total area of each sampled batch.
Table 10. CCR and total area of each sampled batch.
20102020
PlantationsCCRTotal Area (ha)CCRTotal Area (ha)
Smallest 500.4650.7666
Median 500.361240.68141
Largest 500.9619560.923081
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Navarro, R.; Wirkus, L.; Dubovyk, O. Spatio-Temporal Assessment of Olive Orchard Intensification in the Saïss Plain (Morocco) Using k-Means and High-Resolution Satellite Data. Remote Sens. 2023, 15, 50. https://doi.org/10.3390/rs15010050

AMA Style

Navarro R, Wirkus L, Dubovyk O. Spatio-Temporal Assessment of Olive Orchard Intensification in the Saïss Plain (Morocco) Using k-Means and High-Resolution Satellite Data. Remote Sensing. 2023; 15(1):50. https://doi.org/10.3390/rs15010050

Chicago/Turabian Style

Navarro, Rebecca, Lars Wirkus, and Olena Dubovyk. 2023. "Spatio-Temporal Assessment of Olive Orchard Intensification in the Saïss Plain (Morocco) Using k-Means and High-Resolution Satellite Data" Remote Sensing 15, no. 1: 50. https://doi.org/10.3390/rs15010050

APA Style

Navarro, R., Wirkus, L., & Dubovyk, O. (2023). Spatio-Temporal Assessment of Olive Orchard Intensification in the Saïss Plain (Morocco) Using k-Means and High-Resolution Satellite Data. Remote Sensing, 15(1), 50. https://doi.org/10.3390/rs15010050

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop