Linear and Non-Linear Vegetation Trend Analysis throughout Iran Using Two Decades of MODIS NDVI Imagery

: Vegetation is the main component of the terrestrial Earth, and it plays an imperative role in carbon cycle regulation and surface water/energy exchange/balance. The coupled effects of climate change and anthropogenic forcing have undoubtfully impacted the vegetation cover in linear/non-linear manners. Considering the essential beneﬁts of vegetation to the environment, it is vital to investigate the vegetation dynamics through spatially and temporally consistent workﬂows. In this regard, remote sensing, especially Normalized Difference Vegetation Index (NDVI), has offered a reliable data source for vegetation monitoring and trend analysis. In this paper, two decades (2000 to 2020) of Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI datasets (MOD13Q1) were used for vegetation trend analysis throughout Iran. First, the per-pixel annual NDVI dataset was prepared using the Google Earth Engine (GEE) by averaging all available NDVI values within the growing season and was then fed into the PolyTrend algorithm for linear/non-linear trend identiﬁcation. In total, nearly 14 million pixels (44% of Iran) were subjected to trend analysis, and the results indicated a higher rate of greening than browning across the country. Regarding the trend types, linear was the dominant trend type with 14%, followed by concealed (11%), cubic (8%), and quadratic (2%), while 9% of the vegetation area remained stable (no trend). Both positive and negative directions were observed in all trend types, with the slope magnitudes ranging between − 0.048 and 0.047 (NDVI units) per year. Later, precipitation and land cover datasets were employed to further investigate the vegetation dynamics. The correlation coefﬁcient between precipitation and vegetation (NDVI) was 0.54 based on all corresponding observations (n = 1785). The comparison between vegetation and precipitation trends revealed matched trend directions in 60% of cases, suggesting the potential impact of precipitation dynamics on vegetation covers. Further incorporation of land cover data showed that grassland areas experienced signiﬁcant dynamics with the highest proportion compared to other vegetation land cover types. Moreover, forest and cropland had the highest positive and negative trend direction proportions. Finally, independent (from trend analysis) sources were used to examine the vegetation dynamics (greening/browning) from other perspectives, conﬁrming Iran’s greening process and agreeing with the trend analysis results. It is believed that the results could support achieving Sustainable Development Goals (SDGs) by serving as an initial stage study for establishing conservation and restoration practices.


Introduction
Vegetation is a key constituent of the Earth's biosphere, which serves as the interlinkage substance between the soil, water, and atmosphere [1,2]. Vegetation significantly influences the terrestrial carbon cycle, surface water/energy exchange and balance, and magnitudes were also computed for each vegetation pixel throughout Iran. MODIS land cover data were used to examine associations between vegetation trends and vegetation types. Furthermore, a land cover change map of the study area was produced based on MODIS land cover maps to investigate the relationship between land cover transitions and vegetation trend types. Additionally, in-situ precipitation observations for the same period were employed to investigate the relationship between vegetation and precipitation dynamics. Finally, vegetation trend types, in-situ precipitation, and MODIS land cover data were also used to examine associations between vegetation trends, vegetation types, and precipitation. The results would benefit both the scientific community and governors in studying the possible drivers of change and identifying hotspot zones (i.e., highly dynamic regions) for further monitoring. It is believed that the results could support achieving Sustainable Development Goals (SDGs) by serving as an initial stage study for establishing conservation and restoration practices.

Study Area
The study area, Iran, spanning between 24 • -40 • N (latitudes) and 44 • -64 • E (longitudes), is located in the Middle East, the western part of Asia (see Figure 1a). Iran, with an area of approximately 1.7 million km 2 , is recognized as the seventeenth largest country in the world. The study area is a mountainous region with two main mountain ranges of Zagros and Alborz. The former begins in the northwestern parts of Iran and follows the western border of the country to the southern parts, and the latter is located along a belt from west to the east in the northern regions of the country that can be seen in the ALOS Global Digital Surface Elevation Model (see Figure 1b) [45]. These two mountain ranges and their surroundings are among the most vegetated regions in Iran, which include dense and sparse forest, rangeland, and farmland areas (see Figure 1c) [46]. The southern border of Iran is shaped by the Persian Gulf and the Oman Sea, and the Caspian Sea is located in the north. The immense latitudinal range of Iran and its geographical characteristics cause high climate variations from very-humid to hyper-arid climate zones and high precipitation variations [47,48]. These variations are responsible for diverse land cover types in Iran, from dense forests (between the Caspian Sea and the Alborz Mountains) to the uncovered plains and deserts, such as Dasht-e-Kavir and Dasht-e-Lut (central parts).
Remote Sens. 2022, 14, x FOR PEER REVIEW 4 of 22 changes. Meanwhile, the significance of trends, trend directions (positive or negative), and slope magnitudes were also computed for each vegetation pixel throughout Iran. MODIS land cover data were used to examine associations between vegetation trends and vegetation types. Furthermore, a land cover change map of the study area was produced based on MODIS land cover maps to investigate the relationship between land cover transitions and vegetation trend types. Additionally, in-situ precipitation observations for the same period were employed to investigate the relationship between vegetation and precipitation dynamics. Finally, vegetation trend types, in-situ precipitation, and MODIS land cover data were also used to examine associations between vegetation trends, vegetation types, and precipitation. The results would benefit both the scientific community and governors in studying the possible drivers of change and identifying hotspot zones (i.e., highly dynamic regions) for further monitoring. It is believed that the results could support achieving Sustainable Development Goals (SDGs) by serving as an initial stage study for establishing conservation and restoration practices.

Study Area
The study area, Iran, spanning between 24°-40° N (latitudes) and 44°-64° E (longitudes), is located in the Middle East, the western part of Asia (see Figure 1a). Iran, with an area of approximately 1.7 million km 2 , is recognized as the seventeenth largest country in the world. The study area is a mountainous region with two main mountain ranges of Zagros and Alborz. The former begins in the northwestern parts of Iran and follows the western border of the country to the southern parts, and the latter is located along a belt from west to the east in the northern regions of the country that can be seen in the ALOS Global Digital Surface Elevation Model (see Figure 1b) [45]. These two mountain ranges and their surroundings are among the most vegetated regions in Iran, which include dense and sparse forest, rangeland, and farmland areas (see Figure 1c) [46]. The southern border of Iran is shaped by the Persian Gulf and the Oman Sea, and the Caspian Sea is located in the north. The immense latitudinal range of Iran and its geographical characteristics cause high climate variations from very-humid to hyper-arid climate zones and high precipitation variations [47,48]. These variations are responsible for diverse land cover types in Iran, from dense forests (between the Caspian Sea and the Alborz Mountains) to the uncovered plains and deserts, such as Dasht-e-Kavir and Dasht-e-Lut (central parts).  [45], and (c) recent land cover map of Iran [46].  [45], and (c) recent land cover map of Iran [46].

