Assessing Forest Phenology: A Multi-Scale Comparison of Near-Surface (UAV, Spectral Reﬂectance Sensor, PhenoCam) and Satellite (MODIS, Sentinel-2) Remote Sensing

: The monitoring of forest phenology based on observations from near-surface sensors such as Unmanned Aerial Vehicles (UAVs), PhenoCams, and Spectral Reﬂectance Sensors (SRS) over satellite signiﬁcant in the ﬁeld of remote sensing and vegetation exploring different aspects of forest phenology based on observations from these sensors and drawing comparatives from the time series of vegetation indices (VIs) a challenge. Accordingly, this research explores the potential of near-surface sensors to track the temporal dynamics of phenology, cross-compare their results against satellite observations (MODIS, Sentinel-2), and validate satellite-derived phenology. A time series of Normalized Difference Vegetation Index (NDVI), Green Chromatic Coordinate (GCC), and Normalized Difference of Green & Red (VIgreen) indices were extracted from both near-surface and satellite sensor platforms. The regression analysis between time series of NDVI data from different sensors shows the high Pearson’s correlation coefﬁcients (r > 0.75). Despite the good correlations, there was a remarkable offset and signiﬁcant differences in slope during green-up and senescence periods. SRS showed the most distinctive NDVI proﬁle and was different to other sensors. PhenoCam GCC tracked green-up of the canopy better than the other indices, with a well-deﬁned start, end, and peak of the season, and was most closely correlated (r > 0.93) with the satellites, while SRS-based VIgreen accounted for the least correlation (r = 0.58) against Sentinel-2. Phenophase transition dates were estimated and validated against visual inspection of the PhenoCam data. The Start of Spring (SOS) and End of Spring (EOS) could be predicted with an accuracy of <3 days with GCC, while these metrics from VIgreen and NDVI resulted in a slightly higher bias of (3–10) days. The observed agreement between UAV NDVI vs. satellite NDVI and PhenoCam GCC vs. satellite GCC suggests that it is feasible to use PhenoCams and UAVs for satellite data validation and upscaling. Thus, a combination of these near-surface vegetation metrics is promising for a holistic understanding of vegetation phenology from canopy perspective and could serve as a good foundation for analysing the interoperability of different sensors for vegetation dynamics and change analysis.


Introduction
Vegetation phenology describes the life cycle events that occur throughout the year and has been identified as a way of studying ecosystem processes [1]. In particular, the monitoring of phenology dynamics over the years at various ecological scales helps in understanding how the plants are responding to climate change [2][3][4].
The study of phenology has a long history [5,6] and is founded in field observations of phenological events (leaf shooting, flowering, leaf fall, etc.) that respond to internal molecular mechanisms of plants driven by photoperiod, temperature and other factors [7]. Hence, systematic observations of these events provide important evidence of changes

PhenoCam
The study used data from one facing south, mounted on a tower, at a height of 7.25 m above ground level, pointing 45° down from the horizontal, which captures digital imagery of the foreground canopy inside a field-of-view (FOV) of 60°. The PhenoCam images were acquired by a Mobotix 5 MP RGB camera. The camera captures a coloured three-layer R, G and B images at hourly temporal resolution from 5 a.m. to 8 p.m. all year long. The images for the year 2018 were downloaded from the SITES project server [40]. For this study, we used only images between 10 a.m. and 2 p.m., to avoid low solar elevation.

Unmanned Aerial Vehicles (UAVs)
For this study, we used data captured by a Parrot Sequoia multispectral camera (Parrot Drone SAS, Paris, France) [41] mounted on a Solo 3DR UAV (3DR Robotics, California, USA) [42]. The UAV system (UAS) payload capability comprises a Global Navigation Satellite System (GNSS), an Inertial Motion System (IMU), an accelerometer, a compass, and the camera system. The sensor onboard UAV captured data from an area covering around 10 hectares. Figure 2 illustrates the UAV platform and Parrot Sequoia sensor used in the study.

PhenoCam
The study used data from one facing south, mounted on a tower, at a height of 7.25 m above ground level, pointing 45 • down from the horizontal, which captures digital imagery of the foreground canopy inside a field-of-view (FOV) of 60 • . The PhenoCam images were acquired by a Mobotix 5 MP RGB camera. The camera captures a coloured three-layer R, G and B images at hourly temporal resolution from 5 a.m. to 8 p.m. all year long. The images for the year 2018 were downloaded from the SITES project server [40]. For this study, we used only images between 10 a.m. and 2 p.m., to avoid low solar elevation.

Unmanned Aerial Vehicles (UAVs)
For this study, we used data captured by a Parrot Sequoia multispectral camera (Parrot Drone SAS, Paris, France) [41] mounted on a Solo 3DR UAV (3DR Robotics, California, USA) [42]. The UAV system (UAS) payload capability comprises a Global Navigation Satellite System (GNSS), an Inertial Motion System (IMU), an accelerometer, a compass, and the camera system. The sensor onboard UAV captured data from an area covering around 10 hectares. Figure 2 illustrates the UAV platform and Parrot Sequoia sensor used in the study. 3DR quadcopter [42] (right).
The Parrot Sequoia camera is designed to acquire images in four spectral bands: green (G), red (R), Red Edge (RE) and near-infrared (NIR), with wavelengths specified in Table 1. The multispectral camera was installed under the drone facing the forest canopy. The images were captured by 1.2 MP monochrome sensors at a bit depth of 16 bit in raw format and later saved as .tif files. able 1. Spectral band information for the near-surface multispectral sensors (Parrot Sequoia onboard UAV and Spectral eflectance Sensors), satellite spectrometers (MODIS and Sentinel-2) and RGB PhenoCam. Only the spectral bands which ere used in the research are mentioned here. The UAV flight frequency for the study period was roughly one flight per month from April to August 2018. The flights were not conducted in a fixed interval of time mainly because of weather conditions. The UAV was programmed to fly at a height of 80 m above ground level with frontal and side overlap of 80%. Most of the flights were conducted at midday (between 9 a.m. and 4 p.m.) on either clear days or evenly overcast days. The UAV was pre-programmed to follow the waypoints that would cover up the study area in two flights of approximately 10 min each. The camera and image settings were kept identical for all flights, in auto mode (i.e., camera chooses shutter speed and ISO value, while it has a fixed aperture) for exposure compensation. During the missions, the incoming radiation was measured and recorded by the sunshine sensor mounted on top of the UAV. The purpose of using the sunshine sensor was to normalize the images for variations in incoming sunlight, and for calculation of reflectance.

