Exploring the Applicability and Scaling Effects of Satellite-Observed Spring and Autumn Phenology in Complex Terrain Regions Using Four Different Spatial Resolution Products

Abstract

The information on land surface phenology (LSP) was extracted from remote sensing data in many studies. However, few studies have evaluated the impacts of satellite products with different spatial resolutions on LSP extraction over regions with a heterogeneous topography. To bridge this knowledge gap, this study took the Loess Plateau as an example region and employed four types of satellite data with different spatial resolutions (250, 500, and 1000 m MODIS NDVI during the period 2001–2020 and ~10 km GIMMS3g during the period 1982–2015) to investigate the LSP changes that took place. We used the correlation coefficient (r) and root mean square error (RMSE) to evaluate the performances of various satellite products and further analyzed the applicability of the four satellite products. Our results showed that the MODIS-based start of the growing season (SOS) and end of the growing season (EOS) were highly correlated with the ground-observed data with r values of 0.82 and 0.79, respectively (p < 0.01), while the GIMMS3g-based phenology signal performed badly (r < 0.50 and p > 0.05). Spatially, the LSP that was derived from the MODIS products produced more reasonable spatial distributions. The inter-annual averaged MODIS SOS and EOS presented overall advanced and delayed trends during the period 2001–2020, respectively. More than two-thirds of the SOS advances and EOS delays occurred in grasslands, which determined the overall phenological changes across the entire Loess Plateau. However, both inter-annual trends of SOS and EOS derived from the GIMMS3g data were opposite to those seen in the MODIS results. There were no significant differences among the three MODIS datasets (250, 500, and 1000 m) with regard to a bias lower than 2 days, RMSE lower than 1 day, and correlation coefficient greater than 0.95 (p < 0.01). Furthermore, it was found that the phenology that was derived from the data with a 1000 m spatial resolution in the heterogeneous topography regions was feasible. Yet, in forest ecosystems and areas with an accumulated temperature ≥10 °C, the differences in phenological phase between the MODIS products could be amplified.

1. Introduction

Land surface phenology (LSP) has been recognized as one of the most effective indicators of climate change [1,2,3,4] and is closely related to animal migration, gross primary production, and crop productivity [5,6,7]. Methods for measuring phenology include ground observations (i.e., PhenoCam network and phenology network) and satellite observations [8,9,10]. Ground observations usually only reflect the phenological information of local ecological communities [11,12,13]. However, satellite remote sensing has the potential to continuously observe the variation in vegetation phenology at multiple scales [14,15,16,17].

Phenology has been widely monitored in different types of remote sensing data in attempts to understand the interactions between vegetation and climate change during the past few decades [18,19,20,21,22]. A variety of vegetation indexes were used to monitor vegetation phenology in previous studies [23,24,25]. The two most used vegetation indices are the enhanced vegetation index (EVI) and the normalized difference vegetation index (NDVI) [26,27,28,29]. Meanwhile, several sets of freely accessible remote sensing products with different resolutions were released, such as the third generation GIMMS (GIMMS3g) NDVI with a spatial resolution of 1/12 degree that is derived from the Advanced Very-High-Resolution Radiometer (AVHRR) series satellites [30,31]; the Systeme Probatoire d’Observation de la Tarre (SPOT) NDVI with a 1 km spatial resolution [32]; and the moderate-resolution imaging spectroradiometer (MODIS) NDVI with spatial resolutions of 250, 500, or 1000 m [33,34]. Among them, GIMMS3g NDVI is the latest and longest-used product, and MODIS NDVI has several different spatial resolutions. These products have been widely used for studies involving phenology extraction [15]. Recent studies showed that the temporal and spatial variation trends of vegetation phenology that were observed in some areas by GIMMS3g NDVI and MODIS NDVI were consistent, but the conclusions that were reached in other regions were opposite [35,36]. The topography has a significant impact on the phenology of different product identifications [37]. Furthermore, the spatial phenological heterogeneity of data with a different resolution increases with the increase in landscape fragmentation [5]. However, the impacts of satellite products with different spatial resolutions on LSP extraction over regions with a heterogeneous topography have not been well clarified.

