Next Article in Journal
Real-Time Orthophoto Mosaicing on Mobile Devices for Sequential Aerial Images with Low Overlap
Next Article in Special Issue
Green LAI Mapping and Cloud Gap-Filling Using Gaussian Process Regression in Google Earth Engine
Previous Article in Journal
Integrated GNSS/IMU-Gyrocompass with Rotating IMU. Development and Test Results
Previous Article in Special Issue
RUESVMs: An Ensemble Method to Handle the Class Imbalance Problem in Land Cover Mapping Using Google Earth Engine
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Letter

Improved Estimates of Arctic Land Surface Phenology Using Sentinel-2 Time Series

1
CREAF, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain
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
*
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(22), 3738; https://0-doi-org.brum.beds.ac.uk/10.3390/rs12223738
Submission received: 30 September 2020 / Revised: 1 November 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.

Graphical Abstract

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 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 level-2A, where 1 = saturated or defective; 2 = dark area pixels; 3 = cloud shadows; 6 = water; and, 7−10 = clouds and cirrus (Figure 1b,c).

2.3. Phenology Extraction

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 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)):
u = (Vmin − Vmax) × p + Vmin,
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 mid-greendown of the growing season. The threshold metrics estimated with 50% of the amplitude are less affected by biases due to discontinuities in the time series [13]. We applied two variants of the threshold method as follows:
  • 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 exceeded u; then, we applied a linear interpolation between this first observation when the vegetation index was >u and the preceding observation from which to estimate SoS as the value >u (Figure 2c). For the EoS, the linear interpolation was applied between the latest date when the vegetation index was >u, and the subsequent observation, where EoS corresponded to the linearly interpolated value >u.
  • Threshold method after smoothing. The time series data were smoothed prior to the extraction of LSP metrics (Figure 2d), as is common practice in LSP estimation to reduce noise and discontinuities of time series data [21]. The criteria for selection of the processing steps were based on the feasibility of their implementation in GEE, without comprising the recreation of the phenology curve (see GEE code in Supplementary Materials). Excessive smoothing of time series may lead to unrealistic recreations of the growing season. We first applied a moving average window, with an average radius of 10 d, every 20 d (Figure 2d); if a pixel in the 20 d composite window was empty due to a lack of valid observations, the window size was increased to 40 d. Next, we applied a cubic interpolation to convert the 20 d composites to a daily time series. The threshold was estimated from the amplitude of the interpolated time series, rather 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:
GCC   =   Green Blue   +   Green   +   Red   ,
NDVI   =   NIR     Red NIR +   Red   ,
EVI   = 2.5   ( NIR     Red ) NIR +   6   ×   Red 7.5   ×   Blue + 1   ,   and
NDPI   =   NIR   ( alpha   ×   Red + ( 1 alpha )   ×   S W I R 2 ) NIR + ( alpha   ×   Red + ( 1 alpha )   ×   S W I R 2 )
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.

2.5. Implementation of Land Surface Phenology Algorithms in Google Earth Engine

The two variants of the threshold method were vectorized for their correct implementation in GEE. The vectorization is the process of transforming a code so that all components of an array are processed simultaneously [25]. This concept is contrary to the commonly used practice in LSP estimation, in which time series are processed separately, pixel by pixel, in a for loop. The use of for loops are, however, highly discouraged in GEE in preference for the recommended map functions. For instance, the moving average for the 20-day composition was implemented as a function that maps a list of dates (see code in Supplementary Materials). The function takes the dates as the input arguments to filter the Sentinel-2 collection with a window size of 20 days and, then, makes the average out of the selected images. The purpose of the map function is that each element of the array, in this case the dates, is processed separately and, consequently, each 20-day composition is generated at the same time.

2.6. Validation with PhenoCam

