Next Article in Journal
Spatiotemporal Change in Evapotranspiration across the Indus River Basin Detected by Combining GRACE/GRACE-FO and Swarm Observations
Previous Article in Journal
SDRSwin: A Residual Swin Transformer Network with Saliency Detection for Infrared and Visible Image Fusion
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Trend of Changes in Phenological Components of Iran’s Vegetation Using Satellite Observations

by
Hadi Zare Khormizi
1,
Hamid Reza Ghafarian Malamiri
2,
Zahra Kalantari
3 and
Carla Sofia Santos Ferreira
4,5,6,*
1
Faculty of Natural Resources, University of Tehran, Karaj 3158777871, Iran
2
Department of Geography, Yazd University, Yazd 8915818411, Iran
3
Department of Sustainable Development, Environmental Science and Engineering (SEED), KTH Royal Institute of Technology, 10044 Stockholm, Sweden
4
Department of Physical Geography and Bolin Centre for Climate Research, Stockholm University, 11419 Stockholm, Sweden
5
Applied Research Institute, Polytechnic Institute of Coimbra, 3045-601 Coimbra, Portugal
6
Research Centre for Natural Resources, Environment and Society (CERNAS), Polytechnic Institute of Coimbra, 3045-601 Coimbra, Portugal
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(18), 4468; https://doi.org/10.3390/rs15184468
Submission received: 7 August 2023 / Revised: 4 September 2023 / Accepted: 7 September 2023 / Published: 11 September 2023
(This article belongs to the Section Environmental Remote Sensing)

Abstract

:
Investigating vegetation changes, especially plant phenology, can yield valuable information about global warming and climate change. Time series satellite observations and remote sensing methods offer a great source of information on distinctions and changing aspects of vegetation. The current study aimed to determine the trend and rate of changes in some phenological components of Iran’s vegetation. In this regard, the current study employed the daily NDVI (Normalized Difference Vegetation Index) product of the AVHRR sensor with a spatial resolution of 0.05° × 0.05°, named AVH13C1. Then, using the HANTS algorithm, images of amplitude zero, annual amplitude, and annual phase were prepared annually from 1982 to 2019. Using TIMESAT software, the starting, end, and length of time of growing season were calculated for each pixel time series to prepare annual maps. The Mann–Kendall statistical test was used to investigate the significance of changes during the study period. On average in the entire area of Iran, the annual phase was declining with a trend of −0.6° per year, and the time for the start and end of the season was declining by −0.3 and −0.65 days per year, respectively. Major changes were noticed in the northeast, west, and northwest regions of Iran, where the annual phase declined with a trend of −0.9° per year. Since the annual growth cycle of the plant (equivalent to 356 days) was in the form of a sinusoidal signal, and the angular changes in the sine wave were between zero and 360°, each degree of change was equivalent to 1.01 days per year. Therefore, the reduction in the annual phase by −0.9 degrees almost means a change in the time (due to the earlier negative start phase) of the start of the annual growth signal by −0.9 days per year. The time of the start and end of the growing season declined by −0.6 and −1.33 days per year, respectively. The reduction in annual phase and differences in time of the starting season from 1982 to 2019 indicate the acceleration and earlier initiation of various phenological processes in the area.

Graphical Abstract

1. Introduction

Climate change brought on by increasing global warming is one of the most difficult global issues [1]. Global warming has had an impact on the physical and biotic systems on every continent [1,2,3]. Since the start of the industrial revolution (1880–2022), the average global temperature has risen by 0.85 °C, and this trend accelerated after 1950 [4,5]. Additionally, according to [5], the Earth’s surface will be 0.3 to 4.8 °C warmer by the end of the 21st century. One of the most crucial biological factors to consider when examining the effects of global warming on the Earth’s vegetation coverage and ecosystems is the phenological characteristics of vegetation [6], which describe the timing of recurring events in plant life involving both living and non-living factors [7,8]. In other words, many researchers identify vegetation phenology as a sensitive indicator for investigating the effects of global warming and climate change on terrestrial ecosystems. According to [2,9,10], these characteristics point to the flow of energy, water vapor, and carbon between the biosphere and the lowest regions of the atmosphere. Based on [11,12], plant phenology changes will have a negative impact on environmental systems, which will impact forestry, food production, agriculture, human health, and the economics of many countries.
Various biotic and abiotic factors such as temperature, precipitation, hours of sunshine, physical and chemical characteristics of soil, and composition of vegetation coverage affect vegetation phenology [13,14,15,16,17]. Temperature and precipitation are two critical factors associated with vegetation phenology changes, whose effects are different in various regions of the world [18]. In humid temperate locations and high latitudes in particular, temperature is one of the most crucial factors in determining when the season begins and the processes that determine plant phenology [18,19,20,21]. Precipitation is particularly important in phenology changes in tropical and arid regions [22,23]. However, since temperature is thought to be the primary element in plant phenology variations, precipitation will be most effective in modifying plant phenology when the temperature is appropriate for plant growth [9].
There are two different methods for monitoring vegetation phenology [2]. The first approach is based on ground truth data collection about the annual changes of phenological stages in response to environmental parameters and climate change. Although these methods, which concentrate on a small number of land survey locations, are effective and accurate, they are also expensive and, in some cases, impractical. The second strategy was recently created and is based on remote sensing and satellite imagery. Since the plant’s phenology in response to climate variables is different in time and space [15], the remote sensing technology provides the possibility of repeated data collection (time series satellite observations) with a wide area coverage [24]. All phenology studies using remote sensing are based on two methods [18]: one method based on the investigation of trends in the vegetation indices over time, and another method based on extracting and estimating phenological parameters (phenophases). Creating time series of vegetation indices over time requires complete and smooth time series data to overcome missing (gaps), outlier data, as well as errors driven, e.g., by noise induced via cloud cover, mixed pixel effects, sensor failures, and anomalies [25]. The gap-filling procedure could be carried out using statistical methods, curve fitting, and data conversion techniques [26]. In turn, the estimation of phenological parameters or the detection of phenophases is carried out using thresholding, curve-derived, and functional fitting methods [27].
At both regional and global levels, numerous studies have examined how plant cover phenology changes in response to climate change [28]. A study of phenological changes throughout the world from 1982 to 2003 using NDVI time series data, for instance, revealed that the start of the season began 7.9 days earlier while the season’s conclusion was delayed by 9.4 days [29]. In another study, the NDVI time series between 1982 and 2018 were investigated at the global scale and showed that average annual phases decreased by 9 days, indicating an earlier start of the growing season [30]. In the Northern Hemisphere, between the periods of 1982 to 1999 and 2000 to 2008, the start time of the season was 5.2 and 0.2 days earlier, and the end time of the season was delayed by 3.4 and 23 days, respectively [31]. In middle and high latitudes of the Northern Hemisphere, the progress at the start of the season and the delay at the season’s conclusion were slower from 1982 to 2013 than from 1992 to 2013 [32]. In a different study, the start and end of the season were found to have been delayed at a rate of less than 1.5 days per year between 2000 and 2009 using the Enhanced Vegetation Index (EVI), calculated from the MODIS sensor in Northeast China forests [33]. In Northeast China, the time for the start of the growth season significantly occurred earlier and the time for the end of the growth season occurred later from 1982 to 2015 [34]. In alpine grasslands at the north of the Tibetan Plateau, the time for the start and end of the season and the length of the season experienced a trend of delay, progress, and shortening at the rates of 1.32, −2.06, and −3.38 days per decade, respectively [35]. Generally, relatively limited attention has been given to these studies regarding the detailed examination of the effect of temperature and rainfall on changes in phenology components.
In two previous studies in Iran, only two four-year time series of 1982–1985 and 2015–2018 were compared and proved that a significant difference existed between the main phenological factors of plants (e.g., annual phases) [36,37]. In the current study, the complete and continuous time series of the NDVI from 1982 to 2019 was used. Also, in addition to examining the annual changes in various phenological parameters, including the starting of the growing season, the significance of the trend and the slope of the changes were also examined based on Kendall’s test. This study aims to determine the trend and rate of changes in phenological components (i.e., amplitude 0, amplitude 1, phase 1, start of the season, end of the season, and length of the season) of Iran’s vegetation using the time series of the NDVI index based on the Advanced, Very High-Resolution Radiometer (AVHRR) sensor during the period from 1982 to 2019. Besides investigating the trends over a longer period, this study goes beyond previous studies by focusing on complete time series based on high accuracy in the restored data from errors in the time series. Understanding the trend of changes in phenology components during the past few decades can increase our knowledge about the impact of climate change and support adaption for future changes.
The manuscript is structured in five main sections. In the following sections, we present the methodology (Section 2), including a description of the study area, the data series used, and the theoretical background of the algorithms used. Then, we present the results of the trends in the different phenology components (Section 3), discuss the results (Section 4), and present the main conclusions (Section 5).

