Analysis of Differences in Phenology Extracted from the Enhanced Vegetation Index and the Leaf Area Index

Remote-sensing phenology detection can compensate for deficiencies in field observations and has the advantage of capturing the continuous expression of phenology on a large scale. However, there is some variability in the results of remote-sensing phenology detection derived from different vegetation parameters in satellite time-series data. Since the enhanced vegetation index (EVI) and the leaf area index (LAI) are the most widely used vegetation parameters for remote-sensing phenology extraction, this paper aims to assess the differences in phenological information extracted from EVI and LAI time series and to explore whether either index performs well for all vegetation types on a large scale. To this end, a GLASS (Global Land Surface Satellite Product)-LAI-based phenology product (GLP) was generated using the same algorithm as the MODIS (Moderate Resolution Imaging Spectroradiometer)-EVI phenology product (MLCD) over China from 2001 to 2012. The two phenology products were compared in China for different vegetation types and evaluated using ground observations. The results show that the ratio of missing data is 8.3% for the GLP, which is less than the 22.8% for the MLCD. The differences between the GLP and the MLCD become stronger as the latitude decreases, which also vary among different vegetation types. The start of the growing season (SOS) of the GLP is earlier than that of the MLCD in most vegetation types, and the end of the growing season (EOS) of the GLP is generally later than that of the MLCD. Based on ground observations, it can be suggested that the GLP performs better than the MLCD in evergreen needleleaved forests and croplands, while the MLCD performs better than the GLP in shrublands and grasslands.


Introduction
Phenology is the study of the timing of recurring biological events, the causes of their timing with regard to biotic and abiotic forces, and the interactions among phases of the same or different species [1]. Sprouting and flowering of plants in spring, color changes of leaves in fall, bird migration and nesting, insect hatching, and animal hibernation are all examples of phenological events [2]. Phenology is an important indicator of global change and the carbon cycle because of its direct effects on vegetation photosynthesis, carbon sequestration and land-atmosphere water and energy exchange [3]. Vegetation phenology indicates the responses of plants to seasonal and interannual variations of climate, hydrology, soil and anthropogenic factors [4]. Vegetation phenology also affects the climate system by influencing the seasonality of albedo; surface roughness length; canopy conductance; and fluxes of water, energy, CO 2 and biogenic volatile organic compounds [5]. The extraction of those from the MODIS-EVI and -NDVI, the timing of the start, end, and length of season differed by approximately one month at the global scale. The higher uncertainties appear in areas with limited seasonality expressed in the satellite signal and systematic biases due to the differences in the methodologies and datasets. Some other studies have also shown significant discrepancies in vegetation phenological metrics derived from time series of different vegetation parameters [50][51][52]. However, it is not clear which data sources and vegetation parameters are best suited for extracting vegetation phenological metrics accurately.
While EVI and LAI are the most widely used vegetation parameters that can be used for remote-sensing phenological extraction, this paper aims to assess the differences in phenological information extracted using the EVI and LAI time series and to determine whether either EVI or LAI time series performs well for all vegetation types over a large scale. The objective is to compare the phenology detection results in China from the MODIS-EVI and the GLASS-LAI. To achieve this objective, the GLASS-LAI phenology product (GLP) was generated using the same algorithm as the MODIS-EVI phenology product (MLCD) in China from 2001 to 2012. The two phenology products were compared in different vegetation types in China and then evaluated using ground observations from forests, grasslands and croplands. Our work lays the foundation for uniting multisource data and for improving remote-sensing phenology products in the future.