In general, data with a finer spatial resolution possess more information about the seasonality and phenology properties of vegetation [38,39,40]. Meanwhile, data with a fine spatial resolution have the problem of providing a larger amount of data, having a slow computation speed, and being time-consuming [41]. The coarse spatial resolution data is more suitable for monitoring phenology at a landscape scale [42]. Moreover, there will be some differences in vegetation phenology that are estimated from the coarse spatial resolution data to the fine spatial resolution data [5]. Thus, it is important to select remote sensing products with an appropriate spatial resolution in order to investigate vegetation phenology changes.

In this study, we explored the applicability of 250, 500, and 1000 m MODIS NDVI and GIMMS3g NDVI across the entire Loess Plateau, which is a typically ecological fragile region with a heterogeneous topography. The aims of this study were to (1) investigate the spatial and temporal patterns of vegetation phenology in the Loess Plateau, (2) analyze the applicability of different types of satellite data in complex terrain regions, and (3) explore the factors that influence the differences in the phenology of multiple datasets.

2. Materials and Methods

2.1. Study Area

The Loess Plateau lies in the north of China, covering an area of 62.4 × 10⁴ km² (Figure 1). The region is dominated by a continental monsoon climate. The annual accumulated temperature ≥10 °C (annual AT10) increases from ~50 °C·d in the high-elevation western part to ~5800 °C·d in the southern part. Additionally, the annual mean precipitation varies from 50 mm in the northwest to 700 mm in the southeast. The terrain of the Loess Plateau varies significantly, and the altitude ranges from 80 m in the southeast to 5200 m in the west. Land use data with a spatial resolution of 500 m from the MODIS land cover product (MCD12Q1) show that there are seven land cover types in the entire study area; these are grasslands (68.5%), croplands (21.2%), forests (5.1%), barren regions (3.3%), urban and built-up areas (1.8%), and water bodies (0.1%).

2.2. Data Resources and Preprocessing