2. Materials and Methods

2.1. Study Area

This study involves the entire area of Iran. The country has a very diverse climate due to its unique geographical location, high altitude changes (Figure 1), and considerable distance from north to south. The study area extends from the lowlands of the Persian Gulf, Caspian Sea, and Oman Sea coasts to the Alborz and Zagros highlands. The Zagros highlands extend from the northwest to the south, and the Alborz highlands extend from the northwest to the northeast. Therefore, the humidity conditions range from arid lands such as Lut Desert and Dasht Kavir to very humid lands such as the edge of the Caspian Sea.

2.2. Methodology

The methodology is based on a study of the AVHRR-sensor-based annual time series from 1982 to 2019 (Figure 2). The HANTS algorithm [38] was used to rebuild the time series and thus calculate the amplitude zero (mean annual NDVI), amplitude 1 (annual amplitude), and phase 1 (annual phase) parameters for each year. The images of the time for the season’s start (SOS), finish (EOS), and length (LOS) were generated using the TIMESAT software [39], after the reconstructed NDVI time series were loaded. The components of amplitude zero, amplitude 1, phase 1, the beginning of the season, the end of the season, and the duration of the season in the 38-year time series were then examined.
The changes in the mentioned factors during the time period from 1982 to 2019 were displayed at four selected locations in different ecological regions. Location 1 was selected in the Hyrcani ecological zone (forests of northern Iran), locations 2 and 3 in the Zagros ecological zone, and location 4 in the Iran and Turan ecological zones. The locations were selected to represent distinct vegetation. The main feature of the Hyrcani region is dense vegetation with the dominance of deciduous species (such as Fagus orientalis and Carpinus betulus). The ecological zone of Zagros has medium vegetation cover with deciduous plant species (such as Quercus L.) and the vegetation zones of Iran and Turan have weak vegetation cover (central and eastern parts of Iran) with dominant species such as Artemisia sieberi Besser. In the following, the mentioned items are explained in full detail. In the next sections, the theoretical foundations, the investigated parameters, and the working method are fully explained. Each location represents one pixel of the AVHRR image (5 × 5 km).

2.2.1. NDVI Data Series

The daily NDVI result of the AVHRR sensor (AVH13C1), with a spatial resolution of 0.05°, was used in the current investigation from 1982 to 2019 (NOAA CDR Program data). The NOAA CDR (Climate Data Record) software created the NDVI calculation algorithm [40]. One of the most often used indicators for analyzing plant dynamics is the NDVI index. Equation (1) is used to calculate it as follows [41]:
N D V I =   N I R R E D   N I R + R E D
where the NIR and RED represent the spectral reflectance of the red and near-infrared bands, respectively (bands 1 and 2 of the AVHRR sensor).
The time series of NDVI data are intrinsically periodic due to the fact that they have cyclic growth stages from germination to fully mature and then harvesting stage. For periodic phenomena such as NDVI time series, Fourier theory [42] can be used to decompose its signal into significant periodic components with their amplitudes and phases. In the present study, the changing trend in the phenological parameters of NDVI time series was investigated in the period from 1982 to 2019. Specifically, the HATNS algorithm was used to extract the phenological components of amplitude 0 (average of time series), amplitude 1 (first significant periodic component), amplitude 2 (second significant periodic component), phase 1, and phase 2.
Despite many advantages of the NDVI index, it also includes some limitations affecting the interpretation of the results, such as the influence of soil background [41]. However, these possible errors are not affecting the current study, because the phenology components are extracted from vegetation growth cycles. Furthermore, the background soil cannot create a sinusoidal curve over time (due to the absorption and reflection effect of vegetation in the red and near-infrared bands). In this research, the areas with scarce vegetation were not included, which also mitigates possible errors in NDVI index. Nevertheless, this study explores changes over time rather than the exact amount of vegetation cover, which also mitigates the impact of possible errors.
The TIMESAT program was also used to extract the phenological components of the time at the beginning, end, and duration of the season [39]. There are four methods to determine the start and end of the growing season in the TIMESAT software [39]. The first method is based on the seasonal amplitude, defined between the base level and the maximum value for each individual season. In the second method, the start and end of the season occur when the curve reaches an absolute value, defined in the units of the data (e.g., reaching a level of 0.2 NDVI in Figure 3). The third method is based on the relative amplitude for the whole time series. This amplitude is calculated as the difference between the robust mean maximum and the robust mean base level. Unlike method 1 and 2, this method produces an index to determine the start and end of the growing season, which is the relative value of the vegetation index to reach the start and end of growth for one pixel over time, but this value is different for different pixels. In the fourth method, the start and end occur when the curve crosses the STL (Seasonal Trend Decomposition by LOESS) trend line. In this study, considering the variety of vegetation in the study area and the superiority of the third method [36], a relative amplitude value of 0.3 was used to determine the start and end of the growing season [36].
In general, HANTS and TIMESAT require a continuous and consistent time series growth signal without missing (gaps) and outlier data in order to derive phenology components from the time series of vegetation indicators. However, the time series of vegetation indices of satellite images are commonly prone to missing data, outliers, and noise (in spatial and temporal domain) due to atmospheric dust, airborne particles, sensor inefficiency, gases, and especially clouds [37,43]. TIMESAT software has different algorithms such as the Savitzky–Golay, asymmetric Gaussians, and double logistic to fit a curve on the initial data of time series to obtain complete information [39]. Because the HANTS method was more accurate than the techniques used in TIMESAT, it was utilized to rebuild missing data, identify outliers, and exclude them from the growth curve [36,44]. Therefore, at first, the NDVI time series was reconstructed using the HANTS algorithm. Then, these time series were entered into TIMESAT software and different phenology parameters (SOS, EOS, and LOS) were extracted from them. The selection of the HANTS algorithm was based on its reliability in reconstructing time series and because it provides outputs including the phase and range of the time series, useful for understanding plant phenological processes [38,45]. The HANTS algorithm is introduced in Section 2.2.2.

2.2.2. Theoretical Foundations

HANTS Algorithm