Material and Methods
In this section, first, the used datasets (see Table 1), including NDVI, in-situ precipitation, and land cover datasets, are introduced. Afterward, the implemented algorithm for vegetation trend analysis is explained.

NDVI Dataset
The MODIS instrument onboard the Terra platform was launched on 18 December 1999 and became operational in early 2000. MODIS daily captures Earth's surface radiance in 36 spectral bands ranging from 0.4 µm to 14.4 µm of the electromagnetic spectrum. Due to its functional characteristics, MODIS imagery has been extensively used in different applications and for generating many global products [49]. Among various products, MODIS vegetation indices are well-approved products by the scientific community and have been widely incorporated into multiple tasks, such as vegetation dynamics monitoring [11,16,17,38,50]. Here, Collection 6 (version 6) of 16-day NDVI images (MOD13Q1) of the Terra-MODIS product was used. In addition to several advantages of NDVI, its lower sensitivity to topographical variations prioritized its utility over EVI [31,33] since the vegetation areas are mainly located near mountainous regions of the study area. Previous Collection of Terra-MODIS NDVI was adversely affected due to sensor degradation; however, the planned adjustments largely eliminated this effect for Collection 6 NDVI [51]. This NDVI dataset is a level 3 product generated every 16 days at the spatial resolution of 250 m. The MODIS algorithm specifies the best available pixel value within each 16-day interval based on the criteria of low clouds, low view angle (to minimize the BRDF effects), and the highest NDVI value, providing a high-quality NDVI dataset [32,52]. This product and its associated uncertainties have been well quantified through different approaches and achieved the stage 3 validation process [53]. In particular, this product has been validated using AERONET-corrected data, spaceborne and airborne sensors, as well as radiometric field measurements. Based on the validation analysis, the accuracy of NDVI values is within ±0.025, demonstrating high-quality assurance of the MODIS 16-day NDVI dataset [54]. In total, 168 bi-monthly (16-day) NDVI images during the growing season (i.e., 21 April to 21 August) between 2000 and 2020 were employed in the present study. The growing season was determined based on the knowledge of the study area and the commonly considered months in previous relevant studies [10,55,56]. Later, all NDVI images within each year's growing season were averaged to produce yearly NDVI images for trend analysis. It is worth noting that the NDVI dataset was prepared using the Google Earth Engine (GEE) cloud computing platform [57,58] with the snippet of (ee.ImageCollection("MODIS/006/MOD13Q1")). Although the MODIS NDVI dataset includes high-quality NDVI values, the yearly aggregation procedure is expected to reduce the random changes in NDVI values, the possibility of outlier influence, and the temporal autocorrelation in time-series for trend analysis [10,34,37]. Figure 2 shows the per pixel vegetation condition of Iran based on the aggregated yearly NDVI images in 2000 (start), 2020 (end), multi-year NDVI mean, and the NDVI range (i.e., the difference between the largest and lowest NDVI values) over the last two decades. As evident in Figure 2a,b, vegetated areas are mainly located in northern and western parts, along the Zagros and Alborz Mountains, as well as in central southern parts of Iran.

Land Cover Dataset
Different vegetation types may behave multiformly through passing time and thus, might be considered as an impactful criterion for trend analysis. Accordingly, two MODIS land cover maps (MCD12Q1) of 2001 and 2020 with 500 m spatial resolution (relatively closer to the MODIS NDVI dataset) were employed for three purposes. The land cover dataset was downloaded from the GEE data catalog with the snippet of (ee.ImageCollection("MODIS/006/MCD12Q1")). First, the 2001 land cover map was used to explore the dominant land cover type for each linear and non-linear trend pattern (see Section 3.4). Second, a land cover change map was generated by comparing these maps to identify land cover transitions, which then was incorporated to investigate the relationship between land cover transitions and trend patterns (see Section 3.4). Third, it was used to interpret the possible three-way relationship between vegetation trend direction (positive or negative), precipitation trend direction, and land cover types. Different land cover products have been developed using MODIS imagery, and in this study, the land cover product generated based on the International Geosphere-Biosphere Programme (IGBP) scheme was used. This product was generated using nadir Bidirectional Reflectance Distribution Function (BRDF) adjusted surface reflectance data and Random Forest classifier [59]. It contained 17 land cover classes, 12 of which classes are different vegetation covers, and the rest are permanent wetland, urban/built-up, permanent snow/ice, barren, and water bodies.

Land Cover Dataset
Different vegetation types may behave multiformly through passing time and thus, might be considered as an impactful criterion for trend analysis. Accordingly, two MODIS land cover maps (MCD12Q1) of 2001 and 2020 with 500 m spatial resolution (relatively closer to the MODIS NDVI dataset) were employed for three purposes. The land cover dataset was downloaded from the GEE data catalog with the snippet of (ee.ImageCollection ("MODIS/006/MCD12Q1")). First, the 2001 land cover map was used to explore the dominant land cover type for each linear and non-linear trend pattern (see Section 3.4). Second, a land cover change map was generated by comparing these maps to identify land cover transitions, which then was incorporated to investigate the relationship between land cover transitions and trend patterns (see Section 3.4). Third, it was used to interpret the possible three-way relationship between vegetation trend direction (positive or negative), precipitation trend direction, and land cover types. Different land cover products have been developed using MODIS imagery, and in this study, the land cover product generated based on the International Geosphere-Biosphere Programme (IGBP) scheme was used. This product was generated using nadir Bidirectional Reflectance Distribution Function (BRDF) adjusted surface reflectance data and Random Forest classifier [59]. It contained 17 land cover classes, 12 of which classes are different vegetation covers, and the rest are permanent wetland, urban/built-up, permanent snow/ice, barren, and water bodies.

