Next Article in Journal
Automating the Management of 300 Years of Ocean Mapping Effort in Order to Improve the Production of Nautical Cartography and Bathymetric Products: Shom’s Téthys Workflow
Previous Article in Journal
Land Use and Land Cover Change in the Vaal Dam Catchment, South Africa: A Study Based on Remote Sensing and Time Series Analysis
Previous Article in Special Issue
Global Research Trends for Unmanned Aerial Vehicle Remote Sensing Application in Wheat Crop Monitoring
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Google Earth Engine Algorithm to Map Phenological Metrics in Mountain Areas Worldwide with Landsat Collection and Sentinel-2

1
GEO4Agri DISAFA Laboratory, Department of Agricultural, Forest and Food Sciences (DISAFA), Università degli Studi di Torino, Largo Paolo Braccini 2, 10095 Grugliasco, Italy
2
Earth Observation Valle d’Aosta—eoVdA and IN.VA spa, Località L’Île-Blonde, 5, 11020 Brissogne, Italy
3
Istituto Zooprofilattico Sperimentale Piemonte, Liguria e Valle d’Aosta (IZS PLV) SC Aosta—CeRMAS (National Reference Center for Wildlife Diseases), Località Amerique, 7/C, 11020 Quart (AO), Italy
*
Author to whom correspondence should be addressed.
Geomatics 2023, 3(1), 221-238; https://doi.org/10.3390/geomatics3010012
Submission received: 13 December 2022 / Revised: 14 February 2023 / Accepted: 17 February 2023 / Published: 21 February 2023

Abstract

:
Google Earth Engine has deeply changed the way in which Earth observation data are processed, allowing the analysis of wide areas in a faster and more efficient way than ever before. Since its inception, many functions have been implemented by a rapidly expanding community, but none so far has focused on the computation of phenological metrics in mountain areas with high-resolution data. This work aimed to fill this gap by developing an open-source Google Earth Engine algorithm to map phenological metrics (PMs) such as the Start of Season, End of Season, and Length of Season and detect the Peak of Season in mountain areas worldwide using high-resolution free satellite data from the Landsat collection and Sentinel-2. The script was tested considering the entire Alpine chain. The validation was performed by the cross-computation of PMs using the R package greenbrown, which permits land surface phenology and trend analysis, and the Moderate-Resolution Imaging Spectroradiometer (MODIS) in homogeneous quote and land cover alpine landscapes. MAE and RMSE were computed. Therefore, this algorithm permits one to compute with a certain robustness PMs retrieved from higher-resolution free EO data from GEE in mountain areas worldwide.

1. Introduction

Google Earth Engine (GEE) has deeply changed the way in which Earth observation data are processed, allowing the analysis of wide areas in a faster and more efficient way than ever before. Since its inception, many functions and applications have been implemented and developed by a rapidly expanding community [1]. GEE is a cloud-based, global-scale geospatial analysis platform that provides Google’s vast array of solutions to a wide range of serious social issues, including deforestation, drought, disasters, disease, food security, water management, and climate monitoring [2,3,4,5]. Thanks to the power of cloud computing, GEE was designed to support traditional remote-sensing scientists as well as a wider audience who lack the technical ability to leverage traditional supercomputers and large cloud computing resources. As an integration platform, it is unique in this field [6]. This platform has permitted the development of a tremendous amount of geospatially based applications and research [7,8].
A cloud-based platform such as GEE makes it simple to use high-performance computing tools for processing very large amounts of geospatial information without having to deal with the IT issues that are currently plaguing either. Additionally, and in contrast to other supercomputing facilities, Earth Engine was created to make it simple for researchers to share their findings with other scholars, policymakers, non-governmental organizations, field personnel, and even the general public. Users do not need to be experts in application development, web programming, or HTML to deploy interactive apps or produce systematic data products once an algorithm has been developed on Earth Engine [1,9].
It is worth noting that GEE has a high-performance, intrinsically parallel computation service that is co-located with a multi-petabyte data catalog ready for analysis. It can be accessed and managed using an application programming interface (API) that is accessible via the Internet, together with a web-based interactive development environment (IDE) that enables quick prototyping and result visualization [4,10]. The data catalog contains a sizable collection of freely accessible geospatial datasets, such as observations from numerous optical and non-optical satellites and aerial imaging systems, environmental variables, weather and climate forecasts and hindcasts, land cover, and topographic and socioeconomic datasets. All of these data have been preprocessed into a usable but information-preserving format that enables quick access and eliminates many data-management-related obstacles [11,12,13]. Finally, the Earth Engine API offers a library of operators that users may employ to retrieve and analyze both their own private data as well as data from the public catalog. In order to enable high-throughput analysis, these operators are implemented in a sizable parallel processing system that automatically divides and distributes computations. Either a thin client library or a web-based interactive development environment built on top of that client library is employed by users to access the API [14].
Despite the growing number of GEE codes related to several applications, there is still a lack of codes related to phenological metric (PM) estimation [15]. Phenological metrics are numerical parameters that allow the monitoring of the seasonal manifestations of the vegetative process and, therefore, the definition of the different phenological phases [16]. PMs are numerical parameters that can be deduced from the temporal profile of remotely sensed vegetational spectral indices such as NDVI or EVI (for a single pixel) in correspondence with significant moments of the year. Among these, the most commonly used PMs, for whose identification tools and software have been developed in the agronomic, forestry, and environmental fields, are:
(a)
The start of the season (SOS = Start of Season), indicated with reference to the date (often expressed in the form DOY, day of year) on which vegetative activity is observed to begin in a certain position.
(b)
The end of the season (EOS = End of Season), indicated with reference to the date on which the vegetative activity in that position is observed to end.
(c)
The length of the season (LOS = Length of Season), defined as the difference in groups between EOS and SOS.
(d)
The maximum vegetation index value (MAX_VI), identified between SOS and EOS.
(e)
The date on which MAX_VI occurs (MAX_DOY).
(f)
The amplitude of the season (SA = Season Amplitude), which defines the difference between the maximum and minimum values of the index in the considered season.
(g)
Total productivity (TP = Total Productivity), which defines the integral of the interpolated profile over the entire year.
(h)
Seasonal productivity (SP = Seasonal Productivity or SMI = Small Integral), which defines the integral of the interpolated profile of the VI between SOS and EOS.
(i)
The rate of increase at the beginning of the season (Rate of Increase), defined, with reference to the growing part of the phenological bell (left), as the difference between the index values corresponding, respectively, to 80% and 20% SA, divided by the corresponding time interval.
(j)
The rate of decrease at the end of the season (Rate of Decrease), defined, with reference to the decreasing part of the phenological bell (right), as the absolute value of the difference between the index values corresponding, respectively, to 80% and 20% SA, divided by the corresponding time interval.
PM monitoring through remote sensing has garnered great interest in the last few decades due to climate change, both in the forestry and agricultural sectors [15,17,18]. Land surface phenology (LSP) is an important research field in terrestrial remote sensing and has become an indispensable approach in global change research, as evidenced by the many important scientific findings supported by LSP in recent decades [15,16,19]. LSP involves the use of remote sensing to monitor seasonal dynamics in vegetated land surfaces and to retrieve phenological metrics (transition dates, rate of change, annual integrals, etc.). LSP is an essential indicator of global change and has played a pivotal role in shaping our understanding of how terrestrial ecosystems are responding to climate change and human activities. Both regional and global LSP products have been routinely generated and have played prominent roles in modeling crop yield, ecological surveillance, identifying invasive species, modeling the terrestrial biospheric processes, and assessing global change impacts on urban and natural ecosystems [20,21].
Recent advances in field and spaceborne sensor technologies, as well as data fusion techniques, have enabled novel LSP retrieval algorithms that refine LSP retrievals at even higher spatiotemporal resolutions, providing new insights into ecosystem dynamics. Meanwhile, the rigorous assessment of the uncertainties in LSP retrievals is ongoing, and efforts to reduce these uncertainties have also formed an active research field. In addition, open-source software and hardware are being developed and have greatly facilitated the use of LSP metrics by scientists beyond the remote-sensing community [22]. Most PM and LSP studies have involved a temporal scale greater than a decade, with a few using NOAA/AVHRR data for a period longer than three decades [23,24]; despite this, the spatial resolution is not as suitable as that of Himawari to map PMs at the parcel level or in geomorphologically complex areas such as mountains. Very-high-resolution EO data such as PlanetScope data are starting to be used for PM assessment [25]. However, their application remains rare both due to the quality of the data in the Alpine area and because they are free only for research centers upon specific request, which prevents their full and widespread commercial use in developing rural contexts [26].
However, it is worth noting that few studies have evaluated the impacts of satellite products with different spatial resolutions on LSP extraction over regions with a heterogeneous topography. To bridge this knowledge gap, studies such as [24] have employed four types of satellite data with different spatial resolutions (250, 500, and 1000 m MODIS NDVI during the period 2001–2020 and ~10 km GIMMS3g during the period 1982–2015) to investigate the LSP changes taking place in an alpine context such as the Loess Plateau. These studies showed that the MODIS-based start of the growing season (SOS) and end of the growing season (EOS) were highly correlated with the ground-observed data, presenting R values of 0.82 and 0.79, respectively (p < 0.01), while the GIMMS3g-based phenology signal performed badly (R < 0.50 and p > 0.05). Spatially, the LSP that was derived from the MODIS products generated more reasonable spatial distributions. The inter-annual averaged MODIS SOS and EOS presented overall advanced and delayed trends, respectively, during the period 2001–2020. More than two thirds of the SOS advances and EOS delays occurred in grasslands, which determined the overall phenological changes across the entire Loess Plateau. However, both inter-annual trends (SOS and EOS) derived from the GIMMS3g data were opposite to those seen in the MODIS results. There were no significant differences among the three MODIS datasets (250, 500, and 1000 m) with regard to a bias lower than 2 days and an RMSE lower than 1 day. Furthermore, it was found that the phenology derived from the data with a 1000 m spatial resolution in the heterogeneous-topography regions was feasible. Therefore, Landsat and Sentinel EO data should be more suitable in an alpine context [22]. In GEE, the existing algorithms are mainly focused on PM extraction for agricultural crops [5,27,28,29,30,31,32,33,34,35]. However, an application was realized considering arctic areas to monitor PM shifts due to the effects of climate change on vegetation in very sensitive ecosystems [17,20,36,37,38].
Nevertheless, existing algorithms that permit PM computation with an high degree of accuracy in mountain areas have been developed mainly in R with package extensions such as: ndvits [39], greenbrown [40,41,42], phenor [43,44], phenofit [45,46,47], and phenopix [48] related to TIMESAT [49,50]. Currently, no algorithm has been developed in GEE using free high-resolution Earth observation data. Therefore, starting from this evidence, the main aim of this work was to fill this gap by developing an open-source Google Earth Engine algorithm to map phenological metrics (PMs), in particular, Start of Season, End of Season, Length of Season, and Peak of Season, in mountain areas worldwide using high-resolution free satellite data from the Landsat collection and Sentinel-2.
The main objective of this work and its main novelty was the development of a completely open-source algorithm capable of extracting PMs from high-resolution Earth observation data in geomorphologically complex areas [7,8] without having to download any type of satellite data by exploiting the power and performance of cloud processing platforms such as Google Earth Engine [19,43,51,52]. The algorithm developed in GEE should be able to attain the same results as those obtained from PM packages in R in an alpine context. Actually, all PM calculation systems require one to download locally at least one year of satellite data. Conversely, R packages such as rgee [34,45,48] connected to Google Earth Engine permit one to compute PMs more rapidly. However, a higher-performance computer workstation and good skills in R is still necessary, especially when wide areas are considered, such as entire alpine chains.
Rural agriculture and forestry are very important in mountain areas. Therefore, the code adopted and developed in GEE permits one to compute PMs more easily in a cloud-based platform, helping technicians in their planning and management in these areas without the necessity for a higher-performance computer workstation and obtaining results for wide areas with high temporal and spatial resolution data. Moreover, mountain areas are one of the most sensitive proxies of climate change [53], and their monitoring is a key factor for realizing effective adaptation and mitigation policies worldwide [26,54,55,56,57]. In fact, it is worth noting that 24% of the Earth’s landmass can be considered mountainous [58], and concentrating scientific studies and climate adaptation policies in these areas is essential to ensure ecosystem services and global geopolitical stability in relation to conflicts over resources such as water considering future climate projections [59,60].

