remote sensing Letter Improved Estimates of Arctic Land Surface Phenology Using Sentinel-2 Time Series Adrià Descals 1,2, * , Aleixandre Verger 1,2 , Gaofei Yin 1,2,3 and Josep Peñuelas 1,2 1 CREAF, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain; verger@creaf.uab.cat (A.V.); g.yin@creaf.uab.cat (G.Y.); josep.penuelas@uab.cat (J.P.) 2 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, 08193 Barcelona, Catalonia, Spain 3 Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 610031, China * Correspondence: a.descals@creaf.uab.cat  Received: 30 September 2020; Accepted: 11 November 2020; Published: 13 November 2020  Abstract: The high spatial resolution and revisit time of Sentinel-2A/B tandem satellites allow a potentially improved retrieval of land surface phenology (LSP). The biome and regional characteristics, however, greatly constrain the design of the LSP algorithms. In the Arctic, such biome-specific characteristics include prolonged periods of snow cover, persistent cloud cover, and shortness of the growing season. Here, we evaluate the feasibility of Sentinel-2 for deriving high-resolution LSP maps of the Arctic. We extracted the timing of the start and end of season (SoS and EoS, respectively) for the years 2019 and 2020 with a simple implementation of the threshold method in Google Earth Engine (GEE). We found a high level of similarity between Sentinel-2 and PhenoCam metrics; the best results were observed with Sentinel-2 enhanced vegetation index (EVI) (root mean squared error (RMSE) and mean error (ME) of 3.0 d and −0.3 d for the SoS, and 6.5 d and −3.8 d for the EoS, respectively), although other vegetation indices presented similar performances. The phenological maps of Sentinel-2 EVI compared well with the same maps extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) in homogeneous landscapes (RMSE and ME of 9.2 d and 2.9 d for the SoS, and 6.4 and −0.9 d for the EoS, respectively). Unreliable LSP estimates were filtered and a quality flag indicator was activated when the Sentinel-2 time series presented a long period (>40 d) of missing data; discontinuities were lower in spring and early summer (9.2%) than in late summer and autumn (39.4%). The Sentinel-2 high-resolution LSP maps and the GEE phenological extraction method will support vegetation monitoring and contribute to improving the representation of Artic vegetation phenology in land surface models. Keywords: land surface phenology; vegetation monitoring; Sentinel-2; arctic; cloud computing; Google Earth Engine 1. Introduction Studies of the Arctic are becoming increasingly important, particularly in the context of the onset of an expected tipping point in the function of its ecosystems as a result of ongoing climate warming [1] that is lengthening the growing season and increasing vegetation productivity [2,3]. Rising levels of vegetation productivity in the Arctic, characterized as greening, result in reduced albedo that further drives the warming trend [4,5] and has consequences for the regional carbon cycle [6]. For example, permafrost is thawing at an accelerating pace in a process that releases greenhouse gases and, consequently, exacerbates the warming trend in the region [7]. Land surface phenology (LSP), the study of remotely-sensed seasonal patterns in vegetation growth [8], complements the sparse field observations at high latitudes and is essential in the assessment Remote Sens. 2020, 12, 3738; doi:10.3390/rs12223738 www.mdpi.com/journal/remotesensing Remote Sens. 2020, 12, 3738 2 of 13 and monitoring of responses of arctic vegetation to climate warming. Reliable LSP maps of the Arctic could be used to support models that show positive feedbacks between trends in climate warming and a lengthening growing season [3], while studies of Arctic LSP may reveal a link between the advancement of spring onset, drought severity, and an increase in the incidence of artic fires over recent decades [9], and support reports of links between carbon uptake and phenology metrics related to vegetation greenness, such as the amplitude of the vegetation index [10]. Studies of Arctic LSP have tended to use moderate-resolution satellite data, such as the Advanced Very-High-Resolution Radiometer (AVHRR) (1.1 km) and the Moderate Resolution Imaging Spectroradiometer (MODIS) (500 m) [2,3]. At these high latitudes, there are two key challenges for LSP analysis [10]. Firstly, LSP methods tend to erroneously detect the end of a snow period as the start of season (SoS), particularly in methods that extract phenology metrics from greenness indices, such as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI), where distinct changes in reflectance during the snowmelt period are incorrectly detected as the onset of vegetation growth. To account for this problem, vegetation indices that are insensitive to snow have been developed to improve phenology estimation [11,12]. Secondly, a lack of valid satellite observations due to persistent cloud cover hampers LSP estimation, so biome-specific algorithms based on the combination of multi-satellite data and spatio-temporal gap filling methods have been developed [3]. Recent advances in remote sensing technologies present opportunities for the estimation of LSP at greater spatial resolution. Sentinel-2 mission provides decametric images with frequent revisit times (<5 d), allowing the extraction of phenology metrics [13] and study of vegetation dynamics at the canopy scale, while the development of cloud-based platforms, such as Google Earth Engine (GEE) [14], allows processing of large volumes of satellite data for planetary-scale analysis [15,16] and increases the accessibility of high-resolution satellite archive data required for time series analyses. The objectives of this study were (1) evaluating the feasibility of Sentinel-2 for LSP retrieval in the Arctic at a spatial resolution of 10 m, (2) proposing a novel and fast cloud computing implementation in GEE of the widely used threshold phenological extraction method, and (3) assessing the performances of the Arctic Sentinel-2 LSP maps for the SoS and end of season (EoS) for the years 2019 and 2020 based on the comparison with MODIS LSP and PhenoCam ground data. 2. Materials and Methods 2.1. Study Area We generated SoS and EoS maps over regions classified as tundra in the RESOLVE Ecoregions dataset 2017 [17]. High latitudes of the Arctic, where winters are cold and summers are short, are mostly uninhabited by human populations; here, tundra soils contain a layer of permafrost that prevents the growth of trees, but supports the growth of grass and shrub vegetation during the short summer period between June and August, as revealed by vegetation index time series (e.g., NDVI and EVI). 2.2. Data We generated LSP metrics for 2019 and 2020 using Sentinel-2 level-2A data, which provide daily top-of-canopy reflectance at 10, 20, and 60 m of spatial resolution; from these data, we used the 10-m resolution bands 2, 3, 4, and 8, and the 20-m band 12. Sentinel-2A and -2B multispectral satellites were launched in 2015 and 2017, respectively, and have a revisit time of 5 d at the equator that decreases with increasing latitude. The maximum revisit time between 1 May and 30 September 2019 for latitudes between 70 and 75◦ was 1.7 d (average: 0.9 d) (Figure 1a). The LSP metrics estimated with Sentinel-2 were compared with the same metrics estimated from the 500 m MOD09GAv6 product [18]. Maximum discontinuity in the time series from spring and early summer (1 May to 15 July) and later summer and autumn (15 July to 30 September) was plotted after non-valid observations had been filtered using values of the quality band scene classification layer (SCL) provided in the Sentinel-2 Remote Sens. 2020, 12, 3738 3 of 13 Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 13 level-2A, where 2 level-2A, 1 =1 saturated where = saturated defective;2 2==dark orordefective; darkarea pixels;33== cloud areapixels; shadows; 66 ==water; cloud shadows; water;and, and, 7−10 = clouds and cirrus (Figure 1b,c). 7−10 = clouds and cirrus (Figure 1b,c). Figure1.1.Revisit Figure Revisittime timeof ofSentinel-2 Sentinel-2 (a) (a) and and maximum maximum discontinuity in the discontinuity in the Sentinel-2 Sentinel-2 Level Level2A 2Atime time series after cloud masking for spring and early summer (1 May to 15 July) (b) and late summer series after cloud masking for spring and early summer (1 May to 15 July) (b) and late summer and and autumn autumn(15 (15July Julytoto30 30September) September)2019 2019(c). (c). 2.3. Phenology Extraction 2.3. Phenology Extraction We adopted a widely used threshold method [13,19,20] that assigns the SoS and EoS as the first We adopted a widely used threshold method [13,19,20] that assigns the SoS and EoS as the first and last days of the season, respectively, on which a threshold u is exceeded; u may be a constant or and last days of the season, respectively, on which a threshold u is exceeded; u may be a constant or defined dynamically for each pixel [21]. In this study, we estimated u as a dynamic value that depends defined dynamically for each pixel [21]. In this study, we estimated u as a dynamic value that depends on the annual amplitude of the time series (Equation (1)): on the annual amplitude of the time series (Equation (1)): u =u (V=min (Vmin −V −max Vmax) )××pp++ V Vmin min,, (1) (1) where Vmin and Vmax are the minimum and maximum annual values in the time series, respectively, where Vmin and Vmax are the minimum and maximum annual values in the time series, respectively, and p is a given proportion (%) of the amplitude. In this study, we used p = 0.5 as the mid-greenup and p is a given proportion (%) of the amplitude. In this study, we used p = 0.5 as the mid-greenup and and mid-greendown of the growing season. The threshold metrics estimated with 50% of the mid-greendown of the growing season. The threshold metrics estimated with 50% of the amplitude amplitude are less affected by biases due to discontinuities in the time series [13]. We applied two are less affected by biases due to discontinuities in the time series [13]. We applied two variants of the variants of the threshold method as follows: threshold method as follows: 1. Threshold method without smoothing. The threshold method was applied directly to the daily 1. Threshold method without smoothing. The threshold method was applied directly to the daily time series. For the SoS, we searched for the earliest date when the daily vegetation index time series. exceeded For the u; then, weSoS, we searched applied for the earliest a linear interpolation date when between the daily this first vegetation observation whenindex the exceeded u; then, we applied a linear interpolation between this vegetation index was >u and the preceding observation from which to estimate SoS as the first observation when the value vegetation index >u (Figure 2c). the >u Forwas EoS, and thethe preceding linear observation interpolation from which was applied between to estimate the latestSoS aswhen date the value the >u (Figure 2c). For the EoS, the linear interpolation was applied between vegetation index was >u, and the subsequent observation, where EoS corresponded to the linearly the latest date when the vegetation interpolated index value >u.was >u, and the subsequent observation, where EoS corresponded to the linearly interpolated 2. Threshold method after value >u. smoothing. The time series data were smoothed prior to the extraction 2. Threshold method after smoothing. of LSP metrics (Figure 2d), as is common The timepractice series data wereestimation in LSP smoothed prior to reduceto thenoise extraction and of LSP metricsof(Figure discontinuities 2d), data time series as is [21]. common practiceforinselection The criteria LSP estimation to reducesteps of the processing noisewere and discontinuities of time series data [21]. The criteria for selection of the processing based on the feasibility of their implementation in GEE, without comprising the recreation of the steps were based on the phenology feasibility curve (see GEE ofcode theirinimplementation Supplementary in GEE, without Materials). comprising Excessive smoothing theofrecreation time series of the phenology may curve (see lead to unrealistic GEE codeof recreations inthe Supplementary growing season. Materials). We firstExcessive applied asmoothing of time moving average series may window, lead with anto unrealistic average radiusrecreations of 10 d, everyof the 20 growing d (Figure season. We first 2d); if a pixel in theapplied a moving 20 d composite average window, window was empty with duean toaverage a lack ofradius of 10 d, everythe valid observations, 20window d (Figure 2d); size was if aincreased pixel in the to 4020d.d Next, we applied composite window a cubic interpolation was empty due to ato convert lack theobservations, of valid 20 d composites to a daily the window time size wasseries. The increased threshold was we to 40 d. Next, estimated appliedfrom the interpolation a cubic amplitude of to theconvert interpolated the 20time series, rather d composites to a than dailywith time daily series.observations, The thresholdand wasthen the SoS estimated fromandtheEoS were estimated amplitude as the first of the interpolated timeand last rather series, days, respectively, that exceeded the dynamic threshold in the interpolated time series. Remote Sens. 2020, 12, 3738 4 of 13 than with daily observations, and then the SoS and EoS were estimated as the first and last days, respectively, that exceeded the dynamic threshold in the interpolated time series. 2.4. Sentinel-2 Vegetation Indices We extracted the phenology metrics from four vegetation indices: the green chromatic coordinate (GCC) (Equation (2)) [22], NDVI (Equation (3)), EVI (Equation (4)) [23], and normalized difference phenology index (NDPI) (Equation (5)) [11]. These spectral indexes reflect the greenness of vegetation and are calculated from different spectral bands in Senteinl-2: Blue (Band 2), Green (Band 3), Red (Band 4), near-infrared (NIR) (Band 8), and shortwave-infrared (SWIR2) (Band 12). The formulas of these vegetation indices are: Green GCC = , (2) Blue + Green + Red NIR − Red NDVI = , (3) NIR + Red (NIR − Red) EVI = 2.5 , and (4) NIR + 6 × Red − 7.5 × Blue + 1 NIR − (alpha × Red + (1 − alpha) × SWIR2) NDPI = (5) NIR + (alpha × Red + (1 − alpha) × SWIR2) where alpha was originally set to 0.74 for MODIS [11], but re-estimated to 0.51 for Sentinel-2 in the current study (see section Estimation of the optimal alpha in NDPI in Supplementary Materials, Supplementary Figure S1). The NDPI was specifically designed to cope with the snow observations in LSP studies. Reclassification of Snow Observations in Green Chromatic Coordinate (GCC), Normalized Difference Vegetation Index (NDVI), and Enhanced Vegetation Index (EVI) NDVI and EVI show a sharp increase after snowmelt in tundra vegetation (Figure 2a,b), because the previously snow-covered vegetation canopy is predominantly evergreen. The break in the time series attributed to the transition snow-to-vegetation and vegetation-to-snow is problematic in the estimation of LSP metrics. The SoS and EoS can be erroneously assigned to the snow transition dates instead to the actual vegetation dynamics. The reclassification of snow and post-thaw values is a common practice in Arctic LSP studies [13,24] to ensure a consistent threshold value (Equation (1)). We reclassified the snow observations to a fixed minimum value, which is specific for each vegetation index: GCCmin = 0.31, NDVImin = 0.39, EVImin = 0.2 (See histograms in Supplementary Figure S2), and NDPImin = 0.24 (Supplementary Figure S1). This parameter was estimated as the mean value of the first Sentinel-2 snow-free observation of the year, extracted from 400 points randomly distributed. For the NDPI, the NDPImin corresponded to the mean NDPI value of snow observations for the optimal alpha (Supplementary Figure S1). The pixels that presented a maximum vegetation index in the time series below the minimum snow value were masked as non-vegetated pixels and the phenology metrics were not computed. Remote Sens. 2020, 12, 3738 5 of 13 Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 13 Figure 2. Illustration Figure Illustrationofofvegetation vegetationdynamics dynamics forfor a randomly a randomly selected pixel selected in the pixel in Arctic. (a) True the Arctic. (a) color True compositions color of Sentinel-2 compositions timetime of Sentinel-2 series data series during data during2019; 2019;(b) (b)Sentinel-2 Sentinel-2 normalized normalized difference vegetation index (NDVI) vegetation (NDVI) time time series seriesduring during2019, 2019,with withreference referencetototrue truecolor compositions color compositions shown shown in (a);(a); in (c)(c) Sentinel-2 Sentinel-2NDVI NDVItime timeseries seriesafter aftercloud cloudmasking maskingandand snow snow reclassification, reclassification, where start ofof season (SoS) season (SoS) and end of season (EoS) correspond to the the earliest and latest linearly interpolated NDVI values, values, respectively, respectively, above above aa threshold threshold defined defined dynamically dynamically as 50% of the the annual annual amplitude; and, panel panel (d)(d) application application of of the the threshold threshold method methodover overaasmoothed smoothedand andinterpolated interpolatedtime timeseries. series. 2.5. Implementation of 2.5. Implementation of Land Land Surface Surface Phenology Phenology Algorithms Algorithms in in Google Google Earth Earth Engine Engine The The two twovariants variantsofofthe thethreshold threshold method methodwere vectorized were vectorizedfor their correct for their implementation correct implementation in GEE. in The vectorization is the process of transforming a code so that all components of GEE. The vectorization is the process of transforming a code so that all components of an array are an array are processed simultaneously [25]. This concept processed simultaneously [25]. Thisis contrary concept to is the commonly contrary to theused practice used commonly in LSP estimation, practice in in LSP which time series are processed separately, pixel by pixel, in a for loop. The use of for estimation, in which time series are processed separately, pixel by pixel, in a for loop. The use of forloops are, however, highly discouraged loops are, however,inhighly GEE indiscouraged preference for inthe GEE in preferencemap recommended forfunctions. For instance, the recommended mapthe moving functions. average for thethe For instance, 20-day composition moving average was implemented for the as a function 20-day composition that was maps a list ofas implemented dates (see code a function in that Supplementary Materials). The function takes the dates as the input arguments maps a list of dates (see code in Supplementary Materials). The function takes the dates as the inputto filter the Sentinel-2 collection argumentswith a window to filter size of 20 the Sentinel-2 days and, collection withthen, makes the a window average size out ofand, of 20 days the selected images. then, makes the The purpose of the map function is that each element of the array, in this case the average out of the selected images. The purpose of the map function is that each element of the array, dates, is processed separately in this caseand,theconsequently, each 20-day dates, is processed composition separately and, isconsequently, generated at the same each time. composition is 20-day generated at the same time. 2.6. Validation with PhenoCam The SoS and 2.6. Validation EoS with metrics extracted from Sentinel-2 were compared with the same metrics estimated PhenoCam from the near-surface reflectances of the PhenoCam network [26]. PhenoCam provides half-hourly The SoS and EoS metrics extracted from Sentinel-2 were compared with the same metrics images captured from digital cameras for 393 sites across North America and Europe, from which estimated from the near-surface reflectances of the PhenoCam network [26]. PhenoCam provides 15 cameras cover the tundra biome. At the moment of the analysis, GEE provided Sentinel-2 data half-hourly images captured from digital cameras for 393 sites across North America and Europe, from which 15 cameras cover the tundra biome. At the moment of the analysis, GEE provided Sentinel-2 data from 2019 onwards, and only three sites presented available data in the Arctic for the Remote Sens. 2020, 12, 3738 6 of 13 from 2019 onwards, and only three sites presented available data in the Arctic for the years 2019 and 2020. The PhenoCam data used in the study was presented as provisional near-real time and subject to changes. The coordinates of these three PhenoCam cameras, used in the study, are shown in Supplementary Table S1. PhenoCam also provides daily time series of GCC derived from different regions of interest observed by the digital camera [26]. Such regions of interest cover the vegetation types in the site. For the selected sites, however, the camera only observed one vegetation type consisting of tundra grasses. The GCC has been proved a good index for LSP estimation [22] since the time series does not present the sharp increase during the snowmelt. However, we also reclassified the snow values in PhenoCam GCC to the minimum value of the snow-free time series in order to be consistent with the forcing applied in the other Sentinel-2 vegetation indices. The snow observations in GCC were identified and selected by their similar values before and after the growing season (see Supplementary Figures S3–S5). We extracted the LSP metrics with a 50% threshold method from the daily GCC time series for the years 2019 and 2020 and compared them with the same LSP metrics estimated from the four Sentinel-2 vegetation indices (GCC, NDVI, EVI, and NDPI). The statistics that we reported were the mean error (ME) as the bias metric and the root mean squared error (RMSE) as the accuracy metric. Although the selected PhenoCam sites cover a homogeneous area, the location of the sites were manually relocated to ensure consistency with the area observed by the digital camera. The adjustment of the site coordinates was based on the orientation of the camera; the coordinates were relocated to 50 m from the original coordinates in the direction of the observation of the camera. 2.7. Comparison with the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Surface Phenology Our estimates of SoS and EoS using Sentinel-2 time series were compared with estimates based on MODIS time series data, using the same metrics and LSP extraction method; we also applied the reclassification of snow values to the MODIS time series. To compare both satellite datasets, we resized the 10 m Sentinel-2 LSP estimates to the 500 m MODIS projection. The Sentinel-2 was aggregated by averaging the pixels that lie within a 500 m MODIS pixel. The resampled SoS and EoS estimated with Sentinel-2 and the original MODIS LSP metrics, both at 500 m, were compared pixel-wise and the ME and RMSE were reported. To ensure reliable estimates in MODIS, we only considered the MODIS LSP estimates obtained from time series that presented a gap lower than 10 days. This condition did not result in an excessive filtering of pixels since the MODIS time series had a maximum gap that generally did not exceed 20 days (Supplementary Figure S6). To test the reliability of our MODIS LSP estimates, we compared the SoS and EoS with the ‘Mid_greenup’ and ‘Mid_greendown’ layers, respectively, of the MODIS Land Cover Dynamics product (MCD12Q2v6) [27]. The MCD12Q2v6 product was also estimated using the threshold method, with a dynamic threshold of 50% in the ‘Mid_greenup’ and ‘Mid_greendown’ layers [19]. MCD12Q2v6 data for 2019 were unavailable at the time of this analysis, so we compared our MODIS LSP estimates for 2018 with the MCD12Q2v6 layers for 2018. 3. Results The Sentinel-2 vegetation index that showed the best results in the comparison with PhenoCam was the EVI (SoS and EoS ME: −0.3 and −3.8 d, and SoS and EoS RMSE: 3.0 and 6.5 d) for the threshold method without time series smoothing (Figure 3). Results with the threshold method with smoothing were less conclusive (Supplementary Figure S7); RMSE and ME results were uneven between SoS and EoS for the same LSP method and none of the vegetation indices excelled in both ME and RMSE. Overall, the four Sentinel-2 vegetation indices compared well with PhenoCam, and the RMSE and ME were generally below 10 d for the two threshold methods (see time series of PhenoCam and Sentinel-2 in Supplementary Figures S3–S5). Remote Sens. 2020, 12, 3738 7 of 13 Remote Sens. Remote Sens. 2020, 2020, 12, 12, xx FOR FOR PEER PEER REVIEW REVIEW 77 of of 13 13 Figure 3. Figure Figure 3. Comparison Comparison Comparisonof of start ofstart ofofseason startof season season (SoS) (SoS) (SoS) and and and end endendof season of season of season (EoS) dates (EoS) (EoS) between dates dates between between PhenoCam PhenoCam PhenoCam and and four and vegetation four four vegetation vegetation indices indices indices estimated estimated estimated with withwith Sentinel-2. The Sentinel-2.The Sentinel-2. comparison Thecomparison comparisonwas was performed wasperformed performed for for three for three three PhenoCam sites PhenoCam sites in in the the tundra tundra biome biome for for the the years years 2019 2019 the tundra biome for the years 2019 and 2020. and and 2020. 2020. The The vegetation vegetation vegetation indices indices were were the the the green green chromatic chromatic coordinate coordinate (GCC) (GCC) in in PhenoCam PhenoCam and and the the GCC, GCC, normalized normalized green chromatic coordinate (GCC) in PhenoCam and the GCC, normalized difference vegetation difference difference vegetation vegetation index index (NDVI), (NDVI), index (NDVI), enhancedenhanced enhanced vegetation vegetation index indexindex vegetation (EVI), (EVI),(EVI), and normalized normalized and normalized and difference difference phenology phenology difference index phenology index (NDPI) index (NDPI) (NDPI)for for Sentinel-2. Sentinel-2. TheThe for Sentinel-2. The phenology phenology phenology metrics metrics were werewere metrics extracted extracted with with with extracted a 50% a 50%athreshold threshold method 50% threshold method without method without time time series without time series smoothing. series smoothing. smoothing. The bias The The bias between between bias between PhenoCam PhenoCam and Sentinel-2 and Sentinel-2 PhenoCam and Sentinel-2 is reportedis reported is reported with with thewith mean the the mean error mean(ME)error and error (ME) the (ME) and the accuracy accuracy with the with root the mean root mean squared squared error error (RMSE). and the accuracy with the root mean squared error (RMSE). (RMSE). LSP LSP maps LSP maps generated maps generated using generated using Sentinel-2 using Sentinel-2 EVI Sentinel-2 EVI (Figure EVI (Figure 4a,b) (Figure 4a,b) showed 4a,b) showed aaa high showed high level high level of level ofsimilarity of similarity at similarity atthe at the the continental scale continental scale continental compared scale compared compared with with the with the same the same phenology same phenology metrics phenology metrics estimated metrics estimated estimated withwith MODIS with MODIS MODIS time time series time series series (Supplementary (Supplementary Figure S8). However, at the local scale, differences were apparent particularly in the (Supplementary Figure Figure S8). S8). However, However, at at the the local local scale, scale, differences differences were were apparent apparent particularly particularly in in the the EoS; EoS; similarly, similarly, the map the mapof length of of length season of (LoS) season (Figure (LoS) 5a), (Figure which 5a), EoS; similarly, the map of length of season (LoS) (Figure 5a), which represented the difference represented which the represented difference the between difference EoS and SoS, between between EoS showed EoS and SoS,spatial and SoS, showed showed patterns spatial and spatial patterns patternslatitudinal gradients and latitudinal and latitudinal that are that gradients gradients expected that are across the are expected expected region, across across the the with shorter region, with LoS in shorter the LoSsurrounding in the areas surrounding of the areas Kara of Sea the and Kara Arctic Sea and region, with shorter LoS in the surrounding areas of the Kara Sea and Arctic Archipelago (Figure 5b),Archipelago Arctic (Figure Archipelago 5b), (Figureand in 5b), elevated and in areas, elevated such areas,as the such Central as the Siberian Central Plateau Siberian (Figure Plateau 5c). (Figure At the 5c). and in elevated areas, such as the Central Siberian Plateau (Figure 5c). At the local level, the SoS and local At the level, local the level, SoSthe and SoS EoS and maps EoS derived EoS maps maps fromfrom derived derived Sentinel-2 from showed Sentinel-2 Sentinel-2 a high showed showed degree aa high of detail high degree degree that moderate of detail of detail that moderate that scalesscales moderate of resolution scales were of resolution of resolution unable were to capture. unable to For capture.instance, For Figure instance, 5d shows Figure 5d the positive shows the relationship positive were unable to capture. For instance, Figure 5d shows the positive relationship between LoS and between relationship LoS and between elevation LoS and in the Scandinavian elevation elevation in the in mountains.mountains. the Scandinavian Scandinavian The altitudinal mountains. Thegradient The altitudinal altitudinalis observed gradientat gradient isthe is canopyat observed observed scale, at which canopyallows the canopy the scale, scale, the analysis which allows of LSP the metrics analysis of by LSP vegetation metrics by type. Figure vegetation 5e type. shows Figure another which allows the analysis of LSP metrics by vegetation type. Figure 5e shows another example of the5e shows example another of the example effect of of the topography effect of on the topography LoS onin the the Ural LoS inmountains. the Ural In this mountains. second In thisexample, second the effect of topography on the LoS in the Ural mountains. In this second example, the LoS is different in LoS example, is different the LoS in is both differentsites in of thesites both both mountain sites of therange of the mountaindepending mountain range on the presence range depending depending on the on theofpresence glaciers. of presence Supplementary of Figure S9 shows glaciers. Supplementary glaciers. Supplementary Figurethe Figure S9 S9 different shows spring the growth different onset spring depending growth onset on the depending land cover on the type, land shows the different spring growth onset depending on the land cover type, tundra shrublands andtundra cover shrublands type, tundra and herbaceous shrublands and cover in the cover herbaceous herbaceous delta coverof inthe in theLena the deltaRiver. delta of the of the Lena Lena River. River. Figure4. Figure Figure 4. Maps 4. Maps of Maps of (a) of (a) start (a) start of start ofseason of season (SoS) season (SoS) and (SoS) and(b) and (b)end (b) endof end ofseason of season(EoS) season (EoS)in (EoS) inthe in theArctic the Arcticestimated Arctic estimatedusing estimated using using Sentinel-2for Sentinel-2 Sentinel-2 for2019. for 2019. The 2019. The phenology The phenology metrics phenology metrics were metrics were extracted were extracted with extracted with the with the 50% the 50% threshold 50% thresholdmethod threshold methodapplied method appliedto applied to to the the non-smoothed non-smoothed Sentinel-2 Sentinel-2 enhanced enhanced vegetation vegetation index index (EVI) (EVI) time time the non-smoothed Sentinel-2 enhanced vegetation index (EVI) time series. series. series. Remote Sens. 2020, 12, 3738 8 of 13 Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 13 Figure 5. Figure 5. Map Map ofof length length of of season season (LoS) (LoS) in in the the Arctic Arcticestimated estimatedusingusingSentinel-2 Sentinel-2enhanced enhancedvegetation vegetation index (EVI) time series for 2019 (a). EVI time series and start of season (SoS) andend index (EVI) time series for 2019 (a). EVI time series and start of season (SoS) and endofofseason season(EoS) (EoS) fortwo for two randomly randomly selected selected pixels pixels (b,c); (b,c); note, note, the the time timeseries seriesinin(b) (b)show showaashorter shorterLoS LoSthan thaninin(c) (c)due due to the to the latitudinal latitudinal gradient. gradient. Arctic Arctic digital digital elevation elevation model modelmapmap(ArcticDEM) (ArcticDEM)and andLoS LoSexemplify exemplifythe the altitudinal gradient of phenology in the Scandinavian mountains and the Ural altitudinal gradient of phenology in the Scandinavian mountains and the Ural mountains (d,e) andmountains (d,e) and comparison of comparison of LoS LoS estimated estimated with with Sentinel-2 Sentinel-2 and and Moderate Moderate Resolution ResolutionImaging ImagingSpectroradiometer Spectroradiometer (MODIS) (d). The phenology metrics were extracted with the 50% thresholdmethod (MODIS) (d). The phenology metrics were extracted with the 50% threshold methodapplied appliedtotothethe non-smoothed Sentinel-2 non-smoothed Sentinel-2 EVI EVI time time series. series. Thecomparison The comparisonbetween betweenthe theLSPLSPmetrics metricsestimated estimatedusing usingMODIS MODISandandSentinel-2 Sentinel-2resized resizedtoto500 500m, m,which in in which MODIS MODIS pixels pixels werewere filtered filtered whenvariance when varianceof ofthe the Sentinel-2 Sentinel-2 phenology phenologyestimates estimateswithin within the500 the 500mmpixel were >5 pixel were >5 d22,, showed showedaahigh highlevel levelof ofsimilarity similarityininhomogeneous homogeneouslandscapes landscapes(Figure (Figure6),6), regardless of regardless of data data smoothing smoothing (Figure 6a,b), 6a,b), and and greater greater similarity similarityforforSoS SoSthan thanEoS EoS(e.g., (e.g.,threshold threshold method with method with smoothing smoothing SoS SoS andand EoSEoS RMSE: RMSE: 8.1 8.1 and and 10.4 10.4 d, d, respectively). respectively). The Thecomparison comparisonofofour our MODIS estimates and the MCD12Q2v6 product (Supplementary Figure S10) MODIS estimates and the MCD12Q2v6 product (Supplementary Figure S10) showed a slight bias showed a slight bias towards EoS towards EoS (threshold (threshold method method with with smoothing smoothing ME:ME: 4.9 4.9 and and6.4 6.4d, d,respectively) respectively)and andsimilarity similaritywith with Sentinel-2 for SoS and EoS (threshold method without smoothing RMSE = 10.2 and Sentinel-2 for SoS and EoS (threshold method without smoothing RMSE = 10.2 and 9.8 d for SoS and 9.8 d for SoS and EoS,respectively). EoS, respectively). The phenology maps of the Arctic presented some unreliable LSP estimates, particularly when there were continuous gaps in the time series, due to the presence of clouds. The proportion of vegetated land pixels in the study area with a gap length that exceeded 40 days was 9.2% for the spring and early summer period and 39.4% for the late summer and autumn period. These percentages refer only to the pixels that presented an EVI value higher than 0.2 and, thus, were considered as vegetated land in the study. The percentage of pixels flagged as non-vegetated areas (EVI values lower than 0.2 during the entire time series) was 15.4%. Figure 7 illustrates the flags associated with the SoS estimates. The regions with low revisit times (Figure 1a) were prone to discontinuities in the time series. Remote Sens. 2020, 12, 3738 9 of 13 Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 13 Figure 6. Comparison of the phenology metrics, start of season (SoS) and end of season (EoS), estimated with Sentinel-2 and MODIS (MOD09GAv6) time series using the threshold method without (a) and with (b) data smoothing. The bias between PhenoCam and Sentinel-2 is reported with the mean error (ME) and the accuracy with the root mean squared error (RMSE). The phenology maps of the Arctic presented some unreliable LSP estimates, particularly when there were continuous gaps in the time series, due to the presence of clouds. The proportion of vegetated land pixels in the study area with a gap length that exceeded 40 days was 9.2% for the spring and early summer period and 39.4% for the late summer and autumn period. These percentages refer only to the pixels that presented an EVI value higher than 0.2 and, thus, were considered Figure asComparison Figure 6.6. vegetated of Comparison land of in phenology thethe the study. phenology Thestart metrics, percentage metrics, start of of pixels of season season (SoS) flagged (SoS) and end and as non-vegetated end of season of season (EoS), (EoS),areas estimated (EVIwith values estimated lower than 0.2 and with Sentinel-2 Sentinel-2 and MODIS during MODISthe(MOD09GAv6) (MOD09GAv6)entire time time series) time series was series using 15.4%. theusing the Figure threshold threshold method 7 illustrates (a) the method without without and flags associated (a) andwith with the (b) SoS data estimates. smoothing. The The regions bias with between low PhenoCam revisit and times (Figure Sentinel-2 is 1a) reported with (b) data smoothing. The bias between PhenoCam and Sentinel-2 is reported with the mean error were with prone the to discontinuities meanand (ME) error in the the(ME) time and the accuracy series. accuracy with the rootwith meanthesquared root mean errorsquared (RMSE).error (RMSE). The phenology maps of the Arctic presented some unreliable LSP estimates, particularly when there were continuous gaps in the time series, due to the presence of clouds. The proportion of vegetated land pixels in the study area with a gap length that exceeded 40 days was 9.2% for the spring and early summer period and 39.4% for the late summer and autumn period. These percentages refer only to the pixels that presented an EVI value higher than 0.2 and, thus, were considered as vegetated land in the study. The percentage of pixels flagged as non-vegetated areas (EVI values lower than 0.2 during the entire time series) was 15.4%. Figure 7 illustrates the flags associated with the SoS estimates. The regions with low revisit times (Figure 1a) were prone to discontinuities in the time series. Figure 7.7. Example Figure Example ofof flags flags associated associated with with gaps gaps in in the the Sentinel-2 Sentinel-2 time time series series in in the the Scandinavian Scandinavian mountains. mountains. (a)(a) Revisit Revisit time; time; (b) (b) maximum maximum enhanced enhanced vegetation vegetation index index (EVI) (EVI) for for 2019; 2019; (c) (c) flags flags raised raised during the estimation of end-of-season (EoS); (d) EoS estimation. Red: non-vegetated pixels during the estimation of end-of-season (EoS); (d) EoS estimation. Red: non-vegetated pixels (EVI<0.2);(EVI < 0.2); blue: pixels with >40-d gap in time series; EoS was extracted with the 50% threshold method blue: pixels with >40-d gap in time series; EoS was extracted with the 50% threshold method applied applied to to the thenon-smoothed non-smoothedSentinel-2 Sentinel-2EVI EVItime timeseries. series. 4. Discussion The high agreement between the phenology metrics estimated with Sentinel-2 and MODIS corroborates the feasibility of LSP estimation owing to the low revisit time of the Sentinel-2 time series [13] in a region prone to cloud coverage. We found that only 9.2% of pixels showed a discontinuity of >40 d in 7. Figure theExample time series during of flags the green-up associated with gapsperiod. In contrast, in the Sentinel-2 timethe cloud series coverage in the was more Scandinavian persistent in late summer and autumn (39.4% of pixels with >40-d gap), and this high level mountains. (a) Revisit time; (b) maximum enhanced vegetation index (EVI) for 2019; (c) flags raised of cloud occurrence may explain the lower similarity between EoS estimated using Sentinel-2 and during the estimation of end-of-season (EoS); (d) EoS estimation. Red: non-vegetated pixels (EVI<0.2); MODIS. Such blue: discontinuities in time pixels with >40-d gap series in timemake series; the EoS use of gap-filling was extracted techniques with the and method 50% threshold robust applied smoothing techniques necessary for the accurate estimation to the non-smoothed Sentinel-2 EVI time series. of LSP metrics [13]. Moreover, when time series present continuous gaps, the combined use of Sentinel-2 and Landsat-8 is recommended as, overall, the combination of sensors provides improved LSP estimates [28]. Remote Sens. 2020, 12, 3738 10 of 13 The comparison with PhenoCam further proves the feasibility of LSP estimation with Sentinel-2 with a variety of vegetation indices. Results showed that the EVI performed the best. It was shown by [28] that EVI was more suited than NDVI in temperate deciduous forests and [13] used EVI2 (two-band EVI) for continental-scale analysis, both using a combination of Landsat-8 and Sentinel-2. Despite this, an extensive analysis including more in situ observations is required to further justify the use of EVI in Arctic phenology studies. Moreover, the correct reclassification of snow values seems more determinant than the variable selection in high-latitudes. The procedure for flagging pixels with discontinuities in the time series may be used to identify good-quality pixels for further analysis in phenology studies; however, this flag only works when clouds are successfully filtered. Non-valid observations that were not filtered, mainly cloud-contaminated pixels and cloud shadows that were not reflected in the quality band, are included in the time series and may underestimate the long gaps in the Sentinel-2 time series. Furthermore, such contaminated values may change the growing season curve in the vegetation index time series and lead to erroneous estimates of the phenology. We found high levels of similarity between the phenology metrics estimated using Sentinel-2 and MODIS and PhenoCam, even though the comparison reflects changes in vegetation greenness that do not necessarily correspond to vegetation phenophases or vegetation productivity [8]. In situ databases, such as the National Phenology Network (NPN) [29], the Pan European Phenology database (PEP725) [30], and the FLUXNET network [31] are essential to ground-truth remotely sensed phenological changes in vegetation. Near-surface LSP, such as the PhenoCam network, supports the validation of satellite LSP estimates, as observed in the current study. Phenophase observations provide information on the timing of relevant stages in the annual life cycle of vegetation, such as the leaf-out, which may differ from the LSP metrics estimated from a satellite [20]. Similarly, the onset and the end of the photosynthetic activity, estimated from time series of carbon fluxes in FLUXNET, do not match with the LSP metrics depending on the forest ecosystem [32]. Phenophase, carbon fluxes, and near-surface LSP measurements are scarce in the Arctic, with the distribution of digital cameras generally restricted to northern areas of Alaska, so data are lacking for adequate validation of high-latitude LSP maps. The use of cloud platforms, such as GEE, has made possible the estimation of LSP using Sentinel-2 time series for extensive regions, such as the Arctic. In this study, the generation of the phenology maps for the Arctic used approximately 133,000 Sentinel-2 images, solely for the year 2019; this type of platform is a useful tool that allows the scientific community to inspect such dense, high-spatial resolution time series. Furthermore, GEE allows sharing the code and data, and researchers can easily reproduce the algorithm and inspect the time series locally (see GEE code in Supplementary Materials). The resulting phenology maps should be taken with caution when data availability is scarce, which occurs particularly during the EoS (39.4% of pixels with >40 d gap). The combination of Sentinel-2 and Landsat-8 is recommended when continuous gaps are present in the time series, and more elaborated phenological retrieval methods including gap-filling and robust smoothing techniques should be used instead. Despite this, sensor harmonization and temporal smoothing may improve the LSP retrievals but at the expenses of adding complexity to the preprocessing, which impedes the implementation and fast computation in GEE. Our SoS and EoS maps may benefit future Arctic phenology studies that aim to analyze spatial variability in vegetation dynamics and require a high degree of detail at the canopy level; test effects of spatial scaling on temporal changes in phenology; and identify homogeneous landscapes in which phenology dynamics are similar. 5. Conclusions Here, we present the first high-spatial resolution maps (10 m) of SoS and EoS that fully cover the Arctic for the most recent years (2019–2020). We prove the feasibility of phenology metric extraction in the region solely with Sentinel-2 and using basic implementations of the threshold method in GEE; the high revisit time at high latitudes allows dense time series to be obtained under cloud-free Remote Sens. 2020, 12, 3738 11 of 13 conditions in 1–3 days and only 9.2% of pixels showed a discontinuity of >40 d in the time series during the green-up period. We propose a set of forcing values for the reclassification of snow observations in Sentinel-2 for three common vegetation indices used in LSP estimation: GCC, NDVI, and EVI; and also provide a re-adjusted parameterization of the NDPI specifically for Sentinel-2. Future work may use the adjusted parameters and forcing values along with the implementation of the threshold method in GEE for large-scale estimation of phenology metrics. The retrieved Sentinel-2 high-resolution LSP maps and the proposed GEE phenological extraction method will support monitoring vegetation changes at high-spatial resolution and are expected to contribute to the representation of Artic vegetation phenology in land surface models. Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/22/3738/s1, Supplementary Materials contain one table, ten figures, the section Estimation of the optimal alpha in NDPI, and the GitHub link to the Google Earth Engine code. Table S1: PhenoCam sites used in the study, Figure S1: Estimation of the optimal alpha for the normalized difference phenology index (NDPI), Figure S2: Histogram of Sentinel-2 vegetation indices after snowmelt, Figure S3: Time series in the PhenoCam site imcrkrtussock, Figure S4: Time series in the PhenoCam site NEON.D18.TOOL.DP1.0003, Figure S5: Time series in the PhenoCam site NEON.D19.HEAL.DP1, Figure S6: Maximum discontinuity in the Moderate Resolution Imaging Spectroradiometer (MODIS) time series after cloud masking, Figure S7: Comparison of PhenoCam and Sentinel-2 phenology metrics estimated with time series smoothing, Figure S8: Difference map of phenology metrics between MODIS and Sentinel-2, Figure S9: Differences in the vegetation dynamics in the delta of the Lena river, Figure S10: Comparison between the phenology metrics estimated with MODIS and the MODIS Land Cover Dynamics. Author Contributions: Conceptualization, A.D., A.V., G.Y., and J.P.; methodology, A.D. and A.V.; software, A.D.; formal analysis, A.D., A.V., G.Y., and J.P.; supervision, A.V. and J.P.; writing—original draft preparation, A.D.; writing—review and editing, A.D., A.V., G.Y., and J.P. All authors have read and agreed to the published version of the manuscript. Funding: We thank the funding of this research by the European Research Council ERC SyG-610028 2013-Synergy-grant, IMBALANCE-P, the Marie Skłodowska-Curie grant of European Union’s Horizon 2020 research and innovation programme (835541), the Spanish project PID2019-110521GB-I00, the Ministry of Science and Innovation of Spain FPI grant BES–2017-080197, and the Catalan grant SGR 2017-1005. Conflicts of Interest: The authors declare no conflict of interest. References 1. Post, E.; Alley, R.B.; Christensen, T.R.; Macias-Fauria, M.; Forbes, B.C.; Gooseff, M.N.; Iler, A.; Kerby, J.T.; Laidre, K.L.; Mann, M.E.; et al. The polar regions in a 2◦ C warmer world. Sci. Adv. 2019, 5, eaaw9883. [CrossRef] 2. Park, T.; Ganguly, S.; Tømmervik, H.; Euskirchen, E.S.; Høgda, K.-A.; Karlsen, S.R.; Brovkin, V.; Nemani, R.R.; Myneni, R.B. Changes in growing season duration and productivity of northern vegetation inferred from long-term remote sensing data. Environ. Res. Lett. 2016, 11, 084001. [CrossRef] 3. Zeng, H.; Jia, G.; Epstein, H. Recent changes in phenology over the northern high latitudes detected from multi-satellite data. Environ. Res. Lett. 2011, 6, 045508. [CrossRef] 4. Chapin, F.S.; Sturm, M.; Serreze, M.C.; McFadden, J.P.; Key, J.; Lloyd, A.H.; McGuire, A.; Rupp, T.S.; Lynch, A.H.; Schimel, J.P.; et al. Role of land-surface changes in Arctic summer warming. Science 2005, 310, 657–660. [CrossRef] 5. Peñuelas, J.; Filella, I. Phenology feedbacks on climate change. Science 2009, 324, 887–888. [CrossRef] 6. Schuur, E.A.; Vogel, J.G.; Crummer, K.G.; Lee, H.; Sickman, J.O.; Osterkamp, T. The effect of permafrost thaw on old carbon release and net carbon exchange from tundra. Nature 2009, 459, 556–559. [CrossRef] 7. Schuur, E.A.; McGuire, A.D.; Schädel, C.; Grosse, G.; Harden, J.; Hayes, D.J.; Hugelius, G.; Koven, C.D.; Kuhry, P.; Lawrence, D.M.; et al. Climate change and the permafrost carbon feedback. Nature 2015, 520, 171–179. [CrossRef] 8. Helman, D. Land surface phenology: What do we really ‘see’from space? Sci. Total Environ. 2018, 618, 665–673. [CrossRef] 9. Witze, A. The Arctic is burning like never before-and that’s bad news for climate change. Nature 2020, 585, 336–337. [CrossRef] Remote Sens. 2020, 12, 3738 12 of 13 10. Myers-Smith, I.H.; Kerby, J.T.; Phoenix, G.K.; Bjerke, J.W.; Epstein, H.E.; Assmann, J.J.; John, C.; Andreu-Hayles, L.; Angers-Blondin, S.; Beck, P.S.; et al. Complexity revealed in the greening of the Arctic. Nat. Clim. Change 2020, 10, 106–117. [CrossRef] 11. Wang, C.; Chen, J.; Wu, J.; Tang, Y.; Shi, P.; Black, T.A.; Zhu, K. A snow-free vegetation index for improved monitoring of vegetation spring green-up date in deciduous ecosystems. Remote Sens. Environ. 2017, 196, 1–12. [CrossRef] 12. Jin, H.; Eklundh, L. A physically based vegetation index for improved monitoring of plant phenology. Remote Sens. Environ. 2014, 152, 512–525. [CrossRef] 13. Bolton, D.K.; Gray, J.M.; Melaas, E.K.; Moon, M.; Eklundh, L.; Friedl, M.A. Continental-scale land surface phenology from harmonized Landsat 8 and Sentinel-2 imagery. Remote Sens. Environ. 2020, 240, 111685. [CrossRef] 14. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 2017, 202, 18–27. [CrossRef] 15. Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.; Goetz, S.J.; Loveland, T.R.; et al. High-resolution global maps of 21st-century forest cover change. Science 2013, 342, 850–853. [CrossRef] 16. Pekel, J.-F.; Cottam, A.; Gorelick, N.; Belward, A.S. High-resolution mapping of global surface water and its long-term changes. Nature 2016, 540, 418–422. [CrossRef] 17. Dinerstein, E.; Olson, D.; Joshi, A.; Vynne, C.; Burgess, N.D.; Wikramanayake, E.; Hahn, N.; Palminteri, S.; Hedao, P.; Noss, R.; et al. An ecoregion-based approach to protecting half the terrestrial realm. BioScience 2017, 67, 534–545. [CrossRef] 18. Vermote, E.; Wolfe, R. MOD09GA MODIS/Terra Surface Reflectance Daily L2G Global 1 km and 500 m SIN Grid V006. Available online: http://doi.org/10.5067/MODIS/MOD09GA.006 (accessed on 11 November 2020). 19. Gray, J.; Sulla-Menashe, D.; Friedl, M.A. User Guide to Collection 6 MODIS Land Cover Dynamics (MCD12Q2) Product. Available online: https://modis-land.gsfc.nasa.gov/pdf/MCD12Q2_Collection6_UserGuide.pdf (accessed on 11 November 2020). 20. Bórnez, K.; Descals, A.; Verger, A.; Peñuelas, J. Land surface phenology from VEGETATION and PROBA-V data. Assessment over deciduous forests. Int. J. Appl. Earth Obs. Geoinf. 2020, 84, 101974. [CrossRef] 21. Zeng, L.; Wardlow, B.D.; Xiang, D.; Hu, S.; Li, D. A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Remote Sens. Environ. 2020, 237, 111511. [CrossRef] 22. Sonnentag, O.; Hufkens, K.; Teshera-Sterne, C.; Young, A.M.; Friedl, M.; Braswell, B.H.; Milliman, T.; O’Keefe, J.; Richardson, A.D. Digital repeat photography for phenological research in forest ecosystems. Agric. For. Meteorol. 2012, 152, 159–177. [CrossRef] 23. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 2002, 83, 195–213. [CrossRef] 24. Beck, P.; Jönsson, P.; Høgda, K.-A.; Karlsen, S.; Eklundh, L.; Skidmore, A. A ground-validated NDVI dataset for monitoring vegetation dynamics and mapping phenology in Fennoscandia and the Kola peninsula. Int. J. Remote Sens. 2007, 28, 4311–4330. [CrossRef] 25. Van der Walt, S.; Colbert, S.C.; Varoquaux, G. The NumPy array: A structure for efficient numerical computation. Comput. Sci. Eng. 2011, 13, 22–30. [CrossRef] 26. Richardson, A.D.; Hufkens, K.; Milliman, T.; Aubrecht, D.M.; Chen, M.; Gray, J.M.; Johnston, M.R.; Keenan, T.F.; Klosterman, S.T.; Kosmala, M.; et al. Tracking vegetation phenology across diverse North American biomes using PhenoCam imagery. Sci. Data 2018, 5, 180028. [CrossRef] [PubMed] 27. Friedl, M.; Gray, J.; Sulla-Menashe, D. MCD12Q2 MODIS/Terra+ Aqua Land Cover Dynamics Yearly L3 Global 500m SIN Grid V006. Available online: http://doi.org/10.5067/MODIS/MCD12Q2.006 (accessed on 11 November 2020). 28. Kowalski, K.; Senf, C.; Hostert, P.; Pflugmacher, D. Characterizing spring phenology of temperate broadleaf forests using Landsat and Sentinel-2 time series. Int. J. Appl. Earth Obs. Geoinf. 2020, 92, 102172. [CrossRef] 29. Schwartz, M.D.; Betancourt, J.L.; Weltzin, J.F. From Caprio’s lilacs to the USA National Phenology Network. Front. Ecol. Environ. 2012, 10, 324–327. [CrossRef] Remote Sens. 2020, 12, 3738 13 of 13 30. Templ, B.; Koch, E.; Bolmgren, K.; Ungersböck, M.; Paul, A.; Scheifinger, H.; Busto, M.; Chmielewski, F.-M.; Hájková, L.; Hodzić, S.; et al. Pan European Phenological database (PEP725): A single point of access for European data. Int. J. Biometeorol. 2018, 62, 1109–1113. [CrossRef] 31. Baldocchi, D.; Falge, E.; Gu, L.; Olson, R.; Hollinger, D.; Running, S.; Anthoni, P.; Bernhofer, C.; Davis, K.; Evans, R.; et al. FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. Bull. Am. Meteorol. Soc. 2001, 82, 2415–2434. [CrossRef] 32. D’Odorico, P.; Gonsamo, A.; Gough, C.M.; Bohrer, G.; Morison, J.; Wilkinson, M.; Hanson, P.J.; Gianelle, D.; Fuentes, J.D.; Buchmann, N. The match and mismatch between photosynthesis and land surface phenology of deciduous forests. Agric. For. Meteorol. 2015, 214, 25–38. [CrossRef] Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).