In order to deconstruct time series of satellite data, the Fast Fourier Transform (FFT) is utilized [38,45]. This is the foundation of the HANTS method. Here, the technique is first presented, and then the necessary parameters are listed in order to create a trustworthy model. Equation (2) can be used to characterize the time series y i as a Fourier series if it has N observations (i = 1, 2, …, N).
y i = a 0 + j = 1 M a j c o s ( w j t i φ j )
where t i is the moment at which the ith sample is taken and w j is the jth frequency of the harmonic period in the Fourier series. M is the Fourier series’ frequency number (M = N), while a j and φ j are the harmonic period’s jth harmonic’s amplitude and phase, respectively. Zero frequency ( a 0 ) has an amplitude equal to the mean of all N observations of the variable y since zero frequency has no phase. Equation (3) states that harmonic frequencies ( w j ) are created by multiplying a fundamental frequency (wj = 2π/N, for instance) by an integer i (i = 1, 2, …, N).
w j = ( 2 π N ) × i i = 1 , 2 , 3 , , N
The amplitude ( a j ) and the value of the phase (φj), which are obtained by curve fitting the time series of observations using the least squares approach, are the unknown parameters of the Fourier series in the HANTS algorithm after the number of frequencies (M) and frequency ( w j ) have been established. To use the HANTS algorithm to build a trustworthy model of a time series, the following parameters must be defined:
  • Valid data range: The permissible range of observed values. By giving zero weight to values outside of this range, they are eliminated at the first stage.
  • Period: the total number of observations in the Fourier series for each periodic component.
  • Number of frequency (NOF): The number of frequencies affects how much information is used to reconstruct the signal. Compared to a high frequency number, a low frequency number generates a signal with less clarity.
  • Direction of outlier: This shows the direction of an outlier (outliers) regarding the current curve model. For example, cloud cover leads to a decrease in NDVI values, so the direction of an outlier can be selected low in the algorithm.
  • Fit Error Tolerance (FET): The fit error tolerance (FET) parameter specifies the maximum deviation from the curve’s current form that is still allowable in the selected direction. If an observation’s deviation exceeds the FET, it is labeled as an outlier and removed from the calculations by being given 0 weight after each iteration.
  • Degree of Over-determinedness (DOD): This statistic illustrates the absolute minimal number of additional data points that should be used to fit a curve. The number of valid observations should always be greater than the number of parameters required to fully describe the signal (2NOF-1).
The parameters utilized by HANTS to reconstruct each one-year time series from 1982 to 2019 are displayed in Table 1. In accordance with Table 1, the permissible NDVI range for vegetation was taken to be between 0 and 1, and the base period was taken to be 365 days per year, or 365 daily NDVI photos. NDVI values outside the 0–1 range were not included in this study, since they represent other areas such as those covered with water or snow. According to periodic components in a year [23,46], the number of frequencies (NOF) was taken to be 3 and the FET values were taken to be 0.1. This FET value allows the noisy NDVI data (positive values) to be smoothed. Due to the fact that cloud cover causes NDVI values to drop, the direction of an outlier was deemed to be low. Figure 3 shows a one-year time series that was recreated using HANTS and the TIMESAT notion of the starting time, end of the season, and season length.
Due to the lack of growth curves and sinusoidal trends, it is not possible to extract reliable phenology components (such as SOS, EOS, and LOS) in areas with weak vegetation. Therefore, in the present study, areas with an average annual NDVI of less than 0.1 were ignored in the final maps. This value should not be confused with the valid data range in HANTS algorithm for the reconstruction of NDVI time series.

Fast Fourier Transform (FFT) and the Amplitude and Phase Concepts

To break down a complicated signal into sine or cosine components, one can apply Fourier series analysis [42]. The matrix representation of Equation (2) is given by Equation (4):
y 1 y N = f 1 ( t 1 ) f M ( t 1 ) f 1 ( t N ) f M ( t N ) a 1 a M
In Equation (4), observations yi, with i = 1, …, N, in MN time basis functions fj(ti) with j = 1, …, M are expressed. Suppose the time series is given by a data vector y, and the amplitudes are represented by a vector a, then one can write (Equation (5))
y = F a
where F is a matrix of which the elements depend only on the value of N. When the elements of y, F, and a are considered as complex numbers, it becomes particularly easy to find the solution. Equation (6) is the result of multiplying the two sides of the relationship by the matrix F’s transposition:
F T y = F T F a or a = ( F T F ) 1 F T y
Then, we estimate the amplitude and phase of each component of vector (a) using the least squares approach. The Fast Fourier Transform (FFT) is a method that multiplies the matrices in Equation (6) using the fewest possible computations. The amplitudes and phases of all frequencies (N phases and N frequencies) are included in the FFT results obtained for a time series with (N) number of data points. Any alternating function can be created by adding a number of sine or cosine waves with various frequencies, according to the Fourier series expansion [42]. As a result, the non-sinusoidal wave that results can be divided into components that have sinusoids with various frequencies. If these sine waves are combined together, the primary wave (signal) is produced [42]. For instance, the major wave (black line) in Figure 4A represents the change trend of the NDVI index during a 12-month period. This signal was produced from three different signals with distinct amplitude, frequency, and phase (Figure 4A, blue, green, and red lines, respectively). The main signal (black line) is divided into three cosine waves with various amplitudes and frequencies using the HANTS algorithm and the proper settings (three frequencies). The cosine waveform with the same frequency as the original signal is the first harmonic, followed by the second harmonic if the frequency is doubled and the nth harmonic if the frequency is n times larger. The amplitude, phase, and frequency of each sine or cosine wave are unique. This information is quite effective in recognizing and classifying phenology and determining the changes and dynamics of plants [30,37]. Figure 4B shows the concepts of amplitude, phase, and the phase difference between the two waves. The change process of the plant growth curve (i.e., NDVI) is almost like a sine wave. In the trigonometric circle, the angular changes of the sine are between zero and 360 degrees. Therefore, each complete cycle of a wave is 360 degrees. On the other hand, the annual growth changes of plants and the repetition of the same event in the following year last for one year or 365 days. Therefore, in the HANTS algorithm, a change of one degree in the annual phase is approximately equal to one day in the growth cycle of plants. For example, a decrease of one degree in the annual phase means a change of the signal to the negative side of one degree (or one day).

2.2.3. Analysis of Changing Trend in Phenology Components

The current study used the Mann–Kendall test [47,48] to determine how different vegetation phenology components were changing over time in Iran. The data series’ randomness and lack of a trend are the test’s null hypotheses, and if there is a trend in the data series, the alternative hypothesis is proposed. The following steps are used to calculate the test’s statistics.
(A)
Computing the difference between observations, applying the sign function, and extracting the parameters based on Equation (7):
s = k = 1 n 1 j = k + 1 n sgn ( x j x k )
where x j   and x k are the jth and kth data point in the time series, respectively, and n is the number of observations. Equation (8) can also be used to compute the sign function:
sgn ( x ) = + 1 i f ( x j x k ) > 0 0 i f ( x j x k ) = 0 1 i f ( x j x k ) < 0
(B)
Variance calculation based on Equations (9) and (10):
var ( s ) = n ( n 1 ) ( 2 n + 5 ) i 1 m + ( t 1 ) ( 2 t + 5 ) 18 : n > 10
var ( s ) = n ( n 1 ) ( 2 n + 5 ) 18 : n 10
where t is the frequency of data points with the same value, n is the overall number of data points, and m is the total number of series with at least one repeated data point.
(C)
Extracting the z statistic using Equation (11):
z = s 1 var ( s ) i f s > 0 0 i f s = 0 s + 1 var ( s ) i f s < 0
If Equation (12) holds true, the null hypothesis will be accepted in a two-tailed test to determine the trend of the data series:
z z a / 2
z a statistics is the standard normal distribution at the significance level of a, which is utilized as a/2 due to the two-tailed test, and a is the significant level taken into consideration for the test. The upward trend of the data series is taken into consideration if the z-statistic is positive, and the downward trend is taken into consideration if it is negative. Sen’s slope (Sen’s) of the shifting trend line is a crucial signal in the Mann–Kendall test. Here, the slope between each pair of values in the time series is calculated following the approach described by [49,50], and the median of the resulting slopes is extracted (Equation (13)):
S e n , s = m e d i a n ( x t x s t s )
where t is a time unit after s, and x t and x s   are observational data at times t and s, respectively. Sen’s values that are positive or negative represent an upward or downward trend, respectively.