Phenology Extraction Based on the EVI and LAI
The MLCD product was used as the vegetation phenology extracted from the EVI for analysis in the paper. The MLCD provides yearly 500-m resolution global estimates of vegetation phenology [47] and is available for the public at https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_ table/mcd12q2. The MLCD identifies phenophase transition dates based on the logistic functions fit-to-time series of the EVI, which were calculated from composite 8-day normalized BRDF (Bidirectional Reflectance Distribution Function)-adjusted reflectance data [26]. Specifically, a time series of the EVI is assembled for each pixel, the data undergo a gap-filling and smoothing process, periods of sustained EVI increase or decrease are identified, logistic models are fit to the time series, and transition dates are identified as local maxima and minima in the rate of curvature change of the fitted logistic function [22]. The MLCD has four layers corresponding to key phenological stages: onset greenness increase, onset greenness maximum, onset greenness decrease, and onset greenness minimum. MLCD images over China from 2001 to 2012 were acquired for the study. Due to its different spatial resolution from the GLP, the MLCD was upscaled to 1000 m by averaging the value from every four pixels for statistics of different values between GLP and MLCD. But in other analyses, MLCD was used with its original spatial resolution.
The vegetation phenology extracted from the LAI was generated from the GLASS-LAI (http://glass-product.bnu.edu.cn/) that provides yearly 1000-m resolution global estimates of vegetation phenology, which is called the GLASS-LAI-based phenology product (GLP). GLASS-LAI was estimated from MODIS and CYCLOPES LAI time series using general regression neural networks, with 8-day temporal resolution and 1-km spatial resolution. Based on multitemporal data, GLASS-LAI are spatially and temporally continuous with no gaps in time series and no high-frequency temporal variability. As a result of cloud, atmospheric and snow contamination, Savitzky-Golay filtering was used to smooth and gap-fill the CYCLOPES LAI. The MODIS reflectance was reprocessed to remove remaining effects of cloud contamination and other factors [53]. The GLP applied the same algorithm as MLCD (the piecewise logistic model) to fit the LAI time series, which is a flexible, repeatable, and realistic means of monitoring seasonal and interannual vegetation dynamics from remote-sensing data at large scales. The piecewise logistic equation is represented as Equation (1): where t is time in days, y(t) is the LAI value at time t, a and b are fitting parameters, c + d is the maximum LAI value, and d is the initial background LAI value. GLP products with 1-km resolution were generated over China from 2001 to 2012. For both the MLCD and GLP products, the vegetation dynamics can be divided into four stages: greenup, maturity, senescence and dormancy [22], as shown in the Figure 1. The key transition dates define the main phases of vegetation dynamics, including the start of the growing season (SOS), the start of the peak (PS), the end of the peak (PE) and the end of the growing season (EOS): (1) SOS is the date of onset of photosynthetic activity; (2) PS is the date when the greenness reaches its maximum; (3) PE represents the date when the photosynthetic activity and the greenness begin to decrease rapidly; and (4) EOS is the date when the physiological activity approaches zero [14,54].
Sensors 2017, 17, 1982 4 of 16 where t is time in days, y(t) is the LAI value at time t, a and b are fitting parameters, c + d is the maximum LAI value, and d is the initial background LAI value. GLP products with 1-km resolution were generated over China from 2001 to 2012. For both the MLCD and GLP products, the vegetation dynamics can be divided into four stages: greenup, maturity, senescence and dormancy [22], as shown in the Figure 1. The key transition dates define the main phases of vegetation dynamics, including the start of the growing season (SOS), the start of the peak (PS), the end of the peak (PE) and the end of the growing season (EOS): (1) SOS is the date of onset of photosynthetic activity; (2) PS is the date when the greenness reaches its maximum; (3) PE represents the date when the photosynthetic activity and the greenness begin to decrease rapidly; and (4) EOS is the date when the physiological activity approaches zero [14,54].

Remote-Sensing Auxiliary Datasets
The nadir BRDF-adjusted reflectance (NBAR) data from MODIS (MCD43A4) were used to compute the EVI, as shown in Equation (2). The MODIS-EVI products and the GLASS-LAI products were used to compare the difference of vegetation growing states reflected by the LAI and EVI time series. Besides, the land cover type (MOD12Q1) was downloaded from NASA's Earth Observing System Data Center website. The MOD12Q1 was used to identify the land cover type over pixels, as shown in Figure 2.
where , and are the reflectances of near infrared, red and blue bands respectively.

Ground Phenological Observations
The ground phenological observations were taken as the true values to evaluate the accuracy of the phenology products. However, the ground-based phenology was observed in a number of individual plants, while the remote-sensing phenology represented the integrated phenological characteristics of a plant community in one pixel. Ground validation of remote-sensing measurements with coarse resolution entails considerable difficulties. To improve the reliability of the statistical analysis based on ground observations, ground phenological observations meeting the following requirements were selected from the China Meteorological Administration (CMA) and the Chinese Ecosystem Research Network (CERN): (1) Spatial representation. The poor relationship between ground and satellite phenology due to data-scale issues is a drawback of satellite phenology because of the small chance of a singlepoint ground observation being representative of an entire area at the remote-sensing scale (typically ≥1 km in remote-sensing phenological studies) [4,35]. Consequently, the phenological homogeneity and subdued topography of field sites must be ensured in comparison with the remote-sensing data [41]. The phenological homogeneity requires the phenophases of dominant

Remote-Sensing Auxiliary Datasets
The nadir BRDF-adjusted reflectance (NBAR) data from MODIS (MCD43A4) were used to compute the EVI, as shown in Equation (2). The MODIS-EVI products and the GLASS-LAI products were used to compare the difference of vegetation growing states reflected by the LAI and EVI time series. Besides, the land cover type (MOD12Q1) was downloaded from NASA's Earth Observing System Data Center website. The MOD12Q1 was used to identify the land cover type over pixels, as shown in Figure 2.
where ρ N IR , ρ RED and ρ BLUE are the reflectances of near infrared, red and blue bands respectively.

Ground Phenological Observations
The ground phenological observations were taken as the true values to evaluate the accuracy of the phenology products. However, the ground-based phenology was observed in a number of individual plants, while the remote-sensing phenology represented the integrated phenological characteristics of a plant community in one pixel. Ground validation of remote-sensing measurements with coarse resolution entails considerable difficulties. To improve the reliability of the statistical analysis based on ground observations, ground phenological observations meeting the following requirements were selected from the China Meteorological Administration (CMA) and the Chinese Ecosystem Research Network (CERN): (1) Spatial representation. The poor relationship between ground and satellite phenology due to data-scale issues is a drawback of satellite phenology because of the small chance of a single-point ground observation being representative of an entire area at the remote-sensing scale (typically ≥1 km in remote-sensing phenological studies) [4,35]. Consequently, the phenological homogeneity and subdued topography of field sites must be ensured in comparison with the remote-sensing data [41]. The phenological homogeneity requires the phenophases of dominant species at one site be similar (less than 30 days), and the topography is checked by Google Earth to avoid sites in mountains as far as possible. The dominant species were selected after considering the distribution and quantity of the community based on the instruction files from CERN and CMA. (2) Data integrity. The selected ground sites should have phenological phase continuity and few missing records.
Sensors 2017, 17, 1982 5 of 16 species at one site be similar (less than 30 days), and the topography is checked by Google Earth to avoid sites in mountains as far as possible. The dominant species were selected after considering the distribution and quantity of the community based on the instruction files from CERN and CMA. (2) Data integrity. The selected ground sites should have phenological phase continuity and few missing records. The ground phenological observations were selected according to the above requirements, as shown in Table 1. Each site corresponds to a typical vegetation type, including evergreen needleleaved forest (ENF), evergreen broadleaved forest (EBF), deciduous needleleaved forest (DNF), deciduous broadleaved forest (DBF), shrub, grass and crop. Rice was chosen as a representative crop with two sites, because the SY (Shenyang) site from CERN can provide harvest and JT (Jiutai) site from CMA can provide grain-filling. The mean phenological dates of the dominant species were calculated and treated representative of this region. The distribution of the ground phenological observation sites is shown in Figure 2.  The ground phenological observations were selected according to the above requirements, as shown in Table 1. Each site corresponds to a typical vegetation type, including evergreen needleleaved forest (ENF), evergreen broadleaved forest (EBF), deciduous needleleaved forest (DNF), deciduous broadleaved forest (DBF), shrub, grass and crop. Rice was chosen as a representative crop with two sites, because the SY (Shenyang) site from CERN can provide harvest and JT (Jiutai) site from CMA can provide grain-filling. The mean phenological dates of the dominant species were calculated and treated representative of this region. The distribution of the ground phenological observation sites is shown in Figure 2.

Methods for Evaluating the Phenology Products
To analyze the differences in the GLP and the MLCD, the mean phenological metrics for each pixel in each year were calculated as the spatial pattern of China. Furthermore, the mean and standard deviation of the difference values between the two products were calculated based on land cover type data. In addition, the mean bias (MB) and mean absolute bias (MAB) were used to evaluate the remote-sensing phenological metrics based on ground-station observations [55], as shown in Equations (3) and (4).
where rs is the remote-sensing phenological date; obs is the ground observation; and N is the number of rs or obs.

Methods for Evaluating the Phenology Products
To analyze the differences in the GLP and the MLCD, the mean phenological metrics for each pixel in each year were calculated as the spatial pattern of China. Furthermore, the mean and standard deviation of the difference values between the two products were calculated based on land cover type data. In addition, the mean bias (MB) and mean absolute bias (MAB) were used to evaluate the remote-sensing phenological metrics based on ground-station observations [55], as shown in Equations (3) and (4).
where rs is the remote-sensing phenological date; obs is the ground observation; and N is the number of rs or obs.

Missing Data
In general, the ratio of missing data is 8.3% for the GLP, which is obviously less than the 22.8% for the MLCD, as shown in Figure 4. The missing ratios of the MLCD and the GLP differ obviously over different vegetation types. For the MLCD product, the missing ratio is highest in EBF (53.6%), followed by ENF (48.73%), which is much higher than that of GLP (6.2% and 4.0% respectively). Nevertheless, the missing ratio of shrub is the highest for GLP (55.3%), which is higher than that of MLCD (14.9%). In addition, the missing ratio of MLCD in crop is relatively high, especially for the double-cycle crop in the North China Plain ( Figure 2); invalid results are also treated as missing data.

Missing Data
In general, the ratio of missing data is 8.3% for the GLP, which is obviously less than the 22.8% for the MLCD, as shown in Figure 4. The missing ratios of the MLCD and the GLP differ obviously over different vegetation types. For the MLCD product, the missing ratio is highest in EBF (53.6%), followed by ENF (48.73%), which is much higher than that of GLP (6.2% and 4.0% respectively). Nevertheless, the missing ratio of shrub is the highest for GLP (55.3%), which is higher than that of MLCD (14.9%). In addition, the missing ratio of MLCD in crop is relatively high, especially for the double-cycle crop in the North China Plain ( Figure 2); invalid results are also treated as missing data. As shown in Figure 5, the patterns of the missing ratios depend strongly on the latitude between 22° N and 53° N in China. In general, the missing ratios of the MLCD are higher than those of the GLP along each latitude zone except from 37° N to 43° N. The missing ratio of the MLCD gradually decreases as the latitude decreases, especially in midlatitude zones (30° N-53° N), and remained inordinately high between 22° N and 30° N (approximately 50%), where much data are missing in the Sichuan Basin (as shown in Figure 2). However, the missing ratio of GLP largely depends on desert vegetation. For the GLP, the latitudinal distribution of missing data is consistent with the distribution of grass or shrub, as shown in Figures 5 and 6.   In general, the ratio of missing data is 8.3% for the GLP, which is obviously less than the 22.8% for the MLCD, as shown in Figure 4. The missing ratios of the MLCD and the GLP differ obviously over different vegetation types. For the MLCD product, the missing ratio is highest in EBF (53.6%), followed by ENF (48.73%), which is much higher than that of GLP (6.2% and 4.0% respectively). Nevertheless, the missing ratio of shrub is the highest for GLP (55.3%), which is higher than that of MLCD (14.9%). In addition, the missing ratio of MLCD in crop is relatively high, especially for the double-cycle crop in the North China Plain ( Figure 2); invalid results are also treated as missing data. As shown in Figure 5, the patterns of the missing ratios depend strongly on the latitude between 22° N and 53° N in China. In general, the missing ratios of the MLCD are higher than those of the GLP along each latitude zone except from 37° N to 43° N. The missing ratio of the MLCD gradually decreases as the latitude decreases, especially in midlatitude zones (30° N-53° N), and remained inordinately high between 22° N and 30° N (approximately 50%), where much data are missing in the Sichuan Basin (as shown in Figure 2). However, the missing ratio of GLP largely depends on desert vegetation. For the GLP, the latitudinal distribution of missing data is consistent with the distribution of grass or shrub, as shown in Figures 5 and 6.   In general, the ratio of missing data is 8.3% for the GLP, which is obviously less than the 22.8% for the MLCD, as shown in Figure 4. The missing ratios of the MLCD and the GLP differ obviously over different vegetation types. For the MLCD product, the missing ratio is highest in EBF (53.6%), followed by ENF (48.73%), which is much higher than that of GLP (6.2% and 4.0% respectively). Nevertheless, the missing ratio of shrub is the highest for GLP (55.3%), which is higher than that of MLCD (14.9%). In addition, the missing ratio of MLCD in crop is relatively high, especially for the double-cycle crop in the North China Plain ( Figure 2); invalid results are also treated as missing data. As shown in Figure 5, the patterns of the missing ratios depend strongly on the latitude between 22° N and 53° N in China. In general, the missing ratios of the MLCD are higher than those of the GLP along each latitude zone except from 37° N to 43° N. The missing ratio of the MLCD gradually decreases as the latitude decreases, especially in midlatitude zones (30° N-53° N), and remained inordinately high between 22° N and 30° N (approximately 50%), where much data are missing in the Sichuan Basin (as shown in Figure 2). However, the missing ratio of GLP largely depends on desert vegetation. For the GLP, the latitudinal distribution of missing data is consistent with the distribution of grass or shrub, as shown in Figures 5 and 6.   There are several factors causing the missing data. The first is incomplete of time-series source data caused by clouds or snow. Cloud cover severely compromises the quality of NBAR EVI time series in extensive areas of tropics and subtropics [56], as a consequence, many phenological metrics are missing in these areas. In addition, the seasonal change characteristics of evergreen forest are not obvious, so it is hard to extract the phenology information of EBF and ENF. While GLASS-LAI exhibits spatial-temporal continuity [53], there is no gaps between GLASS-LAI time series. This is the main reason why the missing ratio of GLP is less than that of MLCD generally, but the missing ratio of GLP in desert vegetation (e.g., shrub) is relatively high because the range of the seasonal variation in the LAI time series is small. Therefore, fitting by the logistics method fails easily.

Difference Comparison between the GLP and the MLCD
Specifically, the difference between the GLP and the MLCD varies among different vegetation types, as shown in Figure 7. Overall, the SOS of the GLP is earlier than that retrieved by the MLCD in most vegetation types, while the EOS of the GLP is later, including shrub. Among the vegetation types, the difference between the GLP and the MLCD in deciduous broadleaved forest is the smallest, among which the difference in SOS is smaller than that of EOS (the means and standard deviations are −7.2 ± 17.7 and 5.5 ± 36.5, respectively). The greatest difference occurs in evergreen broadleaved forest, among which the difference in SOS is larger than that of EOS (the means and standard deviations are −84.3 ± 76.5 and 38.1 ± 73.3, respectively). In addition, the difference in SOS or EOS between the GLP and the MLCD in croplands is relatively large (the means and standard deviations are −14.5 ± 53.6 and 38.8 ± 61.9, respectively), as double-cycle crops are poorly detected in MLCD. There are several factors causing the missing data. The first is incomplete of time-series source data caused by clouds or snow. Cloud cover severely compromises the quality of NBAR EVI time series in extensive areas of tropics and subtropics [56], as a consequence, many phenological metrics are missing in these areas. In addition, the seasonal change characteristics of evergreen forest are not obvious, so it is hard to extract the phenology information of EBF and ENF. While GLASS-LAI exhibits spatial-temporal continuity [53], there is no gaps between GLASS-LAI time series. This is the main reason why the missing ratio of GLP is less than that of MLCD generally, but the missing ratio of GLP in desert vegetation (e.g., shrub) is relatively high because the range of the seasonal variation in the LAI time series is small. Therefore, fitting by the logistics method fails easily.

Difference Comparison between the GLP and the MLCD
Specifically, the difference between the GLP and the MLCD varies among different vegetation types, as shown in Figure 7. Overall, the SOS of the GLP is earlier than that retrieved by the MLCD in most vegetation types, while the EOS of the GLP is later, including shrub. Among the vegetation types, the difference between the GLP and the MLCD in deciduous broadleaved forest is the smallest, among which the difference in SOS is smaller than that of EOS (the means and standard deviations are −7.2 ± 17.7 and 5.5 ± 36.5, respectively). The greatest difference occurs in evergreen broadleaved forest, among which the difference in SOS is larger than that of EOS (the means and standard deviations are −84.3 ± 76.5 and 38.1 ± 73.3, respectively). In addition, the difference in SOS or EOS between the GLP and the MLCD in croplands is relatively large (the means and standard deviations are −14.5 ± 53.6 and 38.8 ± 61.9, respectively), as double-cycle crops are poorly detected in MLCD. Obviously, the differences between the GLP and the MLCD become stronger as the latitude decreases, as shown in Figures 2 and 8. The difference between the GLP and the MLCD largely depends on different vegetation types. The greatest difference between the GLP and the MLCD is found below 26° N, where EBF is mainly distributed. As the ratio of EBF decreases, the difference decreases. Because DBF is mainly distributed above 40° N, the difference between the GLP and the MLCD become small above 40° N. The difference in the SOS below 40° N is larger than that in the EOS, and the differences in the SOS and the EOS are similar above 40° N. Obviously, the differences between the GLP and the MLCD become stronger as the latitude decreases, as shown in Figures 2 and 8. The difference between the GLP and the MLCD largely depends on different vegetation types. The greatest difference between the GLP and the MLCD is found below 26 • N, where EBF is mainly distributed. As the ratio of EBF decreases, the difference decreases. Because DBF is mainly distributed above 40 • N, the difference between the GLP and the MLCD become small above 40 • N. The difference in the SOS below 40 • N is larger than that in the EOS, and the differences in the SOS and the EOS are similar above 40 • N.

Accuracy Assessment Based on Field Data
To analyze the variances and relationships between ground observations and the phenological information extracted from the MODIS-EVI and GLASS-LAI products, the EVI and LAI time series at representative sites were chosen as reference to compare the phenological results, as shown in Figure 9. The information from ground phenology observation sites is presented in Table 1.

Accuracy Assessment Based on Field Data
To analyze the variances and relationships between ground observations and the phenological information extracted from the MODIS-EVI and GLASS-LAI products, the EVI and LAI time series at representative sites were chosen as reference to compare the phenological results, as shown in Figure 9. The information from ground phenology observation sites is presented in Table 1.

Accuracy Assessment Based on Field Data
To analyze the variances and relationships between ground observations and the phenological information extracted from the MODIS-EVI and GLASS-LAI products, the EVI and LAI time series at representative sites were chosen as reference to compare the phenological results, as shown in Figure 9. The information from ground phenology observation sites is presented in Table 1.
(e) (f) For evergreen forests, significant discrepancies occur in the LAI and EVI time series (Figure 9a,b). Obviously, irregular fluctuations exist in the EVI of ENF and EBF, which are more likely due to a pattern of noise caused by clouds in the EVI data structure rather than to vegetation growing changes. Thus, the phenological metrics extracted by MODIS-EVI may be discontinuous, or even missing, especially in the tropical climate zone. Apart from noise, in many areas with evergreen vegetation, the seasonal change characteristics in EVI is too subtle to retrieve phenology.
For deciduous forests, the LAI time series agree with the EVI time series overall, but the LAI time series have fewer fluctuations than the EVI time series in the peak (Figure 8c,d). The differences between LAI and EVI time series in DNF is larger than that in DBF. In DNF, although the LAI start to increase earlier than the EVI, the EVI is more sensitive than the LAI during the greenup. Therefore, the EVI generally reaches maturity faster than the LAI. And also, the EVI start to decrease later than the LAI, but decreases rapidly to the minimum.
The LAI and the EVI of shrublands are relatively lower. Although the range of seasonal variation of the EVI time series is smaller than that of the LAI, the precision of the current LAI products is 0.1, which is easier for presenting the "ladder" time series than the EVI product, whose precision is 0.01. Thus, the EVI shows the subtle dynamic changes of grasslands more clearly. As shown in Figure 9e, the inflection point of the EVI time series is close to bud burst and leaf defoliation of shrub. When using the LAI, it is difficult to show the subtle dynamic changes, causing a time offset in shrub phenology. Especially for shrubs in desert, the ranges of the seasonal variations in EVI or LAI are so small that the logistics model is difficult to fit to extract phenological information.
For croplands, although the LAI start to increase earlier than the EVI, the EVI is more sensitive to vegetation growth than the LAI during the greenup. Therefore, he EVI generally reaches maturity faster than the LAI, but the peak of the LAI time series is closer to the ground observations, as shown in Figure 9g,h.
In order to take the ground phenology observations as the truth values to evaluate the accuracy of the phenology products, the phenophases observed from ground observations corresponding to remote-sensing phenological metrics are determined, as shown in Table 2. All the observed ground phenology stages are the beginning of the phenophases in both CERN and CMA, so the phenophases observed from ground observations that were closest to the remote-sensing phenological metrics were chosen based on Figure 9. According to previous studies, the key remote-sensing phenological metrics might correspond to the four growing stages for crop [57,58]. For herbs and deciduous broadleaf trees, no phenophases correspond to PS and PE. Figure 10 compares the MB and MAB values between the remote-sensing phenological metrics and the ground data at typical vegetation sites.
For evergreen forests, the phenological metrics of MLCD are missing at the HS and DHS sites; only the phenological metrics of GLP are presented, as shown in Figure 10a  For evergreen forests, significant discrepancies occur in the LAI and EVI time series (Figure 9a,b). Obviously, irregular fluctuations exist in the EVI of ENF and EBF, which are more likely due to a pattern of noise caused by clouds in the EVI data structure rather than to vegetation growing changes. Thus, the phenological metrics extracted by MODIS-EVI may be discontinuous, or even missing, especially in the tropical climate zone. Apart from noise, in many areas with evergreen vegetation, the seasonal change characteristics in EVI is too subtle to retrieve phenology.
For deciduous forests, the LAI time series agree with the EVI time series overall, but the LAI time series have fewer fluctuations than the EVI time series in the peak (Figure 8c,d). The differences between LAI and EVI time series in DNF is larger than that in DBF. In DNF, although the LAI start to increase earlier than the EVI, the EVI is more sensitive than the LAI during the greenup. Therefore, the EVI generally reaches maturity faster than the LAI. And also, the EVI start to decrease later than the LAI, but decreases rapidly to the minimum.
The LAI and the EVI of shrublands are relatively lower. Although the range of seasonal variation of the EVI time series is smaller than that of the LAI, the precision of the current LAI products is 0.1, which is easier for presenting the "ladder" time series than the EVI product, whose precision is 0.01. Thus, the EVI shows the subtle dynamic changes of grasslands more clearly. As shown in Figure 9e, the inflection point of the EVI time series is close to bud burst and leaf defoliation of shrub. When using the LAI, it is difficult to show the subtle dynamic changes, causing a time offset in shrub phenology. Especially for shrubs in desert, the ranges of the seasonal variations in EVI or LAI are so small that the logistics model is difficult to fit to extract phenological information.
For croplands, although the LAI start to increase earlier than the EVI, the EVI is more sensitive to vegetation growth than the LAI during the greenup. Therefore, he EVI generally reaches maturity faster than the LAI, but the peak of the LAI time series is closer to the ground observations, as shown in Figure 9g,h.
In order to take the ground phenology observations as the truth values to evaluate the accuracy of the phenology products, the phenophases observed from ground observations corresponding to remote-sensing phenological metrics are determined, as shown in Table 2. All the observed ground phenology stages are the beginning of the phenophases in both CERN and CMA, so the phenophases observed from ground observations that were closest to the remote-sensing phenological metrics were chosen based on Figure 9. According to previous studies, the key remote-sensing phenological metrics might correspond to the four growing stages for crop [57,58]. For herbs and deciduous broadleaf trees, no phenophases correspond to PS and PE. Figure 10 compares the MB and MAB values between the remote-sensing phenological metrics and the ground data at typical vegetation sites.
For evergreen forests, the phenological metrics of MLCD are missing at the HS and DHS sites; only the phenological metrics of GLP are presented, as shown in Figure 10a  For deciduous forests, the accuracy of GLP or MLCD for SOS is higher than that for EOS, as shown in Figure 10c   For deciduous forests, the accuracy of GLP or MLCD for SOS is higher than that for EOS, as shown in Figure 10c,d. The SOS of MLCD in DBF is very close to the 1:1 line, suggesting that the SOS that extracted by EVI is in good agreement with ground observations. Specifically, the MAB of the MLCD for the SOS is 4.1 days (MB = −1.9 days), and the MAB of the GLP for the SOS is 9.0 days (MB = −8.5 days). The EOS of GLP is later than that of MLCD, and their MABs are 22.9 days and 17.3 days respectively.
For shrublands, because LAI is hard to extract the phenology information in desert at SPT (Shapotou) site, only a few phenological metrics of GLP are presented, as shown in Figure 10e

Discussion and Conclusions
This paper assessed the differences in phenological information extracted from the EVI and LAI time series and explored whether the EVI or the LAI time series performs well for all vegetation types over a large scale in China. The GLP and the MLCD were compared in different regions and evaluated based on ground observations. The missing ratio of the MLCD is more than that of the GLP overall. The main reason is incompletion of time-series source data caused by clouds or snow for the MLCD product. While GLP is not affected by clouds or snow, because the GLASS-LAI is provided with spatial-temporal continuity, there are no gaps in time series and no high-frequency temporal variability. However, the missing ratio of GLP in desert vegetation (e.g., shrub) is relatively high. Because the range of seasonal variation of LAI time series is small, fitting by the logistics method easily fails. Besides, in many areas with evergreen vegetation, the annual variation in EVI is too subtle to retrieve phenology.
Accuracy assessment is always an important concern in any remote-sensing-based analysis, especially one performed at a large scale. To improve the reliability of the statistical analysis based on ground observations, we selected sites based on the requirements of spatial representation and data integrity. The accuracy assessment based on representative sites in this paper show that both phenology products can provide the best accuracy with ground observations for the SOS of DBF and the EOS of rice. Because the leaf defoliation period for DBF lasts for one or two months, and the EOS extracted from LAI may be close to the end of leaf defoliation, a lag can occur between the EOS and the start of leaf defoliation from ground observations. Besides, previous research has shown that the end-of-season metrics of the MLCD associated with vegetation senescence and dormancy are highly uncertain [47,56]. In this paper, it is shown that the EOS in croplands has less uncertainties than that in forests. Thus, it can be inferred that the uncertainties of EOS in forest may come from its complex structure and composition, as the composition of cropland is relatively homogeneous. For grasslands or shrublands, the range of seasonal variation is relatively small, the logistics model couldn't fit well with EVI or LAI time series; as a consequence, the accuracy of grasslands or shrublands is worse than that of deciduous forests or croplands. LAI is easier for presenting the "ladder" time series than the EVI and even failing in fitting when the maximum of LAI is less than 0.3. In addition, the growth of grass or shrub in arid regions or desert is sensitive to precipitation, which may result in irregular fluctuations in EVI time series. Our analysis also shows that the GLP performs better than the MLCD in ENF and croplands, especially for the PS and PE of crop. Although GLP has less missing data in EBF, the disagreement exist between vegetation phenological metrics and ground observations. However, the MLCD performs better than the GLP in shrublands and grasslands. The precision of the current LAI product is 0.1, indicating that it does not show subtle dynamic changes as well as the EVI product, whose precision is 0.01.
It is noted that MLCD is derived from 500-m, 8-day composite EVI data; the GLP is derived from 1000-m, 8-day composite LAI data; while ground observations are daily point-based observations. The temporal and spatial scales of ground observations make them difficult to compare directly with remote-sensing retrievals. Although studied sites are the best representation for each station and the surrounding area, mismatches between phenology products and ground-based observations at both temporal and spatial scales are unavoidable [56]. In addition, the monitoring method may lead to a mismatch with the ground observations; therefore, we used a flexible, repeatable, and realistic means (the piecewise logistic models) to monitor seasonal and inter-annual dynamics in vegetation from remote-sensing data at large scales. However, the method may also lead to an offset in some situations when the remote-sensing time series do not match the piecewise logistic models.
Besides, even though based on the same algorithm and mathematical curvature changes for the same spatial/vegetation type, phenological metrics derived from EVI and LAI time series may have different biophysical meanings. It seems that the SOS and EOS that are extracted from LAI are close to bud burst and leaf defoliation respectively, while the SOS and EOS that are extracted from EVI is close to leaf unfolding and leaf coloring respectively at deciduous broadleaved forest site. Thus, the SOS that is extracted from LAI is earlier than that from EVI, while the EOS that is extracted from LAI is later than that from EVI overall according to the analysis results.
In conclusion, both the LAI and the EVI have advantages in different vegetation types. It should be considered that uniting multisource data is an effective way to improve the accuracy and validity of remote-sensing phenology products. However, for evergreen forests and sparse vegetation, extracting phenological metrics is still challenging using the LAI or EVI time series based on piecewise logistic models. New methods or data sources must be adopted for these land cover types. In addition, due to the coarse spatial resolution of current remote-sensing phenology products, other plants are inevitably present within one square kilometer, which may lead to variability in phenological developmental stages and a weak relationship between remote-sensing data and ground observations. These findings are of great value for improving the spatial resolution of remote-sensing phenology products to promote their application and development in the future.