Next Article in Journal
Assessing Sponge Cities Performance at City Scale Using Remotely Sensed LULC Changes: Case Study Nanjing
Next Article in Special Issue
Detecting Recent Crop Phenology Dynamics in Corn and Soybean Cropping Systems of Kentucky
Previous Article in Journal
Radiometric Correction of Multispectral UAS Images: Evaluating the Accuracy of the Parrot Sequoia Camera and Sunshine Sensor
Previous Article in Special Issue
Mapping Spatiotemporal Changes in Vegetation Growth Peak and the Response to Climate and Spring Phenology over Northeast China
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 13
Issue 4
10.3390/rs13040578
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mapping and Characterization of Phenological Changes over Various Farming Systems in an Arid and Semi-Arid Region Using Multitemporal Moderate Spatial Resolution Data

by
Youssef Lebrini
1,2,*,
Abdelghani Boudhar
1,3,
Ahmed Laamrani
3,
Abdelaziz Htitiou
1,2,
Hayat Lionboui
2,
Adil Salhi
4,
Abdelghani Chehbouni
3,5 and
Tarik Benabdelouahab
2
1
Faculty of Sciences and Techniques, Sultan Moulay Slimane University, Beni Mellal 23023, Morocco
2
Department of Environment and Natural Resources, National Institute of Agronomic Research, Rabat 10090, Morocco
3
Center for Remote Sensing Applications, Mohammed VI Polytechnic University, Ben Guerir 43150, Morocco
4
Geography and Development Group, Abdelmalek Essaadi University, Martil 93150, Morocco
5
Centre d’Etudes Spatiales de la Biosphère (CESBIO)/Institut de Recherche pour le Développement (IRD), UMR CNES-CNRS-INRAE-UPS, Université de Toulouse, 31401 Toulouse CEDEX 9, France
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(4), 578; https://doi.org/10.3390/rs13040578
Submission received: 17 December 2020 / Revised: 23 January 2021 / Accepted: 1 February 2021 / Published: 6 February 2021
(This article belongs to the Special Issue Vegetation Phenology from Remote Sensing data: Monitoring, Mapping, and Modelling)
Browse Figures
Versions Notes

Abstract

:
Changing land use patterns is of great importance in environmental studies and critical for land use management decision making over farming systems in arid and semi-arid regions. Unfortunately, ground data scarcity or inadequacy in many regions can cause large uncertainties in the characterization of phenological changes in arid and semi-arid regions, which can hamper tailored decision making towards best agricultural management practices. Alternatively, state-of-the-art methods for phenological metrics’ extraction and long time-series analysis techniques of multispectral remote sensing imagery provide a viable solution. In this context, this study aims to characterize the changes over farming systems through trend analysis. To this end, four farming systems (fallow, rainfed, irrigated annual, and irrigated perennial) in arid areas of Morocco were studied based on four phenological metrics (PhM) (i.e., great integral, start, end, and length of the season). These were derived from large Moderate resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) time-series using both a machine learning algorithm and a pixel-based change analysis method. Results showed that during the last twenty-year period (i.e., 2000–2019), a significant dynamism of the plant cover was linked to the behavior of farmers who tend to cultivate intensively and to invest in high-income crops. More specifically, a relevant variability in fallow and rainfed areas, closely linked to the weather conditions, was found. In addition, significant lag trends of the start (−6 days) and end (+3 days) were found, which indicate that the length of the season was related to the spatiotemporal variability of rainfall. This study has also highlighted the potential of multitemporal moderate spatial resolution data to accurately monitor agriculture and better manage land resources. In the meantime, for operationally implementing the use of such work in the field, we believe that it is essential consider the perceptions, opinions, and mutual benefits of farmers and stakeholders to improve strategies and synergies whilst ensuring food, welfare, and sustainability.

Graphical Abstract

1. Introduction

Accurate monitoring of land use changes in arid landscapes and understanding of these changes drivers are crucial in arid environment research, especially since changes in farming systems affect directly socioeconomic as well as environmental sectors. Improving global food security will need a good understanding of the behavior of farming systems and their responses to climate- and human-induced factors [1,2]. In African arid and semi-arid regions, there is a pressing need for characterizing the behavior and the distribution of farming systems for better monitoring of agricultural lands and, thus, managing of land resources. In this context, accurate management and monitoring of farming systems, particularly at a large scale, can substantially benefit from geo-spatial technologies and remote sensing data [3]. However, in semi-arid regions, phenological data are scarce and difficult to produce due to the lack of sufficient ground information and limited studies [4]. Therefore, the use of remotely sensed phenological metrics (PhM) may provide a viable alternative to costly and time-consuming field sampling [5,6]. In addition, the change analysis of PhM could provide valuable information that would identify current trends as quickly as possible and would be easily understandable by political authorities as well as by land managers to take actions quickly for better control of agricultural lands. In the same context, information on vegetation cover provides an insight on changes occurring and could be a potential indicator of food security and sustainability [7,8,9,10]. To this end, the use of remote sensing has demonstrated a strong potential for understanding and detecting phenological changes because of consistent and frequent coverage [11]. Data from remote sensing satellites provide large and continuous observations that characterize the changes occurring on the earth [11,12]. Indeed, time-series of satellite data are suitable to monitor the spatiotemporal behavior of plant phenology [13,14,15]. These issues motivate studying changes in farming systems to characterize the spatiotemporal variability that has occurred over long periods related to different drivers, i.e., short- and long-term weather events and public policies.
Many studies have proved the existence of a consistent relationship between vegetation indices and weather indices (i.e., temperature and precipitation) [16,17,18,19,20,21]. During recent decades, the scientific community has investigated the impact of land management on specific indices (e.g., normalized difference vegetation index, NDVI) in order to link vegetation changes to anthropogenic actions [22]. With this in mind, many researchers went further in studying changes in vegetation cover, based on satellite-derived products, by developing techniques and tools for extracting and analyzing PhM, especially for monitoring farming systems [5,23,24,25]. They demonstrated the ability of these metrics to monitor and differentiate vegetation cover based on contrasted phenological profiles [5,26,27,28,29]. Consequently, remote sensing-derived PhM have improved agricultural monitoring, whether it is to map farming systems [30,31], to analyze the dependence of vegetation cover on climate events (Cui et al. 2017), or to characterize the behavior of vegetation cover from different data sources [32].
Indeed, the use of Moderate resolution Imaging Spectroradiometer (MODIS) data provides an opportunity to characterize the spatiotemporal variability of vegetation cover at a large scale [33,34]. It allows the construction of time series of vegetation indices (VIs) (e.g., NDVI, EVI…), by taking advantages of the combination of accurate reflectance, frequent coverage, moderate spatial resolution and the relatively long period of data availability since 2000 [35].
In the scope of farming systems monitoring, different algorithms have been used. These include Random Forest (RF) as in Wang et al. [36], Support Vector Machine (SVM) as in Wessel et al. [37], Artificial Neural Networks (ANNs), and K-Nearest Neighbor (KNN) as in Yu et al. [38]. It is worth noting that all of the above-mentioned algorithms are considered as machine learning algorithms, which are based on automatic learning from a dataset to find the relationships between a predictor and a response [39]. Among all of them, the Random Forest (RF) [39] provides powerful classification of farming systems. It has several advantages, including the ability of running efficiently on a large volume of input variables, resisting noise or over-fitting, being relatively robust concerning outliers, and requiring fewer parameters [40,41]. These advantages make this algorithm the best choice for running the classifications of farming systems.
At present, many studies based on limited observations over some irrigated zones in the Oum Er-Rbia (OER) basin have been conducted, but fewer focused on large-scale analysis. Investigating the changes in farming systems at a large scale in the OER basin will have a far-reaching impact on agricultural sector development and how the system is adapting or not to climate changes. In this context, the present research sought to (i) use phenological metrics derived from twenty years of NDVI MODIS datasets (i.e., 2000–2019) to map and monitor changes in selected farming systems over a large arid-to-semi-arid region in Morocco (i.e., OER basin); (ii) investigate trends in selected farming systems at a large scale in the study area to evaluate how the systems are adapting or not to changes; and (iii) provide information on farming systems’ changes for stakeholders to adopt more accurate and efficient strategies.

2. Materials and Methods

2.1. Study Area