3. Results

3.1. Trends in Amplitude Characteristics

The significance map of the trend changes from 1982 to 2019 regarding amplitude 0 (mean annual NDVI) and Sen’s line slope of amplitude 0 are shown in Figure 5. According to Figure 5A, in the inner plain regions and the eastern regions of Iran, the changes in amplitude 0 show a significant negative trend from 1982 to 2019, which indicates a decline. In the coastal areas of the Caspian Sea and northwestern Iran, the changing trend of amplitude 0 is significant and positive, meaning that the average level of vegetation increased over the study period. Figure 5B shows the Sen’s line slope of amplitude 0 changes from 1982 to 2019. Sen’s slope line is positive in the mountainous areas of Alborz and Zagros, and negative in the inner plain regions of Iran. Figure 6 displays the evolving trend of amplitude 0 at four different locations (the locations of the dots are given in Figure 1) that were chosen to represent four diverse biological zones. In locations 1 and 2, the changing trend in amplitude 0 was positive, whereas in locations 3 and 4, no significant trend was observed in amplitude 0. Locations 1 and 2 indicate the Hyrcani and North Zagros ecological zones, respectively. Here, a positive trend can be seen in the level of vegetation considering the higher rainfall and soil moisture than in locations 3 and 4. According to Figure 6, some anomalies such as the sharp decrease in amplitude 0 in location 2 observed in 2008 were caused by drought. Also, the sharp increase in amplitude 0 in location 3, recorded in 1994 and 2006, is expected to be caused by the annual rainfall exceeding the long-term average.
In most parts of Iran, the changes in amplitude 1 (annual amplitude) do not show a significant trend from 1982 to 2019 (Figure 7A). In some areas located in the west of Iran and on the edge of the Caspian Sea, there is a significant positive trend. In the central areas of Iran, negative and significant trends can be seen. The Sen’s line slope of changes in amplitude 1 is shown in Figure 7B. According to the shape of the western, northwest, and northern regions of Iran, Sen’s line slope of the changes in amplitude 1 is positive, whereas in other areas, it is negative. The changing trend of amplitude 1 at four different locations (same locations as in Figure 1) from 1982 to 2019 is shown in Figure 8. In locations 1 and 4, the changing trend of amplitude 1 is significantly negative, whereas in locations 2 and 3, it is positive.

3.2. Trends in Annual Phase

The significance level of the changing trend of phase 1 (annual phase) and the changes in Sen’s line slope of phase 1 are shown in Figure 9. According to Figure 9A, in the west, northwest, and northeast of Iran, the changes in phase 1 significantly show a negative trend, indicating a decrease from 1982 to 2019. Following the reduction in annual phases, the start time of the season was also accelerated as a dependent variable. According to Figure 9B, phase 1 in the west, northwest, and northeast of Iran started earlier by 0.5 to 3 days per year (average of 0.9 ≈ 1 day) during the study period. Also, phase 1 is decreasing on average across the country with a slope of −0.6 (degrees) days per year. For example, in Figure 10, the changing trend of phase 1 is shown in locations 2 and 3, which are located in areas with a negative trend. According to Figure 10, phase 1 is decreasing with a slope of 1 and 0.5 days per year in locations 2 and 3, respectively. Locations 2 and 3 are placed in the Zagros ecological zone. These areas are generally mountainous and the process of starting the growing season and the annual phase are controlled by the temperature factor. The reduction in annual phases or the earlier start of annual growth curves during the years 1982 to 2019 can indicate an increase in temperature or an increase in temperature that triggered the start of the season in these areas.
In areas such as the edge of the Caspian Sea and the central, eastern, and southern regions of Iran, no significant trend was observed in the changes of phase 1. However, it should be noted that in most of these areas, the slope of phase 1 changes is negative (Figure 9B). For example, the non-significant trend of phase 1 changes is shown in locations 1 and 4 (Figure 10).
As in Figure 6, some anomalies can be seen in Figure 10. For example, in location 1 in 2000, the annual phase has decreased. This could have been driven by drought conditions according to the meteorological statistics. It should be noted that apart from rainfall, temperature is also an important factor in determining the start of the growing season and determining the annual phase.

3.3. Trends in the Start of the Season

In Iran’s western, northwestern, northern, and northeastern regions, the time when the season conventionally begins is trending significantly downward (Figure 11A). In these regions, the time for the start of the season often decreases with a slope of −0.6 days per year, as shown in Figure 11B. The slope line of the period for the start of the season displays a small downward trend in all other parts of Iran except the central regions. Figure 12 displays the trend and Sen’s slope of the changes in the time for the season to begin from 1982 to 2019 at four different locations. In locations 1 and 2, the season’s start time was drastically shortened to 0.56 and 0.63 days every year, respectively, whereas locations 3 and 4 showed no discernible trends. As mentioned in Section 3.1, locations 1 (Hyrcanian ecological zone) and 2 (ecological zone of North Zagros) are placed in mountainous areas and in high latitudes in Iran, and so the beginning of the growing season is controlled by temperature. It seems that during the last four decades, the increase in temperature in these areas has accelerated the beginning of the growing season.
The start of the season and the annual phases are directly related and both indicate the start of the growing season. However, the annual phases can be extracted according to the entire annual cycles of the plant growth signal, if the start of the growing season is obtained using thresholding techniques. Therefore, in this context, annual phases can provide more reliable results.

3.4. Trends in the End of the Season

In the western, northwestern, and portions of the northeastern areas of Iran, the time at which the growing season ends has a considerable downward trend from 1982 to 2019 (Figure 13). The end of the season is approaching sooner in these regions, with a Sen’s slope of −1.33 days annually (−0.65 across Iran). The alterations in the conclusion of the season’s timing generally reveal a poor and declining tendency in Iran. Figure 14 displays the varying trend of the time for the season’s end at four different locations from 1982 to 2019. The season’s end time is advancing with a slope of 1 day each year at location 1 (Hyrcanian ecological zone), while at location 2, it is significantly decreasing at 1.7 days per year. In general, the determining factor for the parameter of the end of the growing season in most of Iran (except the Hyrcanian ecological zone) is soil moisture. Therefore, in most parts of Iran, the length of the growing season shows a negative trend (Figure 13B, green-blue-violet areas). One of the reasons for this can be the increase in temperature during the last four decades.
As mentioned, soil moisture in most of the arid and semi-arid regions of Iran determines the end of the growing season. Therefore, the occurrence of droughts is one of the important factors in reducing the time of the growing season. For example, an anomaly can be seen in location 1 in 2000. Meteorological statistics in this area show the occurrence of a drought this year.