The SoS and EoS metrics extracted from Sentinel-2 were compared with the same metrics estimated from the near-surface reflectances of the PhenoCam network [26]. PhenoCam provides 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 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).
LSP maps generated using Sentinel-2 EVI (Figure 4a,b) showed a high level of similarity at the continental scale compared with the same phenology metrics estimated with MODIS time series (Supplementary Figure S8). However, at the local scale, differences were apparent particularly in the EoS; similarly, the map of length of season (LoS) (Figure 5a), which represented the difference between EoS and SoS, showed spatial patterns and latitudinal gradients that are expected across the region, with shorter LoS in the surrounding areas of the Kara Sea and Arctic Archipelago (Figure 5b), and in elevated areas, such as the Central Siberian Plateau (Figure 5c). At the local level, the SoS and EoS maps derived from Sentinel-2 showed a high degree of detail that moderate scales of resolution were unable to capture. For instance, Figure 5d shows the positive relationship between LoS and elevation in the Scandinavian mountains. The altitudinal gradient is observed at the canopy scale, which allows the analysis of LSP metrics by vegetation type. Figure 5e shows another example of the effect of topography on the LoS in the Ural mountains. In this second example, the LoS is different in both sites of the mountain range depending on the presence of glaciers. Supplementary Figure S9 shows the different spring growth onset depending on the land cover type, tundra shrublands and herbaceous cover in the delta of the Lena River.
The comparison between the LSP metrics estimated using MODIS and Sentinel-2 resized to 500 m, in which MODIS pixels were filtered when variance of the Sentinel-2 phenology estimates within the 500 m pixel were >5 d2, showed a high level of similarity in homogeneous landscapes (Figure 6), regardless of data smoothing (Figure 6a,b), and greater similarity for SoS than EoS (e.g., threshold method with smoothing SoS and EoS RMSE: 8.1 and 10.4 d, respectively). The comparison of our MODIS estimates and the MCD12Q2v6 product (Supplementary Figure S10) showed a slight bias towards EoS (threshold method with smoothing ME: 4.9 and 6.4 d, respectively) and similarity with Sentinel-2 for SoS and EoS (threshold method without smoothing RMSE = 10.2 and 9.8 d for SoS and 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.

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 the time series during the green-up period. In contrast, the cloud coverage was more persistent in late summer and autumn (39.4% of pixels with >40-d gap), and this high level of cloud occurrence may explain the lower similarity between EoS estimated using Sentinel-2 and MODIS. Such discontinuities in time series make the use of gap-filling techniques and robust smoothing techniques necessary for the accurate estimation 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].
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 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 https://0-www-mdpi-com.brum.beds.ac.uk/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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef]
  5. Peñuelas, J.; Filella, I. Phenology feedbacks on climate change. Science 2009, 324, 887–888. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef]
  8. Helman, D. Land surface phenology: What do we really ‘see’from space? Sci. Total Environ. 2018, 618, 665–673. [Google Scholar] [CrossRef]
  9. Witze, A. The Arctic is burning like never before-and that’s bad news for climate change. Nature 2020, 585, 336–337. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef]
  12. Jin, H.; Eklundh, L. A physically based vegetation index for improved monitoring of plant phenology. Remote Sens. Environ. 2014, 152, 512–525. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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://0-doi-org.brum.beds.ac.uk/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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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://0-doi-org.brum.beds.ac.uk/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. [Google Scholar] [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. [Google Scholar] [CrossRef]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef]