Spectral Bands (nm) Sensor Blue (B) Green (G) Red (R) Near-infrared (NIR)
Three Spectralon reflectance panels (Labsphere, North Sutton, NH, USA), often defined as ground calibration targets (GCTs), and eight ground control points (GCPs) were placed within the UAV survey's coverage on all flight dates for the purpose of radiometric calibration and georeferencing, respectively. Detailed information on all of the UAV flights carried out during the study period is listed in Table 2. The Parrot Sequoia camera is designed to acquire images in four spectral bands: green (G), red (R), Red Edge (RE) and near-infrared (NIR), with wavelengths specified in Table 1. The multispectral camera was installed under the drone facing the forest canopy. The images were captured by 1.2 MP monochrome sensors at a bit depth of 16 bit in raw format and later saved as .tif files. The UAV flight frequency for the study period was roughly one flight per month from April to August 2018. The flights were not conducted in a fixed interval of time mainly because of weather conditions. The UAV was programmed to fly at a height of 80 m above ground level with frontal and side overlap of 80%. Most of the flights were conducted at midday (between 9 a.m. and 4 p.m.) on either clear days or evenly overcast days. The UAV was pre-programmed to follow the waypoints that would cover up the study area in two flights of approximately 10 min each. The camera and image settings were kept identical for all flights, in auto mode (i.e., camera chooses shutter speed and ISO value, while it has a fixed aperture) for exposure compensation. During the missions, the incoming radiation was measured and recorded by the sunshine sensor mounted on top of the UAV. The purpose of using the sunshine sensor was to normalize the images for variations in incoming sunlight, and for calculation of reflectance.
Three Spectralon reflectance panels (Labsphere, North Sutton, NH, USA), often defined as ground calibration targets (GCTs), and eight ground control points (GCPs) were placed within the UAV survey's coverage on all flight dates for the purpose of radiometric calibration and georeferencing, respectively. Detailed information on all of the UAV flights carried out during the study period is listed in Table 2. Table 2. Information on the UAV missions conducted during 2018 with corresponding number of images and weather conditions. Note '*'represents one flight for which images were not processed due to large variation in irradiance data during the flight.

No.
Flight

Spectral Reflectance Sensors (SRS)
The SRS used for this study were two two-band Decagon radiometers designed to measure the incident and reflected radiation in green, red and near-infrared (NIR) wavelengths suitable for the computation of NDVI (SRS-NDVI, hereafter) and photochemical reflectance index (PRI) (SRS-PRI, hereafter). The SRS-NDVI system is comprised of a pair of sensors: the SRS hemispherical, and the SRS field-stop lens sensor. The first type has a hemispherical 180 • FOV, and is installed looking up to measure incident radiation, whereas the second type has a FOV confined to 36 • for measuring downward reflected radiation from the forest canopy.
The SRS are mounted on the same tower as the PhenoCam, at the same height of 7.25 m, above the canopy of the spruce forest. Upward and downward-looking sensors measure the incident and reflected canopy radiance, respectively, simultaneously every 10 s, at a wavelength band specified in Table 1, with data aggregated to 10-min intervals. The data have been obtained from the SITES data portal [43].

Earth Observation Satellite Data
A time series of MODIS NDVI global product MOD13Q1 V6 [44] was obtained over the study area for the year 2018. This product, available as an image collection in Google Earth Engine (GEE) [45], provided NDVI values on a per-pixel basis generated every 16 days at 250 m spatial resolution. In addition, for greenness change estimates from visible bands, the MOD09A1 V6 [46] product, also available in GEE [47], was used, which provides surface spectral reflectance of MODIS Terra bands 1-7 at 500 m resolution corrected for atmospheric conditions.
In addition, from the GEE archive, a time series of surface reflectance in red and near-Infrared bands were derived from Sentinel-2A (bands 4 and 8, respectively) [48], which were used to compute NDVI at 10 metres spatial resolution. For the MODIS product, VI values of a single pixel were extracted. The comparison of UAV and Sentinel-2A vegetation indices with MODIS were based on a spatial average of all pixels that falls within ground projected footprint of the MODIS pixel. For the case of Sentinel-2A comparison against SRS, it was equivalent to the spatial average of all pixels within the ground projected footprint of the SRS.

Methodology
Our workflow comprised four processing steps: (1) data collection, (2) data extraction/processing, (3) data smoothing, and (4) phenophase transition dates estima-tion. Finally, an inter-comparison of smoothed time series information from all platforms and evaluation of correlations between data types was carried out. The details are described after the methodology chart ( Figure 3). All the stated python script for PhenoCam image processing and Google Earth Engine (GEE) codes for extracting reflectance and different VIs that were written as a part of this research are available in https://github.com/shangharsha2929/ASA.git (accessed on 7 March 2021).

Methodology
Our workflow comprised four processing steps: (1) data collection, (2) data extraction/processing, (3) data smoothing, and (4) phenophase transition dates estimation. Finally, an inter-comparison of smoothed time series information from all platforms and evaluation of correlations between data types was carried out. The details are described after the methodology chart ( Figure 3). All the stated python script for PhenoCam image processing and Google Earth Engine (GEE) codes for extracting reflectance and different VIs that were written as a part of this research are available in https://github.com/shang-harsha2929/ASA.git (accessed on 7 March 2021).

Phenocam Image Processing
First, the PhenoCam images that appeared disoriented, very bright (i.e., affected by solar glare), foggy, blurred, or had stripes covering more than 90% area of the image were removed manually before processing the images. In addition, snow-covered images were also filtered out so that the VI values from those could be computed separately. All the images that passed these quality control filters were used to compute the time series of vegetation indices. As per the experiment conducted by [31], the images captured at dawn

Phenocam Image Processing
First, the PhenoCam images that appeared disoriented, very bright (i.e., affected by solar glare), foggy, blurred, or had stripes covering more than 90% area of the image were removed manually before processing the images. In addition, snow-covered images were also filtered out so that the VI values from those could be computed separately. All the images that passed these quality control filters were used to compute the time series of vegetation indices. As per the experiment conducted by [31], the images captured at dawn and dusk experience low levels of diffuse sunlight, and are inclined to have lower VI values, compared to those recorded at middle day. With this finding in mind, of all valid images, only images between 10 a.m. to 3 p.m. (solar time) were used in the analysis.
Secondly, a rectangular region of interest (ROI) was defined within each digital image representing a sample of the spruce forest, from which the time series of RGB Digital Number (DN) triplets were extracted. The ROI was constant for all images. An effort was made to confirm that the time series was not affected by camera movements, which would result in a different ROI. When camera movements were detected, these images were processed separately by using a revised ROI that fit the same area of spruce forest. Figure 4 depicts the ROI selected for the study. and dusk experience low levels of diffuse sunlight, and are inclined to have lower VI values, compared to those recorded at middle day. With this finding in mind, of all valid images, only images between 10 a.m. to 3 p.m. (solar time) were used in the analysis.
Secondly, a rectangular region of interest (ROI) was defined within each digital image representing a sample of the spruce forest, from which the time series of RGB Digital Number (DN) triplets were extracted. The ROI was constant for all images. An effort was made to confirm that the time series was not affected by camera movements, which would result in a different ROI. When camera movements were detected, these images were processed separately by using a revised ROI that fit the same area of spruce forest. Figure 4 depicts the ROI selected for the study. The DN values distributed across all the pixels in the corresponding ROI were extracted with a python code (involves use of OpenCV and NumPy module) that extracted the mean DN values of RGB bands across the defined ROIs. Then, the extracted mean DN values were fed into Equations (1) and (2) to compute GCC and VIgreen, respectively.
where = mean red DN, = mean green DN and = mean blue DN for the chosen ROI.
With the calculation of these indices, the impact of variations in scene illumination is minimized [2,31,34,49].
For most of the applications, a very high temporal resolution (e.g., sub-hourly resolution, like in our study) seems unnecessary. This is because, in general, the vegetation colour over such short period of time remains relatively constant. With this consideration in mind, we preferred to calculate 3-day averages [49] of all the valid images for the vegetation indices mentioned above. The DN values distributed across all the pixels in the corresponding ROI were extracted with a python code (involves use of OpenCV and NumPy module) that extracted the mean DN values of RGB bands across the defined ROIs. Then, the extracted mean DN values were fed into Equations (1) and (2) to compute GCC and VIgreen, respectively.
where DN R = mean red DN, DN G = mean green DN and DN B = mean blue DN for the chosen ROI.
With the calculation of these indices, the impact of variations in scene illumination is minimized [2,31,34,49].
For most of the applications, a very high temporal resolution (e.g., sub-hourly resolution, like in our study) seems unnecessary. This is because, in general, the vegetation colour over such short period of time remains relatively constant. With this consideration in mind, we preferred to calculate 3-day averages [49] of all the valid images for the vegetation indices mentioned above.