3.5. Trends in the Length of the Season

The changes in the length of the season in a small part of the western region of Iran show a significant negative trend from 1982 to 2019 (Figure 15A). In these regions, the length of the season is decreasing with a Sen’s slope of −1.25 days per year (Figure 15B). The Sen’s slope of changes in the length of the season in large parts of Iran shows a negative trend (Figure 15B). However, as mentioned before, these changes are significant only in a small part of the western regions of Iran. According to Figure 16, at location 1 (Hyrcanian ecological zone), the trend is positive and significant, and the increase in the length of the season can be seen with a slope of 1.38 days per year. At location 2, the length of the season is decreasing with a slope of 1.7 days per year. Also, at locations 3 and 4, no significant trends are observed in the length of the season.
The length of the growing season is extracted from the time between the start and end of the growing season. Therefore, like the time component of the end of the growing season, factors such as the occurrence of severe droughts can be effective in creating anomalies during the growing season. For example, in location 1 in 2000, location 3 in 1986, and location 4 in 1986, 2009, and 2012, the length of the growing season is reduced due to drought.

4. Discussion

The trend of changes in different phenology components in Iran between the years 1982 and 2019 was investigated using AVHRR images. Amplitude 0 or the average level of vegetation showed a markedly downward tendency in Iran’s interior and eastern regions, which are often plains and deserts. The variations in amplitude 0 displayed a favorable and significant trend (p < 0.01) in the Caspian Sea coast and northwest Iran region. Also, the Sen’s slope of the changes in amplitude 1 was negative in the inner and eastern regions of Iran, but positive in the western and northern regions, including the Alborz and Zagros mountain ranges. The reason for the differences in the northern and northeastern regions with the central plains of Iran can be attributed to water availability [36]. In these regions, the greater access of plants to water due to the more humid climate and more abundant water resources led to an overall increase in vegetation. On the other hand, in the inner central regions of Iran, these changes were opposite and different. Other studies have shown that the growth of plants and the ability to produce ecosystems in areas with no limiting factors for plant growth can increase due to global warming [9,51,52].
In the northeastern, western, and northwestern regions of Iran, phase 1 is decreasing from 1982 to 2019 with a trend of −0.9° per year, and the start and the end of the growing season are decreasing by −0.6 and −1.33 days per year, respectively. Therefore, on average, in the entire area of Iran, phase 1 is decreasing at −0.6° per year, and the times for the start and end of the season are decreasing by −0.3 and −0.65, correspondingly. The negative trend of amplitude 0 (annual phase) and the starting time of the season indicate the acceleration and earlier start of various phenological processes in these areas from 1982 to 2019. In the northwestern regions of Iran, which are generally mountainous, the temperature is a more relevant limiting factor for the start of the season than humidity [36]. Increases in temperature lead to an increase in plant growth and ecosystem production ability, and thus increasing and positive trends in amplitude 0 and amplitude 1. On the other hand, it seems that the increase in temperature caused by global warming has led to a significant negative trend (p < 0.01) in the annual phase and the time of the start of the season and the end of the season. This was mainly observed in the northeastern, western, and northwestern regions of Iran, which are generally located in the Alborz and Zagros mountains.
The findings of other studies conducted in Iran support the acceleration of the start of the season and the shortening of its duration [36,37]. In a study to assess the impact of climate change on plant production and phenology, the findings revealed that the start and end of the season occur earlier in Iran’s Tehran Province’s northern parts [53]. According to the findings of another study, the length of the vegetation season in the Iranian provinces of Tehran and Zanjan was shortened by 2.8 and 7.4 days, respectively [54]. Other studies in this field have been developed on various scales, including the global level [29,30], focusing on the temperate zone of the Northern Hemisphere [31], and regions such as the Tibetan Plateau [55], Northeastern China [33,34], Africa [56], Europe [57], and middle and high latitudes in the Northern Hemisphere [32]. Generally, they revealed that the beginning of the season was decreased (earlier), and the time for the end of the season and the length of season were changed (increased or decreased). Acceleration of the growing season and lengthening of the growing season have been attributed to a variety of climate variables, particularly temperature increases due to global warming or water availability [28,58].
Phenological components can also be influenced by human factors, which have not been investigated in this research. However, it should be noted that the area of residential and agricultural land is less than 15% of the total area of Iran (scattered throughout the entire country). Therefore, the effects of human factors in the observed changes, especially in the trends of amplitude 0, annual phase, and time of the beginning and end of the growing season, are not expected to have a relevant impact. Phenological changes can also occur in agricultural lands due to changes in the main crops, but no significant change was noticed over the study period.
The output of agricultural goods and the environmental systems of a nation can suffer from changes in phenological factors. The performance of important crops like wheat, rice, maize, and soybeans has decreased in China as a result of seasonal climate fluctuations between 1980 and 2000 [59]. The negative trend regarding the start of the season from 1982 to 2019 can be an important threat to the productivity of agricultural and horticultural products in Iran. This is because by reducing the time for the start of the season, plants and growing trees would enter the flowering stage earlier. However, sudden cold and frost will cause damage. As a result, it is less likely that many agricultural and horticultural goods will be produced in susceptible locations.

5. Conclusions

In the present study, the trend and rate of changes in phenological components of Iran vegetation were investigated using the NDVI time series from 1982 to 2019. The results of this research showed changes in the various phenological components investigated (i.e., amplitude 0, amplitude 1, phase 1, start of the season, end of the season, and length of the season). The trend of changes, however, differs between regions, notably from the central plains to the northern and western regions of Iran. Phenological changes are greater in the northeastern, western, and northwestern regions of Iran, generally located in the Alborz and Zagros mountain ranges. In these areas, phase 1 and phase 2 decreased with a trend of −0.9 and −1.5° per year, respectively. Also, the time for the start of the season, and the time for the end of the season decreased by −0.6 and −1.33 days per year, respectively. Additionally, the start and end of the seasons are becoming shorter by −0.6 and −1.33 days annually, correspondingly. Significant changes in various phenological components (especially changes in the start of the season and annual phase) in these areas can be a serious threat to the production of agricultural and horticultural products.

Author Contributions