This study was conducted within the Oum Er-Rbia (OER) basin in west-central Morocco, between 31°15′–33°22′ N and 5°08′–8°23′ W (Figure 1). The OER basin covers an area of 38,000 km2, while its administrative area is around 50,000 km2. The OER basin is made up of five geographical units from the Atlas Mountains, the foothill areas, the plain, the phosphate plateau, and the coastal area [4]. A combination of flat and mountain terrains generally characterize the topography of the basin. Elevation ranges between 100 (e.g., in the western and coastal zones) and 3890 m (e.g., in the eastern zone) above sea level [42] (Figure 1). The OER River sources are located in the mountainous upstream zones, and the river covers a distance of 550 km, overpassing the Tadla irrigated perimeter (TIP), the coastal areas, and the northern zone of the Doukkala irrigated perimeter (DIP), and flows into the Atlantic Ocean at Azemmour city [4]. The climate is variable from humid in coastal and mountainous zones to semi-arid in the plains, with cold winters and dry summers [42]. Annual rainfall average varies from 230 to 1000 mm in the plains and the Atlas Mountains, respectively [43]. The agricultural season generally occurs between September and June, while the most important amount of rainfall is received between October and May (i.e., 70% to 80% of the annual rainfall). The annual temperature varies between a maximum of 46 °C in August and a minimum of −3 °C in January. The investigated areas are primarily agricultural, with irrigated crops (i.e., cereals, sugar beet, and alfalfa) and rainfed areas (wheat). The study area is also dominated by tree cultivations, especially pomegranate and citrus trees, which present a permanent activity during the season (Figure 2).
In this study, four farming system types (Figure 2A–C) have been selected to evaluate the effectiveness of phenological metrics in the characterization of changes. The farming systems are named (i) Irrigated Perennial Crop (IPC), (ii) Irrigated Annual Crop (IAC), (iii) Rainfed Area (RA), and (iv) Fallow (FA).
The phenology of the IPC farming system is characterized by permanent vegetative activity during the cropping season and high biomass production (Figure 2). IAC represents high inter-annual variability with a high amplitude and rapid growth and senescence moments. The IAC farming system is described with high vegetative activity with a maximum NDVI generally in March, low dependency on climate conditions and soil type, and a maintained growth cycle from seed to harvest time. The RA and FA farming systems are typical on semi-arid lands, where the start of vegetative activity depends on the climate conditions, especially the first rainfall. Figure 2 gives deep informative illustrations of the four farming systems investigated in this study.

2.2. Multitemporal MODIS and CHIRPS Data Acquisition and Derived NDVI Processing

In this study, we used the MOD13Q1 NDVI product. A dataset of 437 images covering the OER basin was downloaded for the period between 2000 and 2019 via the U.S. Geological Survey (USGS) Land Processes Distributed Active Archive Center (LP-DAAC). The MOD13Q1 NDVI data are within a spatial resolution of 250 m and a temporal frequency of 16 days in sinusoidal projection [4]. Time series data were re-projected to UTM projection to be used further in this study. Rainfall data from Climate Hazards group Infrared Precipitation with Stations (CHIRPS) [44], available from 1981, were summed over the OER basin for each year within the MOD13Q1 product period. The CHIRPS data are delivered at two spatial resolution levels, 0.05 and 0.25°. In our study, CHIRPS data converted to a 16-day period with a spatial resolution of 0.05° were used to extract rainfall data and plotted against NDVI profile.
We performed NDVI time-series filtering and PhM extraction with TIMESAT software [45], a commonly used package used for (1) data filtering and smoothing and (2) phenological metrics extraction (Figure 3). In order to obtain high-quality time series and to fill the gaps in data, TIMESAT implements three different smoothing methods, which are Savitzky–Golay (SG) [46], asymmetric Gaussian (AG) [27,47], and double logistic (DL) [48]. The AG filtering method is less sensitive to noise than the other approaches [45,49]. The study performed by Fu et al. [50] shows the robustness of this algorithm in smoothing time series and preserving the characteristics of NDVI profile, contrary to the DL and SG methods. The AG method can generate a smoothed NDVI profile while describing minor changes in the NDVI sequence data [51]. Therefore, the AG method was chosen as the time series reconstruction method to be used in this study. Table 1 shows the setting used for NDVI time series processing in TIMESAT. More details on the filtering approach used and phenological metrics extraction for this study can be found in Lebrini et al. [4].
The second step was PhM extraction; such measurements are usually calculated using a common method based on value thresholds of the seasonal vegetation amplitude, assuming that a particular phenomenon has started when the NDVI values surpass a given threshold (Figure 3). Our study was, thus, carried out by setting the proportion of the seasonal amplitude to 10% measured from the left and right minimum values, respectively [45]. In general, four PhM were extracted using the TIMESAT program for trend analysis. Definitions of the PhM used in this study are explicitly defined in Table 2.

2.3. Statistical Analysis

2.3.1. Random Forest Classification and Change Analysis

To characterize the main cropping systems, the ground data used in this research were collected through ground truthing exercises over farming systems during the 2018–2019 cropping season. The collected data were reported using a GPS receiver, with a positional error of less than 2 m to generate data for classification training and accuracy assessments for the 2018/2019 dataset. The reference points were collected in a way so as to represent the full variety of farming system elements in the study area. Accuracy assessment of other classification results was performed using similar or near-similar ancillary data from MODIS and Google Earth images and high-resolution aerial photographs. A land cover map originated from Glob Cover was used to mask the agricultural zones over the study area [4,55]. A summary of the ground data is given in Table 3.
To classify farming systems over the OER basin, the supervised Random Forest (RF) classifier was used based on the CARET package within R [56]. RF is a non-parametric machine learning classifier that combines a random selection of training subsets of data with an ensemble of trees. Recent studies reveled the effectiveness of the RF classifier in the remote sensing field, including land phenology mapping [39,57,58,59]. For measuring the accuracy of the classification results using the RF classifier, ground truth data were split into two sets of training (80%) and testing (20%) samples using a spatial cross-validation approach with 5 folds from the CAST package in R [60]. This spatial cross-validation helped to make sure that the ground truth sample of the same field will be either in the test or in the training data in order to avoid over-fitting. The accuracy assessment was performed by calculation of many accuracy metrics for the classifications results, including overall accuracy (OA), Kappa coefficient, producer’s accuracy (PA), user’s accuracy (UA), and F1-score. The classification of 2000/2001 was generated using the “predict” function from the CARET package in R [56]. In order to increase the reliability and validity of our RF classification model, and as an additional check of the resultant information of study area-specific classification accuracy, a second accuracy assessment was performed for the 2000/2001 classification map using the same methods as in the accuracy assessment for the 2018/2019 classification map. In order to update testing polygons based on the 2000/2001 situation, we used the testing polygons from the 2018/2019 data and plotted them against the smoothed NDVI profiles from MODIS data and Google Earth images. From the MODIS time series, we investigated the correspondence between the NDVI profile of each testing polygon and the farming system. From imagery, we added necessary polygons to perform the accuracy assessment. Furthermore, farming system (FS) maps obtained from the 2000/2001 and 2018/2019 classification results were used to map changes that occurred over the study area. The transition between FS classes revealed by comparing the classification results was used to extract unchanged FS during the study period. In order to have significant and robust results of the further trend analysis, we opted to compute trend on unchanged farming systems resulted from the change analysis step. Figure 4 shows a detailed flowchart of the adopted methodology in this study.

2.3.2. Variability and Trend Analysis