2. Materials and Methods

2.1. Development Area

The GEE algorithm was realized and tested considering the entire Alpine chain in Europe (please see Figure 1). The Alpine chain is defined for much of its distance by the watershed between the drainage basin of the Po in Italy on one side and the divide formed by the Rhone, the Rhine, and the Danube on the other side. Further east, the watershed is between the Adige and the Danube, before heading into Austria and draining on both sides into the Danube. For much of its distance, the watershed lies on or close to the Italian border, although there are numerous deviations—notably, the Swiss canton of Ticino, which lies south of the range in the Po river basin. For only a small portion of its total distance does the Alpine divide form a part of the main European watershed, in the central section where the watershed is between the Po and the Rhine. The Alps are generally divided into the Eastern Alps and the Western Alps, split along a line between Lake Como and Lake Constance, following the Rhine valley. The Eastern Alps (with the main ridge being elongated and broad) belong to Austria, Germany, Italy, Slovenia, and Switzerland. The Western Alps are higher, but their central chain is shorter and very curved; they are located in France, Italy, and Switzerland. Piz Bernina (4049 m) is the highest peak of the Eastern Alps, while the highest peak of the Western Alps is Mont Blanc (4810.45 m). An overview of the code development area is reported in [61].

2.2. Earth Observation Data