UAV Image Processing
First, the sunshine sensor data from all UAV flights were retrieved from the images' metadata and then plotted, to check the incoming light conditions for each single band (G, R, RE and NIR).
The UAV flights that presented irregular and complex light variability during the flight ( Figure 5c) were discarded. The data that presented constant light, under sunny ( Figure 5a) or overcast sky (Figure 5b), were further processed. The information regarding the flights that were accepted or rejected for further processing is listed in Table 2. Based on this information, UAV data were characterized in three ways: (1) consistent sunny incoming light, (2) consistent cloudy incoming light, with slight variation, and (3) complex incoming light recorded during the mission. The three different cases are depicted as a, b and c in Figure 5, respectively.

UAV Image Processing
First, the sunshine sensor data from all UAV flights were retrieved from the images' metadata and then plotted, to check the incoming light conditions for each single band (G, R, RE and NIR).
The UAV flights that presented irregular and complex light variability during the flight ( Figure 5c) were discarded. The data that presented constant light, under sunny ( Figure 5a) or overcast sky (Figure 5b), were further processed. The information regarding the flights that were accepted or rejected for further processing is listed in Table 2. Based on this information, UAV data were characterized in three ways: (1) consistent sunny incoming light, (2) consistent cloudy incoming light, with slight variation, and (3) complex incoming light recorded during the mission. The three different cases are depicted as a, b and c in Figure 5, respectively. After categorizing the flights, the individual images of each flight data were subjected to exposure calibration, vignetting correction and irradiance normalization. These corrections were carried out in accordance with a recently developed methodology [50]. To adjust for different exposure settings of images, exposure calibration according to Equation where is pseudo-radiance defined in an arbitrary unit common to all Sequoia cameras, the aperture f-number (f = 2.2 for Sequoia), the digital number of a pixel (DN), exposure time (ε), ISO value (γ) and A, B, C are calibration coefficients provided in the image metadata.
Vignetting is defined as a radial fall-off in pixel values that results in darker areas near the edges of images [17]. The method for minimizing the vignetting effect was based on the use of vignetting polynomial [51] derived from the EXIF/XMP metadata of the Sequoia images. The vignetting polynomial estimation as explained in [50] was based on finding a correction factor for each pixel from a large number of images captured over a Lambertian surface and these pixel-wise correction factors fitted by a polynomial.
Only flight data showing the consistent irradiance throughout the mission ( Figure  6a) was processed for the exposure and vignetting correction. Similarly, the flight data that showed variation in irradiance data that could be handled by normalizing for incoming light condition from sunshine sensor data ( Figure 6b) were normalized using a polynomial trend of degree 'n' or spline trend, based on the nature of the data [50]. The best fit polynomial or splines were selected based on visual inspection by manually fitting to the data. The degree of polynomial was chosen in such a way that it fit the irradiance data with coherence at the start and end of each flight. This is because the GCT images were captured only at the beginning and end of the mission, which were later used for radiometric calibration of the derived orthophotos. The detailed process regarding the exposure, vignetting and irradiance normalization are explained in [50]. After categorizing the flights, the individual images of each flight data were subjected to exposure calibration, vignetting correction and irradiance normalization. These corrections were carried out in accordance with a recently developed methodology [50]. To adjust for different exposure settings of images, exposure calibration according to Equation (3) provided by Parrot Sequoia was performed.
where L p is pseudo-radiance defined in an arbitrary unit common to all Sequoia cameras, the aperture f-number (f = 2.2 for Sequoia), the digital number of a pixel (DN), exposure time (ε), ISO value (γ) and A, B, C are calibration coefficients provided in the image metadata.
Vignetting is defined as a radial fall-off in pixel values that results in darker areas near the edges of images [17]. The method for minimizing the vignetting effect was based on the use of vignetting polynomial [51] derived from the EXIF/XMP metadata of the Sequoia images. The vignetting polynomial estimation as explained in [50] was based on finding a correction factor for each pixel from a large number of images captured over a Lambertian surface and these pixel-wise correction factors fitted by a polynomial.
Only flight data showing the consistent irradiance throughout the mission (Figure 6a) was processed for the exposure and vignetting correction. Similarly, the flight data that showed variation in irradiance data that could be handled by normalizing for incoming light condition from sunshine sensor data ( Figure 6b) were normalized using a polynomial trend of degree 'n' or spline trend, based on the nature of the data [50]. The best fit polynomial or splines were selected based on visual inspection by manually fitting to the data. The degree of polynomial was chosen in such a way that it fit the irradiance data with coherence at the start and end of each flight. This is because the GCT images were captured only at the beginning and end of the mission, which were later used for radiometric calibration of the derived orthophotos. The detailed process regarding the exposure, vignetting and irradiance normalization are explained in [50].