In-Situ Precipitation Dataset
Many studies have found supporting evidence for the impact of climate change on terrestrial vegetation cover [60][61][62][63][64]. Among climatic variables, precipitation, especially in semi-arid and closed climate zones, has been recognized as the most prominent factor in Remote Sens. 2022, 14, 3683 7 of 23 vegetation cover dynamics due to its dominant role in the water supply [65,66]. Accordingly, accurate precipitation observations from synoptic stations were employed to investigate the relationship between vegetation and precipitation inter-annual dynamics. In this regard, 85 synoptic stations (see Figure 1a) that were located on vegetation pixels (based on MODIS NDVI), selected from a total of 244 synoptic stations, were acquired from the Islamic Republic of Iran Meteorological Organization (IRIMO). These synoptic stations contained daily precipitation observations from 2000 to 2020. In particular, only synoptic stations located on vegetation pixels with 21 years of consistent precipitation observations were employed for further analysis. Furthermore, it can be seen in Figure 1a that synoptic stations are reasonably distributed across vegetation areas in Iran. The daily precipitation values during the growing season (i.e., similar to the NDVI dataset) were accumulated to generate a 21-year precipitation dataset for each synoptic station. Subsequently, they were utilized for trend analysis to explore their trend direction (i.e., positive or negative) for further comparisons with vegetation trend results.

Linear and Non-Linear Trend Analysis
It has been recognized that the vegetation covers do not always follow linear dynamics and may change in non-linear manners instead [34,[42][43][44]67,68]. Accordingly, this paper implemented an automated vegetation trend analysis framework developed by Jamali et al. [34] that considers both linear and non-linear trend patterns. This algorithm, named PolyTrend, works with the per-pixel time-series NDVI values and assigns one of the (1) Cubic, (2) Quadratic, (3) Concealed, (4) Linear, (5) No Trend labels to the examined pixel. The concealed label is assigned to pixels with either cubic or quadratic trends that has insignificant net change during the study period (i.e., negligible NDVI difference between the start and end dates). PolyTrend is a 3-stage workflow (see Figure 3), which is initiated by checking the existence of a cubic (ax 3 + bx 2 + cx + d) trend. In this regard, first (stage 1), a cubic model is fitted to the input time-series NDVI values, followed by checking the statistical significance of the cubic coefficient (a) at the 95% confidence interval using the t-test statistics. In case of passing the significance test, the presence of both local maximum and minimum of the time-series in the fitted model is checked (comparing the roots of the fitted model with time-series NDVI values), and if approved, a linear model is fitted to the examined time-series NDVI values. The significance of the linear coefficient will lead to assigning the cubic label to the corresponding pixel; otherwise (i.e., the insignificance of the linear coefficient), the Concealed label is assigned. In case of rejecting both criteria, the algorithm will test the existence of a quadratic trend (bx 2 + cx + d) at the second stage. The procedure is advanced as in the cubic workflow, except for checking the presence of only one local maximum or minimum in the fitted model. If the considered criteria are not met, the possibility of a linear trend for the time-series NDVI values is explored (stage 3). To this end, a linear model is fitted to the time series, and if the model is found to be statistically insignificant, the No Trend label is assigned, and in case of being statistically significant, the corresponding pixel is identified with the Linear pattern. Further details of the PolyTrend are available at Jamali et al. [34]. Here, time-series NDVI values (see Section 3.1) of each pixel were inserted into the PolyTrend. It is worth noting that only pixels with at least one NDVI value of over 0.15 were considered as vegetation cover and employed in trend analysis [69][70][71]. In other words, non-vegetation and very sparsely vegetated pixels with low NDVI values < 0.15 during the study period (21 years) were excluded to reduce the computational burden of trend analysis due to the massive extent of Iran. In addition to trend types, the significance of trends, direction of trends, and the slope magnitude of trends were also computed for each vegetation pixel. These parameters would provide complementary information to the actual trend analysis results for profound investigations. For instance, vegetation pixels with statistically significant trends and high slope magnitude could present the hotspot locations that have been altered considerably. parameters would provide complementary information to the actual trend analysis results for profound investigations. For instance, vegetation pixels with statistically significant trends and high slope magnitude could present the hotspot locations that have been altered considerably.
It is worth noting that PolyTrend was also applied to in-situ precipitation data to determine their trend directions (i.e., positive or negative). The derived trend directions were then used to compare the vegetation and precipitation mutual dynamics.

Results
This section covers the PolyTrend assessment, followed by providing trend analysis results in Iran over the last two decades. Finally, other data sources are included to investigate the associations between vegetation trends, precipitation trends, and land cover types.

Trend Linearity/Non-Linearity Assessment
The PolyTrend has been statistically tested using both synthesized and actual timeseries NDVI data [34]. Furthermore, it has been employed for trend analysis in other studies and is claimed to be accurate in trend pattern identification [72][73][74][75][76]. However, in the first step, the results of the PolyTrend were evaluated through visual inspections for quality assurance of the detected vegetation trends in Iran. To this end, 100 candidate pixels were randomly (stratified random) selected throughout the study area, containing five different classes (see Section 3.4). Later, the trends detected by the PolyTrend for each individual pixel were visually and statistically examined. In general, the assessment step suggested the applicability of PolyTrend for vegetation trend analysis using two decades of MODIS NDVI imagery in Iran. Figure 4 shows the identified trends for 10 pixels derived from the candidate pixels for justification, including two samples for each trend class. Based on Figure 4, the PolyTrend correctly identified trend types with a 95% confidence interval. For instance, it properly identified the cubic and quadratic trends ( Figure  4a-d), and assigned concealed trends to pixels with possible cubic and quadratic trends (Figure 4e,f), while their linear coefficients were insignificant (i.e., low net NDVI difference). Furthermore, the linear labels were given to pixels in which the required criteria for cubic and quadratic trends were not satisfied, but the linear coefficient was significant (Figure 4g,h). Finally, when no significant trend was observed for the input time series, the no trend label was assigned (Figure 4i,j). These pixels contained fluctuating NDVI values within very low boundaries (i.e., low variation) and ended without notable NDVI changes. It is worth noting that PolyTrend was also applied to in-situ precipitation data to determine their trend directions (i.e., positive or negative). The derived trend directions were then used to compare the vegetation and precipitation mutual dynamics.

Results
This section covers the PolyTrend assessment, followed by providing trend analysis results in Iran over the last two decades. Finally, other data sources are included to investigate the associations between vegetation trends, precipitation trends, and land cover types.