Conceptualization, H.R.G.M. and H.Z.K.; methodology, H.Z.K. and H.R.G.M.; validation, H.R.G.M. and C.S.S.F.; formal analysis, H.Z.K. and H.R.G.M.; writing—original draft preparation, H.Z.K. and H.R.G.M.; writing—review and editing, C.S.S.F. and Z.K.; supervision, H.R.G.M. and Z.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Chen, A.; He, B.; Wang, H.; Huang, L.; Zhu, Y.; Lv, A. Notable shifting in the responses of vegetation activity to climate change in China. Phys. Chem. Earth Parts A/B/C 2015, 87, 60–66. [Google Scholar] [CrossRef]
  2. Cheng, M.; Jin, J.; Zhang, J.; Jiang, H.; Wang, R. Effect of climate change on vegetation phenology of different land-cover types on the Tibetan Plateau. Int. J. Remote Sens. 2018, 39, 470–487. [Google Scholar] [CrossRef]
  3. Brown, C.J.; O’Connor, M.I.; Poloczanska, E.S.; Schoeman, D.S.; Buckley, L.B.; Burrows, M.T.; Richardson, A.J. Ecological and methodological drivers of species’ distribution and phenology responses to climate change. Glob. Chang. Biol. 2016, 22, 1548–1560. [Google Scholar] [CrossRef] [PubMed]
  4. Ahmad, S.; Abbas, Q.; Abbas, G.; Fatima, Z.; Atique-ur-Rehman; Naz, S.; Younis, H.; Khan, R.J.; Nasim, W.; Habib ur Rehman, M.; et al. Quantification of climate warming and crop management impacts on cotton phenology. Plants 2017, 6, 7. [Google Scholar] [CrossRef] [PubMed]
  5. Intergovernmental Panel on Climate Change (IPCC). Impacts, Adaptation, and Vulnerability; Fifth Assessment Report on the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, MA, USA, 2014; p. 688. [Google Scholar]
  6. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol. 2013, 169, 156–173. [Google Scholar] [CrossRef]
  7. Shen, M.; Tang, Y.; Chen, J.; Zhu, X.; Zheng, Y. Influences of temperature and precipitation before the growing season on spring phenology in grasslands of the central and eastern Qinghai-Tibetan Plateau. Agric. For. Meteorol. 2011, 151, 1711–1722. [Google Scholar] [CrossRef]
  8. Zheng, Z.; Zhu, W.; Chen, G.; Jiang, N.; Fan, D.; Zhang, D. Continuous but diverse advancement of spring-summer phenology in response to climate warming across the Qinghai-Tibetan Plateau. Agric. For. Meteorol. 2016, 223, 194–202. [Google Scholar] [CrossRef]
  9. Liu, L.; Liang, L.; Schwartz, M.D.; Donnelly, A.; Wang, Z.; Schaaf, C.B.; Liu, L. Evaluating the potential of MODIS satellite data to track temporal dynamics of autumn phenology in a temperate mixed forest. Remote Sens. Environ. 2015, 160, 156–165. [Google Scholar] [CrossRef]
  10. Liu, Q.; Fu, Y.H.; Zhu, Z.; Liu, Y.; Liu, Z.; Huang, M.; Janssens, I.A.; Piao, S. Delayed autumn phenology in the Northern Hemisphere is related to change in both climate and spring phenology. Glob. Chang. Biol. 2016, 22, 3702–3711. [Google Scholar] [CrossRef]
  11. Peñuelas, J.; Filella, I. Responses to a warming world. Science 2001, 294, 793–795. [Google Scholar] [CrossRef]
  12. Guo, L.; Dai, J.; Wang, M.; Xu, J.; Luedeling, E. Responses of spring phenology in temperate zone trees to climate warming: A case study of apricot flowering in China. Agric. For. Meteorol. 2015, 201, 1–7. [Google Scholar] [CrossRef]
  13. Saxe, H.; Cannell, M.G.; Johnsen, Ø.; Ryan, M.G.; Vourlitis, G. Tree and forest functioning in response to global warming. New Phytol. 2001, 149, 369–399. [Google Scholar] [CrossRef] [PubMed]
  14. Pudas, E.; Leppälä, M.; Tolvanen, A.; Poikolainen, J.; Venäläinen, A.; Kubin, E. Trends in phenology of Betula pubescens across the boreal zone in Finland. Int. J. Biometeorol. 2008, 52, 251–259. [Google Scholar] [CrossRef]
  15. Han, H.; Bai, J.; Ma, G.; Yan, J. Vegetation Phenological Changes in Multiple Landforms and Responses to Climate Change. ISPRS Int. J. Geoinf. 2020, 9, 111. [Google Scholar] [CrossRef]
  16. Wang, L.; Chen, S.; Li, D.; Wang, C.; Jiang, H.; Zheng, Q.; Peng, Z. Estimation of paddy rice nitrogen content and accumulation both at leaf and plant levels from UAV hyperspectral imagery. Remote Sens. 2021, 13, 2956. [Google Scholar] [CrossRef]
  17. Zheng, J.; Fu, H.; Li, W.; Wu, W.; Yu, L.; Yuan, S.; Tao, W.Y.W.; Pang, T.K.; Kanniah, K.D. Growing status observation for oil palm trees using Unmanned Aerial Vehicle (UAV) images. ISPRS J. Photogramm. Remote Sens. 2021, 173, 95–121. [Google Scholar] [CrossRef]
  18. Workie, T.G.; Debella, H.J. Climate change and its effects on vegetation phenology across ecoregions of Ethiopia. Glob. Ecol. Conserv. 2017, 13, e00366. [Google Scholar] [CrossRef]
  19. Zare Khormizie, H.; Ghafarian Malamiri, H. Effect of height and temperature on plant phenological processes using harmonic analysis of MODIS NDVI time series (Case study: Shirkouh, Yazd province). Iran. J. Remote Sens. GIS 2020, 12, 1–22. [Google Scholar] [CrossRef]
  20. Pellerin, M.; Delestrade, A.; Mathieu, G.; Rigault, O.; Yoccoz, N.G. Spring tree phenology in the Alps: Effects of air temperature, altitude and local topography. Eur. J. For. Res. 2012, 131, 1957–1965. [Google Scholar] [CrossRef]
  21. Piao, S.; Friedlingstein, P.; Ciais, P.; Viovy, N.; Demarty, J. Growing season extension and its impact on terrestrial carbon cycle in the Northern Hemisphere over the past 2 decades. Glob. Biogeochem. Cycles 2007, 21, GB3018. [Google Scholar] [CrossRef]
  22. Dubovyk, O.; Landmann, T.; Erasmus, B.F.; Tewes, A.; Schellberg, J. Monitoring vegetation dynamics with medium resolution MODIS-EVI time series at sub-regional scale in southern Africa. Int. J. Appl. Earth Obs. Geoinf. 2015, 38, 175–183. [Google Scholar] [CrossRef]
  23. Zhou, J.; Cai, W.; Qin, Y.; Lai, L.; Guan, T.; Zhang, X.; Zheng, Y. Alpine vegetation phenology dynamic over 16 years and its covariation with climate in a semi-arid region of China. Sci. Total Environ. 2016, 572, 119–128. [Google Scholar] [CrossRef] [PubMed]
  24. Zare Khormizi, H.; Hosseini, S.Z.; Mokhtari, M.H.; Ghafarian Malamiri, H.R. Reconstruction of MODIS NDVI Time Series using Harmonic AN analysis of Time Series algorithm (HANTS). J. Spat. Plan. 2017, 21, 221–255. [Google Scholar]
  25. Ghafarian Malamiri, H.R.; Zare, H.; Rousta, I.; Olafsson, H.; Izquierdo Verdiguier, E.; Zhang, H.; Mushore, T.D. Comparison of Harmonic Analysis of Time Series (HANTS) and Multi-Singular Spectrum Analysis (M-SSA) in Reconstruction of Long-Gap Missing Data in NDVI Time Series. Remote Sens. 2020, 12, 2747. [Google Scholar] [CrossRef]
  26. Atkinson, P.M.; Jeganathan, C.; Dash, J.; Atzberger, C. Inter-comparison of four models for smoothing satellite sensor time-series data to estimate vegetation phenology. Remote Sens. Environ. 2012, 123, 400–417. [Google Scholar] [CrossRef]
  27. De Beurs, K.M.; Henebry, G.M. Spatio-temporal statistical methods for modelling land surface phenology. In Phenological Research; Springer: Dordrecht, The Netherland, 2010; pp. 177–208. [Google Scholar]
  28. Tang, H.; Li, Z.; Zhu, Z.; Chen, B.; Zhang, B.; Xin, X. Variability and climate change trend in vegetation phenology of recent decades in the Greater Khingan Mountain area, Northeastern China. Remote Sens. 2015, 7, 11914–11932. [Google Scholar] [CrossRef]
  29. Julien, Y.; Sobrino, J.A. Global land surface phenology trends from GIMMS database. Int. J. Remote Sens. 2009, 30, 3495–3513. [Google Scholar] [CrossRef]
  30. Khormizi, H.Z.; Malamiri, H.R.G.; Alian, S.; Stein, A.; Kalantari, Z.; Ferreira, C.S. Proof of evidence of changes in global terrestrial biomes using historic and recent NDVI time series. Heliyon 2023, 9, e18686. [Google Scholar] [CrossRef]
  31. Jeong, S.J.; Ho, C.H.; Gim, H.J.; Brown, M.E. Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982–2008. Glob. Chang. Biol. 2011, 17, 2385–2399. [Google Scholar] [CrossRef]
  32. Zhao, J.; Zhang, H.; Zhang, Z.; Guo, X.; Li, X.; Chen, C. Spatial and temporal changes in vegetation phenology at middle and high latitudes of the Northern Hemisphere over the past three decades. Remote Sens. 2015, 7, 10973–10995. [Google Scholar] [CrossRef]
  33. Yu, X.; Wang, Q.; Yan, H.; Wang, Y.; Wen, K.; Zhuang, D.; Wang, Q. Forest phenology dynamics and its responses to meteorological variations in Northeast China. Adv. Meteorol. 2014, 2014, 592106. [Google Scholar] [CrossRef]
  34. Yu, L.; Liu, T.; Bu, K.; Yan, F.; Yang, J.; Chang, L.; Zhang, S. Monitoring the long term vegetation phenology change in Northeast China from 1982 to 2015. Sci. Rep. 2018, 7, 14770. [Google Scholar] [CrossRef] [PubMed]
  35. Zhang, X.; Du, X.; Hong, J.; Du, Z.; Lu, X.; Wang, X. Effects of climate change on the growing season of alpine grassland in Northern Tibet, China. Glob. Ecol. Conserv. 2020, 23, e01126. [Google Scholar] [CrossRef]
  36. Zare Khormizi, H.; Ghafarian Malamiri, H. Investigation of phenological components changes of Iranian vegetation in response to climate change using NDVI products of AVHRR sensor from 1982 to 2018. J. RS GIS Nat. Res. 2020, 11, 87–113. [Google Scholar]
  37. Ghafarian Malamiri, H.; Zare Khormizi, H. Investigating vegetation changes in Iran using NDVI time series of NOAA-AVHRR sensor and Harmonic ANalysis of Time Series (HANTS). Sci.-Res. Q. Geogr. Data 2020, 29, 141–158. [Google Scholar]
  38. Verhoef, W. Application of Harmonic Analysis of NDVI Time Series (HANTS). In Fourier Analysis of Temporal NDVI in Southern Africa and America Continent; Azzali, S., Menenti, M., Eds.; Report 108; DLO Winand Staring Centre: Wageningen, The Netherlands, 1996; pp. 19–24. [Google Scholar]
  39. Jönsson, P.; Eklundh, L. TIMESAT—A program for analyzing time-series of satellite sensor data. Comput. Geosci. 2004, 30, 833–845. [Google Scholar] [CrossRef]
  40. Vermote, E.; Justice, C.; Csiszar, I.; Eidenshink, J.; Myneni, R.; Baret, F.; Masuoka, E.D.; Masuoka, E.; Wolfe, R.E.; Claverie, M. NOAA CDR Program: NOAA Climate Data Record (CDR) of Normalized Difference Vegetation Index (NDVI), Version 4; NOAA National Climatic Data Center: Asheville, NC, USA, 2014. [Google Scholar]
  41. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring vegetation systems in the Great Plains with ERTS. In Third Earth Resources Technology Satellite-1 Symposium: Section A-B. Technical Presentations; NASA SP-351 I; NASA: Washington, DC, USA, 1973; pp. 309–317. [Google Scholar]
  42. Fourier, J. Note relative aux vibrations des surfaces e’lastiques et au movement desondes. Bull. Sci. Socie’te’philomatique 1818, 1, 126–136. [Google Scholar]
  43. Ghafarian Malamiri, H.; Zare Khormizie, H. Reconstruction of cloud-free time series satellite observations of land surface temperature (LST) using harmonic analysis of time series algorithm (HANTS). J. RS GIS Nat. Res. 2017, 8, 37–55. [Google Scholar]
  44. Julien, Y.; Sobrino, J.A. Optimizing and comparing gap-filling techniques using simulated NDVI time series from remotely sensed global data. Int. J. Appl. Earth Obs. Geoinf. 2019, 76, 93–111. [Google Scholar] [CrossRef]
  45. Roerink, G.J.; Menenti, M. Reconstructing cloudfree NDVI composites using Fourier analysis of time series. Int. J. Remote Sens. 2000, 21, 1911–1917. [Google Scholar] [CrossRef]
  46. Zare Khormizie, H.; Hosseini, S.Z.; Mokhtari, M.H.; Ghafarian Malamiri, H.R. Analysis of relationship between drought and NDVI variations in different vegetation types (Case study: Southern rangelands of Yazd Province). J. Arid Biome 2017, 7, 85–101. [Google Scholar] [CrossRef]
  47. Mann, H.B. Nonparametric tests against trend. Econom. J. Econ. Soc. 1945, 13, 245–259. [Google Scholar] [CrossRef]
  48. Kendall, M. Rank Correlation Methods; Griffin: London, UK, 1975. [Google Scholar]
  49. Theil, H. A rank invariant method of linear and Polynomial regression analysis, Part 3. Netherlands Akad. Wetensch. Proc. 1950, 53, 1379–1412. [Google Scholar]
  50. Sen, P. Estimates of the regression coefficient based on Kendall’s tau. J. Am. Stat. Assoc. 1968, 63, 1379–1389. [Google Scholar] [CrossRef]
  51. Eastman, J.R.; Sangermano, F.; Machado, E.A.; Rogan, J.; Anyamba, A. Global trends in seasonality of normalized difference vegetation index (NDVI), 1982–2011. Remote Sens. 2013, 5, 4799–4818. [Google Scholar] [CrossRef]
  52. Guay, K.C.; Beck, P.S.; Berner, L.T.; Goetz, S.J.; Baccini, A.; Buermann, W. Vegetation productivity patterns at high northern latitudes: A multi-sensor satellite data assessment. Glob. Chang. Biol. 2014, 20, 3147–3158. [Google Scholar] [CrossRef]
  53. Rayegani, B.; Arzani, H.; Heydari Alamdarloo, E.; Moghadami, M.M. Application of Remote Sens. to assess climate change effects on plant productivity and phenology (Case study area: Tehran Province). J. Rang. 2019, 3, 450–460. [Google Scholar]
  54. Malayeri, F.; Ashourloo, D.; Shakiba, A.; Matkan, A.A.; Aghighi, H. Investigating the Effects of Climate Change on Vegetation Phenology Using AVHRR Time Series Data. J. Agroecol. 2018, 8, 98–117. [Google Scholar]
  55. Zhang, G.; Zhang, Y.; Dong, J.; Xiao, X. Greenup dates in the Tibetan Plateau have continuously advanced from 1982 to 2011. Proc. Natl. Acad. Sci. USA 2013, 110, 4309–4314. [Google Scholar] [CrossRef]
  56. Vrieling, A.; De Leeuw, J.; Said, M.Y. Length of growing period over Africa: Variability and trends from 30 years of NDVI time series. Remote Sens. 2013, 5, 982–1000. [Google Scholar] [CrossRef]
  57. Fu, Y.H.; Piao, S.; Op de Beeck, M.; Cong, N.; Zhao, H.; Zhang, Y.; Menzel, A.; Janssens, I.A. Recent spring phenology shifts in western Central E urope based on multiscale observations. Glob. Ecol. Biogeogr. 2014, 23, 1255–1263. [Google Scholar] [CrossRef]
  58. Forkel, M.; Migliavacca, M.; Thonicke, K.; Reichstein, M.; Schaphoff, S.; Weber, U.; Carvalhais, N. Codominant water control on global interannual variability and trends in land surface phenology and greenness. Glob. Chang. Biol. 2015, 21, 3414–3435. [Google Scholar] [CrossRef] [PubMed]
  59. Zhang, Z.; Song, X.; Tao, F.; Zhang, S.; Shi, W. Climate trends and crop production in China at county scale, 1980 to 2008. Theor. Appl. Climatol. 2016, 123, 291–302. [Google Scholar] [CrossRef]