Satellite Data Processing
An automated approach was developed to extract the NDVI values for the stu area, from the global MODIS NDVI product MOD13Q1 on the GEE platform (Figure The NDVI from Sentinel-2A was not directly available in GEE; therefore, the correspon ing band values required to compute NDVI (R and NIR reflectance from bands 4 and respectively) were used. In addition, all three RGB band values of MODIS (MOD09 product) and Sentinel-2A, for the common area similar to SRS were extracted and la used to compute GCC and VIgreen by using GEE. These GCC and VIgreen values w downloaded as .csv files, which were further used on to compute the seasonality even A significant number of pixels reached the top of the spectral range (saturation) in the images of each flight, mostly in the G band and least in the NIR band. This problem is also mentioned by [50,52]. The saturated pixels were generally over bright surfaces (i.e., rocks or gravel soil) on the ground. As these saturated pixels influence the orthophotos, they were masked out [50].
All the flight images that had undergone the process of these corrections were thereafter imported into Photoscan 1.4.2 (Agisoft, St. Petersburg, Russia), a photogrammetry software package for image orthorectification (Table 3). Five out of eight GCPs were used to produce a georeferenced orthomosaic with a nominal spatial resolution of 0.07 m. The other three GCP markers were used as check points, which resulted in a root mean squared error (RMSE) of 0.02 m. Only one orthomosaic (dated 2018-07-17) was georeferenced in Photoscan and this was then used as a reference to co-register the rest of the orthomosaics (eight in total) by using same GCP markers.
Finally, the orthomosaics were subjected to radiometric calibration using the empirical line method suggested by [53]. The method relies on mean pixel values of three Spectralon reflectance panels (GCTs) present in the images and their standard reflectance values. The standard reflectance values were black (5%), dark grey (20%) and light grey (50%). The standard reflectance value of these panels, together with the mean DN of the same panel extracted from a set of images, were used to convert the orthomosaics into reflectance, by establishing a linear relationship between the measurements.
The radiometrically corrected orthomosaics were then used to compute NDVI maps, for all flights. NDVI, with a dynamic range of −1 to +1 [16] was calculated according to Equation (4): where ρ N IR = near-infrared reflectance, and ρ RED = visual red reflectance. The values around 0 represent bare soil and values of 1 refers to vigorous vegetation. The footprint of the SRS was projected on top of the NDVI maps to extract the values of NDVI of the UAV flights, for comparison with other sensors. Figure 6 shows the spatial variation of NDVI over the study area. Low-lying areas in the NE part of the region represent wet bare lands, which are characterized by intermediate

Satellite Data Processing
An automated approach was developed to extract the NDVI values for the study area, from the global MODIS NDVI product MOD13Q1 on the GEE platform ( Figure 7). The NDVI from Sentinel-2A was not directly available in GEE; therefore, the corresponding band values required to compute NDVI (R and NIR reflectance from bands 4 and 8, respectively) were used. In addition, all three RGB band values of MODIS (MOD09A1 product) and Sentinel-2A, for the common area similar to SRS were extracted and later used to compute GCC and VIgreen by using GEE. These GCC and VIgreen values were downloaded as .csv files, which were further used on to compute the seasonality events. The cloudy pixels, in the case of Sentinel-2 images, were masked out using available cloud mask band 'QA60' within GEE platform.

Spectral Reflectance Sensor Data Processing
The irradiance data from SRS were calibrated using calibration factors of 0.250 (R) and 0.260 (NIR) in the case of SRS-NDVI sensor, and 0.259 (570 nm) for the SRS-PRI sensor following methodology in [54]. A subset of measurements was extracted from the full time series, which covered exactly between the same hours as the period used for the Pheno-Cam (9 a.m. to 3 p.m.). The calculation of SRS-NDVI was done on reflectance, calculated as a ratio from incoming and reflected radiation, by dividing the downward-looking by the upward-looking sensor data. The data was averaged in similar time aggregates as the PhenoCam, i.e., 3-day averaged NDVI values. In addition to the NDVI, VIgreen was computed in a similar fashion to the PhenoCam, combining the red band from the NDVI-SRS and one of the green bands at 570 nm from the SRS-PRI sensor. The footprint of the SRS was computed using the sensor characteristics (i.e., FOV: 36°, height: 7.25 m, off-nadir angle: 45°) and the tower coordinates. The same footprint was used across all sensors to extract time series of VI values for effective comparison.

Curve Fitting to Time Series Data
The original time series data obtained from different platforms (UAV, PhenoCam, satellite and SRS) were processed with the Savitzky-Golay filter [4,55,56] to remove residual irregular changes of VI values in the annual cycle. This process was performed using a first-degree polynomial with a window size of 3 (in the case of UAV-derived NDVI time series) and 5 (in case of the PhenoCam, satellite and SRS-based VI time series). Of all the different curve fitting methods explained in [57], univariate spline interpolation was used to fit the time series data [49] from all studied platforms, as it led to the best fit.
The phenology transition dates estimation was based on the fitted spline curves, which were calculated based on the rate of change in curvature [4,57,58]. The rate of change in curvature (k) is defined in Equation (4) as: where ( ) and, ′′( ) are the first and second order derivatives, respectively of the spline fit.
The phenological transition dates for the green-up phase were defined as the times at which the rate of change in curvature exhibits local minima or maxima. These extremes were used to extract start, end, and middle of the spring season (SOS, MOS, and EOS,

Spectral Reflectance Sensor Data Processing
The irradiance data from SRS were calibrated using calibration factors of 0.250 (R) and 0.260 (NIR) in the case of SRS-NDVI sensor, and 0.259 (570 nm) for the SRS-PRI sensor following methodology in [54]. A subset of measurements was extracted from the full time series, which covered exactly between the same hours as the period used for the PhenoCam (9 a.m. to 3 p.m.). The calculation of SRS-NDVI was done on reflectance, calculated as a ratio from incoming and reflected radiation, by dividing the downward-looking by the upward-looking sensor data. The data was averaged in similar time aggregates as the PhenoCam, i.e., 3-day averaged NDVI values. In addition to the NDVI, VIgreen was computed in a similar fashion to the PhenoCam, combining the red band from the NDVI-SRS and one of the green bands at 570 nm from the SRS-PRI sensor. The footprint of the SRS was computed using the sensor characteristics (i.e., FOV: 36 • , height: 7.25 m, off-nadir angle: 45 • ) and the tower coordinates. The same footprint was used across all sensors to extract time series of VI values for effective comparison.

Curve Fitting to Time Series Data
The original time series data obtained from different platforms (UAV, PhenoCam, satellite and SRS) were processed with the Savitzky-Golay filter [4,55,56] to remove residual irregular changes of VI values in the annual cycle. This process was performed using a first-degree polynomial with a window size of 3 (in the case of UAV-derived NDVI time series) and 5 (in case of the PhenoCam, satellite and SRS-based VI time series). Of all the different curve fitting methods explained in [57], univariate spline interpolation was used to fit the time series data [49] from all studied platforms, as it led to the best fit.
The phenology transition dates estimation was based on the fitted spline curves, which were calculated based on the rate of change in curvature [4,57,58]. The rate of change in curvature (k) is defined in Equation (4) as: where f (t) and, f (t) are the first and second order derivatives, respectively of the spline fit. The phenological transition dates for the green-up phase were defined as the times at which the rate of change in curvature exhibits local minima or maxima. These extremes were used to extract start, end, and middle of the spring season (SOS, MOS, and EOS, respectively), in a similar fashion to [11,57]. The local extreme values approximately correspond to 10, 50 and 90% of amplitude in the spring green-up phase of vegetation growth [57].

Sensor Inter-Comparison
The time series data and phenophase transition dates from the studied platforms were analysed at three levels: (1) graphical comparison of VI time series from the different sensors, (2) regression analyses between the time series information, and (3) estimation of the biases in number of days of the transition dates, compared to visual inspection of PhenoCam photos. All three steps were applied to curve-fitted time series data. For regression analysis, using root mean square deviation (RMSD) and Pearson's correlation coefficient were used.
NDVI was originally planned as the common VI to compare across the different studies sensors, as it has been widely used in many remote sensing studies. However, since NDVI cannot be calculated from PhenoCam data, as it does not have NIR band, GCC was used instead to compare PhenoCam VIs against satellite derived indices. In the same manner, GCC cannot be calculated for SRS, as a blue band is not available. Therefore, VIgreen was used as a common VI to compare data from all platforms, as red and green bands were available for all studied sensors. Since the different VI ranges are expressed in different units and scale, they were normalized by using the min-max normalization [59], as a means to match all the values to be in the range of 0-1, for graphical comparison of VIs across different sensors, except for the comparison of NDVI from different platforms (Section 5.2.). For NDVI, we used the original values (between −1 and 1) and not the normalized ones. For the regression analysis, also original VI values were used.
The calculation of RMSD and correlation coefficients included fitted time series data during the green-up, growing and senescence phases (DOY 70−300) across all systems. This was done to avoid the snow season and focus on the vegetation growing season.
For estimating bias in transition dates, SOS was defined as the DOY when the green sprouts started to be visible on the PhenoCam images, while EOS was defined as the day when sprouts were fully developed and could no longer be differentiated from the rest of the branches. In SOS, the sprouts are easily recognizable as they are bright green, while EOS was defined when the sprouts get a darker colour, similar to the rest of the tree (Figure 8). respectively), in a similar fashion to [11,57]. The local extreme values approximately correspond to 10, 50 and 90% of amplitude in the spring green-up phase of vegetation growth [57].

Sensor Inter-Comparison
The time series data and phenophase transition dates from the studied platforms were analysed at three levels: (1) graphical comparison of VI time series from the different sensors, (2) regression analyses between the time series information, and (3) estimation of the biases in number of days of the transition dates, compared to visual inspection of Phe-noCam photos. All three steps were applied to curve-fitted time series data. For regression analysis, using root mean square deviation (RMSD) and Pearson's correlation coefficient were used.
NDVI was originally planned as the common VI to compare across the different studies sensors, as it has been widely used in many remote sensing studies. However, since NDVI cannot be calculated from PhenoCam data, as it does not have NIR band, GCC was used instead to compare PhenoCam VIs against satellite derived indices. In the same manner, GCC cannot be calculated for SRS, as a blue band is not available. Therefore, VIgreen was used as a common VI to compare data from all platforms, as red and green bands were available for all studied sensors. Since the different VI ranges are expressed in different units and scale, they were normalized by using the min-max normalization [59], as a means to match all the values to be in the range of 0-1, for graphical comparison of VIs across different sensors, except for the comparison of NDVI from different platforms (Section 5.2.). For NDVI, we used the original values (between −1 and 1) and not the normalized ones. For the regression analysis, also original VI values were used.
The calculation of RMSD and correlation coefficients included fitted time series data during the green-up, growing and senescence phases (DOY 70−300) across all systems. This was done to avoid the snow season and focus on the vegetation growing season.
For estimating bias in transition dates, SOS was defined as the DOY when the green sprouts started to be visible on the PhenoCam images, while EOS was defined as the day when sprouts were fully developed and could no longer be differentiated from the rest of the branches. In SOS, the sprouts are easily recognizable as they are bright green, while EOS was defined when the sprouts get a darker colour, similar to the rest of the tree (Figure 8).  Since the studied forest was evergreen and the changes in vegetation could not be differentiated with a human eye, it was not possible to define the dates for other seasonality events (i.e., end of the growing season and beginning of the senescence phase). The agreement between the phenophase dates from all sensors were compared against visual inspection data in the PhenoCam imagery, and the phenophase dates calculated statistically, expressed as bias in number of days.

Spline Fit Time Series of VIs Derived from Different Sensors
The time series of different VIs fitted with the spline interpolation technique are depicted in Figure 9. The snow period at the beginning and end of the year (blue rectangles in the figure based on PhenoCam images) influenced the VI values in almost all of the sensor types. The middle of the season peak is much more pronounced in GCC (Figure 9a-c) and VIgreen (Figure 9d-f), with well-defined seasonal start, peak, and end, and similar in shape, in general, across different sensors. The exception is VIgreen for SRS (Figure 9g), which does not show a clear peak in the middle of the growing season, nor pronounced slopes at the beginning and end of it. The fitted curve for Sentinel-2 NDVI (Figure 9i) is flatter compared to VIgreen and GCC, with less pronounced slopes at the beginning and end of the growing season. Nevertheless, as per VIgreen and GCC, the shape of NDVI is similar across sensors. UAV-based VI values show the spline fit only during the middle part of growing season as the flights were conducted only during that period (Figure 9k). Since the studied forest was evergreen and the changes in vegetation could not be differentiated with a human eye, it was not possible to define the dates for other seasonality events (i.e., end of the growing season and beginning of the senescence phase). The agreement between the phenophase dates from all sensors were compared against visual inspection data in the PhenoCam imagery, and the phenophase dates calculated statistically, expressed as bias in number of days.

Spline Fit Time Series of VIs Derived from Different Sensors
The time series of different VIs fitted with the spline interpolation technique are depicted in Figure 9. The snow period at the beginning and end of the year (blue rectangles in the figure based on PhenoCam images) influenced the VI values in almost all of the sensor types. The middle of the season peak is much more pronounced in GCC (Figure 9 a-c) and VIgreen (Figure 9d-f), with well-defined seasonal start, peak, and end, and similar in shape, in general, across different sensors. The exception is VIgreen for SRS ( Figure  9g), which does not show a clear peak in the middle of the growing season, nor pronounced slopes at the beginning and end of it. The fitted curve for Sentinel-2 NDVI (Figure 9i) is flatter compared to VIgreen and GCC, with less pronounced slopes at the beginning and end of the growing season. Nevertheless, as per VIgreen and GCC, the shape of NDVI is similar across sensors. UAV-based VI values show the spline fit only during the middle part of growing season as the flights were conducted only during that period (Figure 9k,l).
In the following sections, we plot each of the VIs for the different sensors, for a more detailed comparison. In the following sections, we plot each of the VIs for the different sensors, for a more detailed comparison. Figure 10 shows a plot of spline fitted curves from UAV NDVI , SRS NDVI , MODIS NDVI and Sentinel-2 NDVI . The plot indicates that the four datasets follow the similar trend in general throughout the growing season. However, they differ at the start and end of the time series, as they show significant different slopes and offset. The discrepancy between the SRS NDVI , Sentinel-2 NDVI and MODIS NDVI is most pronounced at the beginning of the year until end of February, being the satellite NDVI values lower than the SRS NDVI values, and significantly lower for Sentinel-2 NDVI also at the end of the season, from October to December. These low NDVI values could be due to the snow cover, higher decline in satellite data as they are viewing more snow-covered ground than the SRS NDVI , which also show reflectance leaves and branches. Additionally, there can be a different light scattering at the top of the atmosphere, observed by the satellites, than near the ground, as observed by the SRS NDVI . After May, we observe a stabilization of NDVI values until the end of October (please note that, for the UAV NDVI sensor, the time series ends in September). There is, however, certain offset between NDVI values across sensors, ranging from 0.8 (UAV NDVI ) to 0.9 (SRS NDVI ) at the peak of the season. As NDVI ranges from −1 to 1, this difference accounts for a 5% difference in reflectance between UAV NDVI and SRS NDVI . The green-up phase (March to May) is depicted very differently across sensors. In this time of the year, the differences cannot be attributed to snow. The slope for satellite data is significantly higher than for SRS NDVI , and the offset rises to a 25% between MODIS NDVI (0.4) and SRS NDVI (0.9).

Comparison of NDVI Derived from Different Sensors
Sentinel-2-NDVI, (j) SRS-NDVI, (k) UAV-VIgreen, and (l) UAV-NDVI. Grey circles represent the normalized raw VI values from respective sensors. Phenophase dates are marked by the vertical dashed lines (legend). The SOS and EOS from visual assessments of PhenoCam images are also shown for comparison purposes. Figure 10 shows a plot of spline fitted curves from UAVNDVI, SRSNDVI, MODISNDVI and Sentinel-2NDVI. The plot indicates that the four datasets follow the similar trend in general throughout the growing season. However, they differ at the start and end of the time series, as they show significant different slopes and offset. The discrepancy between the SRSNDVI, Sentinel-2NDVI and MODISNDVI is most pronounced at the beginning of the year until end of February, being the satelliteNDVI values lower than the SRSNDVI values, and significantly lower for Sentinel-2NDVI also at the end of the season, from October to December. These low NDVI values could be due to the snow cover, higher decline in satellite data as they are viewing more snow-covered ground than the SRSNDVI, which also show reflectance leaves and branches. Additionally, there can be a different light scattering at the top of the atmosphere, observed by the satellites, than near the ground, as observed by the SRSNDVI. After May, we observe a stabilization of NDVI values until the end of October (please note that, for the UAVNDVI sensor, the time series ends in September). There is, however, certain offset between NDVI values across sensors, ranging from 0.8 (UAVNDVI) to 0.9 (SRSNDVI) at the peak of the season. As NDVI ranges from −1 to 1, this difference accounts for a 5% difference in reflectance between UAVNDVI and SRSNDVI. The green-up phase (March to May) is depicted very differently across sensors. In this time of the year, the differences cannot be attributed to snow. The slope for satellite data is significantly higher than for SRSNDVI, and the offset rises to a 25% between MODISNDVI (0.4) and SRSNDVI (0.9).  Figure 6).

