Characterising the Land Surface Phenology of Middle Eastern Countries Using Moderate Resolution Landsat Data

Abstract

Global change impacts including climate change, increased CO₂ and nitrogen deposition can be determined through a more precise characterisation of Land Surface Phenology (LSP) parameters. In addition, accurate estimation of LSP dates is being increasingly used in applications such as mapping vegetation types, yield forecasting, and irrigation management. However, there has not been any attempt to characterise Middle East vegetation phenology at the fine spatial resolution appropriate for such applications. Remote-sensing based approaches have proved to be a useful tool in such regions since access is restricted in some areas due to security issues and their inter-annual vegetation phenology parameters vary considerably because of high uncertainty in rainfall. This study aims to establish for the first time a comprehensive characterisation of the vegetation phenological characteristics of the major vegetation types in the Middle East at a fine spatial resolution of 30 m using Landsat Normalized Difference Vegetation Index (NDVI) time series data over a temporal range of 20 years (2000–2020). Overall, a progressive pattern in phenophases was observed from low to high latitude. The earliest start of the season was concentrated in the central and east of the region associated mainly with grassland and cultivated land, while the significantly delayed end of the season was mainly distributed in northern Turkey and Iran corresponding to the forest, resulting in the prolonged length of the season in the study area. There was a significant positive correlation between LSP parameters and latitude, which indicates a delay in the start of the season of 4.83 days (R² = 0.86, p < 0.001) and a delay in the end of the season of 6.54 days (R² = 0.83, p < 0.001) per degree of latitude increase. In addition, we have discussed the advantages of fine resolution LSP parameters over the available coarse datasets and showed how such outputs can improve many applications in the region. This study shows the potential of Landsat data to quantify the LSP of major land cover types in heterogeneous landscapes of the Middle East which enhances our understanding of the spatial-temporal dynamics of vegetation dynamics in arid and semi-arid settings in the world.

1. Introduction

Vegetation phenology, the annual life cycle of plant species that occurs through specific biological phases (cycles) under the influence of environmental factors as seasons proceed [1], is considered one of the global climate change indicators [2,3,4,5,6]. The biological phases include, but are not limited to, spring phenophase (i.e., the onset of flowering associated with ‘greenness’) and autumn phenophase (leaf colouring associated with senescence’) [1]. Local and global environmental factors particularly climatic conditions, for example, temperature and precipitation, depending on the spatial distributions and types of the vegetation, are primary controls of vegetation phenology [6,7]. As such, in the temperate eco-region’s temperature is one of the key climate conditions that influences the phenological cycles of vegetation [8,9]. On the other hand, in arid and semi-arid regions (e.g., the Middle East), precipitation is one of the key determinants regulating the phenological cycles of plants [10].

Apart from the significant role of plant phenology in tracking climate change, it also plays a key role in monitoring and managing agricultural practices and understanding ecosystem processes (e.g., functions, gas and energy fluxes, structure, and heath) [2,11,12,13,14]. As such, the dynamics of vegetation phenology (i.e., spatial pattern variability) over spaces and time provide a unique opportunity for estimating the changes in the ecosystem processes and functions based on which current and future investigations can be assessed.