Figure 1. Iran study area along with Digital Elevation Model (SRTM- 90 m), and location of the four locations selected to show the trends of different phenology components.
Figure 1. Iran study area along with Digital Elevation Model (SRTM- 90 m), and location of the four locations selected to show the trends of different phenology components.
Remotesensing 15 04468 g001
Figure 2. Methodological framework used in this research.
Figure 2. Methodological framework used in this research.
Remotesensing 15 04468 g002
Figure 3. The curve that the HANTS algorithm fit to an NDVI three-year daily time series with missing, outlier, and noise data shown in red, including representation of several phenology variables, such as the beginning, end, length of the season, and seasonal amplitude (A, B, C, and D, respectively). This is a time series of deciduous Hyrcanian forests in northern Iran (one pixel).
Figure 3. The curve that the HANTS algorithm fit to an NDVI three-year daily time series with missing, outlier, and noise data shown in red, including representation of several phenology variables, such as the beginning, end, length of the season, and seasonal amplitude (A, B, C, and D, respectively). This is a time series of deciduous Hyrcanian forests in northern Iran (one pixel).
Remotesensing 15 04468 g003
Figure 4. Hypothetical wave and its harmonic components (A) and different concepts of amplitude, phase, and phase difference (B).
Figure 4. Hypothetical wave and its harmonic components (A) and different concepts of amplitude, phase, and phase difference (B).
Remotesensing 15 04468 g004
Figure 5. Significance level of the changing trend of amplitude 0 (A) and the Sen’s slope of amplitude 0 changes (B).
Figure 5. Significance level of the changing trend of amplitude 0 (A) and the Sen’s slope of amplitude 0 changes (B).
Remotesensing 15 04468 g005
Figure 6. The changing trend of amplitude 0 at four different locations identified in Figure 1, from 1982 to 2019.
Figure 6. The changing trend of amplitude 0 at four different locations identified in Figure 1, from 1982 to 2019.
Remotesensing 15 04468 g006
Figure 7. Significance level of the changing trend of amplitude 1 (A) and changes in the Sen’s slope of amplitude 1 (B).
Figure 7. Significance level of the changing trend of amplitude 1 (A) and changes in the Sen’s slope of amplitude 1 (B).
Remotesensing 15 04468 g007
Figure 8. The changing trend of amplitude 1 at four different locations identified in Figure 1, from 1982 to 2019.
Figure 8. The changing trend of amplitude 1 at four different locations identified in Figure 1, from 1982 to 2019.
Remotesensing 15 04468 g008
Figure 9. Significance level of the changing trend of phase 1 (A) and changes in the Sen’s slope of phase 1 (B).
Figure 9. Significance level of the changing trend of phase 1 (A) and changes in the Sen’s slope of phase 1 (B).
Remotesensing 15 04468 g009
Figure 10. The changing trend of phase 1 at four different locations (see location in Figure 1) from 1982 to 2019.
Figure 10. The changing trend of phase 1 at four different locations (see location in Figure 1) from 1982 to 2019.
Remotesensing 15 04468 g010
Figure 11. Significance level of the changing trend of the time for the start of the season (A) and changes in the Sen’s slope regarding the time for the start of the season (B).
Figure 11. Significance level of the changing trend of the time for the start of the season (A) and changes in the Sen’s slope regarding the time for the start of the season (B).
Remotesensing 15 04468 g011
Figure 12. Changing trend on the time for the start of the season at four different locations (see location in Figure 1) from 1982 to 2019.
Figure 12. Changing trend on the time for the start of the season at four different locations (see location in Figure 1) from 1982 to 2019.
Remotesensing 15 04468 g012
Figure 13. Significance level of the changing trend of the time for the end of the season (A) and changes in the Sen’s slope of the time for the end of the season (B).
Figure 13. Significance level of the changing trend of the time for the end of the season (A) and changes in the Sen’s slope of the time for the end of the season (B).
Remotesensing 15 04468 g013
Figure 14. Changing trend of time for the end of the season at four different locations (see location in Figure 1) from 1982 to 2019.
Figure 14. Changing trend of time for the end of the season at four different locations (see location in Figure 1) from 1982 to 2019.
Remotesensing 15 04468 g014
Figure 15. Significance level of the changing trend of the length of the season (A) and the Sen’s slope of the changes in the length of the season (B).
Figure 15. Significance level of the changing trend of the length of the season (A) and the Sen’s slope of the changes in the length of the season (B).
Remotesensing 15 04468 g015
Figure 16. Changing trend of the length of the season at four different locations (see location in Figure 1) from 1982 to 2019.
Figure 16. Changing trend of the length of the season at four different locations (see location in Figure 1) from 1982 to 2019.
Remotesensing 15 04468 g016
Table 1. Parameters that the HANTS algorithm uses to rebuild NDVI images.
Table 1. Parameters that the HANTS algorithm uses to rebuild NDVI images.
ParametersRate
Valid data range0–1
Base period365 images
Number of Frequencies (NOF)3
Fit Error Tolerance (FET)0.1
Direction of outliers Low
Degree of Over-Determinedness (DOD)10
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Khormizi, H.Z.; Malamiri, H.R.G.; Kalantari, Z.; Ferreira, C.S.S. Trend of Changes in Phenological Components of Iran’s Vegetation Using Satellite Observations. Remote Sens. 2023, 15, 4468. https://doi.org/10.3390/rs15184468

AMA Style

Khormizi HZ, Malamiri HRG, Kalantari Z, Ferreira CSS. Trend of Changes in Phenological Components of Iran’s Vegetation Using Satellite Observations. Remote Sensing. 2023; 15(18):4468. https://doi.org/10.3390/rs15184468

Chicago/Turabian Style

Khormizi, Hadi Zare, Hamid Reza Ghafarian Malamiri, Zahra Kalantari, and Carla Sofia Santos Ferreira. 2023. "Trend of Changes in Phenological Components of Iran’s Vegetation Using Satellite Observations" Remote Sensing 15, no. 18: 4468. https://doi.org/10.3390/rs15184468

APA Style

Khormizi, H. Z., Malamiri, H. R. G., Kalantari, Z., & Ferreira, C. S. S. (2023). Trend of Changes in Phenological Components of Iran’s Vegetation Using Satellite Observations. Remote Sensing, 15(18), 4468. https://doi.org/10.3390/rs15184468

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

Article Metrics

Back to TopTop