Comparison of NDVI Derived from Different Sensors
According to the Pearson's regression analysis, the NDVI values from all sensors were positively correlated with each other ( Table 4). The correlation between UAVNDVI and MODISNDVI was the lowest (r = 0.75; RMSD = 0.08), while the highest correlation was found between MODISNDVI and Sentinel-2NDVI (r = 0.92; RMSE = 0.08), indicating a better fit between the two satellite datasets. Figure 10. NDVI values comparing seasonal trends from Sentinel-2 NDVI , MODIS NDVI , SRS NDVI and UAV NDVI . All the NDVI values are computed for the ground projected footprint of the SRS (dashed red polygon in Figure 6).
According to the Pearson's regression analysis, the NDVI values from all sensors were positively correlated with each other ( Table 4). The correlation between UAV NDVI and MODIS NDVI was the lowest (r = 0.75; RMSD = 0.08), while the highest correlation was found between MODIS NDVI and Sentinel-2 NDVI (r = 0.92; RMSE = 0.08), indicating a better fit between the two satellite datasets.

Comparison of GCC from Different Sensor Types
The PhenoCam GCC , MODIS GCC and Sentinel-2 GCC followed quite similar seasonal pattern ( Figure 11). The plot shows that the canopy GCC signal followed a similar green-up response, while the GCC trajectories separate a little after reaching the peak of the growing season (PhenoCam GCC vs. satellite GCC ) and during the senescence phase (MODIS GCC against the rest). Additionally, the peak of PhenoCam GCC is separated by a vertical offset of around 0.2 units from the satellite GCC (0.78-0.9). Moreover, the peak of the season is slightly different for PhenoCam GCC (around June 1st) than for satellite GCC (around 24 June).