In this study, we used four types of satellite data: MODIS NDVI with spatial resolutions of 250, 500, and 1000 m and GIMMS3g NDVI with a spatial resolution of ~10 km. These data have been widely applied for vegetation phenology extraction at regional to global scales [43,44,45]. MODIS data for the period 2001–2020 were acquired from NASA (ftp://ladsweb.nascom.nasa.gov/allData/6/, accessed on 15 April 2021). MODIS NDVI with a temporal resolution of 16 days is a gridded level 3 product. GIMMS3g NDVI was obtained from the NASA Earth Exchange platform for the period 1982–2015 (https://nex.nasa.gov/nex/, accessed on 15 April 2021). The spatial and temporal resolution of the GIMMS3g NDVI was 15 days and 1/12 degree (~10 km), respectively. Both the MODIS and GIMMS3g time-series datasets were used to identify the start of the growing season (SOS) and the end of the growing season (EOS).

The daily mean air temperature with a spatial resolution of 0.25 degrees for the period 1982–2020 was obtained from the Copernicus Climate Change Service Climate Data Store (CDS) (https://cds.climate.copernicus.eu/#!/home, accessed on 19 June 2021). The daily air temperature was used to obtain AT10 in order to analyze the interaction between phenology and temperature.

The annual MODIS land cover type product at a 500 m spatial resolution was used to examine the influence of vegetation type on the LSP. The classification scheme that is used for the product is the International Geosphere–Biosphere Program (IGBP). Additionally, the year used for the latest MODIS land cover type data was 2019. The land cover type product for the period 2001–2019 was obtained from NASA’s Land Processes Distributed Active Archive Center (LP DAAC) (https://lpdaac.usgs.gov/products/mcd12q1v006/, accessed on 19 June 2021). The spatial resolution of the product was resampled to 250 m, 1000 m, and 1/12 degree using nearest-neighbor interpolation. The Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) dataset was achieved at the CGIAR (https://srtm.csi.cgiar.org/, accessed on 11 July 2021). The DEM dataset was used to understand the distribution of the topography in the Loess Plateau. The terrain of the central and western parts of the Loess Plateau is heterogeneous, while the terrain of the Weihe Plain is flat [46].

The 271 detailed ground-observed records in the Loess Plateau were obtained from 14 agro-meteorological stations from 2001 to 2013 (a small set of data with low quality was rejected). The distributions and information of the ground-observed sites can be found in Figure 1a and

Table 1. The description of the ground-observed sites.

Site Name	Longitude (°E)	Latitude (°N)	Altitude (m)	Data Range
Fengxiang	107.38	34.51	779	2001–2013
Yongshou	108.15	34.70	1006	2001–2013
Wugong	108.22	34.25	429	2001–2013
Xianyang	108.71	34.40	473	2001–2013
Changan	108.92	34.15	435	2001–2013
Lintong	109.23	34.40	418	2001–2013
Weinan	109.46	34.50	357	2001–2013
Baishui	109.58	34.95	482	2001–2013
Hancheng	110.45	35.46	446	2001–2013
Ruicheng	110.71	34.70	503	2001–2013
Wanrong	110.83	35.40	609	2001–2013
Yuncheng	111.02	35.03	380	2001–2013
Linfen	111.50	36.06	450	2001–2013
Jincheng	112.83	35.51	726	2001–2013

. The observed records were used to validate the vegetation phenology that was derived from MODIS data and more information could be used as prior knowledge of the research area [47].

2.3. Methods

The quality of the NDVI time series was first examined based on the QA information. The LSP was not produced if three serial periods of NDVI data were contaminated by clouds. Second, to reduce the impacts of noise from cloud contamination or other poor atmospheric conditions, the MODIS NDVI time series were smoothed using the modified Savitsky–Golay algorithm (mSG) with the help of a specific MODIS data layer named “composite day of the year” [19,48]. Similarly, the GIMMS3g NDVI data were also smoothed. However, the GIMMS3g data lacked the layer of “composite day of the year”; we thus regarded the 1st and 16th days of each month as the “day of year (DOY)” for each image. The mSG algorithm is a simple but robust method that is based on the Savitsky–Golay algorithm [48,49]. Finally, the smoothed NDVI growth curve was used to estimate the SOS and EOS with the following logistic model [34]:

$$y\left(t\right)=\frac{c}{1+{\mathrm{e}}^{a+bt}}+d,$$

(1)

where t is the DOY, y(t) represents the NDVI value at time t, a and b are the fitting parameters, c is the difference between the maximum and minimum NDVI values, and d is the initial background vegetation index value. Next, the SOS and EOS were produced from the rate of change in curvature:

$$K=-\frac{b^{2}cz\left(1-z\right){\left(1+z\right)}^3}{{[{\left(1+z\right)}^4+{\left(bcz\right)}^2]}^{3/2}}$$

(2)

$$K^{\prime }=b^{3}cz\left\{ \frac{3z\left(1-z\right){\left(1+z\right)}^3[2{\left(1+z\right)}^3+b^2{c}^{2}z]}{{[{\left(1+z\right)}^4+{\left(bcz\right)}^2]}^{5/2}}-\frac{{\left(1+z\right)}^2\left(1+2z-5z^2\right)}{{[{\left(1+z\right)}^4+{\left(bcz\right)}^2]}^{\frac{3}{2}}} \right\}$$

(3)

where K represents the curvature, $z=e^{a+bt}$, and K′ is the rate of change of K.

In order to acquire the deviation and the correlation characteristics between different products, the root mean square error (RMSE) and correlation coefficient (r) were calculated using Equations (4) and (5), respectively:

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(X_i-\overline{X}\right)}^2}{n}}$$