Trend Linearity/Non-Linearity Assessment
The PolyTrend has been statistically tested using both synthesized and actual timeseries NDVI data [34]. Furthermore, it has been employed for trend analysis in other studies and is claimed to be accurate in trend pattern identification [72][73][74][75][76]. However, in the first step, the results of the PolyTrend were evaluated through visual inspections for quality assurance of the detected vegetation trends in Iran. To this end, 100 candidate pixels were randomly (stratified random) selected throughout the study area, containing five different classes (see Section 3.4). Later, the trends detected by the PolyTrend for each individual pixel were visually and statistically examined. In general, the assessment step suggested the applicability of PolyTrend for vegetation trend analysis using two decades of MODIS NDVI imagery in Iran. Figure 4 shows the identified trends for 10 pixels derived from the candidate pixels for justification, including two samples for each trend class. Based on Figure 4, the PolyTrend correctly identified trend types with a 95% confidence interval. For instance, it properly identified the cubic and quadratic trends (Figure 4a-d), and assigned concealed trends to pixels with possible cubic and quadratic trends (Figure 4e,f), while their linear coefficients were insignificant (i.e., low net NDVI difference). Furthermore, the linear labels were given to pixels in which the required criteria for cubic and quadratic trends were not satisfied, but the linear coefficient was significant (Figure 4g,h). Finally, when no significant trend was observed for the input time series, the no trend label was assigned (Figure 4i,j). These pixels contained fluctuating NDVI values within very low boundaries (i.e., low variation) and ended without notable NDVI changes.  Figure 5 shows the per pixel trend map derived by using two decades of MODIS NDVI imagery and the PolyTrend algorithm. After applying the NDVI threshold value (see Section 3.4) for selecting vegetation pixels, nearly 14 million pixels were analyzed for trend identification, accounting for approximately 44% of Iran. Figure 5 illustrates that the vegetation covers were mainly distributed in the northern and western parts of the country. In general, no consistent spatial patterns in terms of different trend types were observed throughout Iran. The concealed class was observed primarily in northeastern and southwestern parts of the country. According to earlier descriptions (see Section 3.4), this class includes both cubic and quadratic trends without a significant net change in their NDVI values. In fact, the corresponding pixels were subjected to a substantial decrease and then increase (or increase and then decrease), which resulted in low NDVI changes between the start and end of the study period. Pixels identified as No Trend were mostly scattered from the west to the east in the northern parts of the country. Pixels with linear trends were mainly distributed from the northwestern to the southwestern parts of the  Figure 5 shows the per pixel trend map derived by using two decades of MODIS NDVI imagery and the PolyTrend algorithm. After applying the NDVI threshold value (see Section 3.4) for selecting vegetation pixels, nearly 14 million pixels were analyzed for trend identification, accounting for approximately 44% of Iran. Figure 5 illustrates that the vegetation covers were mainly distributed in the northern and western parts of the country. In general, no consistent spatial patterns in terms of different trend types were observed throughout Iran. The concealed class was observed primarily in northeastern and southwestern parts of the country. According to earlier descriptions (see Section 3.4), this class includes both cubic and quadratic trends without a significant net change in their NDVI values. In fact, the corresponding pixels were subjected to a substantial decrease and then increase (or increase and then decrease), which resulted in low NDVI changes between the start and end of the study period. Pixels identified as No Trend were mostly scattered from the west to the east in the northern parts of the country. Pixels with linear trends were mainly distributed from the northwestern to the southwestern parts of the country along the Zagros Mountains. However, several other regions of the country in the northeast and central-south had dominant linear trend pixels. Further visual inspections implied that the quadratic trends occurred in fewer locations and were seen in western and southern parts of Iran. Finally, as the dominant non-linear trend type, pixels with the cubic trend were detected in northern, southern, and northeastern parts of Iran. These pixels experienced significant downward and upward trends with one local maximum and minimum (see Figure 4a,b). As a sample, Figure 4b shows a pixel with a cubic trend, in which its time-series NDVI values started at 0.35 and had a downward tendency to nearly 0.12 in 2015, and then reached over 0.22 at the end of the study period.