The code was able to work with NASA USGS Landsat products—in particular, Level-2 collections (surface reflectance) and the ESA Copernicus Sentinel-2 mission collection within GEE. Therefore, the code could be run using the following collections:
Sentinel-2 MSI data from the ee.ImageCollection (“COPERNICUS/S2_SR_HARMONIZED”), which include Sentinel-2A and S2B images from 28 March 2017 to present day with a general temporal resolution of 5 days and a spatial resolution between 10 and 60 m in the spectral band considered. Sentinel-2 is a wide-swath, high-resolution, multi-spectral imaging mission supporting Copernicus Land Monitoring studies, including the monitoring of vegetation, soil, and water cover, as well as the observation of inland waterways and coastal areas. The Sentinel-2 L2 data were downloaded from Sci-hub, computed by running sen2cor, and ingested in GEE. This collection contains 12 spectral bands representing surface reflectance scaled by 10,000. Moreover, three QA bands are present, with one (QA60) being a bitmask band with cloud mask information. The developed algorithm adopted the SCL masks because they presented the highest spatial resolution, equal to 20 m. These masks are provided by ESA and available in the GEE collection, with snow and ice, clouds, shadows, and saturated or defective pixels having been masked in each image. This collection was adopted instead of COPERNICUS/S2_SR because, after 25 January 2022, Sentinel-2 scenes with a processing baseline ‘04.00’ or above had their DN (value) range shifted by 1000. The harmonized collection shifted the data in newer scenes to be in the same range as in older scenes.
Landsat 4 MSS data from the ee.ImageCollection (“LANDSAT/LT04/C02/T1_L2”), which include the Landsat 4 collection from 22 August 1982 to 24 June 1993 with a general temporal resolution of 16 days and a geometrical resolution of 30 m (except for the thermal band, which was not considered in this work). This dataset contains atmospherically corrected surface reflectance images. These images contain 4 visible and near-infrared (VNIR) bands and 2 short-wave infrared (SWIR) bands processed to generate orthorectified surface reflectance. It also contains pixel quality bands, QA bands, and radiometric saturation bands with the same spatial resolution adopted to mask cloud, shadows, snow and ice, and defective or saturated pixels. The Landsat 4 collection was created with the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (version 3.4.0) and then uploaded in GEE.
Landsat 5 MSS data from the ee.ImageCollection (“LANDSAT/LT05/C02/T1_L2”), which include the Landsat 5 collection from 16 March 1984 to 5 May 2012 with a general temporal resolution of 16 days and a geometrical resolution of 30 m (except for the thermal band, which was not considered in this work). This dataset contains atmospherically corrected surface reflectance images. These images contain 4 visible and near-infrared (VNIR) bands and 2 short-wave infrared (SWIR) bands processed to generate orthorectified surface reflectance. It also contains pixel quality bands, QA bands, and radiometric saturation bands with the same spatial resolution adopted to mask cloud, shadows, snow and ice, and defective or saturated pixels. The Landsat 5 collection was created with the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (version 3.4.0) and then uploaded in GEE.
Landsat 7 MSS data from the ee.ImageCollection (“LANDSAT/LE07/C02/T1_L2”), which include the Landsat 7 collection from 28 May 1999 to 31 May 2003 (when the scan line corrector (SLC) failed). After this, the products have data gaps, but they are still useful and maintain the same radiometric and geometric corrections as the data collected prior to the SLC failure (though only in those areas of the image not affected by the failure). Thus, to compute PMs, it is strongly recommended to adopt only images without SCL failure. These data have a general temporal resolution of 16 days and a geometrical resolution of 30 m (except for the thermal band, which was not considered in this work). This dataset contains atmospherically corrected surface reflectance images. These images contain 4 visible and near-infrared (VNIR) bands and 2 short-wave infrared (SWIR) bands processed to generate orthorectified surface reflectance. It also contains pixel quality bands, QA bands, and radiometric saturation bands with the same spatial resolution adopted to mask cloud, shadows, snow and ice, and defective or saturated pixels. The Landsat 7 collection was created with the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (version 3.4.0) and then uploaded in GEE.
Landsat 8 OLI data from the ee.ImageCollection (“LANDSAT/LC08/C02/T1_L2”), which include the Landsat 8 collection from 18 March 2013 to present day with a general temporal resolution of 16 days and a geometrical resolution of 30 m (except for the thermal band, which was not considered in this work). This dataset contains atmospherically corrected surface reflectance images. These images contain 5 visible and near-infrared (VNIR) bands and 2 short-wave infrared (SWIR) bands processed to generate orthorectified surface reflectance. It also contains pixel quality bands, QA bands, and radiometric saturation bands with the same spatial resolution adopted to mask cloud, shadows, snow and ice, and defective or saturated pixels. The Landsat 8 products were created with the Land Surface Reflectance Code (LaSRC) and then uploaded in GEE.
Landsat 9 MSS data from the ee.ImageCollection (“LANDSAT/LC09/C02/T1_L2”), which include the Landsat 9 collection from 31 October 2021 to present day with a general temporal resolution of 16 days and a geometrical resolution of 30 m (except for the thermal band, which was not considered in this work). This dataset contains atmospherically corrected surface reflectance images. These images contain 5 visible and near-infrared (VNIR) bands and 2 short-wave infrared (SWIR) bands processed to generate orthorectified surface reflectance. It also contains pixel quality bands, QA bands, and radiometric saturation bands with the same spatial resolution adopted to mask cloud, shadows, snow and ice, and defective or saturated pixels. The Landsat 8 products were created with the Land Surface Reflectance Code (LaSRC) and then uploaded in GEE.
Terra and Aqua MODIS from the ee.ImageCollection (“MODIS/006/MCD12Q2”). This collection, informally known as MODIS Global Vegetation Phenology products, provides estimates of the timing of vegetation phenology at a global scale. Additionally, it provides information related to the range and summation of the enhanced vegetation index (EVI) computed from MODIS surface reflectance data at each pixel. It identifies the onset of greenness, greenup midpoint, maturity, peak greenness, senescence, greendown midpoint, dormancy, EVI2 minimum, EVI2 amplitude, integrated EVI2 over a vegetation cycle, and overall and phenology-metric-specific quality information. The MCD12Q2 Version 6 data product was derived from time series of the 2-band enhanced vegetation index (EVI2) calculated from MODIS nadir bidirectional reflectance distribution function (BRDF)-adjusted reflectance (NBAR). Vegetation phenology metrics were identified for up to two detected growing cycles per year. For pixels with more than two valid vegetation cycles, the data represent the two cycles with the largest NBAR-EVI2 amplitudes. This dataset was adopted to assess the quality of the results obtained compared to the metrics retrieved by the greenbrown R package.
Finally, the Shuttle Radar Topography Mission (SRTM) DTM available in GEE. The Shuttle Radar Topography Mission digital elevation data resulted from an international research effort to obtain digital elevation models on a near-global scale. This SRTM V4 product (SRTM Plus) was provided by NASA JPL at a resolution of 1 arc-second (approximately 30 m).

2.3. Phenological Metrics Computation