Comparison of GCC from Different Sensor Types
The PhenoCamGCC, MODISGCC and Sentinel-2GCC followed quite similar seasonal pattern ( Figure 11). The plot shows that the canopy GCC signal followed a similar green-up response, while the GCC trajectories separate a little after reaching the peak of the growing season (PhenoCamGCC vs. satelliteGCC) and during the senescence phase (MODISGCC against the rest). Additionally, the peak of PhenoCamGCC is separated by a vertical offset of around 0.2 units from the satelliteGCC (0.78-0.9). Moreover, the peak of the season is slightly different for PhenoCamGCC (around June 1st) than for satelliteGCC (around 24 June). The Pearson's correlation coefficients computed between the GCC values from these sensors are presented in Table 5. PhenoCamGCC was in good agreement with MODISGCC and Sentinel-2GCC as reflected in high positive correlation coefficients (0.93 and 0.97, respectively) with approximately 10% deviation on an average. The Pearson's correlation coefficients computed between the GCC values from these sensors are presented in Table 5. PhenoCam GCC was in good agreement with MODIS GCC and Sentinel-2 GCC as reflected in high positive correlation coefficients (0.93 and 0.97, respectively) with approximately 10% deviation on an average.

Comparison of VIgreen from Different Sensor Types
Plots of VIgreen for all the sensor types are shown in Figure 12. The time series obtained from PhenoCam VIgreen , UAV VIgreen , MODIS VIgreen , and Sentinel-2 VIgreen appeared to follow a similar seasonal pattern, compared to the one obtained from SRS VIgreen , where the time series are flatter along the growing season. A detailed comparison shows that for PhenoCam VIgreen the green-up phase starts later (DOY~125) and the slope is almost vertical. Another significant difference is observed in UAV VIgreen , which peaks later than PhenoCam VIgreen and satellite VIgreen data, around DOY = 200 vs. DOY~150. Additionally, the VIgreen values differ in the senescence phase; the slope is more pronounced and starts earlier for PhenoCam VIgreen ; it is a gradual slope for MODIS VIgreen , and also a steep slope and later slope for Sentinel-2 VIgreen and UAV VIgreen .