From the annual PhM previously retrieved, we calculated their per-pixel temporal mean and coefficient of variation (CV). In order to assess the spatial distribution of PhM showing improvement (positive change) or degradation (negative change), we employed the non-parametric Mann–Kendall (MK) trend test [61] to determine trends on PhM over the OER basin between 2000 and 2019 [34].
The Z statistic follows the standard normal distribution with zero mean and unit variance under the null hypothesis of no trend. A positive Z value indicates an upward trend whereas a negative value indicates a downward trend. Probability (p) represents the measure for evidence to reject the null hypothesis, and positive p values show a negative trend whereas positive p shows a positive trend. We calculated the MK trend test separately for every pixel for the start and end of season, the length of season, and the great integral over the OER basin. These parameters are crucial indicators of seasonal productivity. The Mann–Kendall statistical trend Z was determined as follows:
S =   ∑ k = 1 n − 1 ∑ j = k + 1 n s i g n ( X j − X k )   ,      w h e r e   j > k
where n is the length of time series data, Xk and Xj are the observations at k and j time, respectively, and
s i g n ( X j − X k ) = { 1 i f   ( X j − X k ) > 0 0 i f   ( X j − X k ) = 0 − 1   i f   ( X j − X k ) < 0
The probability linked to the Mann–Kendall statistic “S” and the selected n-data were determined to quantify the level of significance of the trend. Var(S) was calculated, and then the normalized test statistic Z was computed using the following equations:
V a r ( S ) = 1 18 ( n ( n − 1 ) ( 2 n + 5 ) − ∑ t t f t ( f t − 1 ) ( 2 f t + 5 ) )
where t varies over the set of tied ranks and ft is the number of times that the rank t appears (i.e., frequency). The equation used to calculate the Mann–Kendall significance Z-score is as follows:
Z = { S − 1 V a r ( S )    f o r   S > 0 0    f o r   S = 0 S + 1 V a r ( S )    f o r   S < 0
where Var(S) is the variance of the dataset and n is the number of data points.
The equation used to estimate the Theil–Sen (TS) median trend is:
T S = M e d i a n ( X j − X k t j + t k )
This robust non-parametric trend operator is highly recommended for assessing the rate of change in time-series data. It is calculated by determining the slope between every pairwise combination and then finding the median value. Using the Mann–Kendall test and Theil–Sen median trend analysis, the trends in PhM were described accordingly, and the significance level of the changes in NDVI trends was determined using the Z-score at p-value below 0.1 significance level.

3. Results

3.1. NDVI and Rainfall Time Series Analysis

Given the importance of the NDVI in monitoring vegetation cover, the spatiotemporal variation of this index was assessed based on reference data collected from fieldwork over farming system types between the 2000/2001 and 2018/2019 cropping seasons and the average rainfall over the OER basin for the same period (Figure 5).
A visual analysis of the temporal NDVI values shows different patterns associated with each farming system type. These patterns were characterized by a specific range of NDVI values and represent the seasonality and the growth cycle of each farming system. The NDVI values for IPC range between 0.45 and 0.9, which is an indicator of the permanent photosynthetic activity of perennial systems (Figure 5). The studied irrigated tree crops are generally carried out in intensive systems. Owing to their received quantities of water and their long life cycle, this farming system shows a unique characteristic of persistence during the growing season from 2000 to 2019. The IAC class has some identical patterns to the IPC farming system such as the water supply and the high value of production. The IAC farming system shows NDVI values that range from 0.3 to 0.8. The NDVI values start increasing in September and decrease in May of each agricultural season. Their variability in length of the growing season is mainly related to sowing timing, the management of agricultural land, and other plant physiological disease problems.
For the RA farming system, the NDVI values varied from 0.2 to 0.7 for wet seasons and between 0.2 and 0.4 for dry seasons (i.e., 2006/2007 and 2008/2009 cropping seasons). In general, their NDVI response to rainfall is systematic. Indeed, agricultural lands in the RA farming system are entirely dependent on rainfall, which varies in amount from one year to the next and directly affects agricultural productivity. Indeed, in semi-arid regions, this relationship between rainfall and NDVI has been demonstrated in numerous research studies [62,63,64]. Thereby, the RA farming system shows high variability in terms of the start and end of season and the length of season. This variability is linked to the strong dependency of this farming system on the climatic conditions, especially the rainfall amount. The vegetative activity begins earlier in the case of a wet season and later for a dry season, and the same pattern could be observed for the end of the vegetative activity. A similar pattern could be observed in the FA farming system, where the climatic conditions have more effects on the vegetative activity. These effects of rainfall could be well spotted in the 2009/2010 cropping season over FA and RA farming systems, and we can precisely observe the start of the vegetative activity after receiving an important amount of rainfall (Figure 5). Through the coming sections, the paper explores more the behavior of farming systems using trend analysis techniques.

3.2. Spatial Patterns of Phenological Metrics

To mutually characterize PhM behavior (i.e., TSOS, TEOS, LOS, and GINT) and their variability, we combined them with their coefficients of variation (CVs), which are used as a measure of reliability, into the bivariate maps shown in Figure 6. These bivariate maps simultaneously provide a spatial representation of (1) where high or low variability in PhM is expressed for each farming system and (2) the risk involved in the future management and seeding in these farming systems. This combination of factors can assist in developing knowledge about which farming systems can be cultivated in an average year while considering that in such regions, farming system choices are largely based on the potential for attaining good yields and the risk of season failure [65].
Figure 6D shows that areas located inside the irrigated perimeters with a start of the season (TSOS) that occurred between February and June generally have a CV below 0.1. Encouraged by the availability of irrigation water throughout the year, farmers in these areas are diversifying and intensifying their agriculture, thus putting a larger area under crops with high added value, whatever their water consumption [66]. Outside the irrigated perimeters, especially the center of the study area, strong variability is apparent (i.e., variability above 0.2). This high variability is more accentuated by moving towards the mountainous zones. Farmers prefer not to invest excessively in these zones in order to avoid the risks associated with drought, related mainly to the lower stability of TSOS and, therefore, the rainfall amount and distribution. Areas outside the irrigated perimeters with TSOS between August and September have CV values between 0.1 and 0.3 for practically all locations. Nonetheless, the declining availability of water during the studied period causes a decrease in productivity and leads farmers to opt for crops that require less water but offer good margins in order to maximize their profits and avoid the risk of season failure. Other regions outside the irrigated perimeters with a TSOS between days 242 and 365 of the year have a CV below 0.15. Since these areas may present a low-risk investment opportunity, they could be a good choice for the implementation of new farming systems. Towards the western part of the basin, low variability is apparent (Figure 6D).
On the other hand, areas located outside the irrigated perimeters with an end of season (TEOS) that occurred between the day of year (DOY) 33 and 177 express two degrees of variability, the first with a CV below 0.15 and the second with a CV above 0.35. The majority of the area that has a TEOS between DOY 178 and 241 has low variability, with a CV below 0.15. Inside the irrigated perimeters, low variability is expressed, with a CV below 0.15 for areas having the TEOS between DOY 178 and 241 and between 242 and 273. The high variability in terms of CV observed for the TSOS and TEOS could be related to the dependence of the start of a cropping season to climatic conditions—in rainfed farming systems, the season could not start without the first rainfall. Adversely, the end of the cropping season depends on the climatic conditions but also depends on the physiological properties of the plant and its persistence (Figure 6C).
Considering the phenological metric LOS, the results of this research show that the variability is inversely proportional to the value of the length of the season, without regard to particulars or exceptions. Hence, the areas with an LOS value below 160 days generally have a CV value above 0.35, and the areas with an LOS value between 160 and 192 days have CV values between 0.15 and 0.35, especially in the center of the study area. However, these areas may express a deficit in agricultural production during years with low rainfall amounts. The low variability is mostly located inside the irrigated perimeters where the LOS is above 224 days. These zones shows CV values of below 0.15. Various areas in the western part of the basin also show high CV values which are above 0.28 (Figure 6B).
Finally, by combining the fourth phenological metric explored during this study with its CV values (GINT), our results can help in the characterization of the biomass variability of farming systems encountered in this study area (Figure 6A). The areas having GINT values below four have CV values above 0.35. These zones are generally rainfed areas and fallow, where climatic conditions have an effect. Generally, these areas are not recommended for the seeding of plants with high values of production. Instead, areas with GINT values above eight have CV values between 0.15 and 0.25. This moderate variability in biomass productivity could be considered as normal, considering the differences in crop types over the irrigated parameters (i.e., alfalfa, sugar beet, wheat, etc.). Other zones in the OER basin remain with low to moderate biomass productivity with high variability of the GINT metric.

3.3. Determination of Unchanged Farming Systems’ Area

RF classification offers a powerful algorithm to classify the spatial patterns of farming systems in the study area. In this study, the classification of farming systems was produced from the implementation of the Random Forest (RF) classifier based on PhM. The classifier was applied and the accuracy assessment results are summarized and provided in the Supplementary Materials (Figures S5–S7). In terms of the individual accuracy, the results present higher overall accuracy and much more balanced producer/user accuracy. The overall accuracy for RF was 97% and 93% with a kappa value of 0.95 and 0.91 for the 2018/2019 and 2000/2001 cropping seasons, respectively. In general, all farming system classes achieved over 90% user accuracy. The RF classifier also produced over 90% producer accuracy for most farming system classes. The change detection results from 2000 to 2019 reveal an important dynamic between the different FSs in the studied region. The transition between FSs from 2000/2001 to 2018/2019 shows that there has been a marked change in FS classes over the study area (Figure 7). To evaluate trends over the study area, we only considered the unchanged classes during the studied period.

3.4. Trend Analysis Results

This research further analyzed the inter-annual trends for the four PhM (i.e., SOST, EOST, LOS, and GINT). The Mann–Kendall test was used to identify significant trends at 90% confidence level and Sen’s slope estimator was applied to compute the magnitude of trends for the period between 2000 and 2019, alternating between increasing and decreasing trends over the study area. For the start of season (TSOS), Figure 8 displays its trends as obtained using Mann–Kendall and Sen’s slope tests. Positive and negative values refer to delayed and advanced dates of TSOS, respectively.
Most areas in the OER basin did not show a significant trend during the past twenty years for the TSOS, LOS, and GINT metrics (p > 0.1). On the contrary, the TEOS shows a significant trend over most parts of the study area (p < 0.1). In this section, only significant trend pixels are considered for description and analysis. A significant negative TSOS trend was discovered across small parts of the Doukkala irrigated perimeter, the Tadla irrigated perimeter, and in parts of the rainfed area. For the pixels with negative trends, over irrigated perimeters, Sen’s slope shows an average decrease in TSOS of approximately −0.2 day/year; outside the irrigated perimeters, the TSOS shows a decrease between −0.2 and −0.4 day/year. The early start of the agricultural season in these zones could explain these results, especially the irrigated annual crops where we found different crop types and irrigation systems. Hence, this implies that over the studied period of 20 years, the TSOS has been delayed by more than 6 days (−6 day/20 years). Other regions of the OER basin showed a non-significant trend. These regions are not interpreted since the significance level is above 0.1.
On the other hand, we found a significant positive trend for the TEOS across large parts of the OER basin and in a part of the Tadla irrigated perimeter. However, we found a significant negative trend in the rainfed and fallow farming systems (Figure 9).
For the pixels with positive trends, Sen’s slope shows an average increase in TEOS above 0.2 day/year in the TIP and between 0 and 0.2 in the rainfed area. The TEOS shows a decrease between 0 and −0.2 days/year in numerous parts of the OER basin, especially in the southern part. These results could be explained by the satisfaction of plants requirements in the irrigated zones, whereas in other areas, the negative trend is essentially linked to climatic conditions and poor land management. This signifies that over the studied period of 20 years, the TEOS has changed by more than 3 days/20 years in irrigated areas and by −3 days/20 years in rainfed and fallow areas.
After the trend analysis of the TEOS metric, we directed our studies towards the analysis of trends in the length of season (LOS) (Figure 10) to investigate the link between TEOS and TSOS with the LOS, since the LOS is the result of the difference between the TEOS and the TSOS. The results show a significant positive trend in LOS over the south parts of the Doukkala irrigated perimeter and in the north parts of the Tadla irrigated perimeter. In addition, we found a significant negative trend in the rainfed and fallow farming systems.
For the pixels with positive trends, Sen’s slope gives an average increase in LOS above 0.1 day/year in the Tadla irrigated perimeter and approximately between 0 and 0.1 days/year in the DIP. The LOS shows a decrease between 0 and −0.1 day/year in numerous parts of the OER basin. This signifies that over the studied period of 20 years, the LOS has changed by more than 1.6 days/20 years in irrigated areas and by −1.6 days/20 years in rainfed and fallow areas.
Comparing these results with trends in TSOS and TEOS indicates that increases in TEOS and decreases in TSOS dates are mostly responsible for the positive LOS trends in irrigated perimeters of Doukkala and Tadla. Negative LOS trends in the rainfed area and fallow farming systems are related to a delay in TSOS dates and an early TEOS.
Finally, considering the GINT metric, the trend of this index was obtained using Mann–Kendall and Theil–Sen median trend tests (Figure 11). We found significant positive trends over large farming systems. These areas are specially located in the irrigated perimeters of Tadla and Doukkala, in addition to some zones with improved agricultural management practices.
Generally, significant positive trends are located in the western and eastern parts of the OER basin. The central region is dominated by insignificant negative trends, while stable farming systems are dispersed over the study area. Irrigated perimeters show an increase in GINT by 0.2/year, while the rainfed area shows an increase between 0 and 0.1/year. Generally, the trend results reveal an important variation between the different FSs in the studied region.

4. Discussion

In this study, we have demonstrated that the phenological pattern of vegetation cover across different farming systems and across different regions over the last 20 years has an important implication in vegetation dynamics and climate change. It provides an insight of the vegetation cover status of a region affected by climate risks. To our knowledge, this study is the first to demonstrate the ability of annual phenological metrics (PhM) and RF modeling to map and characterize changes over the major farming systems in Morocco. The use of TIMESAT for NDVI time-series processing and analysis has provided a more comprehensive investigation of the NDVI behavior from a phenological perspective. As far as we know, this study is the only one that has explored phenological farming systems’ variability over North Africa at 250-m spatial resolution.
Overall, the accuracies obtained in this study for the classification of farming systems were reliable and consistent with those revealed for other MODIS-based land cover classifications [30,31]. For comparison, the resulting overall accuracy of the RF classifications obtained over the OER basin scale (i.e., 90–96%) is similar to that obtained using MODIS time series for mapping cropping intensity in China at a regional scale (i.e., 92%; [31]). On the side of classification performance, the confusion between rainfed area and fallow, which was the main source of error, reflects the strong similarity between these two farming systems, especially when the climatic conditions are critical. Consequently, our findings have confirmed the advantage of using a combination of PhM from MODIS time series and an RF classifier to discriminate between farming systems [67]. In addition to the influence of climate conditions on phenological changes, changes in crop type were also an important determinant of phenological changes and trends in the OER basin. For this reason, in this study, we used a mask layer to mask only unchanged areas and ejected the effect of annual land use change over farming systems.
In terms of PhM variability and trends, our study revealed substantial variability over the studied PhM. Generally, TSOS and TEOS show high and low variability inside and outside the irrigated perimeters, respectively (Figure 6). The TSOS inside irrigated perimeters displayed earlier onset trends than other zones, which could be due to the sowing date and irrigation (Figure 7). On the other hand, delayed TEOS could be explained by the variability of harvest time and biological factors and the climatic conditions that occurred during the cropping seasons. These delays and early onsets in time of occurrence between TSOS and TEOS were translated to a longer cropping season over the irrigated perimeters (Figure 9). The remarkable shift expressed by GINT and LOS over the irrigated perimeters is the result of an intensification plan and a change of agricultural practices with the encouragement of farmers to plant perennial crops. The question raised here is how sustainable are the water and soil resources supporting this mode of production under the growing pressure applied?
Our study has also revealed that annual variability outside the irrigated perimeters is strong, especially in the center of the study area, where fallow and rainfed farming systems extend. This high variability can be explained by the changes in rainfall amount and distribution received in these regions during the cropping season and the changes to the agricultural production potential. As identified by Ouatiki et al. [68], climatic stations located in these regions experienced a significant decreasing trend and received a low rainfall amount during the season. Other factors leading to the high variability of PhM in these regions were the depletion of groundwater and the quality of agricultural management, since these regions are considered as rainfed and fallow areas where the accessibility to irrigation supplies is absent [4]. A similar situation was observed in the variability of the GINT metric (Figure 6A) as the same regions experienced high variability and low values of GINT, which signify low biomass production and an unbalanced field for agricultural planning. Regarding the trend results, GINT increased during the studied period by 0.1/year in most areas of the basin, while for certain regions of FA and RA, it decreased by −0.1/year (Figure 11). Therefore, we believe that all of these factors affect the productivity of these two farming systems, which are translated in the short cropping season and the low biomass production (Figure 6A,B, respectively).
The wide spatial distribution of annual variability of PhM found in this study is reliable and comparable with the outcomes from previous studies [65,69,70]. For instance, Vrieling et al. [65] identified that a higher variability of Length of Growing Period (LGP)—represented in our study by the length of season metric (LOS)—is generally found in arid and semi-arid areas, with coefficients of variation above 0.25 ( i.e., LOS in our study varied from 0.15 to 0.35). This slight difference between our studies could be due to the higher spatial resolution of MODIS data (i.e., 250 m). Contrarily, in the irrigated perimeters, where land receives significant amounts of water and fertilizer, the variability of PhM remained lower. Furthermore, the variability in crop types may have had some influence on the variability of PhM. However, this hypothesis is uncertain as our study focused only on unchanged farming systems of the OER basin.
Weather extremes and consecutive droughts in this region strongly affected the vegetation cover dynamics and resulted in adaptations of farming system management in response to climatic variation. Socio-economic policies and improving farm practices were also dominant drivers of farming system changes in addition to climate conditions. These factors could all have affected the variability of the start, the end, and the productivity of the cropping season, resulting in phenological changes over time. A good example of agricultural planning that made great changes in the agricultural sector in Morocco is the Green Morocco Plan strategy, which aimed to encourage tree-growing extensions that led to a high variability over some parts of the perennial farming system.
Notably, our findings showed that the variability of PhM is varied and related to farming system types. The results from this study are expected to constitute a background to other research about drought monitoring and desertification studies in arid regions. The PhM datasets and the trend results could be integrated with climate data to perform an estimation of crop water requirements and offer a tool for managers and stakeholders to make decisions for the extension of agricultural areas according to the available water resources in a context of water stress.

5. Conclusions

Finally, variation in phenological metrics over FSs was estimated at the OER basin level during 2000–2019 seasons using MODIS NDVI data and trend analysis tests. Our study findings are the following:
(1) Over-irrigated perimeters’ (TIP and DIP) mean LOS, GINT, and TEOS values showed low variability. On the other hand, moderate variability was observed for the mean TSOS values during the studied period. Contrary to the irrigated zones, PhM over the rainfed and fallow farming systems showed high variability. This variability over RA and FA is justified by the irregularity of rainfall amounts received over these farming system areas.
(2) Trends over farming systems are not uniform at the OER basin level. Most areas in the OER basin did not show a significant trend during the past 20 years for the TSOS, LOS, and GINT metrics (p > 0.1). Contrary to the TEOS, where a significant trend was observed (p < 0.1), TSOS shows early onset over the IPC and IAC (i.e., 0.2 days/year), while over RA and FA, it was delayed by −0.2 days/year, especially in the center of the basin. Other regions of the FA and RA showed extended TSOS by 0.2 days/year. TEOS shows early onset (i.e., −0.4 days/year) over the FA, RA, and part of the IPC. Other regions of the basin showed a TEOS extended by 0.2 days/year. LOS generally slightly increased over the farming systems, except in particular zones of the FA, and did not advance markedly during the study period. GINT increased during the studied period by 0.1/year in most areas of the basin, while for certain regions of FA and RA, it decreased by −0.1/year.
The investigation of changes in PhM over twenty years (19 cropping seasons) proved the ability of these metrics to characterize the spatiotemporal changes over the OER basin. Nevertheless, the need to take into account the perceptions and opinions of local populations is essential in order to reduce the process of vegetation cover degradation and to better manage natural resources.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-4292/13/4/578/s1. Figure S1: Bivariate map showing, simultaneously, the start of the season time and its coefficient of variation. Figure S2: Bivariate map showing, simultaneously, the end of the season time and its coefficient of variation. Figure S3: Bivariate map showing, simultaneously, the length of season and its coefficient of variation. Figure S4: Bivariate map showing, simultaneously, the great integral metric and its coefficient of variation. Figure S5: Classification result of the RF classifier for 2000/2001 cropping season. Figure S6: Classification result of the RF classifier for 2018/2019 cropping season. Table S7: Accuracy assessment of the 2000/2001 and 2018/2019 classification results.

Author Contributions

All the authors of this manuscript have contributed substantially to the work reported. Conceptualization, Y.L., A.B. and T.B.; methodology, Y.L., A.B. and T.B.; software, Y.L.; validation, Y.L., A.B. and T.B.; formal analysis, Y.L., A.B. and T.B.; investigation, Y.L., A.B. and T.B.; resources, Y.L., A.B., A.L., A.C. and T.B.; data curation, Y.L.; writing—original draft preparation, Y.L., A.L., A.B.; writing—review and editing, A.B., A.L., A.H., H.L., A.S., A.C. and T.B.; visualization, Y.L.; supervision, A.B. and T.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Acknowledgments

Support was provided by the National Center of Scientific and Technical Research (CNRST, Morocco) for Youssef Lebrini as a Ph.D. scholarship (Scholarship N°186USMS). The authors are grateful to NASA’s team for creating and making the vegetation index product MOD13Q1 freely available. A first version of this work was presented at a seminar on “Change Detection from satellite data” at the Center for Remote Sensing Applications (CRSA), Mohammed VI Polytechnic University (UM6P). We also thank the academic editor and three anonymous reviewers for their helpful comments which greatly improved earlier versions of the manuscript. We would like to acknowledge the Department of Strategy and Statistics (DSS) from the Ministry of Agriculture in Morocco for the financial support as part of the agreement between National Institute of Agronomic Research (INRA) and DSS.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lionboui, H.; Benabdelouahab, T.; Htitiou, A.; Lebrini, Y.; Boudhar, A.; Hadria, R.; Elame, F.; AbdelGhani, B. Spatial assessment of losses in wheat production value: A need for an innovative approach to guide risk management policies. Remote. Sens. Appl. Soc. Environ. 2020, 18, 100300. [Google Scholar] [CrossRef]
  2. Benabdelouahab, T.; Lebrini, Y.; Boudhar, A.; Hadria, R.; Htitiou, A.; Lionboui, H. Monitoring spatial variability and trends of wheat grain yield over the main cereal regions in Morocco: A remote-based tool for planning and adjusting policies. Geocarto Int. 2019. [Google Scholar] [CrossRef]
  3. Sishodia, R.P.; Ray, R.L.; Singh, S.K. Applications of Remote Sensing in Precision Agriculture: A Review. Remote. Sens. 2020, 12, 3136. [Google Scholar] [CrossRef]
  4. Lebrini, Y.; Boudhar, A.; Htitiou, A.; Hadria, R.; Lionboui, H.; Bounoua, L.; Benabdelouahab, T. Remote monitoring of agricultural systems using NDVI time series and machine learning methods: A tool for an adaptive agricultural policy. Arab. J. Geosci. 2020, 13, 1–14. [Google Scholar] [CrossRef]
  5. Atzberger, C.; Klisch, A.; Mattiuzzi, M.; Vuolo, F. Phenological Metrics Derived over the European Continent from NDVI3g Data and MODIS Time Series. Remote. Sens. 2013, 6, 257–284. [Google Scholar] [CrossRef] [Green Version]
  6. Schreier, J.; Ghazaryan, G.; Dubovyk, O. Crop-specific phenomapping by fusing Landsat and Sentinel data with MODIS time series. Eur. J. Remote. Sens. 2020, 1, 1–12. [Google Scholar] [CrossRef]
  7. Benabdelouahab, T.; Balaghi, R.; Hadria, R.; Lionboui, H.; Minet, J.; Tychon, B. Monitoring surface water content using visible and short-wave infrared SPOT-5 data of wheat plots in irrigated semi-arid regions. Int. J. Remote. Sens. 2015, 36, 4018–4036. [Google Scholar] [CrossRef]
  8. Benabdelouahab, T.; Dominique, D.; Hayat, L.; Rachid, H.; Bernard, T.; Abdelghani, B.; Riad, B.; Youssef, L.; Hamid, M.; Christian, B. Using SAR Data to Detect Wheat Irrigation Supply in an Irrigated Semi-arid Area. J. Agric. Sci. 2018, 11, 1916–9760. [Google Scholar] [CrossRef]
  9. Hentze, K.; Thonfeld, F.; Menz, G. Evaluating Crop Area Mapping from MODIS Time-Series as an Assessment Tool for Zimbabwe’s “Fast Track Land Reform Programme”. PLoS ONE 2016, 11, e0156630. [Google Scholar] [CrossRef]
  10. Li, L.; Friedl, M.A.; Xin, Q.; Gray, J.; Pan, Y.; Frolking, S. Mapping Crop Cycles in China Using MODIS-EVI Time Series. Remote. Sens. 2014, 6, 2473–2493. [Google Scholar] [CrossRef] [Green Version]
  11. Lu, D.; Mausel, P.; Brondízio, E.; Moran, E. Change detection techniques. Int. J. Remote. Sens. 2004, 25, 2365–2401. [Google Scholar] [CrossRef]
  12. Coppin, P.J.I.; Nackaerts, K.; Muys, B.; Lambin, E. Digital change detection methods in ecosystem monitoring: A review. Int. J. Remote Sens. 2004, 25, 1565–1596. [Google Scholar] [CrossRef]
  13. Salhi, A.; Benabdelouahab, T.; Martin-Vide, J.; Okacha, A.; El Hasnaoui, Y.; El Mousaoui, M.; El Morabit, A.; Himi, M.; Benabdelouahab, S.; Lebrini, Y.; et al. Bridging the gap of perception is the only way to align soil protection actions. Sci. Total. Environ. 2020, 718, 137421. [Google Scholar] [CrossRef] [PubMed]
  14. Htitiou, A.; Boudhar, A.; Lebrini, Y.; Hadria, R.; Lionboui, H.; Elmansouri, L.; Tychon, B.; Benabdelouahab, T. The Performance of Random Forest Classification Based on Phenological Metrics Derived from Sentinel-2 and Landsat 8 to Map Crop Cover in an Irrigated Semi-arid Region. Remote. Sens. Earth Syst. Sci. 2019, 2, 208–224. [Google Scholar] [CrossRef]
  15. 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, 1–24. [Google Scholar] [CrossRef]
  16. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting plant phenology in response to global change. Trends Ecol. Evol. 2007, 22, 357–365. [Google Scholar] [CrossRef]
  17. Moulin, S.; Kergoat, L.; Viovy, N.; Dedieu, G. Global-Scale Assessment of Vegetation Phenology Using NOAA/AVHRR Satellite Measurements. J. Clim. 1997, 10, 1154–1170. [Google Scholar] [CrossRef]
  18. Zhou, L.; Tucker, C.; Kaufmann, R.K.; Slayback, D.; Shabanov, N.V.; Myneni, R. Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999. J. Geophys. Res. Atmos. 2001, 106, 20069–20083. [Google Scholar] [CrossRef]
  19. Castro, C.L.; Beltrán-Przekurat, A.B.; Pielke, R.A. Spatiotemporal Variability of Precipitation, Modeled Soil Moisture, and Vegetation Greenness in North America within the Recent Observational Record. J. Hydrometeorol. 2009, 10, 1355–1378. [Google Scholar] [CrossRef] [Green Version]
  20. Hanamean, J.R., Jr.; Pielke, R.A., Sr.; Castro, C.L.; Ojima, D.S.; Reed, B.C.; Gao, Z. Vegetation greenness impacts on maximum and minimum temperatures in northeast Colorado. Meteorol. Appl. 2003, 10, 203–215. [Google Scholar] [CrossRef] [Green Version]
  21. White, M.A.; De Beurs, K.M.; Didan, K.; Inouye, D.W.; Richardson, A.D.; Jensen, O.P.; O’Keefe, J.; Zhang, G.; Nemani, R.R.; Van Leeuwen, W.J.D.; et al. Intercomparison, interpretation, and assessment of spring phenology in North America estimated from remote sensing for 1982–2006. Glob. Chang. Biol. 2009, 15, 2335–2359. [Google Scholar] [CrossRef]
  22. DeFries, R.S.; Field, C.B.; Fung, I.; Collatz, G.J.; Bounoua, L. Combining satellite data and biogeochemical models to estimate global effects of human-induced land cover change on carbon emissions and primary productivity. Glob. Biogeochem. Cycles 1999, 13, 803–815. [Google Scholar] [CrossRef]
  23. Adole, T.; Dash, J.; Atkinson, P.M. Characterising the land surface phenology of Africa using 500 m MODIS EVI. Appl. Geogr. 2018, 90, 187–199. [Google Scholar] [CrossRef] [Green Version]
  24. Chandola, V.; Hui, D.; Gu, L.; Bhaduri, B.; Vatsavai, R.R. Using Time Series Segmentation for Deriving Vegetation Phenology Indices from MODIS NDVI Data. In Proceedings of the IEEE International Conference on Data Mining Workshops, Sydney, NSW, Australia, 13 December 2010; pp. 202–208. [Google Scholar]
  25. Lebrini, Y.; Boudhar, A.; Hadria, R.; Lionboui, H.; Elmansouri, L.; Arrach, R.; Ceccato, P.; Benabdelouahab, T. Identifying Agricultural Systems Using SVM Classification Approach Based on Phenological Metrics in a Semi-arid Region of Morocco. Earth Syst. Environ. 2019, 3, 277–288. [Google Scholar] [CrossRef]
  26. Wessels, K.J.; Bachoo, A.; Archibald, S. Influence of composite period and date of observation on phenological metrics ex-tracted from MODIS data. In Proceedings of the 33rd International Symposium on Remote Sensing of Environment: Sustaining the Millennium Development Goals, Stresa, Lago Magglore, Italy, 4–8 May 2009. [Google Scholar]
  27. Bachoo, A.; Archibald, S. Influence of Using Date-Specific Values when Extracting Phenological Metrics from 8-day Composite NDVI Data. In Proceedings of the 2007 International Workshop on the Analysis of Multi-temporal Remote Sensing Images, Leuven, Belgium, 18–20 July 2007; pp. 1–4. [Google Scholar] [CrossRef] [Green Version]
  28. Archibald, S.; Scholes, R.J. Leaf green-up in a semi-arid African savanna-separating tree and grass responses to environ-mental cues. J. Veg. Sci. 2007, 18, 583–594. [Google Scholar]
  29. Schwartz, M.D. Phenology: An Integrative Environmental Science; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2003; p. 564. [Google Scholar]
  30. Alcantara, C.; Kuemmerle, T.; Prishchepov, A.V.; Radeloff, V.C. Mapping abandoned agriculture with multi-temporal MODIS satellite data. Remote. Sens. Environ. 2012, 124, 334–347. [Google Scholar] [CrossRef]
  31. Qiu, B.; Lu, D.; Tang, Z.; Song, D.; Zeng, Y.; Wang, Z.; Chen, C.; Chen, N.; Huang, H.; Xu, W. Mapping cropping intensity trends in China during 1982–2013. Appl. Geogr. 2017, 79, 212–222. [Google Scholar] [CrossRef]
  32. Jönsson, P.; Cai, Z.; Melaas, E.; Friedl, M.A.; Eklundh, L. A Method for Robust Estimation of Vegetation Seasonality from Landsat and Sentinel-2 Time Series Data. Remote. Sens. 2018, 10, 635. [Google Scholar] [CrossRef] [Green Version]
  33. Bai, Y.; Yang, Y.; Jiang, H. Intercomparison of AVHRR GIMMS3g, Terra MODIS, and SPOT-VGT NDVI Products over the Mongolian Plateau. Remote. Sens. 2019, 11, 2030. [Google Scholar] [CrossRef] [Green Version]
  34. Lebrini, Y.; Benabdelouahab, T.; Boudhar, A.; Htitiou, A.; Hadria, R.; Lionboui, H. Farming systems monitoring using machine learning and trend analysis methods based on fitted NDVI time series data in a semi-arid region of Morocco. In Proceedings of the Remote Sensing for Agriculture, Ecosystems, and Hydrology XXI, Strasbourg, France, 9–11 September 2019; Volume 11149. [Google Scholar] [CrossRef]
  35. Suepa, T.; Qi, J.; Lawawirojwong, S.; Messina, J.P. Understanding spatio-temporal variation of vegetation phenology and rainfall seasonality in the monsoon Southeast Asia. Environ. Res. 2016, 147, 621–629. [Google Scholar] [CrossRef] [Green Version]
  36. Wang, D.; Wan, B.; Qiu, P.; Su, Y.; Guo, Q.; Wang, R.; Sun, F.; Wu, X. Evaluating the Performance of Sentinel-2, Landsat 8 and Pléiades-1 in Mapping Mangrove Extent and Species. Remote. Sens. 2018, 10, 1468. [Google Scholar] [CrossRef] [Green Version]
  37. Wessel, M.; Brandmeier, M.; Tiede, D. Evaluation of Different Machine Learning Algorithms for Scalable Classification of Tree Types and Tree Species Based on Sentinel-2 Data. Remote. Sens. 2018, 10, 1419. [Google Scholar] [CrossRef] [Green Version]
  38. Yu, L.; Su, J.; Li, C.; Wang, L.; Luo, Z.; Yan, B. Improvement of Moderate Resolution Land Use and Land Cover Classification by Introducing Adjacent Region Features. Remote. Sens. 2018, 10, 414. [Google Scholar] [CrossRef] [Green Version]
  39. Breiman, L. Random Forests. Mach. Learn. 2001, 45, 5–32. [Google Scholar] [CrossRef] [Green Version]
  40. Jin, Y.; Sung, S.; Lee, D.K.; Biging, G.S.; Jeong, S. Mapping Deforestation in North Korea Using Phenology-Based Multi-Index and Random Forest. Remote. Sens. 2016, 8, 997. [Google Scholar] [CrossRef] [Green Version]
  41. James, G.; Witten, D.; Hastie, T.; Tibshirani, R. An Introduction to Statistical Learning: With Applications in R; Springer: New York, NY, USA, 2014. [Google Scholar]
  42. Ouatiki, H.; Boudhar, A.; Ouhinou, A.; Beljadid, A.; Leblanc, M.J.; Chehbouni, A. Sensitivity and Interdependency Analysis of the HBV Conceptual Model Parameters in a Semi-Arid Mountainous Watershed. Water 2020, 12, 2440. [Google Scholar] [CrossRef]
  43. Boudhar, A.; Ouatiki, H.; Bouamri, H.; Lebrini, Y.; Karaoui, I.; Hssaisoune, M.; Arioua, A.; Benabdelouahab, T. Hydro-logical Response to Snow Cover Changes Using Remote Sensing over the Oum Er Rbia Upstream Basin, Morocco; Rebai, N., Mastere, M., Eds.; Springer: Berlin/Heidelberg, Germany, 2020; pp. 95–102. [Google Scholar]
  44. CHIRPS. Available online: https://www.chc.ucsb.edu/data/chirps (accessed on 13 May 2020).
  45. Jönsson, P.; Eklundh, L. TIMESAT—A program for analyzing time-series of satellite sensor data. Comput. Geosci. 2004, 30, 833–845. [Google Scholar] [CrossRef] [Green Version]
  46. Chen, J.; Jönsson, P.; Tamura, M.; Gu, Z.; Matsushita, B.; Eklundh, L. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. Remote Sens. Environ. 2004, 91, 332–344. [Google Scholar] [CrossRef]
  47. Chen, J.; Wilson, W.S.; Seo, K.-W. Optimized smoothing of Gravity Recovery and Climate Experiment (GRACE) time-variable gravity observations. J. Geophys. Res. Space Phys. 2006, 111. [Google Scholar] [CrossRef] [Green Version]
  48. Geng, L.; Ma, M.; Wang, X.; Yu, W.; Jia, S.; Wang, H. Comparison of Eight Techniques for Reconstructing Multi-Satellite Sensor Time-Series NDVI Data Sets in the Heihe River Basin, China. Remote. Sens. 2014, 6, 2024–2049. [Google Scholar] [CrossRef] [Green Version]
  49. Jönsson, P.; Eklundh, L. Seasonality extraction by function fitting to time-series of satellite sensor data. IEEE Trans. Geosci. Remote. Sens. 2002, 40, 1824–1832. [Google Scholar] [CrossRef]
  50. Fu, Y.; He, H.; Zhao, J.-J.; Larsen, D.; Zhang, H.-Y.; Sunde, M.G.; Duan, S. Climate and Spring Phenology Effects on Autumn Phenology in the Greater Khingan Mountains, Northeastern China. Remote. Sens. 2018, 10, 449. [Google Scholar] [CrossRef] [Green Version]
  51. Gao, F.; Morisette, J.T.; Wolfe, R.E.; Ederer, G.; Pedelty, J.; Masuoka, E.; Myneni, R.; Tan, B.; Nightingale, J. An Algorithm to Produce Temporally and Spatially Continuous MODIS-LAI Time Series. IEEE Geosci. Remote. Sens. Lett. 2008, 5, 60–64. [Google Scholar] [CrossRef]
  52. Eklundh, L.; Jönsson, P. Timesat 3.2 Software Manual; Lund and Malmö University: Lund, Sweden, 2015. [Google Scholar]
  53. Eklundh, L.; Jönsson, P. TIMESAT for Processing Time-Series Data from Satellite Sensors for Land Surface Monitoring. In Radar Remote Sensing of Urban Areas; Springer: Berlin/Heidelberg, Germany, 2016; pp. 177–194. [Google Scholar]
  54. Reed, B.C.; Brown Jesslyn, F.; VanderZee, D.; Loveland, T.R.; Merchant, J.W.; Ohlen, D.O. Measuring phenological varia-bility from satellite imagery. J. Veg. Sci. 1994, 5, 703–714. [Google Scholar] [CrossRef]
  55. Hadria, R.; Benabdelouahab, T.; Mahyou, H.; Balaghi, R.; Bydekerke, L.; El Hairech, T.; Ceccato, P. Relationships between the three components of air temperature and remotely sensed land surface temperature of agricultural areas in Morocco. Int. J. Remote. Sens. 2018, 39, 356–373. [Google Scholar] [CrossRef]
  56. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2017. [Google Scholar]
  57. Hao, P.; Zhan, Y.; Wang, L.; Niu, Z.; Shakir, M. Feature Selection of Time Series MODIS Data for Early Crop Classification Using Random Forest: A Case Study in Kansas, USA. Remote. Sens. 2015, 7, 5347–5369. [Google Scholar] [CrossRef] [Green Version]
  58. Belgiu, M.; Drăguţ, L. Random forest in remote sensing: A review of applications and future directions. ISPRS J. Photogramm. Remote. Sens. 2016, 114, 24–31. [Google Scholar] [CrossRef]
  59. Immitzer, M.; Atzberger, C.; Koukal, T. Tree Species Classification with Random Forest Using Very High Spatial Resolution 8-Band WorldView-2 Satellite Data. Remote. Sens. 2012, 4, 2661–2693. [Google Scholar] [CrossRef] [Green Version]
  60. Meyer, H.; Reudenbach, C.; Hengl, T.; Katurji, M.; Nauss, T. Improving performance of spatio-temporal machine learning models using forward feature selection and target-oriented validation. Environ. Model. Softw. 2018, 101, 1–9. [Google Scholar] [CrossRef]
  61. Kendall, M.G. Rank Correlation Methods; Charles Griffin; Oxford University Press: London, UK, 1975; Volume 35. [Google Scholar]
  62. Davenport, M.L.; Nicholson, S.E. On the relation between rainfall and the Normalized Difference Vegetation Index for diverse vegetation types in East Africa. Int. J. Remote. Sens. 1993, 14, 2369–2389. [Google Scholar] [CrossRef]
  63. Barbosa, H.A.; Kumar, T.L.; Paredes, F.; Elliott, S.; Rejas, J.G. Assessment of Caatinga response to drought using Meteosat-SEVIRI Normalized Difference Vegetation Index (2008–2016). ISPRS J. Photogramm. Remote. Sens. 2019, 148, 235–252. [Google Scholar] [CrossRef]
  64. Lotsch, A.; Friedl, M.A.; Anderson, B.T.; Tucker, C.J. Coupled vegetation-precipitation variability observed from satellite and climate records. Geophys. Res. Lett. 2003, 30. [Google Scholar] [CrossRef]
  65. Vrieling, A.; De Leeuw, J.; Said, M.Y. Length of Growing Period over Africa: Variability and Trends from 30 Years of NDVI Time Series. Remote. Sens. 2013, 5, 982–1000. [Google Scholar] [CrossRef] [Green Version]
  66. Lionboui, H.; Benabdelouahab, T.; Hasib, A.; Boulli, A. Analysis of Farms Performance Using Different Sources of Irrigation Water: A Case Study in a Semi-Arid Area. Int. J. Agric. Manag. Dev. 2016, 6, 145–154. [Google Scholar]
  67. Nitze, I.; Barrett, B.; Cawkwell, F. Temporal optimisation of image acquisition for land cover classification with Random Forest and MODIS time-series. Int. J. Appl. Earth Obs. Geoinf. 2015, 34, 136–146. [Google Scholar] [CrossRef] [Green Version]
  68. Ouatiki, H.; Boudhar, A.; Ouhinou, A.; Arioua, A.; Hssaisoune, M.; Bouamri, H.; Benabdelouahab, T. Trend analysis of rainfall and drought over the Oum Er-Rbia River Basin in Morocco during 1970–2010. Arab. J. Geosci. 2019, 12, 128. [Google Scholar] [CrossRef]
  69. Luo, Z.; Yu, S. Spatiotemporal Variability of Land Surface Phenology in China from 2001–2014. Remote. Sens. 2017, 9, 65. [Google Scholar] [CrossRef] [Green Version]
  70. Yuan, M.; Wang, L.; Lin, A.; Liu, Z.; Qu, S. Variations in land surface phenology and their response to climate change in Yangtze River basin during 1982–2015. Theor. Appl. Clim. 2018, 137, 1659–1674. [Google Scholar] [CrossRef]
Remotesensing 13 00578 g001 550
Figure 1. Map of the study area including the major irrigated areas (Tadla and Doukkala irrigated perimeters in red and orange, respectively) and the reference data localization: (A) fallow, (B) rainfed area (C) irrigated perennial crop and (D) irrigated annual crop. Colored points correspond to data used for training and validation of the Random Forest (RF) model.
Figure 1. Map of the study area including the major irrigated areas (Tadla and Doukkala irrigated perimeters in red and orange, respectively) and the reference data localization: (A) fallow, (B) rainfed area (C) irrigated perennial crop and (D) irrigated annual crop. Colored points correspond to data used for training and validation of the Random Forest (RF) model.
Remotesensing 13 00578 g001
Remotesensing 13 00578 g002 550
Figure 2. Farming systems encountered in the study area: (A) typical normalized difference vegetation index (NDVI) profile, (B) on-ground view, and (C) farming systems illustration.
Figure 2. Farming systems encountered in the study area: (A) typical normalized difference vegetation index (NDVI) profile, (B) on-ground view, and (C) farming systems illustration.
Remotesensing 13 00578 g002
Remotesensing 13 00578 g003 550
Figure 3. Phenological metrics illustration from the NDVI profile for one cropping season.
Figure 3. Phenological metrics illustration from the NDVI profile for one cropping season.
Remotesensing 13 00578 g003
Remotesensing 13 00578 g004 550
Figure 4. Overview of the methodology adopted in this study.
Figure 4. Overview of the methodology adopted in this study.
Remotesensing 13 00578 g004
Remotesensing 13 00578 g005 550
Figure 5. NDVI variability over the four farming systems over the Oum Er-Rbia (OER) basin between 2000 and 2019. NDVI presents the average values extracted from reference data; rainfall data were retrieved from Climate Hazards group Infrared Precipitation with Stations (CHIRPS) [44] over the OER basin.
Figure 5. NDVI variability over the four farming systems over the Oum Er-Rbia (OER) basin between 2000 and 2019. NDVI presents the average values extracted from reference data; rainfall data were retrieved from Climate Hazards group Infrared Precipitation with Stations (CHIRPS) [44] over the OER basin.
Remotesensing 13 00578 g005
Remotesensing 13 00578 g006 550
Figure 6. Bivariate maps showing, simultaneously, the great integral (A), length (B), end (C), and the start (D) of season and their coefficients of variation. High-resolution corresponding figures are provided with the Supplementary Materials.
Figure 6. Bivariate maps showing, simultaneously, the great integral (A), length (B), end (C), and the start (D) of season and their coefficients of variation. High-resolution corresponding figures are provided with the Supplementary Materials.
Remotesensing 13 00578 g006
Remotesensing 13 00578 g007 550
Figure 7. Change map of farming systems (FSs) over the OER basin from 2000 to 2019. The diagonal of the square represents the unchanged classes of farming systems.
Figure 7. Change map of farming systems (FSs) over the OER basin from 2000 to 2019. The diagonal of the square represents the unchanged classes of farming systems.
Remotesensing 13 00578 g007
Remotesensing 13 00578 g008 550
Figure 8. Start of season time (TSOS) trend: Sen’s slope values (A); Mann–Kendall Z (B); and p-value intervals (C).
Figure 8. Start of season time (TSOS) trend: Sen’s slope values (A); Mann–Kendall Z (B); and p-value intervals (C).
Remotesensing 13 00578 g008
Remotesensing 13 00578 g009 550
Figure 9. Spatial distribution of end of cropping season time (TEOS) trend. Maps (A), (B), and (C) represent Sen’s slope values, Kendall Z, and p-value intervals, respectively.
Figure 9. Spatial distribution of end of cropping season time (TEOS) trend. Maps (A), (B), and (C) represent Sen’s slope values, Kendall Z, and p-value intervals, respectively.
Remotesensing 13 00578 g009
Remotesensing 13 00578 g010 550
Figure 10. Spatial distribution of length of season (LOS) trend. Maps (A), (B), and (C) represent Sen’s slope values, Kendall Z, and p-value intervals, respectively.
Figure 10. Spatial distribution of length of season (LOS) trend. Maps (A), (B), and (C) represent Sen’s slope values, Kendall Z, and p-value intervals, respectively.
Remotesensing 13 00578 g010
Remotesensing 13 00578 g011 550
Figure 11. Spatial distribution of great integral (GINT) trend. Maps (A), (B), and (C) represent Sen’s slope values, Kendall Z, and p-value intervals, respectively.
Figure 11. Spatial distribution of great integral (GINT) trend. Maps (A), (B), and (C) represent Sen’s slope values, Kendall Z, and p-value intervals, respectively.
Remotesensing 13 00578 g011
Table
Table 1. Input setting for NDVI time series processing in TIMESAT used in this study, as recommended by [52].
Table 1. Input setting for NDVI time series processing in TIMESAT used in this study, as recommended by [52].
ParameterDescriptionValue
Spike methodTwo values: 1 for median filter and 2 for decomposition by Loess1
Spike valueDegree of spike removal1.8
Amplitude cutoff valueData with amplitude below this value are masked.0.1
Valid data rangeData range of time series to be processed0–1
Season parameterThe study area contains one cropping season1
Number of envelope iterationsNumber of iterations for envelope adaptation2
Adaptation strengthStrength of the envelope adaptation3
Table
Table 2. Definition of the computed phenological metrics (PhM) [52,53,54].
Table 2. Definition of the computed phenological metrics (PhM) [52,53,54].
Phenological MetricPhenological Definition (for Cropping Season)Unit
(1) Start of season—time (TSOS)Beginning date of photosynthesis activity in the vegetation canopydays
(2) End of season—time (TEOS)End date of photosynthesis activity in the vegetation canopydays
(3) Length of season (LOS)Length of photosynthetic activity during the cropping seasondays
(4) Great integral (GINT)Canopy photosynthetic activity across the entire growing season-
Table
Table 3. Ground truthing data.
Table 3. Ground truthing data.
Farming System ClassNumber of PolygonsTraining Area (Ha)
Irrigated Perennial Crop (IPC)801626.71
Irrigated Annual Crop (IAC)1051889.12
Rainfed Area (RA)1603617.76
Fallow (FA)1252661.06
Total4709794.65
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lebrini, Y.; Boudhar, A.; Laamrani, A.; Htitiou, A.; Lionboui, H.; Salhi, A.; Chehbouni, A.; Benabdelouahab, T. Mapping and Characterization of Phenological Changes over Various Farming Systems in an Arid and Semi-Arid Region Using Multitemporal Moderate Spatial Resolution Data. Remote Sens. 2021, 13, 578. https://doi.org/10.3390/rs13040578

AMA Style

Lebrini Y, Boudhar A, Laamrani A, Htitiou A, Lionboui H, Salhi A, Chehbouni A, Benabdelouahab T. Mapping and Characterization of Phenological Changes over Various Farming Systems in an Arid and Semi-Arid Region Using Multitemporal Moderate Spatial Resolution Data. Remote Sensing. 2021; 13(4):578. https://doi.org/10.3390/rs13040578

Chicago/Turabian Style

Lebrini, Youssef, Abdelghani Boudhar, Ahmed Laamrani, Abdelaziz Htitiou, Hayat Lionboui, Adil Salhi, Abdelghani Chehbouni, and Tarik Benabdelouahab. 2021. "Mapping and Characterization of Phenological Changes over Various Farming Systems in an Arid and Semi-Arid Region Using Multitemporal Moderate Spatial Resolution Data" Remote Sensing 13, no. 4: 578. https://doi.org/10.3390/rs13040578

APA Style

Lebrini, Y., Boudhar, A., Laamrani, A., Htitiou, A., Lionboui, H., Salhi, A., Chehbouni, A., & Benabdelouahab, T. (2021). Mapping and Characterization of Phenological Changes over Various Farming Systems in an Arid and Semi-Arid Region Using Multitemporal Moderate Spatial Resolution Data. Remote Sensing, 13(4), 578. https://doi.org/10.3390/rs13040578

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