(4)

where $\overline{X}$ is the mean value of X and $\left(X_i-\overline{X}\right)$ represents the deviation value, that is, the bias;

$$\mathrm{r}=\frac{{{\displaystyle }}_{i=1}^n\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(X_i-\overline{X}\right)}^2{{\displaystyle \sum }}_{i=1}^n{\left(Y_i-\overline{Y}\right)}^2}}$$

(5)

where $X_i$ and $Y_i$ correspond to two different datasets.

Based on the geographical location information of agro-meteorological stations, we first extracted the SOS and EOS pixels (e.g., 3 × 3 homogeneous pixels) around the center of each station, and then averaged these values to obtain the mean SOS and EOS of each station [50]. Finally, we employed RMSE and r to investigate the correlation between the satellite-based and ground-observed phenology.

3. Results

3.1. The Performances of Satellite-Based SOS and EOS

This study first verified the performances of the SOS and EOS that were produced from the 250 m MODIS NDVI and GIMMS3g NDVI against the ground-observed data from 2001 to 2013. The results showed that a good agreement was observed between the MODIS-derived phenology and ground-observed data, where the r values of the MODIS SOS and EOS were 0.83 and 0.79, respectively (p < 0.01, Figure 2). However, the phenology that was estimated from GIMMS3g NDVI showed poor consistency. The r values of the GIMMS3g SOS and EOS were only 0.49 and 0.21, and both the p-values were more than 0.05.

3.2. Spatial Patterns of Vegetation Phenology

Figure 3 presents the SOS and EOS of the Loess Plateau in 2001. The spatial patterns of the vegetation phenology that were produced from the GIMMS3g data were quite different from those of the MODIS data. However, three sets of MODIS phenology results with different resolutions gave a similar spatial distribution (Figure A1). The spring phenology of croplands in the Weihe Plain was the earliest, with an average SOS of DOY 46. Most SOSs in the grassland region of the central Loess Plateau were later than DOY 170. In particular, the SOS that was derived from MODIS NDVI in the Mu Us Sandy Land region was earlier than that in the surrounding areas. The early EOS in the Loess Plateau was mainly concentrated in the southern region, while the late EOS was mainly concentrated in the western region and the Lvliang mountains. Among the land cover types, the EOS of almost 80% of the croplands was between DOY 260 and 290, the EOS date of the forests was the earliest, and that of the grasslands was the latest.

Compared with the MODIS-derived phenology, the spatial distribution of the phenology period that was identified by the GIMMS3g data was more concentrated. Moreover, the SOS values that were derived from the GIMMS3g product were largely consistent with the phenology that was derived from the MODIS data, which was mainly distributed in the Weihe Plain. The SOS that was produced by the GIMMS3g NDVI data in the northern part of the Loess Plateau was concentrated in DOY 90–110, and the spatial details were greatly neglected. Additionally, there was a significant difference between the GIMMS3g EOS and the MODIS EOS in the south–central region of the Loess Plateau. The difference ranged from 20 days to more than 60 days.

Figure 4a,b presents the raw and smoothed NDVI in heterogeneous areas, where the land use types are grassland and forest, respectively. Figure 4c shows that the time series of NDVI in the Weihe Plain, where the land use type was cropland. Based on the time-series data, we found that no matter what year it was, the inflection point of the time series of the original GIMMS3g data was always concentrated in the seventh or eighth period’s data in heterogeneous areas. However, this phenomenon did not occur in flat areas. This indicates that the problem of GIMMS3g data is one of the important reasons for the spatially aggregated distribution of GIMMS3g phenology.