Phenological studies are conducted using various observation approaches, namely, ground-observed phenology (e.g., traditional human visual surveillance via established phenology networks [15,16], near-surface observations via fixed phenology cameras (PhenoCams) [17,18], drones (unmanned aerial vehicles (UAVs) [19,20], and satellite-estimated phenology [21,22]). Measurements obtained from satellite sensors are better known as Land Surface phenology (LSP) because the field view of the sensors may capture a mixture of spectral profiles (phenological cycles) from a variety of vegetation types [23]. In other words, LSP measures the dynamics of cyclic patterns of land surface vegetation. Ground and near surface approaches (i.e., human surveillance, PhenoCams, and drones) are rather limited to a certain range of applications, with inconsistencies in their spatial and temporal coverage [24]. Satellite-estimated measurements from polar-orbiting and geostationary satellite sensors provide consistent datasets (products) on the dynamics of LSP at larger spatial and temporal extents [25,26]. In the last decades, various satellite sensor-based products have been developed for estimating and monitoring LSP across a wide range of landscapes and regions (e.g., Europe). The most widely used sensors include (but are not limited to) Landsat [27,28], Advanced Very High-Resolution Radiometer (AVHRR) [29,30], Moderate Resolution Imaging Spectroradiometer (MODIS) [31,32], Sentinel-2 [33,34], MERIS (Medium Resolution Imaging Spectrometer) [35,36,37], and Satellite Pour l’Observation de la Terre (SPOT) VGT (Vegetation) [38,39,40]. These sensor products, depending on their spatio-temporal characteristics offer various opportunities and challenges in LSP explorations. For example, products from AVHRR and MODIS sensors provide high temporal resolutions (i.e., capable of capturing changes in plant development) but are insufficient to capture specific phenological cycles (characteristics) due to their low spatial resolutions [41]. Alternatively, medium and high spatial resolution sensor datasets with consistent multi-temporal coverage (such as Landsat), particularly in areas where cloud cover is not an issue (e.g., Middle East), could provide a unique opportunity in LSP investigation.

Satellite-estimated measurements in LSP studies usually undergo a transformation from normal spectral bands into so-called vegetation indices (VIs), based on which phenological cycles (or chemometrics, e.g., the onset of greenness, season length, and end of senescence) can be characterised [42,43]. So far, the most common indices used in LSP studies include Normalized Difference Vegetation Index (NDVI) [44], EVI (Enhanced Vegetation Index) [45], Difference Vegetation Index (DVI) [46], and Terrestrial Chlorophyll Index (MTCI) [47], and Normalized Difference Phenology Index (NDPI) [48].

In the Middle East (ME) region as a whole (including seventeen countries), the vegetation and other land cover types are under immense pressure from both natural and anthropogenic drivers, for example, extreme drought episodes [49,50,51,52] and political instabilities (war, social and economic unrest in multiple countries, e.g., Yemen, Iraq, Syria) [53]. A recent example can be seen in Egypt where population growth has already caused environmental disturbances by exhausting water reservoirs for food production and agricultural practices [54]. In addition, in the ME until the end of the 20th century, rainfed agricultural land areas have decreased by more than 170,000 km² due to climate change (from drought and infrequent precipitation) [55]. The negative impact of these drivers could exacerbate in the wake of climate change. Sufficient information on spatiotemporal variation in vegetation dynamics and phenology can be crucial to surrogate the impact of natural and man-made disasters.

In the ME as a whole, studies on LSP are limited or non-existent and are often concentrated on a single country as case-study investigations. For example, in Iraq despite the long ongoing political, social, and economic instabilities (and this is true for many of the ME countries) several studies have been conducted on spatio-temporal variations of phenological characteristics in relation to elevation [31], forecasting crop (e.g., wheat and barley) productions [56], climate change [57], phenometric-based classification of vegetation [58], and cultivated land change estimation based on crop phenology [59]. Similar studies have also been performed within the boundaries of the ME, such as in Iran [60,61,62], Turkey [63,64,65,66], Egypt [67,68], Jordan [69,70], Israel [71,72], Saudi Arabia [73], and northern Fertile Crescent [74]. The majority of these studies have employed coarse spatial resolution datasets and they are at the country level. Therefore, at the ME regional level, there is a need for generating detailed phenological maps that depict spatio-temporal dynamics of vegetation changes from moderate-resolution satellite datasets (e.g., Landsat). Finer spatial resolution (30 m) phenological information provides better insight into the local variations of the vegetation characteristics. In addition, detailed LSP characterisation help scientists and policy makers to derive accurate crop type classifications, land-use change measurements, yield prediction and yield gap quantification in the Middle Eastern smallholder farm system. Furthermore, such studies benefit agricultural practices (e.g., applying fertilisers efficiently, and establishing efficient irrigation systems) and further provide useful insights into the influence of climate change on the dynamics of LSP in the ME region.

2. Data and Methods

2.1. Study Area

The term “Middle East” (ME) as used in this paper refers to the areas covered by 17 countries (Figure 1a). According to World Atlas, ME consists of 17 countries (Iran, Iraq, Bahrain, Egypt, Cyprus, Jordan, Israel, Kuwait, Oman, Lebanon, Palestine, Saudi Arabia, Qatar, Turkey, The Syrian Arab Republic, Yemen, and the United Arab Emirates) [75]. The ME is a predominantly semiarid region that contains a strong north to south precipitation gradient [76]. Annual precipitation may reach more than 1000 mm in humid regions of Turkey and Transcaucasia, while the deserts south of the Euphrates River may receive up to 100 mm/year [76]. According to the recent GlobeLand30, cultivated land, grassland, shrubland, and forest are occupying 11.11%, 9.47%, 2.69%, and 2.69% of ME land area in 2020 (Figure 1b) [77]. Bare land covers by far the largest proportion (72%), while around 1.05% of the ME area is covered by artificial surfaces (Figure 1b) [77]. The contributions of other landcover types such as Wetland, Waterbody, and Permanent snow are less than 1% (Figure 1b) [77].

2.2. Dataset

2.2.1. GlobeLand30

To represent the land cover of ME, the GlobeLand30, which is a 30 m resolution global land cover data product, was developed by the National Geomatics Center of China, and was obtained from the GlobeLand30 website [77]. The 2010 version was released in 2014 and its 2020 version was released recently [77,79]. To cover the entire ME, 48 tiles were downloaded and mosaiced. The GlobeLand30 datasets comprise 10 land cover classes in total, including grassland, cultivated land, shrubland, forest, wetland, water bodies, artificial surface, tundra, ice, perennial snow, and bare land [79].

2.2.2. Landsat Land Surface Reflectance Data

We use country boundaries provided by DIVA-GIS [78] to define the area which is later segmented into 22 processing tiles (see Figure 2) in order to comply with the batch processing quota defined by the Google Earth Engine which is the computing platform employed for data acquisition and processing [80]. For each tile, we acquired surface reflectance level2 images from Landsat 5 and Landsat 7 with timestamps ranging from 1999 to 2021 in order to extract phenology parameters for the years 2000 to 2020. The combination of both Landsat 5 and 7 increases the probability of obtaining cloud-free data covering the area of interest [81].

Landsat 5 carried Thematic Mapper (TM) to capture the Earth’s surface reflectance in visible and infrared bands covering 0.45 to 2.35-micron wavelengths [82]. This satellite operated effectively from March 1984 to May 2012 collecting images with 30 m resolution, 185 km swath width and an average revisit time of 16 days. Carrying the improved instrument named Enhanced Thematic Mapper (ETM+), Landsat 7 has been in operation since May 1999. The responsivities of both TM and ETM+ are consistent with each with a difference of less than 2% in red and near-infrared bands [83]. Thus, data from those instruments can be utilised in tandem without any adjustment. However, due to the permanent failure of the Scan-Line Correlator (SLC) since 2003, defects are present in the Landsat 7 data products [84]. There are gaps in the images (22% unread area per scene) produced by Landsat 7 while some pixels have invalid values due to duplicate measurements. Some algorithms have been developed to overcome the problem (e.g., [85]), but we do not implement any spatial gap-filling algorithm to correct the defects. Images from Landsat 5 are expected to be the complement. Instead, we identified and corrected invalid pixel values due to duplicate measurements based on the variation in the temporal domain.

Time variability of the normalized difference vegetation index (NDVI, [86]) becomes the key variable to extracting the phenology in this study. This index is calculated using the following formula:

$${\mathrm{NDVI}}_{\mathrm{raw}}=\frac{\mathrm{B}4-\mathrm{B}3}{\mathrm{B}4+\mathrm{B}3}$$

(1)

where B3 and B4 are images/data acquired in red (0.63–0.69 micron) and near-infrared (0.77–0.90 micron) bands. We use surface reflectance data which were atmospherically corrected using LEDAPS (Landsat Ecosystem Disturbance Adaptive Processing System) software [87]. Quality assessment bands wherein (pixel_qa band) include cloud and shadow masks produced using the CFMASK algorithm [88] were also considered to mask out the pixels covered by high confidence clouds or cloud shadows. We masked pixels which are flagged as cloud shadow (bit 3) and cloud (bit 5) with high confidence indication (bit 7).

2.2.3. Phenology Extraction

For every year, we collected and processed around 2-year long data ranging from July year y − 1 to February year y + 1. By using a wider time window, more data points can be used in the process to establish the time series and to capture growing seasons that cross a calendar year. The process itself consists of three main stages, namely, outlier correction, harmonic extraction, and the extraction of phenological parameters.

The first stage aims to correct the unreasonably high NDVI mainly caused by the SLC failure in Landsat 7 products. The outliers are identified as data points that deviate more than 0.2 from the moving average values calculated along with a 40-day window. Such outliers are replaced with the moving average value. When the time series data is sufficiently dense, a moving average can be obtained such that this approach successfully corrects the outliers. Conversely, outliers in sparse time series data remain uncorrected. To relieve this issue, we relied on the two-year persistent pattern to capture the variability of NDVI. In other words, two-year data were averaged to identify the seasonal variability in a specific year.

After outlier correction, we fitted the following harmonic function [89,90]:

$${\mathrm{NDVI}}_{\mathrm{smoothed}}=\mathrm{f}\left(\mathrm{t}\right)=\mathrm{A}+\mathrm{Bt}+{\displaystyle \sum }_{\mathrm{n}=1}^{\mathrm{n}={\mathrm{n}}_{\mathrm{H}}}{\mathrm{C}}_{\mathrm{n}}\sin \left(2\mathsf{\pi }\frac{\mathrm{nt}}{\mathrm{T}}\right)+{\displaystyle \sum }_{\mathrm{n}=1}^{\mathrm{n}={\mathrm{n}}_{\mathrm{H}}}{\mathrm{D}}_{\mathrm{n}}\cos \left(2\mathsf{\pi }\frac{\mathrm{nt}}{\mathrm{T}}\right)$$

(2)

where t is time expressed in year fraction, T is the harmonic period, while the A, B, C_n, and D_n are the regression coefficients. Basically, this method is similar to the Fourier analysis where periodic signals are identified through Fourier transform and the smoothed signal is generated by inverse transform after eliminating high-frequency variations. In Google Earth Engine, the fitting was performed using multivariate linear regression with variables associated with t, namely t, $\sin \left(2\mathsf{\pi }\frac{\mathrm{t}}{\mathrm{T}}\right), \cos \left(2\mathsf{\pi }\frac{\mathrm{t}}{\mathrm{T}}\right), \sin \left(4\mathsf{\pi }\frac{\mathrm{t}}{\mathrm{T}}\right),$ etc. To capture the variability more accurately, we used T = 2 years and the number of harmonics n_H = 6. In principle, this parameter choice enables us to identify seasons with an overall duration of four months or more while interannual variability (e.g., growing in one year, no growing in the next year) can still be accounted for. We implemented an upper envelope fitting scheme in order to correct the sudden drops in NDVI due to cloud cover or some other factors. In this scheme, we applied more weights to data points with higher NDVI during the iterative least square regression of the harmonic function. In the first iteration, regression was performed to the raw NDVI, and the drops were replaced with the values from the regression solution. These adjusted values were then used for the next iteration. Three iterations were performed to obtain a reasonable upper fit. Examples of the fitting results are presented in Figure 3. The implementation of upper envelope fitting successfully corrects the drops, but overestimation of local minima becomes an inevitable consequence, especially while dealing with the narrow valley of NDVI between two seasons as depicted in the bottom panel of Figure 3.

From Figure 3, It is evident that the intervals between two consecutive data points are not homogeneous such that identification of phenological parameters can be problematic. To overcome this issue, the fitted harmonic function was calculated for a list of dates with 10-day intervals starting from the start of July year y − 1 to February y + 1. These limits were selected according to the typical crop calendar in the ME countries. The temporal resolution of 10 days is approximately half of the average revisit time of Landsat 5 and 7. Derivative-based phenology extraction (see [2] for review) was performed to this considerably homogeneous time series data. Time derivative, which is the difference of NDVI on two consecutive dates, was computed and the local minima (valley) and maxima (peak) were identified according to the sign change of the derivatives. This method is also known as the inflection point method [35] where the valley is identified as the point where the derivative changes from negative to positive. Conversely, the peak is the inflection point where the derivative changes from positive to negative. Valid start and end of the season should be valleys with NDVI less than 80% of the associated peak while the valid peak should be more than 60% of the global maximum. Peaks with a maximum NDVI of less than 0.2 were ignored. The start, peak, and end of the season for a specific year were recorded where peak time should be in that year while the start and end can be outside. If there is more than one season observed in the data, then we also recorded the phenology for the second season. The maximum value of NDVI in a specific season was also recorded. The outputs for these processes are stored as 8-bit raster images with specific multiplication factors. Covering an area of about 10 million square kilometres, we produced at least five 10-gigapixel rasters for each season (see

Table 1. Extracted phenological parameters.

Parameter	Description
Start of season (SOS)	Valid valley before the peak of season
Peak of season (tmax)	Time at which NDVI reaches maximum in a single year
End of season (EOS)	Valid valley after the peak of season
Length of season (LOS)	Time difference between EOS and SOS
Maximum NDVI	Maximum value of NDVI between SOS and EOS

) each year from 2000 to 2020 summing to ~200 Gigabyte’s data. For the purpose of this study, we have only employed the major phenological parameters including the start of the season (SOS), end of the season (EOS) and length of the season (LOS).

It is noteworthy that the coverages of Landsat 5 and 7 are not homogeneous and the SLC failure leaves stripes in the NDVI maps and subsequent outputs. Some errors/outliers were corrected in the above-mentioned processes, but some others remain. To alleviate this, we masked out pixels with incomplete seasons (e.g., pixels without EOS) and pixels with intermittent phenology detection (e.g., there are only detections for four years or less). Such erroneous detection is likely caused by uncorrected pixel values.

2.2.4. Statistical Analysis

To demonstrate the significant impact of the Latitudinal gradient on LSP parameters, the majority of LSP parameters (SOS and EOS) were estimated per degree increase in latitude between 15° to 45°. The simple linear regressions were fitted for different latitude degrees and the slope and coefficient of determination (R²) were calculated as well as the significance of the model.

To demonstrate the inter-annual variability of LSP parameters from 2000 to 2020 over the entire ME, the temporal standard deviation (STD) values for each LSP parameter in each pixel were computed. A large magnitude of STD can highlight areas that have unstable seasons in ME. In addition, to determine the spatial distribution of each of the phenological parameters for major land cover types, including cultivated land, shrubland, grassland and forest, boxplots were generated.

3. Results

3.1. Spatiotemporal Variation in Phenological Parameters

Figure 4 shows the median start of the season (SOS), end of the season (EOS), and length of the season for the study period 2000–2020 across ME. SOS dates are widely distributed from 9 (September) to 16 (April), later SOS dates were seen at higher latitudes and vice versa (Figure 4a). In general, a large amount of variation in the SOS dates across the region was observed (Figure 4a) since the majority of the region directly relies on rainfall to plant or indirectly recharge the surface water. The earliest SOS (September–November) was observed in Egypt, Northwest Syria, West Turkey, middle and south of Iraq, north and west of Iran, north Jordan, Palestine, and Lebanon. These areas are predominantly irrigated cropland. However, in areas associated mainly with rainfed croplands, slightly later SOS (November–December) was observed (some parts of Iraq, Iran, and Syria). The SOS of the majority of shrubland cover types in the west of Iran seems to range between January and February, while in Yemen shrubland cover type has slightly late SOS dates (February–March). Grassland and forest cover types are mainly located in high altitude parts of ME (Turkey, north of Iraq, and north of Iran) and their SOS dates fall approximately between January and mid-March. For some areas in Yemen and other parts of ME, the SOS was delayed to May–July. This is an indication of growing summer crops in these areas.

The EOS dates range from 13 (January) to 24 (December) and their spatial variation is smaller than SOS (Figure 4b). This might be due to a lack of moisture and facing a hot summer across the region. The earlier SOS is generally a coincidence with an earlier EOS and vice versa. Areas with a later EOS exhibit a longer growing season than areas with an earlier SOS (Figure 4b). However, this might not be applicable to the majority of Grassland cover types as it has slightly late SOS and EOS compared to cultivated land, but its growing season is shorter.

LOS result is presented in Figure 4c. Shorter growing seasons (approximately 3 to 4 months) are associated with shrublands in west Yemen and west Iran. Maximum LOS values (approximately 10 months to a year) can be seen in the north and west of Turkey and north of Iran associated with forests.

3.2. Variability in LSP Parameters

Figure 5 shows the standard deviation of LSP estimates to highlight the inter-annual variability of SOS, EOS, and LOS across ME over the 20-year period (Figure 5). In general, greater variability values were observed in SOS compared to EOS. This is because the majority of the ME region is reliant on rainfall for starting the growing season directly in rainfed areas or indirectly to recharge the source of water for irrigation. Some areas in Egypt, Iraq, Syria, Jordan, Palestine, and Iran, mainly cultivated lands, produced small STDs for SOS (<45 days) (Figure 5a). For the same areas, a smaller STD was observed for EOS (<30 days) (Figure 5b) due to a lack of moisture as the area consistently approaches hot summer. However, the largest EOS values up to 60 days were observed in north Turkey and north Iran, mainly associated with forest land cover type. Similar high STD values can be seen in the small proposition of cultivated lands due to mainly heterogeneous agricultural management and crop practising.

3.3. Latitudinal Variation of LSP Parameters

The trends and variations in LSP indicate that LSP is dependent on latitude (Figure 6). As we can see from the results, latitude had more influence on SOS than EOS. Approximately 86% of SOS dates and 83% of EOS dates can be explained by latitude (p < 0.0001). The relationship between SOS and latitude was significant. It indicates a delay in SOS of 4.83 days (R² = 0.86, p < 0.001) per degree of Latitude increase. A significant relationship was also found between EOS and latitude. A degree increase in latitude will result in an approximately 6.54 days delay in EOS dates (R² = 0.83, p < 0.001) (see

Table 2. The coefficient of determination (R²), slopes, and y-intercept of the linear regression between LSP parameters and latitude during.

LSP Parameters	R2	Slop	y-Intercept	p (Sig.)
SOS	0.86	4.83	188.6	<0.0001
EOS	0.83	6.54	325.89	<0.0001

).

3.4. Characterising of the LSP of Major Land Cover Types in ME

The spatiotemporal variation of vegetation phenological patterns in ME is largely driven by different climatic factors (rainfall and temperature). Different patterns were seen in the LSP parameters across the four types of land cover in ME (Figure 7). In general, cultivated lands and grassland present the earliest SOS dates, starting mainly around September, followed by shrublands for which SOS starts from early November. The latest SOS dates were recorded for the forest which start around December. The majority of SOS dates for cultivated lands, grassland, shrublands, and forests are falling into relatively the same window for which SOS dates start around October until mid-January. There is greater variability in SOS dates among the landcover types compared to EOS dates except for forest. The earliest EOS were recorded for shrublands from mid-May until mid-June. Sightly later EOS dates were observed for cultivated lands starting from late May until early July, followed by grasslands for which EOS dates start from June until late July. The latest EOS dates were recorded for the forest for which EOS dates start from mid-June until mid-November. Regarding the LOS, forests recorded the longest growing season (between ~5 to 10 months, followed by cultivated lands (~5.2 to 6.3 months), grasslands (between ~5 to 6.3 months), and the shortest growing season was recorded for shrublands (between ~5 to 5.6 months). The above records were based on the majority of the observation. However, due to the lack of accuracy of the used land cover types in some areas, a possible mixture of different land cover types might have occurred which resulted in inaccurate LSP estimates.

4. Discussion

Accurate characterisation of LSP parameters can be used as a useful tool to investigate global change impacts including climate change, increased CO₂, and nitrogen deposition [91,92]. In addition, the precise estimation of LSP parameters could be crucial for other applications including mapping of vegetation types, yield forecasting and irrigation management [56,58]. The current research was based on the estimation of the date of different LSP parameters using Landsat 5 and 7 time-series data, providing a description of how SOS, EOS, and LOS changes for major land cover types across ME. These results provide a description of Middle Eastern vegetation phenology for the last two decades and help better understand future investigations of vegetation dynamics of various biomes, particularly in relation to climate change. It is worth noting that this is the first attempt to map and characterise vegetation phenology at the ME regional level, using data from Landsat 5 and 7 from 2000 to 2020, at the spatial resolution of 30 m. This work only used Landsat 5 and 7 data despite the availability of Landsat 8 over the study period. One reason for not including Landsat 8 is that some studies have highlighted the differences in the values of average reflectance for each band between the two datasets [93,94]. Therefore, only Landsat 5 and 7 were employed to keep the consistency over the study period. Meanwhile, a detailed comparison between LSP parameters derived from Landsat 5 and 7 and Landsat 8 in arid and semi-arid regions will be investigated in future work.

Although, coarse resolution LSP estimates such as MODIS are available at the global scale [95], our finer resolution LSP estimates can offer much more granular information which can be used for many policies making decisions. In the Middle Eastern complex smallholder landscape, coarse-resolution data may not be able to capture the local variation among farmlands (Figure 8(a1) vs. Figure 8(a2)), whereas these variations were captured well in the Landsat based estimates (Figure 8(a2)). In addition, individual MODIS-500 m pixel covers several agricultural fields, while 30 spatial resolution datasets can better represent small agriculture field sizes in ME (Figure 8(a3)). Furthermore, considerable areas have been left without LSP estimates in the coarse-resolution dataset due to insufficient spatial resolution compared to medium-resolution data (Figure 8(b1,b2)). Finally, besides its significant improvement compared to coarse LSP estimates, our outputs can be vital to enhance many applications in the region, including (i) crop type classification (Figure 8(a2)), (ii) urban expansion over time and urban phenology (Figure 8(c1,c2)), (iii) land-use change, (iv) monitoring crop rotation and map cropland fallows which are widely practised in some areas in ME (Figure 8(d2)), and (v) crop yield estimation and yield gap modelling at the farm level.

In arid to semi-arid regions, spatial and temporal variations in LSP parameters are largely driven by changes in water availability, climate condition, vegetation composition, and soil type [96,97,98,99]. Using traditional approaches to map vegetation phenology at the regional scale is a great challenge due to their inconsistencies, lack of adequate spatial information, cost-inefficiency, and being time-consuming [100,101]. However, remote sensing monitoring systems are faster, can provide information over a large spatial scale on a systematic basis, and can monitor their changes over time at a low cost [21,100,102,103]. In addition, these remotely sensed monitoring systems are highly important in the context of ME where access is restricted in some areas due to insecurity issues. This factor among other limitations has made the validation of LSP parameters derived from remotely sensed data in this region even harder. In addition, based on our knowledge, ground phenological stations do not exist across the region. Furthermore, there has not been any attempt to characterise ME’s vegetation phenology at fine spatial resolution and, generally, a very limited number of studies are available at the local scale.

The overall findings from this study are in line with the expected vegetation phenology patterns across ME. However, there were difficulties and limitations that might have impacted the overall accuracies of estimating the LSP parameters. This study employed the GlobeLand30 for 2020 [77] to discriminate major land cover types, assuming that there were no changes in all major land cover types during the study period. In addition, all the cultivated lands are classified as one class and there is no discrimination between irrigated and rainfed croplands, whereas the previous study showed that irrigated and rainfed croplands had different growing seasons and their phenological patterns are relatively different in Iraq [31]. Regarding the accuracy of the dataset, it has been claimed that the overall accuracy of GlobeLand30 2020 is 85% and the Kappa coefficient is 0.78 [77], however, the product might not be able to provide accurate mapping at the regional level for some land cover types. Although we did not find any studies that assessed the accuracy of the recent version (2020) of the product, the older version of the dataset was assessed in a similar environment. For instance, six study sites around larger cities representing dissimilar eco-regions spanning rural and urban areas were selected in Iran to assess the accuracy of GlobeLand30 and the result showed an overall accuracy of 77% [104]. In central Asia, classification accuracy for GlobeLand30-2010 was quantified using error matrix and Kapa coefficient and the outputs represented an overall accuracy of 46% [105]. The lack of classification accuracy may result in mixing different land cover types which can lead to inaccurate LSP characterisation for individual land cover types presented in this study.

The STD was used to determine the most unstable locations over the last two decades in terms of LSP parameters including SOS, EOS, and LOS. Based on the observation, the most unstable land cover type in terms of LSP parameters over the study period was a forest. This variation was greater in EOS compared to the other LSP parameters. Possible explanations for this result include (i) the occurrence of deforestation and reforestation over the study period (2000–2020) [106], (ii) natural and man-made disasters such as fire [107], (iii) limitation of the method to capture the phenological parameters correctly due to lack of data, and noises that might come from cloud coverage and snow, and (iv) difficulties in the extraction of autumn LSP parameters due to their high uncertainties and greater complexity. In addition, previous studies indicated that vegetation cover and heterogeneity are also vital factors that could cause difficulties in LSP retrieval in arid and semi-arid regions [108]. However, our findings are consistent with the results that are shown in [31]. For instance, The STD of SOS for most of Iraq ranged between 0 and 72 days, while the STD of EOS for the majority of the country ranged between 0 and 32 days [31]. In some cultivated lands, higher STD values were observed due to mainly human interactions such as the application of heterogeneous agricultural management, practising different crop growing systems over the years, and implementation of the traditional biennial system to recharge the depleted soil [109,110].

As it can be seen in Figure 6, the LSP parameters along with latitudinal gradients correlated significantly with each other (p < 0.01). To our knowledge, this is the first study that investigates the latitudinal variation of vegetation phenology in the ME using moderate resolution satellite data. Previous studies at the country scale have evaluated the altitudinal variation of vegetation in Iraq [31] and the study found that an increase in elevation delayed the SOS and EOS. Similar patterns were recorded in our study since elevation generally increases from the south to the north of the region with increasing latitude. The availability of the water at the beginning of the season (irrigated cropland) in the middle and south of the region advanced the SOS dates, while rainfall and temperature constraints as a function of variation in altitude delayed the SOS in the northern parts of the region. In addition, the temperature has long been considered to be the major controlling factor of vegetation phenology at the middle and high altitudes. Other studies have also achieved similar results in a similar environment. For example, in North Africa, Adole et al. [111] reported that a one-degree increase in latitude will result in an approximately 5 days delay in SOS and 5 days advance in EOS dates.

Comparisons of the LSP parameter dates in this work with those from previous studies could improve our understanding of the spatial patterns of vegetation phenology across ME. Although several studies have used different approaches to investigate the vegetation phenology in different parts of ME, most of them are related to phenology trends or changes across time, rather than the characterisation of LSP parameters (73, 72, 63, 70, 65, 2012, 59, 64, 58, 74, 61, 68, 66, 56, 69, 60, 57, 62, 67, 71). Therefore, a direct comparison between the results was not possible, but the overall phenological patterns presented in these studies were consistent with our findings.

5. Conclusions

This study shows the capacity of time series Landsat data to map and characterise the phenology of major land cover types along a latitudinal gradient in ME. For the first time, this research established the most comprehensive characterisation of the LSP of the major land cover types from GlobeLand30 and NDVI derived from Landsat 5 and 7 at a moderate resolution of 30 m over a long temporal range of 20 years (2000–2020). The generated maps of LSP parameters in this study represent the most up to date, the finest spatial resolution and most detailed characterisation of the vegetation phenology in ME. In addition, for the first time, the inter-annual variability of all LSP parameters for all of ME was quantified.

The proposed method performed well for deriving LSP parameters for major land cover types in ME. Different phenology patterns were captured in derived phenology maps despite heterogeneous vegetation types in the study area. Based on the median LSP parameters of the time series from 2000 to 2020, we observed inter-annual variation in vegetation phenology in the study area fluctuated to a different degree. Spatially, the SOS shows a delaying trend from south to north. With an increase in latitude, this trend was more evident (high altitude). A stronger delaying trend in SOS was observed alongside the latitudinal gradient. Temporarily, the inter-annual changes in LSP parameters were quantified and the results can be used to investigate land cover change over the study period in ME.

Future research should be dedicated to validating the LSP parameters when ground data become available. In addition, more research will be needed to understand the major driver of vegetation phenology in ME.