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

: 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 efﬁcient 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 ﬁll 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.


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 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 higherperformance 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 cloudbased 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].

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].

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_HARMO-NIZED"), 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

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_HARMO-NIZED"), 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).

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: Green chromatic coordinate (GCC) [69] 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: where τ is the dynamic value that depends on the annual amplitude of the time series; VI min and VI max 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: GCC min = 0.31, NDVI min = 0.39, EVI min = 0.2, and NDPI min = 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 NDPI min 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.

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: where p i 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).

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.

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). Geomatics 2023, 3, FOR PEER REVIEW 10 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).  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).
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). 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).
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 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).

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 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).

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 goodquality 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 gapfilling 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 remotesensing 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.

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.