Comparison of VIgreen from Different Sensor Types
Plots of VIgreen for all the sensor types are shown in Figure 12. The time series obtained from PhenoCamVIgreen, UAVVIgreen, MODISVIgreen, and Sentinel-2VIgreen appeared to follow a similar seasonal pattern, compared to the one obtained from SRSVIgreen, where the time series are flatter along the growing season. A detailed comparison shows that for PhenoCamVIgreen the green-up phase starts later (DOY~ 125) and the slope is almost vertical. Another significant difference is observed in UAVVIgreen, which peaks later than Phe-noCamVIgreen and satelliteVIgreen data, around DOY = 200 vs. DOY~150. Additionally, the VIgreen values differ in the senescence phase; the slope is more pronounced and starts earlier for PhenoCamVIgreen; it is a gradual slope for MODISVIgreen, and also a steep slope and later slope for Sentinel-2VIgreen and UAVVIgreen. The Pearson's correlation coefficients computed for VIgreen from these sensors are presented in Table 6. PhenoCamVIgreen, MODISVIgreen and Sentinel-2VIgreen were not so strongly correlated (r = 0.63, 0.64 and 0.58, respectively) with SRSVIgreen, while the rest of sensor types are moderately correlated among each other (r > 0.7). The Pearson's correlation coefficients computed for VIgreen from these sensors are presented in Table 6. PhenoCam VIgreen , MODIS VIgreen and Sentinel-2 VIgreen were not so strongly correlated (r = 0.63, 0.64 and 0.58, respectively) with SRS VIgreen , while the rest of sensor types are moderately correlated among each other (r > 0.7).

Evaluating Phenophase Transition Dates from All Platforms
Phenological transition dates derived from visual inspection were compared with those from spline fitted time VIs from the different platforms to detect the different phases of the growing season in the studied spruce forest. The bias (in number of days) between the reference start and end of season (based on visual inspection) compared to the phenophase transition dates calculated from VIs, for all the platforms except the SRS, are shown in Table 7 and Figure 9. SRS-based phenophase dates are not mentioned in the table and are not considered in the comparison with phenophase from other platforms. Irregular patterns of SRS-based VIs did not allow us to extract meaningful phenophase dates from the spline fitted curves, due to the flat nature of the SRS-VI profiles and the method employed for phenophase extraction (change in curvature). According to the visual inspection of PhenoCam images, the Start of the Spring (SOS) happened in DOY: 90 (31 March) and ended in DOY: 150 (30 May). Seasonality parameters, mainly the SOS, MOS, and EOS extracted from all the platforms using the spline fitted time series of GCC, VIgreen and NDVI were observed to be consistently similar, except for the VIgreen from Sentinel-2, which showed very early SOS when compared to visual inspection data. Visually assessed dates (SOS and EOS) compared against dates estimated from other sensors showed a bias ranging from 3-10 days. The negative bias refers to earlier SOS and EOS. The maximum bias in SOS and EOS was observed in seasonality extracted from MODIS and Sentinel derived NDVI values, which accounted for 8 days in case of MODIS and 10 days in case of Sentinel-2. The minimum bias across all sensor types were the ones extracted from GCC with bias ranging between 3 and 6 days, while in case of VIgreen, the biases in SOS and EOS ranged between 3 and 9 days. No UAV flights were conducted as early as the SOS; hence, it was not possible to compute the bias.