Trend Analysis
x FOR PEER REVIEW 10 of 22 country along the Zagros Mountains. However, several other regions of the country in the northeast and central-south had dominant linear trend pixels. Further visual inspections implied that the quadratic trends occurred in fewer locations and were seen in western and southern parts of Iran. Finally, as the dominant non-linear trend type, pixels with the cubic trend were detected in northern, southern, and northeastern parts of Iran. These pixels experienced significant downward and upward trends with one local maximum and minimum (see Figure 4a,b). As a sample, Figure 4b shows a pixel with a cubic trend, in which its time-series NDVI values started at 0.35 and had a downward tendency to nearly 0.12 in 2015, and then reached over 0.22 at the end of the study period. In addition to the trend types map ( Figure 5), three other complementary maps, including trend directions, trend significance, and trend slopes, were also generated by Pol-yTrend, and are provided in Figure 6. Figure 6a shows that the positive trends are dominant over Iran, accounting for nearly 90% of vegetation pixels. This indicates that most parts of the country were subjected to the greening process. However, several hotspot locations can be seen with negative trend directions. Furthermore, Figure 6b illustrates that nearly 56% of vegetation pixels had significant trends over the past two decades that contained linear, quadratic, and cubic trend types. On the other hand, 46% of vegetation pixels were insignificantly changed and were assigned canceled and no trend labels. Figure 6c presents the slope magnitude (linear coefficient) of vegetation dynamics, derived by fitting an independent linear model to the time-series NDVI values, which are classified into eight classes for better illustrations. As anticipated, most vegetation pixels obtained positive slope values, and overall, slope magnitudes mainly fell in two classes with minor slope values between −0.002 and 0.002 (NDVI units) per year. However, visual inspections revealed that in some regions, homogeneous patches existed with notable (either positive or negative) slope values. These areas could be considered hotspot zones, In addition to the trend types map (Figure 5), three other complementary maps, including trend directions, trend significance, and trend slopes, were also generated by PolyTrend, and are provided in Figure 6. Figure 6a shows that the positive trends are dominant over Iran, accounting for nearly 90% of vegetation pixels. This indicates that most parts of the country were subjected to the greening process. However, several hotspot locations can be seen with negative trend directions. Furthermore, Figure 6b illustrates that nearly 56% of vegetation pixels had significant trends over the past two decades that contained linear, quadratic, and cubic trend types. On the other hand, 46% of vegetation pixels were insignificantly changed and were assigned canceled and no trend labels. Figure 6c presents the slope magnitude (linear coefficient) of vegetation dynamics, derived by fitting an independent linear model to the time-series NDVI values, which are classified into eight classes for better illustrations. As anticipated, most vegetation pixels obtained positive slope values, and overall, slope magnitudes mainly fell in two classes with minor slope values between −0.002 and 0.002 (NDVI units) per year. However, visual inspections revealed that in some regions, homogeneous patches existed with notable (either positive or negative) slope values. These areas could be considered hotspot zones, being the candidates with higher priority for further examinations and management practices.
change per year), with the highest slope values for cubic and linear trends. Based on the numerical results, concealed, no trend, linear, quadratic, and cubic patterns covered positive (negative) areas of about 179,181 (33,608) km 2 , 156,042 (28,273) km 2 , 258,521 (6096) km 2 , 39,849 (3228) km 2 , and 143,756 (9575) km 2 , respectively. Based on the obtained results, over 68% of vegetation pixels experienced minor dynamics over the past two decades, with slope magnitude ranges of −0.002 and 0.002. It was also observed that more pixels (30%) were subjected to vegetation gain with higher slope magnitudes (0.002-0.05), accounting for approximately 260,640 km 2 , while nearly 3% (27,827 km 2 ) experienced higher negative slopes (−0.05-−0.002).   Further numerical results are provided in Tables 2 and 3. According to Table 1, the linear trend covering 14% of Iran was the most dominant trend pattern, followed by concealed (11%), cubic (8%), and quadratic (2%) trends, while 9% of the vegetation area remained stable (no trend). Slope magnitudes ranged between −0.048-0.047 (NDVI change per year), with the highest slope values for cubic and linear trends. Based on the numerical results, concealed, no trend, linear, quadratic, and cubic patterns covered positive (negative) areas of about 179,181 (33,608) km 2 , 156,042 (28,273) km 2 , 258,521 (6096) km 2 , 39,849 (3228) km 2 , and 143,756 (9575) km 2 , respectively. Based on the obtained results, over 68% of vegetation pixels experienced minor dynamics over the past two decades, with slope magnitude ranges of −0.002 and 0.002. It was also observed that more pixels (30%) were subjected to vegetation gain with higher slope magnitudes (0.002-0.05), accounting for approximately 260,640 km 2 , while nearly 3% (27,827 km 2 ) experienced higher negative slopes (−0.05-−0.002).  The trend analysis results are valid at a 95% confidence interval, and the area values are computed by multiplying the pixel area (250 × 250 m 2 ) by the number of pixels in each trend type. The percentage values of Positive and Negative directions were computed as the ratio to the total number of vegetation pixels in the study area. The slope magnitudes are calculated by fitting a line to 21 years of NDVI values for each pixel, and the average standard error of slope magnitudes (linear coefficient) was 4.1 × 10 −6 .
Further investigations using visual interpretation of satellite images (true-color Landsat Archive images with 30 m spatial resolution at the start and end of the study period) were carried out over three hotspot regions (see black-orange stars in Figure 6c) with high NDVI slope magnitudes. Figure 7a,b presents two satellite images captured in 2000 and 2020 over a region with high NDVI gain (lowermost star in Figure 6c). The trend analysis revealed that this region had a dominant positive trend with the highest slope magnitude range (class 7, 8). It was observed that high NDVI gains were mostly associated with agricultural expansion and vegetated wetland formation in this area. On the other hand, Figure 7c,d illustrates a location (uppermost star in Figure 6c) with a relatively intense negative slope magnitude, which was found to be linked to mountainous vegetation loss. Likewise, Figure 7e,f presents a region with agricultural expansion approximate to the waterbody areas, justifying the obtained greening process (positive trend) with high NDVI slope magnitudes using the PolyTrend. magnitude range (class 7, 8). It was observed that high NDVI gains were mostly associated with agricultural expansion and vegetated wetland formation in this area. On the other hand, Figure 7c,d illustrates a location (uppermost star in Figure 6c) with a relatively intense negative slope magnitude, which was found to be linked to mountainous vegetation loss. Likewise, Figure 7e,f presents a region with agricultural expansion approximate to the waterbody areas, justifying the obtained greening process (positive trend) with high NDVI slope magnitudes using the PolyTrend.  row (a,c,e)) and 2020 (second row (b,d,f)) over three candidate regions (black-orange stars in Figure 6c) with high NDVI slope magnitude changes for further visual interpretations. Figure 8 shows the distribution of trend types and trend directions among different vegetation covers from the MODIS land cover dataset (see Section 3.2). Figure 8 shows the proportional amounts, meaning the number of pixels divided by the total number of pixels in the corresponding land cover class. To extract such information, first, the land cover dataset was resampled to 250 m spatial resolution to match the trend analysis results  row (a,c,e)) and 2020 (second row (b,d,f)) over three candidate regions (black-orange stars in Figure 6c) with high NDVI slope magnitude changes for further visual interpretations. Figure 8 shows the distribution of trend types and trend directions among different vegetation covers from the MODIS land cover dataset (see Section 3.2). Figure 8 shows the proportional amounts, meaning the number of pixels divided by the total number of pixels in the corresponding land cover class. To extract such information, first, the land cover dataset was resampled to 250 m spatial resolution to match the trend analysis results based on the nearest neighbor method. Afterward, a reclassification step was adopted to lower the land cover types into five more general vegetation land covers, namely (1) Forest, (2) Shrubland, (3) Savannas, (4) Grassland, and (5) Cropland. In addition to inherent misclassification errors in land cover maps, the spatial mismatch between original MODIS land cover and NDVI (used for trend analysis) datasets, as well as the employed threshold for vegetation pixel selection, several inconsistencies occurred. As a consequence, several vegetation pixels (determined through NDVI threshold > 0.15) in the trend analysis fell in non-vegetation classes in the MODIS land cover data. Therefore, they were excluded for further analysis, and only vegetation pixels identified using both datasets were used in this section. The results show the dominant trend types were no trend, concealed, no trend, linear, and no trend for land cover classes of forest, shrubland, savannas, grassland, and cropland, respectively, indicating the relative dominance of no trend for the considered vegetation types. Linear had the highest number and portion in all vegetation classes, considering the significant trend types, followed by the cubic trend pattern. The results also illustrate that the vegetation cover as grassland experienced the highest statistically significant dynamics, which were mainly located along the Zagros Mountains. Furthermore, Figure 8b shows that all considered vegetation covers were subjected to NDVI gain (i.e., positive trend directions), with the two highest portions for forest (99.63%) and grassland (97.65%). However, the amounts of negative trend direction for savannas (13.82%) and cropland (17.10%) were much higher than forest (0.37%), shrubland (3.45%), and grassland (2.35%).