In order to compare the differences between the GIMMS3g-derived and MODIS-derived phenology in flat areas, we calculated the differences of the two datasets in the Weihe Plain during the period 2001–2015 (Figure 5). The findings showed that the differences in SOS (GIMMS3g SOS–250 m MODIS SOS) were mainly less than 5 days and 10 days in 43.49% and 70.66%. The GIMMS3g data performed well at monitoring SOS over the flat areas. In addition, we found that the frequency with which the differences in SOS were greater than 20 days decreased significantly and became close to 0. In the results showing the SOS differences greater than 25 days, the value of the GIMMS3g SOS was always greater than that of the 250 m MODIS SOS. This shows that GIMMS3g SOS tended to be later than MODIS SOS. However, differences in the EOS (GIMMS3g EOS–250 m MODIS EOS) that were less than 5 days or 10 days were only 13.85% and 29.88%. In addition, for the entire Loess Plateau, the proportion of differences in EOS between the GIMMS3g and 250 m MODIS data that were less than 5 days was still less than 20%. This indicated that the consistency of the GIMMS3g EOS and the MODIS EOS was poor, even in flat areas.

3.3. Temporal Variation in Vegetation Phenology

The inter-annual trends of the SOS and EOS during the period 1982–2020 are presented in Figure 6. The results showed that the inter-annual trends of the SOS and EOS from GIMMS3g were the reverse of those of the MODIS data. The SOS showed a trend of postponing and the EOS presented a trend of advancing from 1982 to 2015 based on the GIMMS3g data. In addition, the SOS (EOS) of the GIMMS3g data delay (advance) trend during the period 2001–2015 was more significant, and the trend lines K of the SOS and EOS were 0.34 and −0.14, respectively. In contrast, the SOS (EOS) showed an advanced (delayed) trend based on the MODIS data, and the trend line K of the MODIS SOS and EOS was −0.63 and 0.19, respectively, during the period 2001–2015.

In the comparison of MODIS products, our findings showed that the phenological periods that were derived from MODIS products with various spatial resolutions gave only small differences. The average difference between the 500 m MODIS SOS (EOS) and 250 m MODIS SOS (EOS) was only 1.2 (0.3) days. The correlation coefficient and RMSE between the 500 m MODIS and 250 m MODIS results were greater than 0.99 and less than 0.60, respectively (

Table 2. The correlation coefficient and RMSE that were used for the comparison of the vegetation phenology between different MODIS products.

	500 m MODIS–250 m MODIS			1000 m MODIS–250 m MODIS
	Bias	Correlation Coefficient	RMSE	Bias	Correlation Coefficient	RMSE
SOS	−1.2	0.9977 **	0.5665	−1.7	0.9957 **	0.9101
EOS	0.3	0.9952 **	0.3176	1.4	0.9775 **	0.8241

The symbol ** indicates significance at the 0.01 level.

). Moreover, the average difference between the 1000 m MODIS SOS (EOS) and 250 m MODIS slightly increased to 1.7 (1.4) days. The correlation coefficient and RMSE between the 1000 m MODIS and 250 m MODIS results were greater than 0.95 and approximately equal to 1.0, respectively. This demonstrated that there was little difference between the 1000 m MODIS NDVI and the 250 m MODIS NDVI. Therefore, the 1000 m MODIS NDVI could be used to monitor the LSP of the Loess Plateau. However, the GIMMS3g product may not be able to accurately monitor heterogeneous areas, such as the Loess Plateau.

Figure 7 shows the spatial trend of the phenology over the Loess Plateau that was calculated from the MODIS data during the period 2001–2020, which passed the significance test of α = 0.05. The area with a significantly advanced SOS was about 16.7 × 10⁴ km², which was about a quarter of the area of the Loess Plateau. Areas with a significantly advanced SOS were mainly in the central and northeastern regions of the Loess Plateau. The area of the delayed SOS was only one-third that of the advanced SOS. Meanwhile, the delayed SOS was mainly distributed across the croplands of the Weihe and Hetao Plains. Moreover, the area with a significantly delayed EOS was about 9.3 × 10⁴ km², which was more than three times that of the area with an advanced EOS. More than two-thirds of the advanced SOS and delayed EOS occurred in grasslands, which determined the overall phenological changes in the Loess Plateau. Although the land cover types in the Loess Plateau changed dramatically, the phenological changes were still dominated by the areas where the land use did not change. Only about 20% of the SOS or EOS significant trends occurred in areas with changes in the land cover type.