PMs were computed considering as a reference the year 2021. Specifically, Sentinel-2 level-2A data and Landsat 8 OLI were adopted. PMs estimated with these EO collections were compared with the same metrics extracted from MCD12Q2 V6 Land Cover Dynamics products with a 500 spatial resolution. Maximum discontinuity in the time series was plotted after non-valid observations had been filtered using values of the quality band scene classification layer (SCL) provided in Sentinel-2 and QA in Landsat.
In the algorithm, the PMs could be computed from the following vegetation spectral indexes:
Normalized-difference vegetation index (NDVI) [62,63,64,65,66]
NDVI = NIR RED NIR + RED
Enhanced vegetation index (EVI) [67,68]
EVI = 2.5 × NIR RED NIR + 6 × RED 7.5 × BLUE + 1  
Green chromatic coordinate (GCC) [69]
GCC = GREEN RED + GREEN + BLUE
Normalized-difference phenology index (NDPI) [70]
NDPI = NIR ( α × RED + 1 α × SWIR 2 NIR + ( α × RED + 1 α × SWIR 2
where α was set to 0.74 for MODIS [71] but re-defined to 0.51 for Sentinel-2 and 0.56 for Landsat. The normalized-difference phenology index was adopted to overcome the snow observations in the collection that normally affect mountain areas.
Specifically, a commonly used threshold approach was adopted [19,72]. This approach determines the SOS and EOS as the first and final days of the season, respectively, on which a threshold τ is overtaken. This τ value can be assigned dynamically or as a constant for each pixel [73]. The altimetric gradient strongly influences phenology; therefore, in this work τ was retrieved as a dynamic value related to the yearly time-series amplitude:
τ   =   φ   VI min   VI max +   VI min
where τ is the dynamic value that depends on the annual amplitude of the time series; VImin and VImax are the minimum and maximum vegetation index yearly values in the time series, respectively; and φ is a given proportion (%) of the amplitude. A value of φ = 0.5 was adopted as the mid-greenup and mid-greendown of the growing season. Biases resulting from time-series discontinuities are less of an issue for threshold metrics computed with 50% of the amplitude [19].
In order to eliminate noise and discontinuities, the time-series data were filtered and smoothed before the extraction of PMs [73]. Unrealistic recreations of the growing season may result from overly smoothing time series; therefore, a moving average window with a mean radius of 10 days was realized every 20 days. In cases where a pixel in the 20-day composite window was empty because of a lack of valid observations, the window size was increased to 40 days. Afterwards, a cubic spline interpolation to convert the 20-day composites to a daily time series was applied. It is worth noting that the threshold was established using the interpolated time series amplitude rather than daily observations, and the SOS and EOS were calculated as the first and last days, respectively, in the interpolated time series that overtook the dynamic threshold. Due to the primarily evergreen nature of the formerly snow-covered vegetation canopy, the NDVI and EVI for alpine vegetation exhibit a substantial increase after snow melting. Therefore, the PM estimation was complicated by the time-series break related to the change from snow to vegetation and from vegetation to snow. It is possible to incorrectly attribute the SOS and EOS to the snow transition dates rather than to the real vegetation dynamics. In order to maintain a constant threshold value, considering the alpine environment and as suggested by investigations in high-latitude areas [19,40], we reclassified the snow and post-thaw readings. Each vegetation index’s unique minimum value was used to reclassify the snow observations: GCCmin = 0.31, NDVImin = 0.39, EVImin = 0.2, and NDPImin = 0.24.
The mean value of the first Sentinel-2 and Landsat snow-free observation of the year, taken from 400 randomly dispersed sites, was used to estimate these parameters. When using the best alpha for the NDPI, the NDPImin matched the mean NDPI value of the snow observations. The PMs were not computed for the pixels that displayed a maximum vegetation index in the time series below the lowest snow value.
In order to ensure proper functioning in GEE, the threshold approach was vectorized. The vectorization technique involves changing a code so that every element of an array is processed at once [74]. This idea runs counter to the standard method for PM estimation, which involves processing each time series separately, pixel by pixel, in a for loop. However, GEE strongly discourages the use of for loops in favor of the suggested map functions—for instance, a function that converts a series of dates to the moving average for the 20-day composition (see Algorithm S1 in Supplementary Materials). The function uses dates as input arguments to filter the Earth observation collection with a 20-day window size before averaging the chosen pictures. Each element of the array, in this case the dates, is processed separately by the map function in order to construct each 20-day composition simultaneously.

2.4. Algorithm Validation

The GEE algorithm was tested by performing PM extraction from other satellite missions and R packages. Therefore, the mean absolute error (MAE) and root mean square error (RMSE) were computed as follows:
MAE = i = 1 n p i o i n        
RMSE = i = 1 n p i o i 2 n
where pi is the prediction (assumed as the PMs computed with the GEE script code); oi is the observed true value (assumed as the PMs computed through MODIS and the greenbrown R package with a higher satellite resolution); and n is the number of samples (see Table 1).

3. Results

The algorithm developed permitted us to compute and extract PMs directly from GEE without the necessity to download EO data locally (in this case, using Sentinel-2 and the Landsat collection for the reference year of interest). The mean processing time necessary to compute and download PMs (SOS, EOS, LOS, and POS) into the Google Drive was around 20–25 min for the entire Alpine chain with a spatial resolution of 20 m and 30 m, respectively, for yearly Sentinel-2 and Landsat data. The only mandatory conditions were the possession of an internet connection, a GEE account, and basic knowledge of GEE (to set the code developed according to one’s needs). The code worked well, and the maps shown below (please see Figure 2) were computed for SOS, EOS, and LOS.
In order to perform a strong validation, 500 points were randomly generated from the vegetated areas detected by a threshold of NDVI > 0.3 assumed as vegetated in SAGA GIS v.8.3.0 over the entire area of study (AOI), in this case the Alpine chain. A snapshot of the AOI is shown below (please see Figure 3).
The PMs, namely SOS, EOS, LOS, and POS, were computed using the developed GEE code. At the same time, by adopting the greenbrown package in R Studio, the metrics were computed starting from Sentinel-2 and Landsat VI time series. MODIS/006/MCD12Q2 and greenbrown R PMs (from Landsat and Sentinel) were considered as the validation PM set. For validation, in addition to MODIS, we also considered Sentinel and Landsat data by computing PMs in R due to the fact that the spatial and temporal resolution was different for these sensors and we wished to check the consistency by preserving these aspects. The MAE and RMSE are reported in Table 1.
Both Sentinel and Landsat data seemed to accurately depict PMs, though Sentinel-2 produced better results, probably due to its higher temporal resolution (please see Table 2).
It is worth noting that the wide standard deviations are due to the fact that we summarized the PMs extracted from 500 points. This variability could be explained by the altimetry conditioning the phenological seasons.
The results were also plotted with their relative standard deviations (please see Figure 4 below).
In order to evaluate the consistency of the results, we adopted SRTM DTM to assess the PM errors for all 500 points at the same quote, and an overall RMSE < 2 was found.
The comparison of Start of Season (SOS) and End of Season (EOS) dates computed, respectively, from Sentinel-2 and Landsat retrieved in R with the greenbrown package and from the GEE algorithm from the same satellite showed a strong similarity. Specifically, four vegetation indices were computed: GCC, NDVI, EVI, and NDPI. The phenology metrics were extracted with a 50% threshold method without time-series smoothing. The bias between the R package greenbrown and the GEE code was less the 5%, according to the MAE and RMSE results. Therefore, the GEE algorithm was proven suitable for correctly extracting PMs. Below, the graph is depicted (Figure 5).
The comparison between the PMs estimated using the GEE code from Sentinel-2 resized to 500 m and MODIS phenology products (MCD12Q2), in which MODIS pixels were filtered when the variance of the Sentinel-2 phenology estimates within the 500 m pixel were >5 days, showed a high level of similarity in homogeneous landscapes, regardless of data smoothing, and a greater similarity for SOS than EOS (see Table 1). The comparison of the GEE PMs retrieved from Sentinel and Landsat and MODIS MCD12Q2 phenology products showed a slight bias towards SOS and LOS and similarity for EOS and POS (see Table 1). The graph is depicted below (Figure 6).
The GEE code seemed to be consistent according to the results obtained. The bias between PMs computed in R and GEE, respectively, as well as MODIS, could be considered acceptable, as demonstrated in Table 1 by the mean absolute error (MAE) and the accuracy according to the root mean squared error (RMSE).

4. Discussion

The GEE scripts were able to accurately compute the PMs in geomorphologically complex areas, as suggested by the validation performed. Therefore, the script realized may be adopted to define PMs in alpine areas. Nevertheless, it is worth noting that the code could certainly be improved; in particular, a ground truth validation based on worldwide mountain phenology detection would represent a stronger and more rigorous validation approach. Therefore, we hope that the open-source code will be implemented and a stronger validation performed with ground samples, as previously stated. The results obtained suggest a good degree of confidence in the use of this GEE code to retrieve PMs in alpine areas. Despite the comparison reflecting changes in vegetation greenness that may not always correspond to vegetation pheno-phases or plant productivity, we discovered significant similarities between the PMs computed using Sentinel-2, Landsat, and MODIS. The results obtained were compared with PMs retrieved from the R package greenbrown and phenofit [45,46,47,75], and a difference below 5% was found for each of the PMs considered (SOS, EOS, LOS, and POS) for each collection method, in this case the Sentinel and Landsat missions. Moreover, the method designed for identifying good-quality pixels for additional analysis in phenology research may be used to flag pixels with discontinuities in the time series; however, this flagging method only functions when clouds are effectively filtered. The Earth observation time series may have been overestimated by non-valid observations that were not filtered, namely cloud-contaminated pixels and cloud shadows that were not reflected in the quality band. Additionally, such polluted data might have altered the vegetation index (VI) time series growing season curve and resulted in inaccurate phenology estimations. However, it is important to note that coupling remote sensing with ground data will always be welcome and necessary. In fact, remote sensing and the algorithm developed could be combined to generalize PMs and not follow each phenological phase of each grassland species, for example. The limiting factors of phenological remote sensing are represented on the one hand by the resolution of the images, which do not allow the identification of early phenological phases that began much earlier than a few decades ago, such as the swelling and bursting of leaf buds on trees [76]. On the other hand, the limitations are also related to identifying phenological phases in rare species. It is important to remember that individual species can respond differently to meteorological conditions, and their phenological phases are not necessarily linked. Moreover, as suggested by [77,78,79], studies using satellite imagery in some cases seem to be unable to determine the exact growing season of trees. This is particularly evident in ring-porous trees, for which some research [77,78,79] has shown that there is no clear relationship between the timing of the onset of wood formation and leaf development. Despite this, EO data can be useful to obtain information on phenology at the landscape level worldwide. Therefore, for a complete understanding of phenology, it would be beneficial to use different approaches that complement each other, and the code developed could be one of these approaches to support ground studies and help enhance databases using citizen-science ground observations.
As previously discussed, the code permits one to retrieve PMs from different spectral indexes (such as EVI, NDVI, NDPI, and GCC). It is worth noting that the developed GEE code was tested across the whole study area using the NDVI, because this is the most frequently adopted index and provides simple estimation; easy availability at different spatial and temporal resolutions; and the cancellation of noise caused by the solar angle, topographic illumination, clouds, and atmospheric conditions, especially in alpine areas [80,81]. Despite these several advantages, the NDVI is more saturated at higher biomass levels due to leaf canopy variations [82]. Therefore, as an alternative, we encourage the computation of the EVI to retrieve PMs in order to minimize these errors, as this index improves the estimation of the biomass level under saturation conditions [83]. In addition, the EVI range is more extensive and dynamic and allows the capture of more variations than the NDVI, as it includes the coefficient of resistance term, which corrects the influence of aerosols. Although the EVI has been used as a vegetation proxy in many studies due to its improved performance in many regions across the world [84], the use of the EVI (as well as GCC and NDPI) in mountainous terrain remains limited, and NDVI is generally preferred. It is necessary to keep in mind that in the case of the alpine study of PMs in dense forests, the EVI is strongly advised, while for studies at a higher latitude with the sporadic presence of snow, the NDPI is more suitable, as demonstrated by [70]. Additionally, GCC has similar results to the NDVI [69].
Earth observation optic multispectral higher-resolution data such as Landsat and Sentinel-2 time series may now be used to estimate PMs for wide mountain range regions, such as the Alps, Himalayas, Andes, Urals, Rocky Mountains, Great Dividing Range, Great Escarpment, Apennines, Alaska Range, Scandinavian Alps, Japanese Alps, Hindu Kush, Altai Mountain, Western Ghats, Drakensberg, Aravalli Range, Appalachian Mountains, Jura, and Pyrenees, thanks to the utilization of cloud platforms such as GEE. Approximately 91,000 EO data were utilized in this work to determine the PMs for the European Alps; this kind of platform is a helpful tool that enables the scientific community to examine such rich, high-spatial-resolution time series. Additionally, GEE permits data and code exchange, making it simple for researchers to replicate algorithms, locally examine time series, and therefore afford a real technological transfer. It is worth noting that, when data are limited, it is especially important to exercise caution when using the generated PMs maps (especially when considering SOS and EOS). When there are continuous gaps in the time series, it would be better to combine Sentinel-2 with Landsat-8; otherwise, more sophisticated phenological retrieval methods, such as robust smoothing and gap-filling techniques, should be utilized. In spite of this, sensor harmonization and temporal smoothing may enhance PM retrievals, but at the expense of making preprocessing more difficult and slowing down GEE implementation.
Regarding future research and studies on mountain phenological metrics worldwide or in a given area, such as the area of the Alps considered herein, it would be interesting to use the developed GEE script to analyze spatial variability in vegetation dynamics and acquire a high level of detail at the canopy level, or to test the effects of spatial scaling on temporal changes in phenology and identify homogeneous landscapes where the phenology dynamics are similar. Recent research using multi-annual optical time-series data for vegetation phenology assessment, together with the functions of newly developed packages, has improved our ability to extract PMs from EO data [33,48,75]. Nevertheless, until now, no research has been focused on PM GEE codes in alpine ecosystems, except for the R package Phenopix [48]. However, this R package (quite similar to greenbrown due to its specificity in geomorphologically complex areas) remains applicable only outside of cloud processing environments such as GEE. The computation of PMs using this R package therefore requires the downloading of locally processed satellite time-series data or the use of R packages such as rgee that depend on GEE but require a high-performance workstation. Therefore, to overcome all these issues, the algorithm developed in this work was directly written in the GEE front-end code editor, avoiding all the above-mentioned matters.
One of the possible future developments will certainly consist in the implementation of a function for the direct calculation of phenological metrics in the dashboard of Google Earth Engine, which will facilitate use by all users and be scalable to any type of multispectral optical data. Particular attention will have to be paid to mountainous areas where remote-sensing applications are increasingly complex due to the topographical and morphological characteristics of these places. However, we hope this work will help achieve this goal. In fact, we must not forget that most of the mountain areas in the world not only play a fundamental role from the point of view of biodiversity and a key role in monitoring climate change (and therefore deserve continuous attention, management, safeguarding, and investment), but also represent important centers from a social and economic point of view. This includes ecosystem services that are difficult to monetize to simpler forms of agriculture, animal husbandry, forest management, and subsistence that are often more sustainable and deserving of funding towards sustainable development.
In conclusion, from the results obtained, when using time-series data not polluted by a high number of clouds during the year, the developed script seemed to be a good tool for monitoring and mapping phenological metrics at the field scale in mountain areas across the world.
It is worth noting that very-high-resolution EO data are not currently available in an open-source format. Moreover, this kind of data is still not available in cloud processing platforms such as Google Earth Engine or Planetary Microsoft (except for PlanetScope data, which are only available for CONUS areas). Certainly, it would be interesting if, in the near future, very-high-resolution data are made available in cloud processing platforms for free for research purposes, strengthening our knowledge of mountain areas and PMs at an accurate level.

5. Conclusions

The work carried out made it possible to compute the main phenological metrics starting from high-level satellite data, such as those gathered from the Sentinel-2 and Landsat missions. The added value is represented by the fact that the Javascript code works on Google Earth Engine, thus avoiding the need to download data by operating in the cloud and making operations much faster and scalable to any area of the globe, although the algorithm was tested and designed to work in areas where remote sensing is more limited, i.e., mountainous areas. The results obtained seem to suggest the suitability of the implemented code, as indicated by the low MAE and RMSE values. However, it should be emphasized that the validation was carried out in two ways: (a) first with coarser remote-sensing data, and (b) second through the computation of phenological metrics with the greenbrown package available in R software. The only limiting factor of the results obtained is represented by the lack of ground validation data. Since the code is open-source, we hope in the future to validate it with ground-based phenological data collected in various locations around the world. While acknowledging that much still needs to be improved in terms of both validation and the implementation of other metrics (as well as the direct implementation of this code as a server-side function in GEE), we hope that this work can help agronomists, farmers, climatologists, foresters, veterinarians, and, more generally [75,85,86], territorial planners in forming direct management and planning strategies and territorial policies in the Alpine area [87,88], with a view to adaptation, the mitigation of climate change, and sustainability.
In conclusion, we hope that this open-source algorithm permits users to compute with a certain robustness PMs retrieved from higher-resolution free EO data in GEE, helping them face the new challenges in mountain areas worldwide. The algorithm developed in GEE seemed to be capable of attaining the same robust results that can be obtained from PM packages in R in alpine contexts.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/geomatics3010012/s1, Algorithm S1: The GEE code algorithm developed and the data adopted can be found in this section.

Author Contributions

Conceptualization, T.O.; methodology, T.O.; software, T.O.; validation, T.O. and E.B.M.; formal analysis, T.O.; investigation, T.O., A.V. and D.C.; resources, T.O.; data curation, T.O. and A.V.; writing—original draft preparation, T.O.; writing—review and editing, T.O.; visualization, T.O., A.V. and D.C.; supervision, E.B.M.; project administration, T.O. and E.B.M.; funding acquisition, T.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Google Earth Engine.

Acknowledgments

Special thanks to our GEO4Agri DISAFA Lab colleagues.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. 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]
  2. Amani, M.; Ghorbanian, A.; Ahmadi, S.A.; Kakooei, M.; Moghimi, A.; Mirmazloumi, S.M.; Moghaddam, S.H.A.; Mahdavi, S.; Ghahremanloo, M.; Parsian, S.; et al. Google Earth Engine Cloud Computing Platform for Remote Sensing Big Data Applications: A Comprehensive Review. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2020, 13, 5326–5350. [Google Scholar] [CrossRef]
  3. Kumar, L.; Mutanga, O. Google Earth Engine Applications since Inception: Usage, Trends, and Potential. Remote Sens. 2018, 10, 1509. [Google Scholar] [CrossRef] [Green Version]
  4. Aybar, C.; Wu, Q.; Bautista, L.; Yali, R.; Barja, A. Rgee: An R Package for Interacting with Google Earth Engine. J. Open Source Softw. 2020, 5, 2272. [Google Scholar] [CrossRef]
  5. Guo, Y.; Xia, H.; Pan, L.; Zhao, X.; Li, R. Mapping the Northern Limit of Double Cropping Using a Phenology-Based Algorithm and Google Earth Engine. Remote Sens. 2022, 14, 1004. [Google Scholar] [CrossRef]
  6. Mutanga, O.; Kumar, L. Google Earth Engine Applications. Remote Sens. 2019, 11, 591. [Google Scholar] [CrossRef] [Green Version]
  7. Orusa, T.; Cammareri, D.; Borgogno Mondino, E. A Possible Land Cover EAGLE Approach to Overcome Remote Sensing Limitations in the Alps Based on Sentinel-1 and Sentinel-2: The Case of Aosta Valley (NW Italy). Remote Sens. 2022, 15, 178. [Google Scholar] [CrossRef]
  8. Orusa, T.; Cammareri, D.; Borgogno Mondino, E. A Scalable Earth Observation Service to Map Land Cover in Geomorphological Complex Areas beyond the Dynamic World: An Application in Aosta Valley (NW Italy). Appl. Sci. 2022, 13, 390. [Google Scholar] [CrossRef]
  9. Orusa, T.; Mondino, E.B. Landsat 8 Thermal Data to Support Urban Management and Planning in the Climate Change Era: A Case Study in Torino Area, NW Italy. In Proceedings of the Remote Sensing Technologies and Applications in Urban Environments IV, Strasbourg, France, 9–10 September 2019; International Society for Optics and Photonics. Volume 11157, pp. 133–149. [Google Scholar]
  10. Carella, E.; Orusa, T.; Viani, A.; Meloni, D.; Borgogno-Mondino, E.; Orusa, R. An Integrated, Tentative Remote-Sensing Approach Based on NDVI Entropy to Model Canine Distemper Virus in Wildlife and to Prompt Science-Based Management Policies. Animals 2022, 12, 1049. [Google Scholar] [CrossRef]
  11. Zhao, Q.; Yu, L.; Li, X.; Peng, D.; Zhang, Y.; Gong, P. Progress and Trends in the Application of Google Earth and Google Earth Engine. Remote Sens. 2021, 13, 3778. [Google Scholar] [CrossRef]
  12. Kennedy, R.E.; Yang, Z.; Gorelick, N.; Braaten, J.; Cavalcante, L.; Cohen, W.B.; Healey, S. Implementation of the LandTrendr Algorithm on Google Earth Engine. Remote Sens. 2018, 10, 691. [Google Scholar] [CrossRef] [Green Version]
  13. Tamiminia, H.; Salehi, B.; Mahdianpari, M.; Quackenbush, L.; Adeli, S.; Brisco, B. Google Earth Engine for Geo-Big Data Applications: A Meta-Analysis and Systematic Review. ISPRS J. Photogramm. Remote Sens. 2020, 164, 152–170. [Google Scholar] [CrossRef]
  14. Yang, L.; Driscol, J.; Sarigai, S.; Wu, Q.; Chen, H.; Lippitt, C.D. Google Earth Engine and Artificial Intelligence (Ai): A Comprehensive Review. Remote Sens. 2022, 14, 3253. [Google Scholar] [CrossRef]
  15. Orusa, T.; Borgogno Mondino, E. Exploring Short-Term Climate Change Effects on Rangelands and Broad-Leaved Forests by Free Satellite Data in Aosta Valley (Northwest Italy). Climate 2021, 9, 47. [Google Scholar] [CrossRef]
  16. Orusa, T.; Orusa, R.; Viani, A.; Carella, E.; Borgogno Mondino, E. Geomatics and EO Data to Support Wildlife Diseases Assessment at Landscape Level: A Pilot Experience to Map Infectious Keratoconjunctivitis in Chamois and Phenological Trends in Aosta Valley (NW Italy). Remote Sens. 2020, 12, 3542. [Google Scholar] [CrossRef]
  17. Peñuelas, J.; Filella, I. Phenology Feedbacks on Climate Change. Science 2009, 324, 887–888. [Google Scholar] [CrossRef] [Green Version]
  18. Piao, S.; Liu, Q.; Chen, A.; Janssens, I.A.; Fu, Y.; Dai, J.; Liu, L.; Lian, X.; Shen, M.; Zhu, X. Plant Phenology and Global Climate Change: Current Progresses and Challenges. Glob. Change Biol. 2019, 25, 1922–1940. [Google Scholar] [CrossRef]
  19. 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]
  20. Yuan, M.; Wang, L.; Lin, A.; Liu, Z.; Qu, S. Variations in Land Surface Phenology and Their Response to Climate Change in Yangtze River Basin during 1982–2015. Theor. Appl. Climatol. 2019, 137, 1659–1674. [Google Scholar] [CrossRef]
  21. Pourghasemi, H.R.; Pouyan, S.; Heidari, B.; Farajzadeh, Z.; Shamsi, S.R.F.; Babaei, S.; Khosravi, R.; Etemadi, M.; Ghanbarian, G.; Farhadi, A.; et al. Spatial Modeling, Risk Mapping, Change Detection, and Outbreak Trend Analysis of Coronavirus (COVID-19) in Iran (Days between February 19 and June 14, 2020). Int. J. Infect. Dis. 2020, 98, 90–108. [Google Scholar] [CrossRef]
  22. Borgogno-Mondino, E.; Farbo, A.; Novello, V.; de Palma, L. A Fast Regression-Based Approach to Map Water Status of Pomegranate Orchards with Sentinel 2 Data. Horticulturae 2022, 8, 759. [Google Scholar] [CrossRef]
  23. Yang, Y.; Qi, N.; Zhao, J.; Meng, N.; Lu, Z.; Wang, X.; Kang, L.; Wang, B.; Li, R.; Ma, J.; et al. Detecting the Turning Points of Grassland Autumn Phenology on the Qinghai-Tibetan Plateau: Spatial Heterogeneity and Controls. Remote Sens. 2021, 13, 4797. [Google Scholar] [CrossRef]
  24. Chen, F.; Liu, Z.; Zhong, H.; Wang, S. Exploring the Applicability and Scaling Effects of Satellite-Observed Spring and Autumn Phenology in Complex Terrain Regions Using Four Different Spatial Resolution Products. Remote Sens. 2021, 13, 4582. [Google Scholar] [CrossRef]
  25. Vizzari, M. PlanetScope, Sentinel-2, and Sentinel-1 Data Integration for Object-Based Land Cover Classification in Google Earth Engine. Remote Sens. 2022, 14, 2628. [Google Scholar] [CrossRef]
  26. De Marinis, P.; De Petris, S.; Sarvia, F.; Manfron, G.; Momo, E.J.; Orusa, T.; Corvino, G.; Sali, G.; Borgogno, E.M. Supporting Pro-Poor Reforms of Agricultural Systems in Eastern DRC (Africa) with Remotely Sensed Data: A Possible Contribution of Spatial Entropy to Interpret Land Management Practices. Land 2021, 10, 1368. [Google Scholar] [CrossRef]
  27. da Silva Junior, C.A.; Leonel-Junior, A.H.S.; Rossi, F.S.; Correia Filho, W.L.F.; de Barros Santiago, D.; de Oliveira-Júnior, J.F.; Teodoro, P.E.; Lima, M.; Capristo-Silva, G.F. Mapping Soybean Planting Area in Midwest Brazil with Remotely Sensed Images and Phenology-Based Algorithm Using the Google Earth Engine Platform. Comput. Electron. Agric. 2020, 169, 105194. [Google Scholar] [CrossRef]
  28. dela Torre, D.M.G.; Gao, J.; Macinnis-Ng, C.; Shi, Y. Phenology-Based Delineation of Irrigated and Rain-Fed Paddy Fields with Sentinel-2 Imagery in Google Earth Engine. Geo-Spat. Inf. Sci. 2021, 24, 695–710. [Google Scholar] [CrossRef]
  29. Kumar, M.; Phukon, S.N.; Paygude, A.C.; Tyagi, K.; Singh, H. Mapping Phenological Functional Types (PhFT) in the Indian Eastern Himalayas Using Machine Learning Algorithm in Google Earth Engine. Comput. Geosci. 2022, 158, 104982. [Google Scholar] [CrossRef]
  30. Ni, R.; Tian, J.; Li, X.; Yin, D.; Li, J.; Gong, H.; Zhang, J.; Zhu, L.; Wu, D. An Enhanced Pixel-Based Phenological Feature for Accurate Paddy Rice Mapping with Sentinel-2 Imagery in Google Earth Engine. ISPRS J. Photogramm. Remote Sens. 2021, 178, 282–296. [Google Scholar] [CrossRef]
  31. Guo, Y.; Xia, H.; Pan, L.; Zhao, X.; Li, R.; Bian, X.; Wang, R.; Yu, C. Development of a New Phenology Algorithm for Fine Mapping of Cropping Intensity in Complex Planting Areas Using Sentinel-2 and Google Earth Engine. ISPRS Int. J. Geo-Inf. 2021, 10, 587. [Google Scholar] [CrossRef]
  32. Pan, L.; Xia, H.; Zhao, X.; Guo, Y.; Qin, Y. Mapping Winter Crops Using a Phenology Algorithm, Time-Series Sentinel-2 and Landsat-7/8 Images, and Google Earth Engine. Remote Sens. 2021, 13, 2510. [Google Scholar] [CrossRef]
  33. Pan, L.; Xia, H.; Yang, J.; Niu, W.; Wang, R.; Song, H.; Guo, Y.; Qin, Y. Mapping Cropping Intensity in Huaihe Basin Using Phenology Algorithm, All Sentinel-2 and Landsat Images in Google Earth Engine. Int. J. Appl. Earth Obs. Geoinf. 2021, 102, 102376. [Google Scholar] [CrossRef]
  34. Kibret, K.S.; Marohn, C.; Cadisch, G. Use of MODIS EVI to Map Crop Phenology, Identify Cropping Systems, Detect Land Use Change and Drought Risk in Ethiopia–an Application of Google Earth Engine. Eur. J. Remote Sens. 2020, 53, 176–191. [Google Scholar] [CrossRef]
  35. Salinero-Delgado, M.; Estévez, J.; Pipia, L.; Belda, S.; Berger, K.; Paredes Gómez, V.; Verrelst, J. Monitoring Cropland Phenology on Google Earth Engine Using Gaussian Process Regression. Remote Sens. 2021, 14, 146. [Google Scholar] [CrossRef] [PubMed]
  36. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting Plant Phenology in Response to Global Change. Trends Ecol. Evol. 2007, 22, 357–365. [Google Scholar] [CrossRef]
  37. Carter, S.K.; Saenz, D.; Rudolf, V.H. Shifts in Phenological Distributions Reshape Interaction Potential in Natural Communities. Ecol. Lett. 2018, 21, 1143–1151. [Google Scholar] [CrossRef]
  38. Sarvia, F.; De Petris, S.; Borgogno-Mondino, E. Exploring Climate Change Effects on Vegetation Phenology by MOD13Q1 Data: The Piemonte Region Case Study in the Period 2001–2019. Agronomy 2021, 11, 555. [Google Scholar] [CrossRef]
  39. Brown, M.E.; de Beurs, K.; Vrieling, A. The Response of African Land Surface Phenology to Large Scale Climate Oscillations. Remote Sens. Environ. 2010, 114, 2286–2296. [Google Scholar] [CrossRef] [Green Version]
  40. Beck, P.S.; Atzberger, C.; Høgda, K.A.; Johansen, B.; Skidmore, A.K. Improved Monitoring of Vegetation Dynamics at Very High Latitudes: A New Method Using MODIS NDVI. Remote Sens. Environ. 2006, 100, 321–334. [Google Scholar] [CrossRef]
  41. Elmore, A.J.; Guinn, S.M.; Minsley, B.J.; Richardson, A.D. Landscape Controls on the Timing of Spring, Autumn, and Growing Season Length in Mid-A Tlantic Forests. Glob. Change Biol. 2012, 18, 656–674. [Google Scholar] [CrossRef] [Green Version]
  42. White, M.A.; Thornton, P.E.; Running, S.W. A Continental Phenology Model for Monitoring Vegetation Responses to Interannual Climatic Variability. Glob. Biogeochem. Cycles 1997, 11, 217–234. [Google Scholar] [CrossRef]
  43. Hufkens, K.; Basler, D.; Milliman, T.; Melaas, E.K.; Richardson, A.D. An Integrated Phenology Modelling Framework in R. Methods Ecol. Evol. 2018, 9, 1276–1285. [Google Scholar] [CrossRef] [Green Version]
  44. 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]
  45. Kong, D. R Package: A State-of-the-Art Vegetation Phenology Extraction Package, Phenofit Version 0.2 6; R Package: Vienna, Austria, 2020. [Google Scholar]
  46. Kong, D.; Zhang, Y.; Wang, D.; Chen, J.; Gu, X. Photoperiod Explains the Asynchronization between Vegetation Carbon Phenology and Vegetation Greenness Phenology. J. Geophys. Res. Biogeosci. 2020, 125, e2020JG005636. [Google Scholar] [CrossRef]
  47. Kong, D.; Zhang, Y.; Gu, X.; Wang, D. A Robust Method for Reconstructing Global MODIS EVI Time Series on the Google Earth Engine. ISPRS J. Photogramm. Remote Sens. 2019, 155, 13–24. [Google Scholar] [CrossRef]
  48. Filippa, G.; Cremonese, E.; Migliavacca, M.; Galvagno, M.; Forkel, M.; Wingate, L.; Tomelleri, E.; Di Cella, U.M.; Richardson, A.D. Phenopix: AR Package for Image-Based Vegetation Phenology. Agric. For. Meteorol. 2016, 220, 141–150. [Google Scholar] [CrossRef] [Green Version]
  49. Eklundh, L.; Jönsson, P. TIMESAT: A Software Package for Time-Series Processing and Assessment of Vegetation Dynamics. In Remote Sensing Time Series; Springer: Stockholm, Sweden, 2015; pp. 141–158. [Google Scholar]
  50. Jönsson, P.; Eklundh, L. TIMESAT—A Program for Analyzing Time-Series of Satellite Sensor Data. Comput. Geosci. 2004, 30, 833–845. [Google Scholar] [CrossRef] [Green Version]
  51. Araya, S.; Ostendorf, B.; Lyle, G.; Lewis, M. CropPhenology: An R Package for Extracting Crop Phenology from Time Series Remotely Sensed Vegetation Index Imagery. Ecol. Inform. 2018, 46, 45–56. [Google Scholar] [CrossRef]
  52. 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]
  53. Armstrong, A.; Ostle, N.J.; Whitaker, J. Solar Park Microclimate and Vegetation Management Effects on Grassland Carbon Cycling. Environ. Res. Lett. 2016, 11, 074016. [Google Scholar] [CrossRef] [Green Version]
  54. Kohler, T.; Giger, M.; Hurni, H.; Ott, C.; Wiesmann, U.; von Dach, S.W.; Maselli, D. Mountains and Climate Change: A Global Concern. Mt. Res. Dev. 2010, 30, 53–55. [Google Scholar] [CrossRef] [Green Version]
  55. Adler, C.; Huggel, C.; Orlove, B.; Nolin, A. Climate Change in the Mountain Cryosphere: Impacts and Responses. Reg. Environ. Change 2019, 19, 1225–1228. [Google Scholar] [CrossRef] [Green Version]
  56. Sarvia, F.; Petris, S.D.; Orusa, T.; Borgogno-Mondino, E. MAIA S2 Versus Sentinel 2: Spectral Issues and Their Effects in the Precision Farming Context. In Proceedings of the International Conference on Computational Science and Its Applications, Cagliari, Italy, 13–16 September 2021; pp. 63–77. [Google Scholar]
  57. Samuele, D.P.; Filippo, S.; Orusa, T.; Enrico, B.-M. Mapping SAR Geometric Distortions and Their Stability along Time: A New Tool in Google Earth Engine Based on Sentinel-1 Image Time Series. Int. J. Remote Sens. 2021, 42, 9135–9154. [Google Scholar] [CrossRef]
  58. Meybeck, M.; Green, P.; Vörösmarty, C. A New Typology for Mountains and Other Relief Classes. Mt. Res. Dev. 2001, 21, 34–45. [Google Scholar] [CrossRef] [Green Version]
  59. Clarke-Sather, A.; Crow-Miller, B.; Banister, J.M.; Anh Thomas, K.; Norman, E.S.; Stephenson, S.R. The Shifting Geopolitics of Water in the Anthropocene. Geopolitics 2017, 22, 332–359. [Google Scholar] [CrossRef]
  60. Mittermeier, R.A.; Mittermeier, C.G.; Brooks, T.M.; Pilgrim, J.D.; Konstant, W.R.; Da Fonseca, G.A.; Kormos, C. Wilderness and Biodiversity Conservation. Proc. Natl. Acad. Sci. USA 2003, 100, 10309–10313. [Google Scholar] [CrossRef] [Green Version]
  61. Schlagintweit, H.; Schlagintweit, A. On the Physical Geography of the Alps. 353. Am. J. Sci. 1852, 64, 359. [Google Scholar]
  62. Deering, D. Measuring” Forage Production” of Grazing Units from Landsat MSS Data. In Proceedings of the Tenth International Symposium of Remote Sensing of the Envrionment, Ann Arbor, MI, USA, 6–10 October 1975; pp. 1169–1198. [Google Scholar]
  63. Deering, D.W. Rangeland Reflectance Characteristics Measured by Aircraft and Spacecraftsensors; Texas A&M University: College Station, TX, USA, 1978. [Google Scholar]
  64. Rouse, J.; Haas, R.; Schell, J.; Deering, D. Monitoring vegetation systems in the greant plains with ERTS. In Proceedings of the Third ERTS (Earth Resources Technology Satellite) Symposium, NASA SP-351, Washington, DC, USA, 10–14 December 1973; Volume 1, pp. 309–317. [Google Scholar]
  65. Rouse, J., Jr.; Haas, R.H.; Deering, D.; Schell, J.; Harlan, J.C. Monitoring the Vernal Advancement and Retrogradation (Green Wave Effect) of Natural Vegetation; Springer: New York, NY, USA, 1974. [Google Scholar]
  66. Tucker, C.J. Red and Photographic Infrared Linear Combinations for Monitoring Vegetation. Remote Sens. Environ. 1979, 8, 127–150. [Google Scholar] [CrossRef] [Green Version]
  67. Wang, C.; Li, J.; Liu, Q.; Zhong, B.; Wu, S.; Xia, C. Analysis of Differences in Phenology Extracted from the Enhanced Vegetation Index and the Leaf Area Index. Sensors 2017, 17, 1982. [Google Scholar] [CrossRef] [Green Version]
  68. Wu, C.; Chen, J.M.; Huang, N. Predicting Gross Primary Production from the Enhanced Vegetation Index and Photosynthetically Active Radiation: Evaluation and Calibration. Remote Sens. Environ. 2011, 115, 3424–3435. [Google Scholar] [CrossRef]
  69. Liu, Z.; Liu, K.; Zhang, J.; Yan, C.; Lock, T.R.; Kallenbach, R.L.; Yuan, Z. Fractional Coverage Rather than Green Chromatic Coordinate Is a Robust Indicator to Track Grassland Phenology Using Smartphone Photography. Ecol. Inform. 2022, 68, 101544. [Google Scholar] [CrossRef]
  70. Xu, D.; Wang, C.; Chen, J.; Shen, M.; Shen, B.; Yan, R.; Li, Z.; Karnieli, A.; Chen, J.; Yan, Y.; et al. The Superiority of the Normalized Difference Phenology Index (NDPI) for Estimating Grassland Aboveground Fresh Biomass. Remote Sens. Environ. 2021, 264, 112578. [Google Scholar] [CrossRef]
  71. 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]
  72. 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]
  73. 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]
  74. 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]
  75. Kong, D.; McVicar, T.R.; Xiao, M.; Zhang, Y.; Peña-Arancibia, J.L.; Filippa, G.; Xie, Y.; Gu, X. Phenofit: An R Package for Extracting Vegetation Phenology from Time Series Remote Sensing. Methods Ecol. Evol. 2022, 13, 1508–1527. [Google Scholar] [CrossRef]
  76. Stenseth, N.C.; Mysterud, A. Climate, Changing Phenology, and Other Life History Traits: Nonlinearity and Match–Mismatch to the Environment. Proc. Natl. Acad. Sci. USA 2002, 99, 13379–13381. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Puchalka, R.; Klisz, M.; Koniakin, S.; Czortek, P.; Dylewski, L.; Paź-Dyderska, S.; Vítková, M.; Sádlo, J.; Rašomavičius, V.; Čarni, A.; et al. Citizen Science Helps Predictions of Climate Change Impact on Flowering Phenology: A Study on Anemone Nemorosa. Agric. For. Meteorol. 2022, 325, 109133. [Google Scholar] [CrossRef]
  78. Sass-Klaassen, U.; Sabajo, C.R.; den Ouden, J. Vessel Formation in Relation to Leaf Phenology in Pedunculate Oak and European Ash. Dendrochronologia 2011, 29, 171–175. [Google Scholar] [CrossRef]
  79. Puchalka, R.; Koprowski, M.; Gričar, J.; Przybylak, R. Does Tree-Ring Formation Follow Leaf Phenology in Pedunculate Oak (Quercus Robur L.)? Eur. J. For. Res. 2017, 136, 259–268. [Google Scholar] [CrossRef] [Green Version]
  80. Matsushita, B.; Yang, W.; Chen, J.; Onda, Y.; Qiu, G. Sensitivity of the Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) to Topographic Effects: A Case Study in High-Density Cypress Forest. Sensors 2007, 7, 2636–2651. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  81. Kumari, N.; Srivastava, A.; Dumka, U.C. A Long-Term Spatiotemporal Analysis of Vegetation Greenness over the Himalayan Region Using Google Earth Engine. Climate 2021, 9, 109. [Google Scholar] [CrossRef]
  82. Gao, X.; Huete, A.R.; Ni, W.; Miura, T. Optical–Biophysical Relationships of Vegetation Spectra without Background Contamination. Remote Sens. Environ. 2000, 74, 609–620. [Google Scholar] [CrossRef]
  83. 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]
  84. Yan, E.; Wang, G.; Lin, H.; Xia, C.; Sun, H. Phenology-Based Classification of Vegetation Cover Types in Northeast China Using MODIS NDVI and EVI Time Series. Int. J. Remote Sens. 2015, 36, 489–512. [Google Scholar] [CrossRef]
  85. Dong, T.; Liu, J.; Qian, B.; Zhao, T.; Jing, Q.; Geng, X.; Wang, J.; Huffman, T.; Shang, J. Estimating Winter Wheat Biomass by Assimilating Leaf Area Index Derived from Fusion of Landsat-8 and MODIS Data. Int. J. Appl. Earth Obs. Geoinf. 2016, 49, 63–74. [Google Scholar] [CrossRef]
  86. Dong, J.; Xiao, X.; Menarguez, M.A.; Zhang, G.; Qin, Y.; Thau, D.; Biradar, C.; Moore III, B. Mapping Paddy Rice Planting Area in Northeastern Asia with Landsat 8 Images, Phenology-Based Algorithm and Google Earth Engine. Remote Sens. Environ. 2016, 185, 142–154. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Bagliani, M.M.; Caimotto, M.C.; Latini, G.; Orusa, T. Lessico e Nuvole: Le Parole del Cambiamento Climatico; Università degli studi di Torino: Torino, Italy, 2019. [Google Scholar]
  88. Latini, G.; Bagliani, M.; Orusa, T. Lessico e Nuvole: Le Parole del Cambiamento Climatico; Università degli studi di Torino Youcanprint: Torino, Italy, 2021. [Google Scholar]