Vegetation Trend vs. Land Cover and Precipitation
also illustrate that the vegetation cover as grassland experienced the highest statistically significant dynamics, which were mainly located along the Zagros Mountains. Furthermore, Figure 8b shows that all considered vegetation covers were subjected to NDVI gain (i.e., positive trend directions), with the two highest portions for forest (99.63%) and grassland (97.65%). However, the amounts of negative trend direction for savannas (13.82%) and cropland (17.10%) were much higher than forest (0.37%), shrubland (3.45%), and grassland (2.35%).  (9) Water, were compared. The result was a change map containing the class-wise land cover transitions. Later, the generated change map and the trend types map were integrated to explore the relationship between land cover transitions and trend types. In this regard, the land cover transitions for linear, quadratic, and cubic patterns (statistically significant trends) were derived from the integrated map. The proportion of each land cover transition class to the total land cover transitions for each trend type was then calculated. The minor portions were excluded by focusing on land cover transitions with the five highest proportions that occurred more frequently > 5. Table 4 summarizes the calculated proportions for linear, quadratic, and cubic patterns. Accordingly, the land cover transition between Barren (2001) and Shrubland (2020) had the highest proportion for the three trend patterns. This was followed by grassland transition to shrubland for linear and cubic patterns, while the reverse change was observed for quadratic. The third-highest transition for linear, quadratic, and cubic patterns were shrubland to grassland, grassland to shrubland, and barren to grassland, respectively. Linear and quadratic patterns had the same land cover transition (i.e., barren to grassland), while the shrubland to grassland  (9) Water, were compared. The result was a change map containing the class-wise land cover transitions. Later, the generated change map and the trend types map were integrated to explore the relationship between land cover transitions and trend types. In this regard, the land cover transitions for linear, quadratic, and cubic patterns (statistically significant trends) were derived from the integrated map. The proportion of each land cover transition class to the total land cover transitions for each trend type was then calculated. The minor portions were excluded by focusing on land cover transitions with the five highest proportions that occurred more frequently > 5. Table 4 summarizes the calculated proportions for linear, quadratic, and cubic patterns. Accordingly, the land cover transition between Barren (2001) and Shrubland (2020) had the highest proportion for the three trend patterns. This was followed by grassland transition to shrubland for linear and cubic patterns, while the reverse change was observed for quadratic. The third-highest transition for linear, quadratic, and cubic patterns were shrubland to grassland, grassland to shrubland, and barren to grassland, respectively. Linear and quadratic patterns had the same land cover transition (i.e., barren to grassland), while the shrubland to grassland transition was found for cubic as the fourth-highest change. The grassland to cropland transition was also observed for all three trend patterns as the fifth-highest transition. In general, the results did not provide any specific relationship between land cover transitions and trend patterns. However, the results indicated the greening process in Iran since a high portion of Barren was changed to vegetation cover classes. In the next step, the relationship between NDVI and precipitation values was investigated for each synoptic station. In this regard, the time-series NDVI values at each synoptic station location were extracted. Afterward, the correlation coefficient between NDVI and precipitation values was computed similarly to previous studies [33,56,77,78]. It is worth noting that the correlation coefficient values were calculated using all corresponding NDVI and precipitation values throughout Iran (n = 1785) and for each synoptic station individually (n = 21). Figure 9a presents the spatial distribution of obtained correlation coefficient values for each synoptic station, indicating higher values in northern and northwestern parts of the country. The correlation coefficient values ranged from 0.002 to 0.847, with an average value of 0.371. Figure 9b shows the precipitation and NDVI values scatterplot using all corresponding observations. The computed correlation coefficient was 0.54, indicating a relatively strong influence of precipitation on vegetation dynamics.
In the next step, the relationship between NDVI and precipitation values was investigated for each synoptic station. In this regard, the time-series NDVI values at each synoptic station location were extracted. Afterward, the correlation coefficient between NDVI and precipitation values was computed similarly to previous studies [33,56,77,78]. It is worth noting that the correlation coefficient values were calculated using all corresponding NDVI and precipitation values throughout Iran (n = 1785) and for each synoptic station individually (n = 21). Figure 9a presents the spatial distribution of obtained correlation coefficient values for each synoptic station, indicating higher values in northern and northwestern parts of the country. The correlation coefficient values ranged from 0.002 to 0.847, with an average value of 0.371. Figure 9b shows the precipitation and NDVI values scatterplot using all corresponding observations. The computed correlation coefficient was 0.54, indicating a relatively strong influence of precipitation on vegetation dynamics. Furthermore, the trend of time-series NDVI and precipitation observations were also investigated since they might provide complementary information to the previous analysis. Figure 10 presents the results of comparing vegetation and precipitation trend directions (derived by applying PolyTrend). In this regard, the trend matching procedure was implemented based on 4 cases: (1) positive directions in both (i.e., NDVI and precipitation increased), (2) positive/negative directions in vegetation/precipitation (i.e., NDVI increased and precipitation decreased), (3) negative directions in both, and (4) Furthermore, the trend of time-series NDVI and precipitation observations were also investigated since they might provide complementary information to the previous analysis. Figure 10 presents the results of comparing vegetation and precipitation trend directions (derived by applying PolyTrend). In this regard, the trend matching procedure was implemented based on 4 cases: (1) positive directions in both (i.e., NDVI and precipitation increased), (2) positive/negative directions in vegetation/precipitation (i.e., NDVI increased and precipitation decreased), (3) negative directions in both, and (4) negative/positive directions in vegetation/precipitation (see Figure 10a). In general, the results implied that in nearly 60% of the cases (at 49 synoptic stations), both vegetation and precipitation had similar trend directions (either positive or negative), indicating the impact of precipitation dynamics on vegetation cover change. On the other hand, vegetation and precipitation had reverse trend directions in 36 cases, in 22 cases of which vegetation increased and precipitation decreased. Further investigations revealed that in these 22 cases, the synoptic stations are located in regions with a high amount of annual precipitation rates. This probably suggests that although the precipitation rates have declined, the required amount of water supply for vegetation growth still exists. Therefore, other influential factors may exist in these regions, as well as in other parts of Iran. In fact, although precipitation has been recognized as the most impactful criterion for vegetation dynamics, there are other impactful parameters when the water supply is enough to support biological processes [79,80]. Figure 10b illustrates that in the similar trend category (both positive/negative), the linear pattern was the dominant trend, while in the opposite trend, the concealed was dominant. This indicates that significant trend types were the dominant trend patterns in similar directions, meaning that vegetation with statistically significant trends had relatively more similar trend directions to precipitation dynamics. dynamics, there are other impactful parameters when the water supply is enough to support biological processes [79,80]. Figure 10b illustrates that in the similar trend category (both positive/negative), the linear pattern was the dominant trend, while in the opposite trend, the concealed was dominant. This indicates that significant trend types were the dominant trend patterns in similar directions, meaning that vegetation with statistically significant trends had relatively more similar trend directions to precipitation dynamics. The final analysis included land cover data for five vegetation types to investigate the possible three-way associations. It is worth noting that no synoptic stations were located in Forest, and thus, only four vegetation covers were employed to derive subsequent information. In general, the findings were relatively discrepant when considering vegetation and precipitation trends over different vegetation land covers. In fact, the existence of both similar and opposite tend directions among four vegetation covers prevented a robust depiction of any association. In particular, 18(13), 7(4), 8(6), and 1(1) synoptic stations had similar (opposite) trend directions in grassland, shrubland, cropland, and savannas, respectively. Although dominant cases with statistically significant trends (vegetation) in each land cover type had similar trend directions, the high amount of opposite trend directions hindered an authentic conclusion.