3.4. Impact Factors on MODIS Products

3.4.1. Influences of Vegetation on MODIS Products

To examine the possible causes of the difference in LSP from the MODIS data, we investigated the differences in the SOSs and EOSs between MODIS products that were used on different land cover types (Figure 8). Our findings showed that the largest difference between the 250 m MODIS SOS and the 1000 m MODIS SOS (1000 m MODIS SOS–250 m MODIS SOS) in forests was 3.5 days, which was larger than the difference found in grasslands (1.9 days) and croplands (0.6 days). Additionally, the standard deviation of the inter-annual difference between the 250 m and 1000 m MODIS SOS in forests was 1.1 days, which was the largest value that was obtained among all vegetation types. The differences that were obtained between the 250 m MODIS EOS and the 1000 m MODIS EOS were 0.9 days (forests), 1.2 days (grasslands), and 1.1 days (croplands). The standard deviation of the inter-annual difference was the largest in forests. In addition, we found that the differences in the SOS between the 250 m MODIS and 500 m MODIS (500 m MODIS SOS–250 m MODIS SOS) were both less than one day.

3.4.2. Influences of AT10 on the Phenology

Figure 9 presents the variations in the differences that were obtained between multiple sets of the SOS and EOS with an AT10 from 1 January to 30 April and from 1 September to 31 October, respectively. As the AT10 from January to April increased, the difference between each MODIS SOS gradually increased. In particular, the difference between the 1000 m MODIS SOS and the 250 m MODIS SOS (1000 m MODIS SOS–250 m MODIS SOS) was greater than the difference between the 500 m MODIS SOS and the 250 m MODIS SOS (500 m MODIS SOS–250 m MODIS SOS). The standard deviation of the 1000 m MODIS SOS was greater than that of the 500 m MODIS SOS. The relationship between the EOS and AT10 was found to be opposite to that of the SOS and AT10. Our results showed that with the increase in AT10 from September to October, the difference between each MODIS EOS gradually decreased. No matter whether the SOS or EOS were used, the difference in vegetation phenology between the 250 m and 1000 m products was less than 3 days. Additionally, we found that there was almost no difference between the 250 m and 500 m vegetation phenology.

4. Discussion

4.1. Difference between Satellite-Based LSP and Observations

Based on the MODIS NDVI and the layer named “composite day of the year,” we calculated three sets of MODIS phenology with different spatial resolutions from 2001 to 2020. The GIMMS3g data lacked a “composite day of the year” layer, therefore we regarded the 1st and 16th days of each month as the DOY and calculated the GIMMS3g phenology from 1982 to 2015. We used the ground-observed data during the period 2001–2013 to verify the performances of the phenology that was derived from the MODIS and GIMMS3g data. Our results showed that the MODIS-derived phenology and ground-observed data had a good agreement, but that the correlation between the GIMMS3g-derived phenology and ground-observed data was very bad. This indicates that the SOS and EOS that were identified by the MODIS NDVI were more robust, but that the GIMMS3g phenology performed badly.

Moreover, the MODIS SOS that was identified by the logistics model was earlier than the SOS that was observed on the ground, while the MODIS EOS was later than the ground-observed data. Similarly, previous studies also found that in the north of China, the SOS that was identified by the logistics model based on SPOT satellite data was earlier than the SOS observed on the ground [49]. We suggest that this phenomenon was caused by the phenological recognition algorithm and the time resolution of the data used. The finer the temporal resolution of the image is, the more accurate the identified phenological phase will be [51,52,53,54]. However, among the existing vegetation index products, product data with a high time resolution are still limited. Therefore, an alternative suggestion is that using a proper method may make up for the lack of advanced or delayed phenological phases, such as the cumulative NDVI [49,55]. However, it is also worth noting that different methods will cause different problems.