Figure 1. Revisit time of Sentinel-2 (a) and maximum discontinuity in the Sentinel-2 Level 2A time series after cloud masking for spring and early summer (1 May to 15 July) (b) and late summer and autumn (15 July to 30 September) 2019 (c).
Figure 1. Revisit time of Sentinel-2 (a) and maximum discontinuity in the Sentinel-2 Level 2A time series after cloud masking for spring and early summer (1 May to 15 July) (b) and late summer and autumn (15 July to 30 September) 2019 (c).
Remotesensing 12 03738 g001
Figure 2. Illustration of vegetation dynamics for a randomly selected pixel in the Arctic. (a) True color compositions of Sentinel-2 time series data during 2019; (b) Sentinel-2 normalized difference vegetation index (NDVI) time series during 2019, with reference to true color compositions shown in (a); (c) Sentinel-2 NDVI time series after cloud masking and snow reclassification, where start of season (SoS) and end of season (EoS) correspond to the earliest and latest linearly interpolated NDVI values, respectively, above a threshold defined dynamically as 50% of the annual amplitude; and, panel (d) application of the threshold method over a smoothed and interpolated time series.
Figure 2. Illustration of vegetation dynamics for a randomly selected pixel in the Arctic. (a) True color compositions of Sentinel-2 time series data during 2019; (b) Sentinel-2 normalized difference vegetation index (NDVI) time series during 2019, with reference to true color compositions shown in (a); (c) Sentinel-2 NDVI time series after cloud masking and snow reclassification, where start of season (SoS) and end of season (EoS) correspond to the earliest and latest linearly interpolated NDVI values, respectively, above a threshold defined dynamically as 50% of the annual amplitude; and, panel (d) application of the threshold method over a smoothed and interpolated time series.
Remotesensing 12 03738 g002
Figure 3. Comparison of start of season (SoS) and end of season (EoS) dates between PhenoCam and four vegetation indices estimated with Sentinel-2. The comparison was performed for three PhenoCam sites in the tundra biome for the years 2019 and 2020. The vegetation indices were the green chromatic coordinate (GCC) in PhenoCam and the GCC, normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and normalized difference phenology index (NDPI) for Sentinel-2. The phenology metrics were extracted with a 50% threshold method without time series 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).
Figure 3. Comparison of start of season (SoS) and end of season (EoS) dates between PhenoCam and four vegetation indices estimated with Sentinel-2. The comparison was performed for three PhenoCam sites in the tundra biome for the years 2019 and 2020. The vegetation indices were the green chromatic coordinate (GCC) in PhenoCam and the GCC, normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and normalized difference phenology index (NDPI) for Sentinel-2. The phenology metrics were extracted with a 50% threshold method without time series 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).
Remotesensing 12 03738 g003
Figure 4. Maps of (a) start of season (SoS) and (b) end of season (EoS) in the Arctic estimated using Sentinel-2 for 2019. The phenology metrics were extracted with the 50% threshold method applied to the non-smoothed Sentinel-2 enhanced vegetation index (EVI) time series.
Figure 4. Maps of (a) start of season (SoS) and (b) end of season (EoS) in the Arctic estimated using Sentinel-2 for 2019. The phenology metrics were extracted with the 50% threshold method applied to the non-smoothed Sentinel-2 enhanced vegetation index (EVI) time series.
Remotesensing 12 03738 g004
Figure 5. Map of length of season (LoS) in the Arctic estimated using Sentinel-2 enhanced vegetation index (EVI) time series for 2019 (a). EVI time series and start of season (SoS) and end of season (EoS) for two randomly selected pixels (b,c); note, the time series in (b) show a shorter LoS than in (c) due to the latitudinal gradient. Arctic digital elevation model map (ArcticDEM) and LoS exemplify the altitudinal gradient of phenology in the Scandinavian mountains and the Ural mountains (d,e) and comparison of LoS estimated with Sentinel-2 and Moderate Resolution Imaging Spectroradiometer (MODIS) (d). The phenology metrics were extracted with the 50% threshold method applied to the non-smoothed Sentinel-2 EVI time series.
Figure 5. Map of length of season (LoS) in the Arctic estimated using Sentinel-2 enhanced vegetation index (EVI) time series for 2019 (a). EVI time series and start of season (SoS) and end of season (EoS) for two randomly selected pixels (b,c); note, the time series in (b) show a shorter LoS than in (c) due to the latitudinal gradient. Arctic digital elevation model map (ArcticDEM) and LoS exemplify the altitudinal gradient of phenology in the Scandinavian mountains and the Ural mountains (d,e) and comparison of LoS estimated with Sentinel-2 and Moderate Resolution Imaging Spectroradiometer (MODIS) (d). The phenology metrics were extracted with the 50% threshold method applied to the non-smoothed Sentinel-2 EVI time series.
Remotesensing 12 03738 g005
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).
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).
Remotesensing 12 03738 g006
Figure 7. Example of flags associated with gaps in the Sentinel-2 time series in the Scandinavian mountains. (a) Revisit time; (b) maximum enhanced vegetation index (EVI) for 2019; (c) flags raised during the estimation of end-of-season (EoS); (d) EoS estimation. Red: non-vegetated pixels (EVI < 0.2); blue: pixels with >40-d gap in time series; EoS was extracted with the 50% threshold method applied to the non-smoothed Sentinel-2 EVI time series.
Figure 7. Example of flags associated with gaps in the Sentinel-2 time series in the Scandinavian mountains. (a) Revisit time; (b) maximum enhanced vegetation index (EVI) for 2019; (c) flags raised during the estimation of end-of-season (EoS); (d) EoS estimation. Red: non-vegetated pixels (EVI < 0.2); blue: pixels with >40-d gap in time series; EoS was extracted with the 50% threshold method applied to the non-smoothed Sentinel-2 EVI time series.
Remotesensing 12 03738 g007
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Descals, A.; Verger, A.; Yin, G.; Peñuelas, J. Improved Estimates of Arctic Land Surface Phenology Using Sentinel-2 Time Series. Remote Sens. 2020, 12, 3738. https://0-doi-org.brum.beds.ac.uk/10.3390/rs12223738

AMA Style

Descals A, Verger A, Yin G, Peñuelas J. Improved Estimates of Arctic Land Surface Phenology Using Sentinel-2 Time Series. Remote Sensing. 2020; 12(22):3738. https://0-doi-org.brum.beds.ac.uk/10.3390/rs12223738

Chicago/Turabian Style

Descals, Adrià, Aleixandre Verger, Gaofei Yin, and Josep Peñuelas. 2020. "Improved Estimates of Arctic Land Surface Phenology Using Sentinel-2 Time Series" Remote Sensing 12, no. 22: 3738. https://0-doi-org.brum.beds.ac.uk/10.3390/rs12223738

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop