Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Contents lists available at ScienceDirect Int J Appl Earth Obs Geoinformation journal homepage: www.elsevier.com/locate/jag Land surface phenology from VEGETATION and PROBA-V data. Assessment T over deciduous forests Kevin Bórneza,b, , Adrià Descalsa,b, Aleixandre Vergera,b, Josep Peñuelasa,b ⁎ a CREAF, Cerdanyola del, Vallès 08193, Catalonia, Spain b CSIC, Global Ecology Unit, Cerdanyola del Vallès 08193, Catalonia, Spain ARTICLE INFO ABSTRACT Keywords: Land surface phenology has been widely retrieved although no consensus exists on the optimal satellite dataset Land surface phenology and the method to extract phenology metrics. This study is the first comprehensive comparison of vegetation SPOT-VEGETATION variables and methods to retrieve land surface phenology for 1999–2017 time series of Copernicus Global Land PROBA-V products derived from SPOT-VEGETATION and PROBA-V data. We investigated the sensitivity of phenology to Leaf area index (I) the input vegetation variable: normalized difference vegetation index (NDVI), leaf area index (LAI), fraction Ground observations of absorbed photosynthetically active radiation (FAPAR), and fraction of vegetation cover (FCOVER); (II) the smoothing and gap filling method for deriving seasonal trajectories; and (III) the method to extract phenological metrics: thresholds based on a percentile of the annual amplitude of the vegetation variable, autoregressive moving averages, logistic function fitting, and first derivative methods. We validated the derived satellite phenological metrics (start of the season (SoS) and end of the season (EoS)) using available ground observations of Betula pendula, B. alleghaniensis, Acer rubrum, Fagus grandifolia, and Quercus rubra in Europe (Pan-European PEP725 network) and the USA (National Phenology Network, USA-NPN). The threshold-based method applied to the smoothed and gap-filled LAI V2 time series agreed best with the ground phenology, with root mean square errors of ˜10 d and ˜25 d for the timing of SoS and EoS respectively. This research is expected to contribute for the operational retrieval of land surface phenology within the Copernicus Global Land Service. 1. Introduction daily observations that allow improving the characterization of land surface phenology on a global scale (Verger et al., 2016; Atkinson et al., Phenology is the study of the timing of recurrent biological and 2012; Zhang et al., 2004). However, the noise in the data and missing seasonal events and their biotic and abiotic factors (Beaumont et al., observations mainly due to cloud contamination may induce significant 2015). Studies of plant phenology focus on how these events and factors uncertainties in the estimation of phenological metrics (Kandasamy are influenced by seasonal and interannual variations in climate and et al., 2013; Verger et al., 2013). The literature shows a broad variety of how they modulate abundance and diversity (Beaumont et al., 2015). time-series processing methods designed to reconstruct gap-filled ve- Phenology is, moreover, key to control physicochemical and biological getation seasonal trajectories from noisy satellites signals. This includes processes, especially albedo, surface roughness, canopy conductance the best index slope method (Viovy et al., 1992), mean filters (Reed and fluxes of carbon, water and energy (Peñuelas et al., 2009; et al., 1994), moving-window filters (Sweets et al. 1999), asymmetric Richardson et al., 2013). Phenological metrics are thus relevant para- Gaussian functions (Jönsson and Eklundh, 2002), Savitzky–Golay filters meters for modeling land surface processes and the global carbon cycle (Chen et al., 2004) or the Whittaker smoother (Eilers, 2003). However, (Wu et al., 2014). no single method always performs better than others for smoothing Phenological metrics are estimated based on ground observations vegetation time series (Cai et al., 2017) and their performance vary and data derived from satellites. Ground observations provide accurate spatially and temporally with land surface conditions and cloud influ- timing of vegetation phenophases but cannot cover continuously large- ence (Atkinson et al., 2012; Kandasamy and Fernandes, 2015). scale areas (Garrity et al., 2011; Yu et al., 2017). Satellite sensors with A broad variety of statistical methods have been designed to extract moderate spatial resolutions, including AVHRR, MODIS, MERIS, phenological metrics from satellite time series. Metrics typically include SPOT-VEGETATION and PROBA-V, provide long-term time series of the start of the season (SoS), the end of the season (EoS), the timing of ⁎ Corresponding author at: CREAF, Cerdanyola del Vallès, 08193, Catalonia, Spain. E-mail address: k.bornez@creaf.uab.cat (K. Bórnez). https://doi.org/10.1016/j.jag.2019.101974 Received 9 February 2019; Received in revised form 22 August 2019; Accepted 9 September 2019 Available online 03 October 2019 0303-2434/ © 2019 Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/). K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 maximum growth and the length of the growing season (LoS) (Reed We discarded ground sites with less than four yearly measurements et al., 1994; Zhang et al., 2004). de Beurs and Henebry (2010) provided to obtain consistent data records over the time series. We used ground- a comprehensive review of the exiting phenology retrieval approaches site located pixels. The spatial heterogeneity hamper comparing that can be classified in four main categories: thresholds and percentile ground-based phenology for individual plants with satellite phenology based methods (Atzberger and Eilers, 2011; Verger et al., 2016), at a resolution of 1 km. We filtered the ground sites located in agri- moving averages (Reed et al., 1994), first derivatives (White et al., cultural or urban areas using high-resolution images from Google Earth 2009) and fitted models (de Beurs and Henebry, 2005). White et al. (https://earth.google.com/) and the ESA Land Cover Map (CCI-LC) (2009) compared ten different phenology retrieval methods applied to (http://maps.elie.ucl.ac.be/CCI/viewer/index.php). In Europe, we used AVHRR NDVI in North America and found large discrepancies of up to only the points with forest coverages > 5% using a tree cover map for two months in the detection of the SoS. European forests (Brus et al., 2012). In addition to the sensitivity to the smoothing and phenological extraction algorithm, the derived phenological metrics are also de- 2.2. Satellite time series pendent on the sensor, spatial and temporal resolution, processing chain, and satellite data set. The satellite-derived spectral vegetation The time series of satellite imagery used for estimating the pheno- indices (e.g. the Normalized Difference Vegetation Index (NDVI)) vary logical metrics were from the SPOT-VEGETATION (1999–2013) and in their strength of phenological prediction across sites and plant PROBA-V satellites (2014–2017) with spatial resolutions of 1 km and functional types (Wu et al., 2014). Unlike previous studies based on temporal frequencies of 10 d. In particular, we used LAI, FAPAR, vegetation indices, the present study aimed to characterize the phe- FCOVER V1 (Baret et al., 2013) and V2 (Verger et al., 2014) and NDVI nology not only with NDVI but also with biophysical variables: the leaf V2.1 (Toté et al., 2017) products generated within the CGLS (https:// area index (LAI), the fraction of absorbed photosynthetically active land.copernicus.eu/global/themes/vegetation) (Table1). radiation (FAPAR), the fraction of vegetation cover (FCOVER). We used Time series of LAI, FAPAR and FCOVER V1 and NDVI V2.1 con- NDVI version V2.1 (Toté et al., 2017), LAI, FAPAR and FCOVER V1 tained frequent disturbances caused by residual cloud contamination, (Baret et al., 2013) and V2 (Verger et al., 2014) time series derived atmospheric variability, snow and bi-directional effects (Fig. 2). These within the Copernicus Global Land Service (CGLS) from SPOT-VEGE- time series required the application of smoothing and gap filling tech- TATION and PROBA-V data. Verger et al. (2017) showed that the niques to generate consistent and gap filled seasonal trajectories (sec- phenology derived from the interannual climatology of LAI V1 im- tion 2.3) before the extraction of phenological metrics. The LAI and proved other existing products including MODIS-EVI when compared to FCOVER V2 products improved over the V1 products in terms of tem- ground observations for the average date of the SoS and EoS. However, poral consistency and continuity (Fig. 2). The V2 algorithm included their study was limited to the baseline LAI phenology as derived from a multi-step data filtering, smoothing and gap-filling techniques that single extraction method. This paper is a continuation of the previous rendered the products suitable for phenological estimation without paper by Verger et al. (2017) and we address now the interannual additional pre-processing. Filtering of outliers is based on an upper variation of the yearly phenology, the impact of the input vegetation envelope approach, the pixel climatology (interannual mean) is used to variable and the phenological extraction method. Further, we in- fill missing data and a Savitzky-Golay filter is used for the smoothing corporate LAI, FAPAR and FCOVER V2 that improved continuity (no (Verger et al., 2014). missing data in V2) and smoothness as compared to V1. Our study had two main objectives: to select the best biophysical 2.3. Smoothing methods variable or vegetation index for estimating phenological metrics on a global scale within the portfolio of the CGLS vegetation products We tested several smoothing methods for reducing noise and re- (NDVI, LAI, FCOVER or FAPAR) and to define the method that best constructing gap filled seasonal trajectories from CGLS time series matched the ground data. (Eerens and Haesen, 2015): 2. Materials and methods • WHITTAKER smoother (Atzberger and Eilers, 2011): It minimizes a cost function describing the balance between fidelity (quadratic 2.1. Phenological ground observations difference between estimates and actual observations) and rough- ness (quadratic difference between successive estimates). Ground-based phenological data from PEP725 and USA-NPN were • BISE (Best Index Slope Extraction) (Viovy et al., 1992): It retains the examined, focusing on the dates of leaf out and leaf senescence for good observations in a local window and replaces missing or Betula (birch) in Europe (Fig. 1a) and the USA and for Quercus (oak), eliminated suspect values by linear interpolation. Fagus (beech) and Acer (maple) in the USA (Fig. 1b). These genera were • MEAN: A linear interpolation is first applied to fill missing data. A chosen because they are present in both Europe and the USA and have running mean filter with a sliding window of 50 d length is then large numbers of records in the combined data set. applied. The PEP725 Pan-European Phenology database (Templ et al., 2018) • SWETS method (Swets et al., 1999): A linear interpolation is first (www.pep725.eu) has complete records from 1990 to the present. The applied to fill missing data. A weighted linear regression over a local phenophases defined in PEP725 are based on BBCH (Biologische Bun- window is then applied. desanstalt, Bundessortenamt and Chemical industry) code (Meier et al., 2009). We used the phenophases corresponding to the first visible 2.4. Methods for extracting phenological metrics leaves (BBCH 11) as the reference for the timing of SoS and the date corresponding to 50% of leaves with autumn coloration (BBCH 94) for We tested four state of the art methods to extract phenological the timing of EoS. metrics from CGLS time series (Fig. 3): The USA National Phenology Network was established in 2007 to collect, store and share historical and contemporary phenological data • Thresholds based on a pixel percentile value: SoS is defined as the day on a North American scale (Schwartz et al., 2012). The data are freely of the year (DoY) when a vegetation variable exceeds a particular available at https://www.usanpn.org/. This network provides mea- threshold. EoS is defined as the DoY when an index remains below a surements of several phenophases. We used the phenophases “leaves” particular threshold. We established dynamic thresholds per pixel which corresponds to first visible leaves and “increasing leaf size” for based on a percentile of the annual amplitude of the vegetation SoS and “colored leaves” for EoS. variable (Verger et al., 2016). The selected percentiles were 2 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Fig. 1. Location of selected phenological ground observations at (a) PEP725 sites for Silver birch (Betula pendula) and (b) USA-NPN sites for Red maple (Acer rubrum), Yellow birch (Betula alleghaniensis), American beech (Fagus grandifolia) and Northern red oak (Quercus rubra). Table 1 Algorithm principles of NDVI V2.1, and LAI, FAPAR, FCOVER V1 and V2. NDVI LAI, FAPAR, FCOVER Version 1 LAI, FAPAR, FCOVER Version 2 Version 2.1 Inputs Top of the canopy (TOC) reflectances in Nadir normalized TOC reflectances in the red, NIR TOC reflectances in the red, NIR and SWIR the red and near infrared (NIR) spectral and short-wave infrared (SWIR) spectral bands, and spectral bands, and cosine of the 3 angles of sun bands cosine of the sun zenith angle at 10:00 local time and view directions Temporal composition 10 d compositing period. The maximum 30 d compositing period with Gaussian weighting Adaptive compositing within 15 and 60d semi- NDVI value in the composition window (minimum of two valid observations) periods defined by the availability of 6 valid is retained. Nominal dates: 3rd, 13th and 21-24th day of each observations at each side of the date being Starting date of composition:1st, 11th month processed and 21th day of the month Nominal dates: 10th, 20th and last day of the month Temporal smoothing and Not applied Not applied Multi-step filtering, temporal smoothing and gap- gap filling filling determined based on the comparison with available ground mea- from 50 to 150 d and selected a time lag of 100 d based on the surements. We tested the 20th, 30th, 40th and 50th percentiles of comparison with ground measurements. SoS and EoS are then de- the annual amplitude for SoS, and the 30th, 40th, 50th and 60th fined as the DOY when the moving average curves cross the original percentiles for EoS. This method is also the basis of SPIRITS phe- curve of the vegetation variable. nological approach. • First derivative: SoS is defined as the DoY of the maximum increase • Autoregressive moving average: A moving average is first computed at (maximum first derivative) in the curve (Tateishi and Ebata, 2004). a randomly chosen time lag (Ivits et al., 2009). We tested time lags EoS is defined as the DoY of the maximum decrease in the curve. Fig. 2. LAI (V1 and V2) and NDVI (V2.1) time series for the PEP725 site 5449 (50°42′20.49″N, 13°46′59.55″E) representative of birch forest in Europe. 3 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Fig. 3. Schematic representation of SoS (on the left of the peak) and EoS (on the right of the peak) retrieved with the four methods for the PEP725 site 4959 (50°42′20.49″N, 13°46′59.55″E) for 2011. The black circles correspond to the original LAI data at a 10-d frequency, and the green line corresponds to the data interpolated at daily steps, which is used for phenological estimation. Table 2 Statistics of comparisons between the derived phenology metrics from re- constructed LAI V1 time series using the different smoothing methods and PEP725 ground measurements (Betula pendula) for the SoS and EoS. The SPIRITS percentile method was used to extract the phenology using 30% of LAI amplitude for the SoS and 40% for the EoS. * Significant correlations with p < 0.05 (**p < 0.001). Metric Method RMSE (d) Bias (d) Abs. Bias (d) R Slope SoS WHITTAKER 23.80 −8.16 16.76 0.41 0.76 "visible BISE 27.22 −7.52 19.84 0.43* 0.89 leaves" MEAN 23.98 −5.84 18.43 0.39 0.74 (n = 359) SWETS 16.44 −1.38 17.39 0.49** 0.67 EoS WHITTAKER 44.56 13.51 31.65 −0.19 −0.12 Fig. 4. Boxplots of the bias errors for (a) SoS and (b) EoS estimated from the "colored BISE 49.45 19.50 37.25 −0.03 −0.10 LAI, FCOVER, FAPAR V1 and V2 and NDVI V2.1 time series minus the PEP725 leaves" MEAN 44.10 13.18 30.54 −0.08 −0.01 ground measurements. An elongated boxplot indicates a greater dispersion of (n = 359) SWETS 34.28 12.01 28.96 −0.15 −0.19 the average bias. • Logistic function: SoS is defined as the DoY of the first local maximum 2.5. Methodological approach rate of change in the curvature of a logistic function fitted to the time series (Zhang et al., 2003). EoS is defined as the DoY of the first The several satellite-derived vegetation variables, smoothing local minimum rate of change in the curvature. methods and phenological extraction approaches lead to a large Table 3 Statistics of comparisons between the derived phenologies from LAI, FCOVER, FAPAR V1 and V2 and NDVI V2.1 time series and PEP725 ground measurements (Betula pendula) for the start of the season (SoS) and end of the season (EoS). The SPIRITS percentile method was used to extract the phenology with specific thresholds per vegetation variable and phenological metric. *mark indicates significant correlations with p < 0.05 (** indicates p < 0.001). Metric Index Version Definition RMSE (d) BIAS (d) Abs. BIAS (d) R Slope SoS LAI V1 30% 16.44 −1.38 17.39 0.49** 0.67 (n = 359) V2 30% 12.49 1.65 10.22 0.62** 0.78 FCOVER V1 40% 17.25 −4.44 23.54 0.52** 0.74 V2 40% 13.95 −4.46 11.21 0.54** 0.95 FAPAR V1 50% 31.40 −18.35 31.23 0.48** 0.85 V2 50% 23.01 −14.26 13.94 0.57** 0.88 NDVI V2.1 50% 20.70 −13.18 15.38 0.25 0.82 EoS LAI V1 40% 34.28 12.01 28.96 −0.15 −0.19 (n = 359) V2 40% 32.72 −6.18 21.12 0.05 0.11 FCOVER V1 40% 30.69 −14.33 42.25 −0.13 −0.20 V2 40% 25.35 6.42 20.34 0.26 0.38 FAPAR V1 50% 44.06 18.75 45.87 −0.23 −0.65 V2 50% 39.99 18.61 37.97 −0.04 −0.51 NDVI V2.1 50% 48.35 14.75 34.58 0.02 −0.25 4 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Fig. 5. Boxplots of the bias error for the SoS estimated from LAI V2 minus the ground measurements at the USA-NPN (a) and PEP725 (b) sites. An elongated boxplot indicates a greater dispersion of the average bias in each method. Table 4 Statistics of comparisons between LAI V2 derived phenology for the start of the season (SoS) and end of the season (EoS), and the ground measurements for the various methods (percentiles, logistic function, derivative and moving average). *mark indicates significant correlations with p < 0.05 (** indicates p < 0.001). Validation source Metric Method RMSE (d) Bias (d) Abs. Bias (d) R Slope USA-NPN SoS 30th percentile 9.19 1.59 6.77 0.81** 0.73 (n = 462) "Leaves" SPIRITS 30% 13.75 0.58 7.05 0.60** 0.58 Logistic 14.14 −5.25 8.66 0.63** 0.60 Derivative 16.82 9.63 12.01 0.68** 0.82 Moving average 30.06 −24.33 23.95 0.59* 0.82 USA-NPN SoS 30th percentile 8.87 −2.20 6.11 0.83** 0.67 (n = 158) "Increaising SPIRITS 30% 10.40 −3.46 6.18 0.78** 0.67 leaf size" Logistic 12.81 −7.24 6.86 0.74** 0.63 Derivative 10.04 4.22 7.73 0.83** 0.84 Moving average 34.49 −28.07 27.86 0.49* 0.70 PEP725 SoS 30th percentile 11.50 1.69 9.84 0.60** 0.90 (n = 359) "visible SPIRITS 30% 12.49 1.65 10.22 0.62** 0.78 leaves" Logistic 17.96 −9.22 14.21 0.53* 0.71 Derivative 28.31 19.19 20.42 0.49* 0.88 Moving average 56.24 −51.35 54.27 0.50* 1.14 USA-NPN EoS 40th percentile 25.61 6.39 17.60 0.14 0.10 (n = 241) "Colored SPIRITS 40% 30.79 −3.69 20.17 0.06 0.07 leaves" Logistic 30.70 21.80 23.66 0.10 0.15 Derivative 40.40 −10.11 26.70 0.02 0.01 Moving average 58.75 50.35 44.26 0.13 0.01 PEP725 EoS 40th percentile 27.69 −5.15 18.89 0.11 0.15 (n = 359) "Colored SPIRITS 40% 32.72 −6.18 21.12 0.05 0.11 leaves" Logistic 30.15 11.91 24.29 0.03 0.28 Derivative 64.93 −43.96 48.49 0.00 −0.11 Moving average 54.95 37.81 52.55 0.01 0.21 number of combinations. We sequentially investigated the impact of validation we used ground measurements of Betula Pendula, which smoothing, variable and extraction method based on an initial set of showed a greater latitudinal distribution. The analysis 3 was performed modalities. The initial modalities were defined a posteriori based on the at the global scale. For the sensitivity analysis 1 and 2, we used the analysis of all the combinations: Software for the Processing and Interpretation of Remotely Sensed Image Time Series (SPIRITS) (Eerens et al., 2014; Eerens and Haesen, 1 Sensitivity analysis of the smoothing method: We used LAI V1 as 2015). For the sensitivity analysis 3, we used Google Earth Engine input dataset and the percentile phenology method. (GEE) (https://earthengine.google.org) which allowed implementing 2 Sensitivity analysis of the vegetation variable: We used the SWETS dedicated algorithms while only the threshold method is available in smoothing method and the percentile phenology method. SPIRITS. The input 10 d time series were linearly interpolated at daily 3 Sensitivity analysis of the method to extract phenological metrics: steps before phenological retrieval. For SPIRITS, the precision of phe- We used the LAI V2. Note that in this case the application of a nological estimates is limited by the frequency of the input time series smoothing method is not required because LAI V2 is already (10 d in our case) (non interpolation). For pixels with multiple growing smoothed and gap-filled. seasons, we computed the phenological metrics for the growing season having the highest LAI amplitude. The analysis 1 and 2 were carried out in Europe and for the The agreement between metric estimates from satellite imagery and 5 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Fig. 6. Scatterplots between the SoS predicted from LAI V2 by the percentile method (a), SPIRITS (b), logistic function (c), derivative (d) and moving average (e) compared with the ground phenology (USA-NPN “leaves”). Values are given in DoY. Statistics of the comparison are indicated in Table 4. ground-based measurements was quantified using the slope of the The best performances for SoS were obtained using the LAI and linear regression, the Pearson correlation coefficient (R), bias, i.e. the FCOVER V2 time series (Fig. 4a), with RMSEs of ˜12 and 14 d, re- average difference between the satellite-derived phenology and the spectively (Table 3). In contrast, RMSEs were ˜21 and 23 d for NDVI observed date (a positive bias indicated that SoS and EoS occurred later V2.1 and FAPAR V2, respectively (Table 3), and contained many out- than the observed leaf out and autumnal coloring, respectively), the liers (Fig. 4a). SoS was slightly underestimated for all cases, except absolute bias and root mean square error (RMSE) calculated (e.g. for when using LAI (Fig. 4a). EoS had a higher RMSE (25–48 d) and a lower SoS) as: R (< 0.3) (Table 3). The estimates of EoS from the LAI and FCOVER V2 n time series also agreed best with ground data although no significant RMSESoS = 1 (SoSground SoSest.) 2 correlations were found (Table 3, Fig. 4b). n j=1 (1) 3.3. Sensitivity analysis of the method to extract phenological metrics where n is the number of samples The 30th percentile and SPIRITS applied to LAI V2 provided the best 3. Results performances among the different analyzed methods when compared both with USA-NPN and Europe-PEP725 measurements of the SoS 3.1. Sensitivity analysis of the smoothing method (Fig. 5, Table 4). The results for the 30th percentile and SPIRITS were similar because both methods use the same definition of SoS based on The SWETS method performed the best for the reconstruction of the 30% threshold of annual amplitude but the 30th percentile method seasonal trajectories and the estimation of phenological metrics:16 d in slightly improved SPIRITS in terms of precision (c.f. Fig. 6a–b, terms of RMSE and 1 d in terms of bias for the timing of SoS, and 34 d Fig. 7a–b) and accuracy (RMSE of 9 vs 14 d for USA-NPN “leaves”, 9 vs (RMSE) and 12 d (bias) for the timing of EoS (Table 2). 10 d for USA-NPN “increasing leaf size” and 11 vs 12 d for PEP725 “first visible leaves” (Table 4)) because of the daily interpolation applied to 3.2. Sensitivity analysis of the vegetation variable the input 10 d data. The logistic function also performed well (RMSE from 13 to 18 d) but provided slightly advanced SoS as compared to The best agreement with ground measurements for the timing of ground measurements (bias from -5 to -9 d). The derivative method SoS was found for the 30% threshold of LAI amplitude, 40% of provided showed higher RMSE (up to 28 d) and positive bias (up to 19 FCOVER, 50% of FAPAR and 50% of NDVI. The best metric definitions d). The moving-average method performance the worst (RMSE from 30 for the EOS were based on 40% of LAI and FCOVER amplitudes and to 59 d) and systematically advanced the timing of SoS as compared to 50% for FAPAR and NDVI. Phenological metrics derived from V2 time ground measurements (bias from -24 to -51d). series improved over V1 for all variables (compare V1 and V2 statistics The different methods provided poorer performances for the EoS in Table 3). (Table 4, Fig. 8). The best agreement with ground measurements was 6 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Fig. 7. Scatterplots between SoS predicted from LAI V2 by the percentile (a), SPIRITS (b), logistic function (c), derivative (d) and moving average (e) methods and ground phenology (PEP725 “first visible leaves”). Values are given in DoY. Statistics of the comparison are indicated in Table 4. Fig. 8. Boxplots of the bias error for the EoS estimated from LAI V2 minus the ground measurements at the USA-NPN (a) and PEP725 (b) sites. An elongated boxplot indicates a greater dispersion of the average bias in each method. found for the 40th percentile method (RMSE of 25 d for USA-NPN and average retrievals showed a positive delay as compared to ground data 28 d for PEP725 “colored leaves” measurements). The logistic function (bias of 38 d for PEP725 and 50 d for USA-NPN). provided similar performances in terms of RMSE but slightly over- estimated ground measurements (bias from 12 to 22 d). The derivative method showed higher scattering (RMSE from 40 to 65 d) and lower 3.4. Spatial patterns of land surface phenology correlation than other methods (Table 4, Figs. 9 and 10). The moving Fig. 11 shows the average timing of SoS, EoS and LoS phenophases 7 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Fig. 9. Scatterplots between the EoS predicted from LAI V2 by the percentile (a), SPIRITS (b), logistic function (c), derivative (d) and moving average (e) methods and ground phenology (USA-NPN “colored leaves”) for Acer rubrum (diamond), Betula alleghaniensis (circle), Fagus grandifolia (triangle) and Quercus rubra (square). Values are given in DoY. Statistics of the comparison are indicated in Table 4. Fig. 10. Scatterplots between the EoS predicted from LAI V2 by the percentile (a), SPIRITS (b), logistic function (c), derivative (d) and moving average (e) methods and ground phenology (PEP725 “colored leaves”). Triangles for heterogeneous sites (forest cover < 50%) and circles for homogeneous sites (forest cover > 50%). Values are given in DoY. Statistics of the comparison are indicated in Table 4. 8 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Fig. 11. Global Map of average SoS (a), EoS (b) and LoS (c) derived from the LAI V2 time series (1999–2017) and the threshold-based method. The continental areas in white are deserts and evergreen forests with limited seasonality where phenology was not computed. at the global scale by using the percentile method (30th and 40th broad variation in the range of values and their spatial pattern at the percentiles for SoS and EoS, respectively) and time series of V2 LAI middle and high latitudes of the Northern Hemisphere strongly depend (1999–2017). The LoS is estimated as the length of time between the on the thermal and photoperiod latitudinal gradient (Verger et al., SoS and the EoS. The derived maps show consistent spatial patterns of 2016). the seasonality of vegetation at the global scale which is driven by the The LoS (Fig. 11c) in some ecoclimatic and biogeographic regions of distribution of latitudinal climatic patterns, type of vegetation and to- transition such as Sahel shows a broad range of variation from 10 to pographic elements, among other factors (Verger et al., 2015; Zhang 200 d following the positive north-south gradient of rainfall (Verger et al., 2004). The timing of SoS (Fig. 11a) and EoS (Fig. 11b) reflected a et al., 2016). On the contrary, in northern latitudes > 50° the LoS 9 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 showed a limited range of variation from 10 to 60 d with shorter days at for detecting EoS in agreement with Verger et al. (2016). The accuracy higher latitudes following the latitudinal gradient of temperature and for SoS using ground observations produced overall RMSEs of 9 and 11 radiation (Verger et al., 2016). d for the date of increasing leaf size and leaf unfolding, respectively, for the forests in Europe and the USA. The biases were < 2 d (˜10 d of 4. Discussion and conclusion standard deviation). We found poorer performances for EoS, than SoS, with higher RMSE of 25 d (28 d) and a bias of 8 d (-6 d) in USA Land surface phenology (LSP) provides a synoptic view of vegeta- (Europe). The lower performances for the EoS as compared to the SoS is tion dynamics and it can substantially improve our macroecological associated to higher uncertainties of both satellite (atmospheric effects, knowledge and the representation of phenology in earth-system snow and poor illumination conditions) (Delbart et al., 2005) and models. Unlike previous studies limited to NDVI, we used three addi- ground (the timing of leaf colouring is more subjective and difficult to tional biophysical variables: LAI, FAPAR and FCOVER generated within identify than spring phenophases like leaf unfolding) (Estrella and CGLS from SPOT-VEGETATION (1999–2013) and PROBA-V Menzel, 2006) phenology for autumn. Richardson et al. (2009) also (2014–2017) satellite imagery. We found that the phenology derived reported higher variabilities across the canopy of the timing and rate of from LAI (or FCOVER) was more closely related to actual ground ob- foliar development in autumn than spring. servation than the NDVI-derived phenology (Table 3). LAI-phenology is The validation of land surface phenology with ground observations based on leaf development rather than on proxies provided by vege- presented some difficulties, such as the spatial distribution and the tation indices which are not driven solely by the amount of leaves but spatial representability of the data. The ground measurements represent also by the canopy structure and the leaf biochemical properties the phenology for a limited number of individual plants that are not (Richardson et al., 2009). LAI is more sensitive than vegetation indices necessarily the most representative species of the 1-km satellite pixels. such as NDVI to larger amounts of vegetation (Myneni and Williams, Conversely, satellite phenology at 1-km resolution represented an in- 1994). In addition, NDVI V2.1 is affected by variations in solar zenith tegrated response across landscapes with diverse species and phenolo- and viewing angles and surface reflectance bidirectional effects (Tote gical behaviors. The phenology of each species and their characteristics et al., 2017). (sizes, ages, homogeneity), though, influenced the satellite signal, de- In addition to the vegetation variable, the derived phenology was pending on its abundance within the pixel sampling area and on the found to be highly sensitive to the retrieval algorithm and processing timing of their phenophases (Delbart et al., 2015). Statistics of the chain. We found that the retrieved phenology performed the best using comparison between LAI V2 derived EoS using the percentile method LAI V2 (Table 3) due to the improved continuity (no missing data in V2) and ground measurements improved significantly when the analysis is and smoothness as compared to V1. In V1 products and in NDVI time- restricted to Acer rubrum in USA-NPN (RMSE of 10 d, bias of -6 d and series, the noise primarily due to cloud contamination and atmospheric significant correlation (p < 0.001) of 0.8 (c.f. Table 4)) and to effects is an important shortcoming in the study of land surface phe- homogeneous sites (forest cover > 50% based on GEE high resolution nology. To overcome this limitation, smoothing methods were applied imagery) in PEP725 (RMSE of 22 d, bias of 6 d and significant corre- in SPIRITS approach as a pre-processing step to phenological retrieval. lation (p < 0.05) of 0.4 (c.f. Table 4)). The choice of the smoothing method introduced differences of up to The differences in the definition of the ground phenophases and 50% in the performance of the phenology derived from LAI V1 satellite metrics hamper the comparison. The logistic function and, (Table 2). These conclusions agree with previous literature studies specially, the moving-average approaches systematically advanced the highlighting the importance of the temporal reconstruction methods SoS as compared to ground measurements since these methods de- (Kandasamy et al., 2013; Verger et al., 2013). However, the phenology termine the SoS as the timing when the vegetation variable starts to derived from original LAI V2 (smoothing and gap filling was already increase (Table 4). On the contrary, the derivative approach based on included in the retrieval algorithm) outperforms the phenology derived the most rapid increase of the signal introduces positive delays in the from the seasonal trajectories derived after smoothing LAI V1 (Table 3). SoS. The opposite trend is observed for the EoS (Table 4): the logistic We tested four state of the art methods to extract phenological function and moving-average show a delay in the detection of the EoS metrics: thresholds, logistic function, derivative and moving average. while the derivatives advances the timing of EoS. On the other hand, Each method has its own strengths and limitations (de Beurs and ground measurements are subjective and some ambiguity exists in the Henebry, 2010). The threshold approach based on a percentage of the definition of phenophases. In this sense, we found positive bias of 2 d annual amplitude is simple and robust but it is sensitive to the for the SoS retrieved with the 30th percentile of LAI amplitude when minimum and maximum values that may be affected by noise in the comparing with USA-NPN “leaves” phenophase but negative bias of -2 d signal. The logistic function approach has been widely used (e.g. Zhang when comparing to “increasing leaf size”. Further studies will focus on et al., 2003) but it is limited to the performance of the model fitting the comparison of the retrieved land surface phenology with continuous (Beck et al., 2006) and it may fail when the curvature function is too ground observations from PhenoCam (Zhang et al., 2018). This should flat to determine the phenophases (de Beurs and Henebry, 2010). The ultimately lead to propose a standardization in the definition of phe- derivative approach based on the maximum increase and decrease of nological metrics. the vegetation variable is very sensitive to the noise in the signal and This research is expected to contribute for the development of a the temporal smoothing and composition approach and it cannot re- dedicated algorithm for the operational retrieval of land surface phe- present short growing seasons well, especially when the increase and nology within CGLS. Validation using ground observations was limited decrease in the annual time series occur rapidly and abruptly (Beck to deciduous broadleaf forests. Further studies should extend the ana- et al., 2006). The moving average approach is based on the assumption lysis to other vegetation types. The methods may need to be adapted to that vegetation growth follows a well-defined temporal profile and it handle multiple and irregularly occurring vegetation growing cycles. may fail in cases of disturbances and abrupt changes. Further, the se- Finally, forecasting approaches need to be developed for near-real time lection of the time lag is arbitrary. land surface phenology retrieval. We found that the choice of the extraction method introduced dif- ferences > 150% and > 85% in the performance of the SoS and EoS, Acknowledgements respectively, derived from LAI V2 when compared to ground observa- tions (PEP725 and USA-NPN) (Table 4). The percentile method agreed This research was supported by an FPU grant (Formación del the best with ground measurements. The validation over ground ob- Profesorado Universitario) from the Spanish Ministry of Education and servations indicated that the 30% threshold of the LAI amplitude was Professional Training to the first author (FPU15-04798), the Copernicus optimal for detecting SoS but that a 40% threshold was more suitable Global Land Service (CGLOPS-1, 199494-JRC), the Spanish 10 K. Bórnez, et al. Int J Appl Earth Obs Geoinformation 84 (2020) 101974 Government grant CGL2016-79835-P, the Catalan Government grant Schnock, U., Stauß, R., van den Boom, T., Weber, E., Zwerger, P., 2009. The BBCH SGR 2017-1005, and the European Research Council Synergy grant system to coding the phenological growth stages of plants history and publications. Journal für Kulturpflanzen 61, 41–52. ERC-2013-SyG-610028 IMBALANCE-P. Myneni, R.B., Williams, D.L., 1994. On the relationship between FAPAR and NDVI. Remote Sens. Environ. 49, 200–211. References Peñuelas, J., Rutishauser, T., Filella, I., 2009. Phenology feedbacks on climate change. Science 324 (5929), 887–888. Reed, B.C., Brown, J.F., Vander, D., 1994. Measuring phenological variability from sa- Atkinson, P.M., Jeganathan, C., Dash, J., Atzberger, C., 2012. Inter-comparison of four tellite imagery. J. Veg. Sci. 5, 703–714. models for smoothing satellite sensor time-series data to estimate vegetation phe- Richardson, A., Braswell, B., Hollinger, D., Jenkins, J., Ollinger, S., 2009. Near-surface nology. Remote Sens. Environ. 123, 400–417. remote sensing of spatial and temporal variation in canopy phenology. Ecol. Appl. 19, Atzberger, C., Eilers, P.H., 2011. Evaluating the effectiveness of smoothing algorithms in 1417–1428. the absence of ground reference measurements. Int. J. Remote Sens. 32, 3689–3709. Richardson, A., Keenan, T., Migliavacca, M., Ryu, Y., Sonnentag, O., Toomey, M., 2013. Baret, F., Weiss, M., Lacaze, R., Camacho, F., Makhmara, H., Pacholcyzk, P., Smets, B., Climate change, phenology, and phenological control of vegetation feedbacks to the 2013. GEOV1: LAI, FAPAR essential climate variables and FCOVER global time series climate system. Agric. For. Meteorol. 169, 156–173. capitalizing over existing products. Part 1: principles of development and production. Schwartz, M., Betancourt, J., Weltzin, J., 2012. From caprio’s lilacs to the USA national Remote Sens. Environ. 137, 299–309. phenology network front. Ecol. Environ. 10, 324–327. Beaumont, L., Hartenthaler, T., Keatley, M., Chambers, L., 2015. Shifting time: recent Swets, D.L., Reed, B.C., Rowland, J.R., Marko, S.E., 1999. A weighted least-squares ap- changes to the phenology of Australian species. Clim. Chang. Res. Lett. 63, 203–214. proach to temporal smoothing of NDVI. ASPRS Annual Conference, From Image to Beck, P.S.A., Atzberger, C., Hogda, K.A., Johansen, B., Skidmore, A.K., 2006. Improved Information, Portland, Oregon, May 17-21, Proceedings: Bethesda, Maryland. monitoring of vegetation dynamics at very high latitudes: a new method using American Society for Photogrammetry and Remote Sens. MODIS NDVI. Remote Sens. Environ. 100, 321–334. Tateishi, R., Ebata, M., 2004. Analysis of phenological change patterns using 1982-2000 Brus, D., Hengeveld, G., Walvoort, D., Goedhart, P., Heidema, A., Nabuurs, G., Gunia, K., Advanced very High-Resolution Radiometer (AVHRR) data. Remote Sens. 25, 2012. Statistical mapping of tree species over Europe. Eur. J. For. Res. Vol. 131, 2287–2300. 145–157. Templ, B., Koch, E., Bolmgren, K., Ungersböck, M., Paul, A., Scheifinger, H., 2018. Pan Cai, Z., Jönsson, P., Jin, H., Eklundh, L., 2017. Performance of smoothing methods for European Phenological database (PEP725): a single point of access for European data. reconstructing NDVI time-series and estimating vegetation phenology from MODIS Int. J. Biometeoroly 62, 1109–1113. data. Remote Sens. 9, 1271. Toté, C., Swinnen, E., Sterckx, S., Clarijs, D., Quang, C., Maes, R., 2017. Evaluation of the Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. A simple SPOT/VEGETATION collection 3 reprocessed dataset: surface reflectances and NDVI. method for reconstructing a high-quality NDVI time-series data set based on the Remote Sens. Environ. 201, 219–233. Savitzky–Golay filter. Remote Sens. Environ. 91, 332–344. Verger, A., Baret, F., Weiss, M., Kandasamy, S., Vermote, E., 2013. The CACAO method de Beurs, K., Henebry, G., 2005. Land surface phenology and temperature variation in the for smoothing, gap filling and characterizing seasonal anomalies in satellite time International Geosphere–Biosphere Program high-latitude transects. Glob. Change series. IEEE Trans. Geosci. Remote. Sens. 51, 1963–1972. Biol. 11, 779–790. Verger, A., Baret, F., Weiss, M., 2014. Near real time vegetation monitoring at global de Beurs, K.M., Henebry, G.M., 2010. Spatio-temporal statistical methods for modelling scale. IEEE J. Stars 7, 3473–3481. land surface phenology. In: Hudson, Keatley (Eds.), Phenological Research: Methods Verger, A., Bater, M., Weiss, M., Filella, Peñuelas, J., 2015. GEOCLIM A global clima- for Environmental and Climate Change Analysis, pp. 177–208. tology of LAI, FAPAR, and FCOVER from VEGETATION observations for 1999-2010. Delbart, N., Kergoat, L., Le Toan, T., Lhermitte, J., Picard, G., 2005. Determination of phe- Remote Sens. Environ. 166, 126–137. nological dates in boreal regions using normalized difference water index. Remote Verger, A., Filella, I., Baret, F., Peñuelas, J., 2016. Vegetation baseline phenology from Sens. Environ. 97, 26–38. kilometric global LAI satellite products. Remote Sens. Environ. 178, 1–14. Delbart, N., Beaubien, E., Kergoat, L., Le Toan, T., 2015. Comparing land surface phe- Verger, A., Baret, F., Weiss, M., 2017. Algorithm Theoretical Basis Document: LAI, nology with leafing and flowering observations from the PlantWatch citizen network. FAPAR, FCOVER Collection 1km, Version 2, Issue I1.40. Available at: https:// Remote Sens. Environ. 160, 273–280. land.copernicus.eu/global/sites/cgls.vito.be/files/products/ Eerens, H., Haesen, D., 2015. User’s manual of Software for the Processing and GIOGL1_ATBD_FAPAR1kmV2_I1.40.pdf. . Interpretation of Remotely Sensed Image Time Series. Version 1.4.0. Available at. Viovy, N., Arino, O., Belward, A.S., 1992. The best index slope extraction (BISE): a http://spirits.jrc.ec.europa.eu/files/SpiritsManual_140.pdf. method for reducing noise in NDVI time-series. Int. J. Remote Sens. 13, 1585–1590. Eerens, H., Haesen, D., Rembold, F., Urbano, F., Tote, C., Bydekerke, L., 2014. Image time White, M.A., de Beurs, K., Didan, K., Inouye, D., Richardson, A., Jensen, O., Lauenroth, series processing for agriculture monitoring. Environ. Model. Softw. 53, 154–162. W., 2009. Intercomparison, interpretation, and assessment of spring phenology in Eilers, P.H.C., 2003. A perfect smoother. Anal. Chem. 75 (14), 3631–3636. North America estimated from remote sensing for 1982–2006. Glob. Change Biol. 15, Estrella, N., Menzel, A., 2006. Responses of leaf colouring in four deciduous tree species 2335–2359. to climate and weather in Germany. Clim. Res. 32, 253–267. Wu, C., Gonsamo, I., Gough, A., Chen, J., Xu, S., 2014. Modeling growing season phe- Garrity, S., Bohrer, G., Maurer, K., Mueller, K., Vogel, C., Curtis, P., 2011. A comparison nology in North American forests using seasonal mean vegetation indices from of multiple phenology data sources for estimating seasonal transitions in deciduous MODIS. Remote Sens. Environ. 147, 79–88. forest carbon exchange. Agric. For. Meteorol. 151, 1741–1752. Yu, L., Liu, T., Bu, K., Yan, F., Yang, J., Chang, L., Zhang, S., 2017. Monitoring the long- Ivits, E., Cherlet, M., Mehl, W., Sommer, S., 2009. Estimating the ecological status and term vegetation phenology change in northeast China from 1982 to 2015. Sci. Rep. 7, change of riparian zones in Andalusia assessed by multi-temporal AVHHR datasets. 14770. Ecol. Indic. 9, 422–431. Zhang, X., Friedl, M.A., Schaaf, C.B., Strahler, A.H., Hodges, J.C.F., Gao, F., Huete, A., Jönsson, P., Eklundh, L., 2002. Seasonality extraction by function fitting to time-series of 2003. Monitoring vegetation phenology using MODIS. Remote Sens. Environ. 84 471- satellite sensor data. IEEE Trans. Geosci. Remote. Sens. 40, 1824–1832. 475. Kandasamy, S., Baret, F., Verger, A., Neveux, P., Weiss, M., 2013. A comparison of Zhang, X., Friedl, M.A., Schaaf, C.B., 2004. Climate controls on vegetation phenological methods for smoothing and gap filling time series of remote sensing observations patterns in northern mid- and high latitudes inferred from MODIS data. Glob. Change application to MODIS LAI products. Biogeosciences 10, 4055–4071. Biol. 10, 1133–1145. Kandasamy, S., Fernandes, R., 2015. An approach for evaluating the impact of gaps and Zhang, X., Jayavelu, S., Liu, L., Friedl, M.A., Henebry, G.M., Liu, Y., Schaaf, C.B., measurement errors on satellite land surface phenology algorithms: application to Richardson, A.D., Gray, J., 2018. Evaluation of land surface phenology from VIIRS 20year NOAA AVHRR data over Canada. Remote Sens. Environ. 164, 114–129. data using time series of PhenoCam imagery. Agric. For. Meteorol. 256, 137–149. Meier, U., Bleiholder, H., Buhr, L., Feller, C., Hack, H., Heß, M., Lancashire, P.D., 11