The agreement between the EOS and the ground-observed data was not as good as for the SOS, especially GIMMS3g EOS. A previous study found that the NDVI at harvest time will be increased due to the noise-reduction algorithm [56]. The land use type of the Weihe Plain is mainly cropland, and the locations of the agro-meteorological stations are distributed around the Weihe Plain. This may lead to worse accuracy for the EOS in cultivated land than for the SOS. In addition, due to the limitation of data from observation stations, the phenological results that were gained in other areas cannot be verified. Compared with other studies, we found that the MODIS phenology was similar to those of other studies in terms of their spatial patterns [49,57].

4.2. Comparisons of Different Product Data

Effectiveness and using the smallest possible amount of data are matters that must be considered first in experiments. In this study, we compared four sets of remote sensing product data with different resolutions for the extraction of phenology in the Loess Plateau. Notably, the three sets of MODIS results with different resolutions showed good consistency. From the correlation coefficient and RMSE, we found that there was little difference between the 1000 m MODIS and 250 m MODIS results. This indicated that lower-resolution data could achieve the same effect as relatively higher-resolution data. Moreover, computers could process the 1000 m MODIS data much faster. Therefore, the comprehensive performance of the 1000 m MODIS data was better. If errors within the range of 1–3 days are allowed, we can use the 1000 m MODIS NDVI rather than the 250 m MODIS NDVI in future phenological studies.

In contrast, the GIMMS3g products were less effective in the Loess Plateau. From the time-series phenological results of the GIMMS3g products, we observed that the phenology was delayed in spring and advanced in autumn. However, the phenological results based on the MODIS products showed the opposite situation. As we know, many studies have confirmed that global warming advances spring phenology and delays autumn phenology [58,59,60]. This indicates that the phenological results obtained using the GIMMS3g data may contain errors in complex-terrain regions, such as the Loess Plateau. Additionally, we found that the strong spatial homogeneity of the GIMMS3g NDVI was the main reason for these phenological differences. Except for croplands, from the time series of the original data, no matter what year was used, the inflection point of the GIMMS3g data was always concentrated in the seventh or eighth period, and this led to the maximum curvature of the smoothed timing signal focus in this period. Therefore, the SOS that was identified by the GIMMS3g NDVI data in the northern part of the Loess Plateau was concentrated in DOY 90–110 (Figure A2).

However, previous studies showed that coarse-resolution SOS was comparable with finer-resolution SOS in homogeneous areas [5,42]. In this study, we calculated the difference between the MODIS and GIMMS3g results in the Weihe Plain during the period 2001–2015. Similarly, our findings showed that the GIMMS3g data performed well in monitoring the SOS over flat areas. Furthermore, the cropping intensity in a large area of the Weihe Plain during 1982–2013 changed from cropping twice a year to a single cropping taking place each year [61]. Due to the effectiveness of the GIMMS3g image in the Weihe Plain, the GIMMS3g SOS also showed a delayed trend that was the same as the MODIS phenological trend. Except for the flat area, the large heterogeneous areas in the Loess Plateau were affected by the original GIMMS3g data. Therefore, the original GIMMS3g data from the heterogeneous area incorrectly identified the phenological trend in the whole Loess Plateau.

