*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(t) = \frac{c}{1 + e^{a + bt}} + d,\tag{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 c z (1 - z)(1 + z)^3}{\left[ (1 + z)^4 + (b c z)^2 \right]^{3/2}} \tag{2}$$

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

where *K* represents the curvature, *z* = *<sup>e</sup>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{\sum\_{i=1}^{n} \left(X\_i - \overline{X}\right)^2}{n}} \tag{4}$$

where *X* is the mean value of *X* and *Xi* − *X* represents the deviation value, that is, the bias;

$$\mathbf{r} = \frac{\stackrel{\text{\tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \tiny \raisebox{\text{\tiny \tiny \tiny \tiny \raisebox{\text{\tiny \tiny \tiny \tiny \raisebox{\text{\tiny \tiny \tiny \raisebox{1em}{1em}1em}{1em}}}}{\sqrt{\sum\_{i=1}^{n}{\left(\boldsymbol{X}\_{i}-\boldsymbol{X}\right)^{2}\sum\_{i=1}^{n}{\left(\boldsymbol{Y}\_{i}-\boldsymbol{Y}\right)^{2}}}}}{\sqrt{\sum\_{i=1}^{n}{\left(\boldsymbol{X}\_{i}-\boldsymbol{X}\right)^{2}\sum\_{i=1}^{n}{\left(\boldsymbol{Y}\_{i}-\boldsymbol{Y}\right)^{2}}}}} \tag{5}$$

where *Xi* and *Yi* 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.