Discussion
Using remote sensors to study vegetation phenology at different scales opens up the opportunity to gather information on different aspects of the life cycle of plants to better understand the interaction between climate change and the biosphere. While different nearsurface and satellite sensors are available for vegetation phenology monitoring, few studies provide sensor overview and inter-comparisons [1,4,24,25]. In this research, an effort has been made to compare the performance of different remote sensors to characterize the phenology of a forest. Using combined methods of UAV photogrammetry, phenocamera data analysis, and multispectral sensor data analysis (SRS and satellite), we explored the difference between phenological time series obtained from these platforms, by means of spectral vegetation indices such as NDVI, VIgreen and GCC.
In this study, we observed the complementarity of the data from the studied remote sensors for characterizing vegetation phenology in an evergreen coniferous forest. PhenoCams and SRS offer the flexibility of defining regions of interest within very fine spatio-temporal resolution images (hourly, tree level), and thus allows researchers to model the phenology dynamics of the observed ecosystem, for either individual or a group of species within the landscape. On the other hand, satellite observations cover large regions for global observations [5], at coarser spatial (decametres to kilometres) and temporal resolutions (days to months). UAVs have the benefit of a fine pixel size (nominally, a few centimetres) and can be operated frequently and avoid atmospheric effects [4].
Our results show that there are similarities between VIs of the different studied sensors in terms of shape of the phenology curves as indicated by Pearson' correlation and RMSD showing good agreement between the VIs across the studied sensors. However, there are significant differences in slopes and offsets when comparing the time series of different sensors, for a given VI. As expected, the time series of satellite data (i.e., MODIS and Sentinel-2) are closely correlated and follow a similar time series shape. On the other hand, depending on the VI (NDVI, VIgreen or GCC), we observed different behaviours among sensors. For instance, we noted that UAV NDVI data show a similar pattern to MODIS NDVI and Sentinel-2 NDVI during the summertime (April to August). Additionally, there is a good agreement between PhenoCam GCC and satellite GCC for the green-up season (April to May). Regarding VIgreen, there is a relatively good match in determining the peak of the growing season (DOY~150) between PhenoCam VIgreen , MODIS VIgreen and Sentinel-2 VIgreen , while the green-up slope of UAV VIgreen and satellite VIgreen also follow similar patterns. However, these are the only agreements between sensors; apart from these, all sensors differ in slopes on the green-up and senescence, the maximum VI value at the peak of the season, and the phenology transition dates (i.e., beginning, end and length of the growing season). A remarkable difference was observed in the behaviour of SRS sensors; SRS NDVI and SRS VIgreen show significantly higher values than any other sensors, and the length of the growing season is notably longer than for the rest of sensors. Additionally, RMSD related to SRS versus other sensors were generally large.
There are many reasons for the differences between the phenology time series, as described by the different VIs and sensor types. Primarily, it must be noted that the visual seasonal cycle is considerably weaker in evergreen forests than in deciduous forests, and the low annual production of green buds in boreal coniferous forests, in combination with the influence of snow, presents considerable difficulties for detecting the green-up and senescence phases using remote sensing [60]. Further reasons for the different performances of the sensors are the different spatial resolutions, different spectral band configurations, sensor calibrations, different viewing angles, or differences in data processing [14,24,25,32]. In addition, the spectral sensors have different spectral response functions, which could cause systematic deviations in the time series of spectral VIs derived from them [24,61]. One important factor that might cause differences in the VIs is the viewing angle. In particular, the PhenoCam points towards the horizon and views the canopies of the studied trees (i.e., profile view), while the other studied sensors view from the zenith, viewing the canopy of the trees vertically, including the understory vegetation, soil and snow, through the forest clearings. The varying viewing angles of the studied sensors influence the bidirectional reflectance distribution function (BRDF) [16]. This affects the spectral and VI values captured by each sensor type [4]. Moreover, the sun-sensor target viewing geometry is different between the near-ground sensors and the satellites. If we consider the geometry of a typical spruce forest as a collection of conical trees, the shadows projected by the trees influence the reflectance of the forest very differently [62,63], depending on the position of the observer (the sensor). If the sensor and the sun are aligned, the sunlight side of the forests is in the view field, whereas, if the observed target is between the sensor and the sun, the shaded side of the forest is in the view field [62].
Another potential reason for the offset in VIs is the difference in the spectral range and bandwidth of each sensor. The highest difference in bandwidth range is observed in NIR and R (Sequoia: 40 nm; MODIS: 35 nm; Sentinel-2: 106 nm; SRS: 10 nm, for R and NIR). Additionally, the central wavelengths of the NIR bands are different ( Table 2). Even though the central wavelengths of G and R bands were close in all sensors, there are significant differences in bandwidth (Sequoia: 40/40 nm, MODIS: 20/50 nm, Sentinel-2: 36/31 nm, and SRS: 10/10 nm for G and R, respectively). These differences affect the spectral VI values (NDVI, GCC and VIgreen) among sensors. Another reason for deviations in values of the VI time series could be the mixed reflectance characteristics caused by the different spatial resolution of the observed sensors. In addition, the specific days used for image compositing, observation time, and solar elevation angle might add differences in the measured spectral values [64]. In this study, there was a huge difference in spatial resolution of the sensors (MODIS: 250 m and 500 m, Sentinel-2: 10 m, UAV: 8 cm, PhenoCam:~17 m and SRS:~8 m). With the increase or decrease of remote sensing image scale, the observed targets within a pixel at the ground surface, as seen by each of these sensors, are quite different. This will result in generation of mixed pixels [16] causing different signal intensities of mixed objects, when compared with the signal intensity of large-scale images. Despite the relative homogeneity of the analysed spruce forest, it is likely that data from the satellite sensors are affected by mixed targets, e.g., understory vegetation and soil. In addition, even though we selected the same observation time for UAV, SRS and PhenoCam for inter-comparison, there was a different acquisition time for the satellite sensors.
Despite the different performance of the selected VIs and studied platforms, we can conclude that different sensors address different aspects of the vegetation. These platforms are complementary when characterising phenology and they should be used for different purposes. As observed in our study, we found strong correlation and similar curve fitting between UAV NDVI and satellite NDVI . Therefore, UAV can be a solid tool for upscaling from near-ground observations to satellite data. Despite the differences in pixel size, Parrot Sequoia has a similar band configuration as MODIS and Sentinel-2 [41]. That is not the case for UAV VIgreen versus MODIS VIgreen and Sentinel-2 VIgreen , where the time series offset and slopes are very different. Our results coincide with the findings of [59], where VIgreen index from multiple cameras had the lowest correlation to MODIS in comparison to the GCC. Despite VIgreen being a normalized index like NDVI, the NIR band seem to play an important role in differentiating vegetation productivity [65,66]. We also observed a good agreement between PhenoCam GCC and the satellites GCC , especially in the green-up phase, as was also observed by [57,67], despite GCC not using the NIR band, and both sensor types having different viewing angles and pixel sizes [14]. This finding validates MODIS and Sentinel-2 as useful sensors for characterizing vegetation phenology at the landscape scale. In summary, the results suggest that VIgreen is not a good proxy for phenology, while NDVI and GCC are, as was previously observed by [2,4,13].
We observed strongly significant different time series pattern between NDVI and GCC ( Figure 9). NDVI saturates at a high value throughout the growing season, with a short green-up and senescence phases. Additionally, for SRS NDVI there is not a significant increase in NDVI between the green and the non-green seasons. On the contrary, GCC presents a sharper peak in the middle of the growing season, with steep slopes for the senescence phase, and especially for the green-up phase. Therefore, we hypothesize that NDVI, through its sensitivity to red light absorption, is more sensitive to the overall albedo change going from a snowy landscape to the fully developed green canopy. The sensitivity of NDVI is well documented [60,68,69]. On the other hand, the GCC normalizes for albedo and thereby amplifies the green signal of the canopy that relates more directly to the phenology phases. These VIs relate to different aspects of the changing landscape during spring. To note is that both VIs work well at both near-surface and landscape scales, as shown by the agreement between UAV NDVI vs. satellite NDVI and PhenoCam GCC vs. satellite GCC . Therefore, PhenoCams and UAVs have important and complementary roles for satellite data validation and upscaling [26,30,57,59,70].
An important observation of this study is the differences in offset between platforms. As mentioned, the shapes of the time series and correlations among VIs across sensors was strong, especially in summer time. However, we found significant differences in the slopes of the green-up and senescence phases, as well as significant offset in VI values, across sensor types. This suggests that these remote sensors can be interchangeable for qualitative characterization of vegetation phenology, especially of the growing season, but cannot be used indistinctively for quantitative studies or accurate determination of phenology phases [4,59,70]. Another result that supports this is that the statistical calculation of the phenology transition dates (SOS, MOS and EOS) using fitted time series curves from VIs shows certain biases with the real dates, captured in the PhenoCam images by visual inspection. The bias in estimating SOS and EOS was around three days for GCC, while it was up to 8-10 days for NDVI. For VIgreen, the bias ranged between 3 and 8 days. This supports the conclusion that NDVI and VIgreen may be less appropriate estimators of phenology phases than GCC [4,14,59].
A second conclusion of our study is that there is a large uncertainty and differences in VI performance among sensors during the winter season, when there is snow on the ground and trees. The presence of snow beneath the forest canopy is assumed to remain longer compared to the canopy mostly during the spring and can be detected differently across near-surface and satellite sensors [71]. This eventually affects SOS events derived from time series obtained from different platforms [72]. While the VI time series are somehow comparable when there is only vegetation present, the behaviour of the VIs is quite random from sensor to sensor in the presence of snow. Complex light scattering and viewing geometry, as well as the likelihood of snow on up-looking sensors, affect the different VI signals observed for snow [16,73,74].

Conclusions
In the present paper, near-surface and satellite remote sensors at different scales were inter-compared for the purpose of vegetation phenology characterization. We defined phenology time series from multispectral ground sensors (SRS), PhenoCam, UAV, Sentinel-2 and MODIS data, depicted by three different spectral vegetation indices (VIs), i.e., NDVI, GCC and VIgreen. These platforms differ significantly in spatial and temporal resolutions, viewing angles and observation scales. The aim of this study was therefore to explore differences in the time series phenology curves from sensor to sensor and for different Vis, and their performances in determining the transition dates in phenology along the year.
Based on the results, we conclude that the studied remote sensing platforms (SRS, UAV, PhenoCam, MODIS and Sentinel-2) and VIs (NDVI, GCC and VIgreen) all produce useful data for phenology research, although the platforms and indices serve different purposes.
GCC is the recommended VI for characterizing phenology phases and transition dates, whereas NDVI is suitable for detecting when the landscape changes from winter to summer following snow-melt and canopy green-up. The tested satellite data were validated by the near-surface remote sensors as useful tools for phenology and productivity studies, since PhenoCam GCC and satellite GCC correlate well, similarly to UAV NDVI and satellite NDVI , which show good agreement. This proves the potential for upscaling vegetation phenology from near-surface to landscape level by using a combination of sensor systems.