Moreover, although the inflection point did not occur at the end of the vegetation growth of the raw GIMMS3g data, the proportion of the differences in EOSs between the GIMMS3g and MODIS data that were less than 5 days was still less than 20% for the entire Loess Plateau. This means that the GIMMS3g data had a weak ability to predict the SOS and EOS in areas with complex terrain, while they could better monitor the change in the SOS in relatively flat areas. If the GIMMS3g product cannot be replaced in an experiment on phenology production, other phenology estimation methods may be used, such as the dynamic threshold method [62,63,64]. Based on the maximum and minimum NDVI each year, the dynamic threshold method was used to determine the SOS and EOS in spring and autumn with the threshold ratio. This method could reduce the deviation in phenological estimation that is caused by the mutation of the temporal signal to a certain extent [65,66].

4.3. Factors for the Differences from MODIS Products

The differences in vegetation phenology that were determined from MODIS products with different spatial resolutions were mainly due to the land-cover types and temperatures involved. This finding can also be seen in other regions and ecosystems [67,68,69,70]. It is usually the case that when the same time resolution of the image is used, the finer the spatial resolution is, the more accurate the phenology properties will be. For example, in forest ecosystems, the vegetation phenology that is identified by MODIS images with a 1000 m resolution is more variable than that identified using images with a 500 m resolution. The main reason for this may be that the structures and functions of forest ecosystems are more complex than those of other natural ecosystems, such as cropland ecosystems [71,72]. Croplands are homogeneous and mainly affected by crop management. In contrast, forests are usually controlled by multiple environmental factors. Previous studies also suggested that the performance of coarse-resolution images over homogeneous areas is better than that over other regions [42,73,74].

Temperature is the main factor that leads to the advancement of spring phenology and the delay of autumn phenology [75,76,77]. Further, we found that the higher the AT10 during the early growth season was, the greater the differences and variability between the SOSs were. Meanwhile, the lower the AT10 during the late growth season was, the greater the differences and volatility of the EOSs were. Areas with a high AT10 from January to April were mainly distributed south of the Weihe Plain. We suggested that the main reason for this phenomenon was that there are more types of ecosystems in the southeast of the Loess Plateau, while the northwest of the plateau has a relatively homogenous ecosystem. Due to differences in the sensitivity of various types of vegetation to AT10, the SOS in areas with a relatively high AT10 showed greater differences. In addition, areas with a low AT10 from September to October were mainly concentrated in the southwest of the Loess Plateau. The huge elevation fluctuation in this area may be the reason for the phenological differences that were seen in products with different spatial resolutions. Although the AT10 and vegetation type or terrain had an impact on the data, the maximum averaged differences of the SOS and EOS between the 250 m MODIS products and the 1000 m MODIS products were less than three days. Additionally, the phenological difference remained within an acceptable range.

5. Conclusions

Based on the 250, 500, and 1000 m MODIS data during the period 2001–2020 and the ~10 km GIMMS3g data during the period 1982–2015, as well as the rate of change in the curvature of the logistic models, this study investigated the applicability and spatial scaling effects of various remote sensing products with different spatial resolutions on phenology extraction in a complex-terrain region. Our study showed that the MODIS products performed better in phenology analysis across the Loess Plateau, and the phenology results that were derived from the different MODIS products showed only small differences. However, the GIMMS3g-based SOSs that were derived from logistic models had a good performance in the flat region (i.e., the Weihe Plain) but a poor performance in regions with a more heterogeneous topography. Additionally, the performances of the GIMMS3g-based EOSs across the whole Loess Plateau were poor. The phenology results that were derived from the MODIS data presented advanced SOS trends and delayed EOS trends during the period 2001–2020 for the entire Loess Plateau. However, both the SOS and EOS trends that were identified by the GIMMS3g products were the opposite. Our finding emphasized that the 1000 m MODIS product can be used to extract phenology from areas with a complex terrain, such as the Loess Plateau, and almost no difference was found in the phenology extraction between the 500 m MODIS product and the 250 m MODIS products. Furthermore, we also investigated the effects of vegetation and AT10 on the spatiotemporal variability of vegetation phenology, which could help us to understand the driving factors of such phenological changes in the future.