Figure 1. Alpine chain considered as development area for the GEE PMs algorithm. EPSG: 4326.
Figure 1. Alpine chain considered as development area for the GEE PMs algorithm. EPSG: 4326.
Geomatics 03 00012 g001
Figure 2. PMs computed from NDVI at 10 m GSD with the GEE algorithm developed in 2021 using Sentinel-2 for the entire Alpine chain. EPSG: 4326.
Figure 2. PMs computed from NDVI at 10 m GSD with the GEE algorithm developed in 2021 using Sentinel-2 for the entire Alpine chain. EPSG: 4326.
Geomatics 03 00012 g002
Figure 3. Validation points randomly selected over the AOI in vegetated areas with NDVI threshold greater than 0.3. EPSG: 4326.
Figure 3. Validation points randomly selected over the AOI in vegetated areas with NDVI threshold greater than 0.3. EPSG: 4326.
Geomatics 03 00012 g003
Figure 4. Boxplot considering all PMs computed in different software from different satellite sensors.
Figure 4. Boxplot considering all PMs computed in different software from different satellite sensors.
Geomatics 03 00012 g004
Figure 5. Comparison of Start of Season (SOS) and End of Season (EOS) dates computed, respectively, from Sentinel-2 and Landsat retrieved in R with the greenbrown package and from the GEE algorithm.
Figure 5. Comparison of Start of Season (SOS) and End of Season (EOS) dates computed, respectively, from Sentinel-2 and Landsat retrieved in R with the greenbrown package and from the GEE algorithm.
Geomatics 03 00012 g005
Figure 6. PM comparison, Start of Season (SOS), End of Season (EOS), and Length of Season (LOS) estimated with Sentinel-2 in GEE and MODIS phenology products (ready to use, provided by NASA).
Figure 6. PM comparison, Start of Season (SOS), End of Season (EOS), and Length of Season (LOS) estimated with Sentinel-2 in GEE and MODIS phenology products (ready to use, provided by NASA).
Geomatics 03 00012 g006
Table 1. Validation.
Table 1. Validation.
EO Data
and PMs
GEE Algorithm (pi)
μ ± σ (DOY)
MODIS Terra and Aqua (oi)
μ ± σ (DOY)
MAERMSE
Sentinel-2 SOS123 ± 28125 ± 2822
Sentinel-2 EOS 295 ± 34296 ± 3411
Sentinel-2 POS209 ± 18210 ± 1811
Sentinel-2 LOS172 ± 52174 ± 5222
Landsat 8 SOS124 ± 28 125 ± 2822
Landsat 8 EOS296 ± 34296 ± 3411
Landsat 8 POS210 ± 18210 ± 1811
Landsat 8 LOS172 ± 52174 ± 5222
Table 2. Phenological Metrics comparison.
Table 2. Phenological Metrics comparison.
EO Data
and PMs
GEE Algorithm (pi)
μ ± σ (DOY)
R Greenbrown Package (oi)
μ ± σ (DOY)
MAERMSE
Sentinel-2 SOS123 ± 28124 ± 2811
Sentinel-2 EOS295 ± 34296 ± 3411
Sentinel-2 POS209 ± 18210 ± 1811
Sentinel-2 LOS172 ± 52172 ± 5200
Landsat 8 SOS124 ± 28 124 ± 2800
Landsat 8 EOS296 ± 34298 ± 3422
Landsat 8 POS210 ± 18208 ± 1822
Landsat 8 LOS172 ± 52171 ± 5211
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Orusa, T.; Viani, A.; Cammareri, D.; Borgogno Mondino, E. A Google Earth Engine Algorithm to Map Phenological Metrics in Mountain Areas Worldwide with Landsat Collection and Sentinel-2. Geomatics 2023, 3, 221-238. https://doi.org/10.3390/geomatics3010012

AMA Style

Orusa T, Viani A, Cammareri D, Borgogno Mondino E. A Google Earth Engine Algorithm to Map Phenological Metrics in Mountain Areas Worldwide with Landsat Collection and Sentinel-2. Geomatics. 2023; 3(1):221-238. https://doi.org/10.3390/geomatics3010012

Chicago/Turabian Style

Orusa, Tommaso, Annalisa Viani, Duke Cammareri, and Enrico Borgogno Mondino. 2023. "A Google Earth Engine Algorithm to Map Phenological Metrics in Mountain Areas Worldwide with Landsat Collection and Sentinel-2" Geomatics 3, no. 1: 221-238. https://doi.org/10.3390/geomatics3010012

APA Style

Orusa, T., Viani, A., Cammareri, D., & Borgogno Mondino, E. (2023). A Google Earth Engine Algorithm to Map Phenological Metrics in Mountain Areas Worldwide with Landsat Collection and Sentinel-2. Geomatics, 3(1), 221-238. https://doi.org/10.3390/geomatics3010012

Article Metrics

Back to TopTop