Discussion
Long-term vegetation trend analysis is essential for vegetation condition monitoring and identifying hotspot locations for further studies. Such information could be used to identify areas prone to vegetation loss to avoid further land degradation [81]. Not only do areas with negative trends require further analysis, but regions with positive trends, The final analysis included land cover data for five vegetation types to investigate the possible three-way associations. It is worth noting that no synoptic stations were located in Forest, and thus, only four vegetation covers were employed to derive subsequent information. In general, the findings were relatively discrepant when considering vegetation and precipitation trends over different vegetation land covers. In fact, the existence of both similar and opposite tend directions among four vegetation covers prevented a robust depiction of any association. In particular, 18(13), 7(4), 8(6), and 1(1) synoptic stations had similar (opposite) trend directions in grassland, shrubland, cropland, and savannas, respectively. Although dominant cases with statistically significant trends (vegetation) in each land cover type had similar trend directions, the high amount of opposite trend directions hindered an authentic conclusion.

Discussion
Long-term vegetation trend analysis is essential for vegetation condition monitoring and identifying hotspot locations for further studies. Such information could be used to identify areas prone to vegetation loss to avoid further land degradation [81]. Not only do areas with negative trends require further analysis, but regions with positive trends, especially with high slope magnitudes, should be candidate cases for further explorations. This is because such regions may be subjected to invasive vegetation growth or arbitrary agricultural expansion, which could bring adverse environmental issues [82,83]. Therefore, the trend analysis results could benefit scholars and governors in conducting research to determine the driving factors and adopting effective legislation. In this study, the longterm vegetation trend analysis was successfully performed using the PolyTrend algorithm that considers both linear and non-linear patterns, satisfying the main objective of this paper. The results indicated a higher rate of greening than browning across Iran, with a dominance of the linear trend pattern. Furthermore, the results of this study are also provided at province scales (see Figure 11), which could be of importance for sub-national decision-makers. This information, either at national or sub-national scales, could serve as initial information for restoration and conservation practices. Moreover, the relationships between vegetation dynamics and land cover, land cover transition, and precipitation were also investigated as other objectives of the current study. The grassland experienced the highest statistically significant dynamics between 2000 and 2020. Besides, the cropland and forest had the highest statistically significant negative and positive trends over the last two serve as initial information for restoration and conservation practices. Moreover, the relationships between vegetation dynamics and land cover, land cover transition, and precipitation were also investigated as other objectives of the current study. The grassland experienced the highest statistically significant dynamics between 2000 and 2020. Besides, the cropland and forest had the highest statistically significant negative and positive trends over the last two decades. The results suggested the relatively strong impact of precipitation on vegetation dynamics in general, with different correlation coefficient values across the country. In this study, an automatic workflow was implemented for linear and non-linear trend analysis using the MODIS NDVI dataset. To the best of the authors' knowledge, this is the first long-term (two decades) time-series analysis of vegetation trends performed over Iran, which also considered both linear and non-linear trend patterns. The importance of considering non-linear trends is that vegetation may respond non-linearly to different disturbances [42][43][44]. In fact, external disturbances, such as fire, clear-cut, insect outbreaks, and land cover conversion, as well as natural processes and ecological shifts could introduce non-linear alterations. Therefore, considering non-linear trend types could manifest more details and enable further investigations of vegetation dynamics in the surrounding environment. Comparing the obtained results with another previous study that employed linear trend analysis using the Theil-Sen algorithm [41], the results were relatively in agreement. Although the study periods, the implemented algorithms for trend detection, and the further assumptions differed, which could affect the results, the general depictions were in accordance. For instance, both studies (i.e., this study and [41]) observed a greening process (NDVI gain) over the western parts of the country, especially along the Zagros Mountains. Moreover, our results agreed with those reported by Kazemzadeh et al. [10], which performed a linear/non-linear trend analysis over the In this study, an automatic workflow was implemented for linear and non-linear trend analysis using the MODIS NDVI dataset. To the best of the authors' knowledge, this is the first long-term (two decades) time-series analysis of vegetation trends performed over Iran, which also considered both linear and non-linear trend patterns. The importance of considering non-linear trends is that vegetation may respond non-linearly to different disturbances [42][43][44]. In fact, external disturbances, such as fire, clear-cut, insect outbreaks, and land cover conversion, as well as natural processes and ecological shifts, could introduce non-linear alterations. Therefore, considering non-linear trend types could manifest more details and enable further investigations of vegetation dynamics in the surrounding environment. Comparing the obtained results with another previous study that employed linear trend analysis using the Theil-Sen algorithm [41], the results were relatively in agreement. Although the study periods, the implemented algorithms for trend detection, and the further assumptions differed, which could affect the results, the general depictions were in accordance. For instance, both studies (i.e., this study and [41]) observed a greening process (NDVI gain) over the western parts of the country, especially along the Zagros Mountains. Moreover, our results agreed with those reported by Kazemzadeh et al. [10], which performed a linear/non-linear trend analysis over the Taleghan watershed, Alborz, Iran. In particular, the dominant positive trends were observed in both studies.
The trend analysis results demonstrated a greening process in Iran (see Figure 6). Two other data sources were also employed to investigate the greening process in Iran. These datasets were independent of the trend analysis and thus, can affirm the greening process from other perspectives. The first dataset included two MODIS land cover maps at 500 m spatial resolution for two years, near the start (2001) and the end (2020) of the study period [59]. It is worth noting that the land cover map of 2001 was employed due to the unavailability of the 2000 land cover map. In this regard, original maps (with 17 classes) were initially reclassified into two broad classes: vegetation and non-vegetation, and then a change map was generated through image differencing. The results indicated that some regions were subjected to vegetation loss, while vegetation gain was observed in most areas. Considering the total transitions, nearly 36,240 km 2 of vegetation gain was recorded in Iran. It should be pointed out that this result suggested the greening process by having more vegetation pixels. Furthermore, the annual MODIS Net Primary Productivity (NPP) data at 500 m spatial resolution between 2001 and 2020 were used as the second independent data source. This is due to the fact that higher vegetation would lead to higher NPP, and thus NPP changes can be employed to explore the greening process [84]. The findings of comparing the NPP data also confirmed the greening process in Iran, with an NPP increase of about 1.5 × 10 5 Kg * C m 2 for the entire country of Iran. However, the spatial patterns depicted that some regions experienced NPP decrease as expected (i.e., negative direction in NDVI) but pixels with NPP increase were dominant. Therefore, the greening process in Iran was also confirmed with two other datasets that were not used for trend analysis. These could imply the robustness of the implemented workflow for vegetation trend analysis.
Currently, NDVI datasets can be obtained from AVHRR, Landsat−8, Himawari−8, and Sentinel−2. However, the utility of these data is limited due to their spatial resolutions or the unavailability of long-term NDVI for such study. For instance, AVHRR can provide long-term NDVI with a very coarse spatial resolution (~9 km). On the other hand, Landsat−8 could provide long-term NDVI data back to 1984, while its medium spatial resolution (30 m) makes it challenging for country-wide studies. This is due to the computational burden of employing numerous pixels for trend analysis, and thus it has been used mostly for local studies [10,82]. The other sources, e.g., Sentinel−2 and Himawari−8, lack long-term NDVI data since they were launched after 2014. Additionally, Landsat−8 and Sentinel−2 could not provide daily observations, and the existence of cloud cover might lead to a large NDVI data gap [85,86]. These points make MODIS NDVI an appealing choice for long-term and large-scale studies. Therefore, in this study, MODIS NDVI data at 250 m spatial resolution were employed for vegetation trend analysis throughout Iran. The results provided an excellent overview of vegetation dynamics in different parts of Iran. Additionally, it assisted in determining hotspot locations (i.e., high vegetation dynamics) for further studies. Subsequently, the identified hotspot zones could be studied with higher resolution satellite data, such as Landsat−8 and Sentinel−2. For instance, higher resolution NDVI data could be used to monitor vegetation growth in smaller-scale areas, providing detailed information [87]. Additionally, integrating NDVI datasets from different sources could alleviate the drawback of each, leading to a more appropriate NDVI dataset for the desired application [88]. In particular, multi-sensor NDVI fusion enables the production of spatially and temporally consistent NDVI datasets for vegetation cover monitoring [88,89]. This is because higher resolution satellite imagery could provide further details about the occurred changes. Meanwhile, further studies are required to investigate the exact time of abrupt changes in hotspot areas. This could help identify the associate drivers of abrupt changes, classifying them as natural or anthropogenic forcings [90][91][92].
This study employed MODIS land cover maps to investigate the relationship between vegetation trend types, land cover, and land cover transitions between 2001 and 2020. Although these results could be helpful in interpreting the vegetation cover changes for further studies, they included several uncertainties that should be considered. First, the spatial mismatch between original MODIS land cover maps (500 m) and trend analysis results (250 m) could be a source of uncertainty. This is due to the fact that spectral values of pixels with 250 m and 500 m spatial resolution in locations without homogenous land covers may differ, and they might be recognized as two different land covers. Second, the inherent misclassification error of land cover maps would propagate into the final results when producing change maps, which could also be considered another source of error when interpreting the results [24].
In this study, only precipitation data from synoptic stations were employed to study its relationship with vegetation trend directions. Despite their associate errors, these in-situ observations are the most accurate precipitation data due to their intuitive measurements. However, they are point observations, and their spatial inconsistency with satellite data (pixels) may introduce some uncertainties in the comparison procedure. Therefore, it is suggested to use long-term satellite data with higher-spatial resolutions for specified locations to minimize spatial inconsistencies. Although precipitation has been recognized as the main driving force of vegetation dynamics, which was also observed for 60% of cases in this study, and the correlation coefficient was 0.54, other climatic factors such as temperature [79] and solar radiation [80] affect the gradual vegetation dynamics. Therefore, utilizing different climatic parameters, along with topographical variables, could result in obtaining a more robust conclusion about the driving forces of vegetation dynamics in Iran. Consideration of other climatic parameters was beyond the scope of this study and shall be pursued in future studies.

Conclusions
In this study, a time series of MODIS NDVI data collected between 2000 and 2020 was employed to study the annual vegetation trends on a national scale in Iran. The PolyTrend, an automated trend analysis algorithm that considers both linear and nonlinear trend patterns, was implemented. A visual assessment using candidate samples suggested the applicability of the PolyTrend for trend type identification. The results indicated a higher rate of greening in Iran with the dominance of the linear, followed by the cubic as the statistically significant trends. Further information illustrated some locations with more intensified vegetation change (positive or negative) that could be nominated regions in future studies. It is believed that the results could support achieving Sustainable Development Goals (SDGs) by serving as an initial stage study for establishing conservation and restoration practices. The correlation coefficient between precipitation and NDVI data was 0.54 considering all corresponding observations across the country. Furthermore, the vegetation and precipitation dynamics comparison suggested their mutual directions in 60% of cases with more statistically significant changes. However, some observed discrepancies suggested the consideration of other climatic and topographical variables to study the driving forces of vegetation cover more profoundly. Land cover types were also used to study two-way and three-way relationships between vegetation trend-land cover, and vegetation-precipitation-land cover. Besides, an attempt was made to investigate the relationship between trend types and land cover transitions. However, the observed inconsistencies hamper the driving of a robust relationship between these data sources. Accordingly, more accurate data, especially land cover maps, with more spatially consistent resolutions should be considered. Data Availability Statement: The NDVI and land cover datasets supporting the reported results of this paper are publicly available from the Google Earth Engine data catalog.