Phenology of short vegetation cycles in a Kenyan rangeland from PlanetScope and Sentinel-2

.


Introduction
Rangelands are the land use with the largest spatial extent globally, and comprise various land cover types such as savannahs, grasslands, prairies, steppe, and shrubland.In Africa alone, rangelands cover over 60% of the land surface and are an essential resource for wildlife and livestock (Reid et al., 2008).Millions of pastoral households depend on their livestock for milk and meat production (Sayre et al., 2013).Their livelihoods are strongly affected by weather shocks, such as drought (Blackwell, 2010;Little et al., 2008).During the past century, climatic shifts have resulted in increasing rainfall variability in rangelands (Sloat et al., 2018), which in turn affects rangeland productivity (Briske et al., 2015;Knapp et al., 2008) and composition (Scheiter and Higgins, 2009), and consequently the animals and humans that depend on the rangelands for forage and livelihoods.The severity of impacts of recurrent droughts depends on their time and on rangeland type.Besides rainfall variability, factors such as carbon dioxide enrichment (e.g., Morgan et al., 2007), soil properties, herbivory, and fire regime equally influence rangeland status (Briske et al., 2005;House et al., 2003).
Assessing phenological dynamics of vegetation permits us to better understand how climatic variability and management affect different rangeland functional groups.For example, the deployment and retention of leaves differs between grasses and trees (Higgins et al., 2011;Liu et al., 2017;Scholes and Archer, 1997) and this temporal separation in photosynthetic activity may reduce competition for water (Ogle and Reynolds, 2004;Sankaran et al., 2004).Vegetation phenology in turn also influences ecosystem functions (e.g., Gross, 2017;Guuroh et al., 2018), among others due to the timing and abundance of resource availability, affecting migratory species (Cleland et al., 2007;Miles et al., 2017), including desert locusts, which caused a major plague in 2019-2020 in East Africa and South Asia (Meynard et al., 2020).
Time series of optical remote sensing data are a critical input to studying phenology at the landscape level.One option for this can be referred to as near-surface remote sensing, i.e. mounting a camera that frequently takes pictures of vegetation from a fixed position in the field.This digital repeat photography allows to either directly observe phenological events like flowering, or to automatically assess changes in greenness by deriving vegetation indices, such as the Greenness Chromatic Coordinate (GCC) from the photograph series (Richardson et al., 2018).For example, Alberton et al. (2014) used camera-derived GCC series to show interspecies phenological differences for a tropical savannah system in Brazil.While this approach can provide useful insights, particularly if employed for long time periods, satellite imagery is the sole option to synoptically assess phenology for larger areas.Moderate to coarse spatial resolution data (100-1000 m) have been used for this purpose, because their short revisit times increase the probability of obtaining sufficient cloud-free observations throughout the year.Typically these observations are subjected to cloud and shadow screening, temporal compositing, and smoothing of resulting vegetation indices to achieve a reliable description of temporal changes in vegetation activity (Hmimina et al., 2013;Zeng et al., 2020).
Estimates of vegetation phenology have been frequently obtained for the African continent, or parts thereof (for a review see Adole et al., 2016), using vegetation index series derived from various moderate and coarse resolution sensors, including the Moderate Resolution Imaging Spectroradiometer (MODIS), the Satellite Pour l'Observation de la Terre (SPOT) Vegetation, and the Advanced Very High Resolution Radiometer (AVHRR) series.Estimates of vegetation phenological timings have various practical applications.For instance they contribute to drought early warning systems and index insurance schemes by indicating when crops or rangelands are in a productive phase (Rembold et al., 2019;Vrieling et al., 2016).Start-and end-of-season (SOS/EOS) estimates are also inputs to calculating a seasonal accumulation of a spectral vegetation index, an often-used proxy for seasonal biomass production (Diouf et al., 2015;Meroni et al., 2014a;Schucknecht et al., 2017).Trends and interannual variability can be assessed when estimating phenological estimates with satellite series over longer time periods (Brown et al., 2010;Heumann et al., 2007), but attribution of this variability to specific factors (i.e.management, climate) can be cumbersome (Vrieling et al., 2011).Water availability is known to be the major determinant of phenological variability in rangelands.Growth onset of herbaceous vegetation tends to follow precipitation onset, but the green-up of woody vegetation can show a deviating behaviour (Archibald and Scholes, 2007).In some cases, woody vegetation becomes photosynthetically active before rains start, which has been attributed to day length cues (Archibald and Scholes, 2007) and the tree root system's access to deep ground water (Guan et al., 2014).In other cases woody green-up may also be delayed with respect to herbaceous vegetation (Liu et al., 2017).Taking advantage of dissimilar responses of vegetation to water variability, phenological behaviour as captured by satellites can be utilized for differentiating functional groups or vegetation communities, as for example Hüttich et al. (2009) illustrated for the Kalahari in Namibia.Other applications include a prediction of seasonal movements of livestock herds (Butt et al., 2011).While some success was achieved in estimating fractional cover of trees and grasses from MODIS-derived phenology in Kruger National Park (Ibrahim et al., 2018), Hmimina et al. (2013) showed that the large spatial heterogeneity of rangelands cannot properly be represented by MODIS-sized pixels (250 m), resulting in poor overall phenology estimates.
Finer spatial resolution satellite sensors (≤30 m) have potential to more accurately describe the phenology of spatially heterogeneous rangelands.However, the longer revisit times of such sensors, in combination with frequent cloud cover, often makes it hard to obtain sufficient cloud-free observations to effectively describe seasonal changes in greenness.To overcome this, a few main strategies have been proposed: a) Fusion of Landsat with coarser-resolution MODIS or VIIRS (Visible Infrared Imaging Radiometer Suite) series to reconstruct a denser dataset at Landsat resolution, which is then used as the basis for phenology estimation (e.g., Frantz et al., 2016;Walker et al., 2014;Zhang et al., 2020).While promising results have been attained, this approach is less effective for heterogeneous landscapes with rapid changes.b) Combination of multiple years of Landsat acquisitions into a single 'synthetic' year (Elmore et al., 2012;Fisher et al., 2006;Nijland et al., 2016).This has been particularly used for locations where adjacent Landsat orbits overlap and consequently result in higher image density.The resulting multi-annual average phenology can subsequently be adjusted to meet individual year observations (Melaas et al., 2013;Melaas et al., 2016).This strategy only works in case the overall shape of the temporal profile remains relatively similar between years.c) Taking advantage of reduced revisit times of new satellite missions, such as Sentinel-2 (10-60 m resolution), phenological metrics may be directly estimated using imagery from a single season (Pan et al. 2015;Schwieder et al., 2018;Vrieling et al., 2018;Vrieling et al., 2017;White et al., 2014).These attempts mostly focussed on areas with long annual vegetative seasons.d) To increase robustness of single-season phenology retrievals, a combination of (b) and (c) was proposed (Jönsson et al., 2018).In this case, if insufficient observations are available for (part of) the season under consideration, specific parameters that describe the temporal trajectory of vegetation indices are pulled from the multiyear average phenology.This concept was also used by Bolton et al. (2020) who fitted cubic splines to generate daily vegetation index series.They assigned weights to alternate-year observations, whereby more weight was given to observations of years with a similar trajectory as the year under consideration.By applying this approach to the new harmonized Landsat 8 and Sentinel-2 dataset (Claverie et al., 2018) they retrieved a 30-m phenology for the North American continent.
Semi-arid rangelands, particularly in East Africa, are characterized by short vegetation cycles that show large interannual variability both in the timing and magnitude of vegetation productivity.In combination with the high degree of cloud cover during rainy seasons, this puts stronger requirements on the observation frequency as compared to many of the 10-30 m -resolution phenology studies executed until present.This is because the density of observations required to reliably sample the dynamic signal of vegetation is roughly inversely proportional to the length of the season.
The aim of this study is to assess if two fine-resolution satellite systems with short revisit times, i.e.PlanetScope and Sentinel-2, can be effectively applied to retrieve and map phenological metrics for a semi-arid rangeland site in Kenya.For the purpose of this paper we will use "fine-resolution" and "fine-scale" to refer to the 3-10 m spatial resolution of these systems and the resulting phenology estimates.Given the large interannual variability in seasonal characteristics, we followed the above-mentioned strategy "c", i.e. the direct estimation of phenology from single-season imagery.Using a model-fit approach on time series of the normalized difference vegetation index (NDVI), we retrieved two temporal phenological metrics (SOS and EOS) and two metrics related to productivity (maximum NDVI and accumulation of NDVI).We then compared the independent retrievals from PlanetScope and Sentinel-2 data to those obtained from in-situ camera photograph series, and from MODIS 250 m observations, as a commonly-used satellite data source for phenology studies.Our purpose is not to propose a new robust phenology extraction methodology (as attempted in other papers, e.g.Bolton et al., 2020), but rather to process different data sources with a standard curve-fitting approach in order to evaluate whether the data quality and frequency is sufficient to obtain accurate and spatially-consistent estimates of phenological metrics in a challenging environment with short and variable vegetation seasons.

Study area
The study is executed for the Kapiti Farm research station in Machakos County, Kenya (Fig. 1).The Kapiti Farm is property of the International Livestock Research Institute (ILRI) and used to conduct research on sustainable livestock and rangeland management.The farm covers approximately 128 km 2 on which about 2500 heads of beef cattle, 1200 sheep, and 250 goats are maintained.In addition, it serves as an important habitat and ecological corridor for wildlife.The vegetation is predominantly herbaceous with various densities of shrub cover (mostly Acacia drepanolobium).The dominant herbaceous species in Kapiti Farm are Themeda (red oat grass), Panicum (switchgrass), Chloris (finger grass), Pennisetum (fountain grass), Cenchrus (African foxtail grass), and Setaria (bristle grass).On the central higher part of Kapiti, Balanites trees (Balanites aegyptiaca, > 2 m) are dominant.

PlanetScope imagery
PlanetScope is a satellite constellation consisting presently of 130+ CubeSats (4-kg satellites) operated by Planet Labs (Planet Labs Inc, 2020).The majority of these CubeSats are in a sun-synchronous orbit with an equator crossing time between 9:30 and 11:30 (local solar time) (Planet Labs Inc, 2020).PlanetScope imagery contains four spectral bands, i.e. blue (455-515 nm), green (500-590 nm), red (590-670 nm), and near-infrared (NIR, 780-860 nm).We used the Level-3B surface reflectance products that have been atmospherically corrected by Planet Labs using the 6S radiative transfer model with ancillary data from MODIS (Kotchenova and Vermote, 2007;Kotchenova et al., 2006;Planet Labs Inc, 2020).The spatial resolution of these products is 3 m and the RMSE positional error is less than 10 m, and we found that different PlanetScope acquisitions were well co-registered (no apparent shifts).
A total of 760 Level-3B surface reflectance products between March 2017 and February 2019 were downloaded.The average revisit interval for Kapiti Farm is 2.22 days.The range (mean ± standard deviation) of view zenith angle, sun azimuth, and sun altitude for all of these surface products is 1.17 ± 1.44, 95 ± 36, and 54 ± 5 degrees, respectively.Each product contained about 200 km 2 that partially overlapped Kapiti Farm in a non-consistent way.Images acquired on the same day from the same satellite and orbit were mosaicked and clipped to the study area extent, resulting in 313 mosaicked images.
Although Planet Labs has recently improved their cloud and shadow masking algorithm (i.e., the so-called "unusable data mask"), its results were not available for the entire temporal range of this study.We therefore performed cloud and shadow screening ourselves.Existing cloud screening algorithms for Sentinel-2, MODIS, and other satellite imagery mostly use thresholds of spectral bands and derived indices (e.g., Ackerman et al., 1998;Hollstein et al., 2016).Given that spectral information in PlanetScope imagery is limited to four bands, we included texture features using the Gray Level Co-occurrence Matrix (GLCM; Bai et al., 2016;Soh and Tsatsoulis, 1999) to take advantage of the fact that clouds and shadows generally have a smoother texture than the land surface, particularly in fine-resolution imagery (Li et al., 2017;Lu, 2007;Schröder et al., 2002).We generated 120 spectral and texture features from four bands of PlanetScope imagery, including a large set of band algebra combinations and per band 18 GLCM features calculated for 9 × 9 pixel moving windows.We then used a Random Forest classifier to test the relative importance of the four spectral bands plus 120 features in discriminating cloud and cloud shadow.A reference set was made by visually interpreting 23,782 pixels from multiple images that were either clear, or had clouds, or cloud shadows; half of this set (11,981 pixels) were used for training, and the other half for validation.The features of highest relative importance for discriminating clouds were found to be band 4 (NIR) reflectance and its GLCM sum average (b4_savg).For cloud shadows the best performing combination was band 2 (green) and its GLCM sum average (b2_savg).Because shadows reduce the estimated reflectance of the land surface, which in itself varies over time (e.g., due to phenological change), the thresholds for classifying clouds and cloud shadows were set separately for each month in each year.The decision tree (Fig. A1) for the classification of clear, clouds, and cloud shadows combines the monthly thresholds of b2, b2_savg, b4 and b4_savg.A Python script (https://gist.github.com/YanCheng-go/6f692136ab05e1f2345892fd0abb03dc)was developed to automate the cloud and cloud shadow detection process for PlanetScope imagery.
The cloud and cloud shadow detection algorithm was implemented for all 313 PlanetScope mosaics.Visual inspection of classification results revealed a reasonable detection for most images.Based on the 11,891 independent validation pixels, we estimated a true positive rate of 71.2% (actual cloud/shadows correctly identified) and a true negative rate (actual clear pixels correctly identified) of 96.6%.We further eliminated pixels within a three-pixel buffer around the identified clouds and cloud shadows, which increased true positive rate to 72.9% and slightly decreased true negative rate to 96.3%.Fig. A2 illustrates examples of the classification results.
The NDVI was then calculated for each PlanetScope mosaic.The NDVI was selected because 1) it is commonly used in phenology studies, and 2) it has a good dynamic range in low biomass areas and in such conditions does not suffer from saturation (Huete et al., 2002).We further refer to the PlanetScope-derived NDVI time series as NDVI P (the subscript P referring to PlanetScope).

Sentinel-2 imagery
We focused the phenology analysis with Sentinel-2 only on the period September 2017 to February 2019, given that before June 2017 no Sentinel-2B images were available for Kapiti and image density was poor for LR 2017 as a result.We used a total of 89 Sentinel-2 surface reflectance Level-2A products for tile code 37MBU that were atmospherically corrected using the Sen2Cor processor (version 2.5.5).The scene classification file was used to mask pixels classified as cloud shadow, cloud, and thin cirrus.In addition, pixels with blue band reflectance below 0.01 were also filtered out, because they were found to predominantly relate to undetected cloud shadows.We then applied a three-pixel buffer around the masked pixels to reduce poor quality observations near cloud and cloud shadow edges.We calculated NDVI from the 10 m resolution spectral bands 4 and 8.The Sentinel-2-derived NDVI time series are referred to as NDVI S in this manuscript (the subscript S referring to Sentinel-2).

Field camera time series
The first data set used to evaluate the phenology retrievals from PlanetScope and Sentinel-2 was obtained through digital repeat photography.To monitor greenness changes of vegetation on a daily basis, three Bushnell Trophy Cam Essential (model 119,736) trail cameras were installed in Kapiti Farm (Fig. 1) in October 2017.We purposively selected relatively homogenous locations of the main vegetation communities to enable subsequent comparison with satellite image series, including those with coarser resolution like MODIS.Fig. 3 shows sample photos captured by each camera, as well as their footprints.camera B of an Acacia drepanolobium shrubland with herbaceous vegetation, and camera C images a woodland containing ~5 m tall Balanites trees and herbaceous vegetation.Each of these three cameras is set up to take one RGB photograph (JPEG format) every 30 min from 8:00 to 17:30 Eastern Africa Time (EAT).The white balance and exposure is determined automatically and cannot be adjusted for these cameras.
Photographs were collected between October 2017 and February 2019, corresponding to two SR seasons and one LR season.The photographs that were blurred, over-or underexposed, or for which animal presence prevented the view of the vegetation were visually identified and discarded from the database.Subsequently, we calculated the GCC (Gillespie et al., 1987) for each retained photograph as follows: The GCC accounts for the influence of scene illumination on brightness levels (Sonnentag et al., 2012;Woebbecke et al., 1995) and is often used in phenology studies (Migliavacca et al., 2011;Richardson et al., 2007).
GCC values were averaged within the regions of interest (ROI; Fig. 3) to generate representative time series of greenness changes.Each ROI represents a dominant vegetation community.ROI1 of camera A (A1) represents open grassland; B1 represents shrub, and B2 grasses; and C1 represents tree canopy and C2 grasses.To relate the GCC to the satellite phenology retrievals, we additionally identified the ROI B0 for camera B that comprises a combination of the grass and shrub vegetation communities within the field of view.Although for camera C we also have multiple vegetation communities, only the tree canopy and the grasses not hidden by the tree canopy (i.e. the open grasses excluding the understory) are visible to a nadir-viewing satellite.For that reason, we estimated the fraction of trees and open grasses from the WorldView-2 image of February 2017, resulting in approximately 50% tree cover and 50% grass cover in the area within the field of view of camera C. For making this estimate we considered a larger area than the selected Sentinel-2 pixel (Fig. 3c), given possible geolocation uncertainty of the imagery.We consequently calculated the average GCC from C1 and C2 for comparison with the satellite NDVI series.To further reduce the influence of different scene illumination on the ROIaveraged GCC, the 90th percentile of all GCC values within a nonoverlapping three-day window (GCC90) was extracted and assigned to the centre day following Sonnentag et al. (2012).

MODIS imagery
The second dataset to which we compared our fine-resolution phenology retrievals consisted of MODIS imagery, given its widespread use in phenology studies.MOD13Q1 and MYD13Q1 Version 6250 m resolution vegetation index products (Didan, 2015a(Didan, , 2015b) ) from March 2017 to February 2019 were accessed through Google Earth Engine (GEE).These products are generated using a maximum value compositing technique with constrained viewing angle for acquisitions from the Terra and Aqua satellite, respectively.Each product includes a 16-day composite NDVI and Enhanced Vegetation Index (EVI).Due to the eight-day shift between MOD13Q1 and MYD13Q1, their combination results in an average eight-day interval for the combined VI time series.Rather than using daily MODIS acquisitions, we selected these products to have a high-quality dataset, which is used in various phenology studies (e.g.Hmimina et al., 2013;Liu et al., 2015;Peng et al., 2017;Vrieling et al., 2017).We used the pixel-specific NDVI and acquisition date, and masked out pixels that were flagged as poorquality observations due to the presence of clouds, cloud shadows, or large viewing angles in the quality reliability layer.The MODIS-derived NDVI time series are named NDVI M in the following (the subscript M referring to MODIS).

Phenology retrieval from vegetation index time series
Many approaches exist to retrieve phenology from vegetation index time series (de Beurs and Henebry, 2010).Despite recent advances in devising robust approaches for estimating phenology from 10 to 30 m resolution data (e.g., Bolton et al., 2020;Jönsson et al., 2018), there is no consensus on a single best approach for all cases, as this depends on study area, targeted phenological metrics, the shape of the VI series, and the presence of noise (White et al., 2009;Zeng et al., 2020).To reduce the impact of noise in the series, the fitting of mathematical functions to VI series has often been used.This curve fitting can be applied iteratively to use the fact that larger VI values are more reliable (Chen et al., 2004;Zeng et al., 2020).While acknowledging the possibility of applying alternative approaches, here we chose to consistently apply a common curve fitting approach to different data sources with the aim of assessing if data quality and frequency is sufficient to obtain accurate and spatially-consistent estimates of phenological metrics.For this purpose, we selected the double hyperbolic tangent function (Meroni et al., 2014b).A hyperbolic tangent function is a rescaled version of (and thus equivalent to) the logistic function, commonly used in phenology studies (e.g., White et al., 2014;Zhang et al., 2003).Contrary to the six-parameter formulations of double logistic models (e.g., Beck et al., 2006;Jönsson et al., 2018), we used an additional parameter to account for different minimum VI values before and after the season (e.g., which happens in our area due to different lengths of the dry season).We fitted the double hyperbolic tangent function to the VI time series generated from the satellite imagery (NDVI P , NDVI S , and NDVI M ) and from the camera photographs (GCC90).Based on seasonspecific VI series we found that this function could accurately describe the temporal VI variability.The double hyperbolic tangent function can be written as: (2 where t is time (days).The parameters were initialized as in Vrieling et al. (2018) and include: -a 0 : the baseline minimum VI value; -a 1 (a 4 ): the difference between the maximum and minimum VI value in the green-up (senescence) phase; -a 2 (a 5 ): the inflection point of the green-up (senescence) phase; -a 3 (a 6 ): controls the steepness of the green-up (senescence) phases.
The curve fitting was executed separately for each season through a Python script (https://gist.github.com/YanCheng-go/d4e17831f294-199443d0f7682558e608).The Levenberg-Marquardt algorithm (Moré, 1978), as implemented in the lmfit Python package (Newville et al., 2014), was used to find the optimum model parameters that minimize the sum of squared residuals between fitted values and actual values.Because lower VI values, relative to preceding or subsequent observations, can be due to remaining contamination of atmospheric effects or undetected clouds and shadows, the fitted curves were adapted to the upper envelope based on an iterative weighing method as used by Meroni et al. (2014b) and Vrieling et al. (2018).The method follows the weighing idea of Sellers et al. (1994), which was later developed into an iterative scheme by Chen et al. (2004).In short, if observations in the original VI time series are smaller than the corresponding fitted values, their weights will be decreased in the next iteration of the curve fitting.Preliminary tests indicated that for NDVI time series from PlanetScopeand Sentinel-2 the fitting converged in less than five iterations for more than 90% of the cases.To avoid infinite iterations, a maximum of 10 iterations was set.
The March 2017 to February 2019 timeframe of this study spans four vegetation seasons.We fitted each season separately, whereby the season included the dry period before and after the vegetation cycle.In this way, the same observations belonging to a dry season that separates the two cycles are used twice; once as the end period of the earlier season, and once as the initial period of the later one.We acknowledge that this may result in fitted curves that are not overlapping precisely during the shared dry season.Nonetheless, to keep a simple and straightforward approach (as implemented also in e.g., Vrieling et al., 2018), we preferred to optimally tune our models to the VI variability for each season separately.To achieve this, the timeframe was divided into four parts with one-month overlapping periods between successive vegetation seasons: 1 March 2017-1 October 2017 (LR2017), 1 September 2017-1 March 2018 (SR2017), 1 February 2018-1 October 2018 (LR2018), 1 September 2018-1 March 2019 (SR2018).We decided on a different starting point for the two LR seasons, given that LR2017 started in April, i.e. one month later than for LR2018 (Fig. 2), which was reflected in the VI series.For field camera series, the curve fitting was applied for both landscape-level ROI and community-level ROIs in the field of view of the three cameras (Fig. 3).For satellite imagery, the curve fitting was applied for all pixels contained within Kapiti Farm.
After fitting curves for satellite-and camera-derived VI time series (i.e., GCC90 and NDVI), two temporal metrics, the start of season (SOS) and end of season (EOS) dates, were estimated from the fitted curves by taking the 50% amplitude threshold, which roughly indicates the time of fastest green-up or senescence (White et al., 1997).The amplitude in the green-up and decay phases is computed as the difference between the maximum of the fitted function minus the minimum fitted initial and final values, respectively.For clarity, these minima are the fitted values within the time frames considered (e.g., for LR2017 the fitted values for 1 March 2017 and 1 October 2017) rather than the the parameter a 0 of Eq. (2).
Apart from SOS and EOS, two NDVI-related phenological metrics were retrieved from satellite-derived NDVI time series.These are maxNDVI, which relates to maximum biomass levels and cumNDVI, which is often used as a proxy for seasonal productivity (Bailey et al., 2004;Heumann et al., 2007).This results in the following four metrics: -SOS 50 : the start of the season; the date when VI first reaches 50% of the difference between the modelled maximum and minimum VI in green-up phase; -EOS 50 : the end of the season; the date when VI first reaches 50% of the difference between the modelled maximum and minimum VI in the senescence phase; -maxNDVI: the modelled maximum NDVI; -cumNDVI 50 : the cumulative value of the modelled NDVI between SOS 50 and EOS 50 .

Analysis of results
Based on the phenological metrics estimated from NDVI series of Sentinel-2, PlanetScope, MODIS, and camera-derived GCC90 time series, we performed a number of analyses to 1) evaluate the relationship between satellite-with camera-derived phenological metrics; 2) visually determine if the resulting spatial patterns are realistic and/or present artefacts; and 3) compare phenological metrics derived from Sentinel-2, PlanetScope and MODIS.With artefacts, we here refer to spatial patterns that are induced by factors that do not relate to the land surface dynamics, but rather to variability in image availability (e.g.due to cloud cover).Different statistical measures were used for these analyses, including: 1) Root Mean Squared Deviation (RMSD) for quantifying the difference across datasets, 2) Mean Signed Deviation (MSD) for assessing the bias, and 3) coefficient of determination (R 2 ) or Pearson's correlation coefficient (r) for evaluating the relationship between two datasets.
For the purpose of comparing satellite-with camera-derived phenological metrics, NDVI P series were aggregated to 10 m resolution to more closely match the field of view of cameras, to mitigate geometric errors, and to be compatible with Sentinel-2 series.This aggregation was achieved by taking the average NDVI P of all 3 m resolution pixels for which the centre is located within the geometry of a 10 m Sentinel-2 grid cell.The SOS 50 and EOS 50 retrieved from the aggregated PlanetScope, the Sentinel-2, and the MODIS NDVI time series for the pixels corresponding to the three camera locations were then compared against the SOS 50 and EOS 50 from the camera-derived GCC90 time series.We only used the landscape-level camera ROIs for this comparison, given that the satellite measures the combined reflectance from multiple vegetation communities within a pixel.The RMSD and MSD between the satellite-derived and camera-derived phenological metric were used to quantify their difference.
We visually analyzed the resulting phenology maps to assess a) if the retrieved dates are sensible; b) if the spatial patterns match expectations; and c) if there are spatial artefacts that are indicative of lacking or erroneous observations, hence resulting in improper fits and retrieved phenological parameters.For clarity, this visual analysis was not intended to eliminate data points, but rather to evaluate the effectiveness of the phenology retrieval for the different satellite sources across the spatial extent of the study area.
To quantitatively compare phenological metrics between PlanetScope, Sentinel-2, and MODIS, we resampled PlanetScope-derived phenology maps to 10 m and 250 m resolution by averaging the values of all pixels whose centre falls within the Sentinel-2 and MODIS grid cells.Subsequently, all seasons were combined by calculating for each phenological metric the deviation from the season-specific mean, as obtained from the PlanetScope retrievals.After that, density scatterplots were generated and the R 2 , RMSD, and MSD were calculated.Due to the anomaly calculation, the plots and statistics illustrate how well the spatial variability of the metrics compare over the different seasons.We note that if instead we would plot the non-normalized metrics, the resulting R 2 would be significantly higher due to the strong interseasonal variability of phenology and NDVI levels.
We also evaluated to what extent the fitted curves correspond with the original data series.For this purpose, r, RMSD, and MSD were calculated for each fitted VI time series as compared to the original VI time series (Vrieling et al., 2017).Lower values of r, and higher values of RMSD or MSD do not necessarily imply a worse fitting performance because the iterative fitting may account for erroneously low VI values by adapting to the upper envelope.Nevertheless, these statistics can provide a combined quantification of model fit and noise level differences between data sources.

Phenological retrievals: camera versus satellite
Fig. 4 compares the GCC90 time series and phenological metrics of different vegetation communities in the field of view of camera A, B, and C for SR2017, LR2018, and SR2018.The most striking phenological difference between vegetation communities was for EOS 50 , which was substantially later for shrubs and trees as compared to grasses for all three seasons considered.On average, after the grass entered the senescence phase, the shrubs and trees maintained the green leaf canopy for about one more month.The time lag was larger with increasing total seasonal precipitation.The SOS 50 of grass was similar to shrubs while earlier than trees in two out of the three seasons (LR2018 and SR2018, while for SR2017 SOS 50 was similar).Overall, grass had a shorter growing season and responded more rapidly (i.e.faster increase) to precipitation than trees.It is worth noting that the trees and grass in the field of view of camera C had a short secondary green-up phase after the main green-up phase in SR2017, caused by an approximate 43-day dry period in late November and December 2017 followed with substantial rainfall on 2 and 3 January 2018, and in LR2018 caused by a small rainfall event on 31 July 2018 (Fig. 4).However, this secondary greenup is less apparent for the ROIs of camera A and B.Moreover, the speed of green-up for the grass in the field of view of camera C is faster than the grass in the field of view of cameras A and B. Possibly also the terrain position plays a role here (Fig. 1), given that the camera C location may be receiving larger amounts of runoff from upslope areas and potentially more rain due to a small orographic effect.In short, the camera-based analysis indicates that vegetation phenology in Kapiti Farm is spatially and temporally heterogeneous, likely as a function of the complex composition of vegetation, rainfall variability, and terrain characteristics.
Fig. 5 compares satellite-derived NDVI time series with cameraderived GCC90 time series for the three camera locations (using ROIs A1, B0, and C0).The seasonal signal is clear in both satellite-and camera-derived vegetation index time series.Because the same MODIS pixel covers the field of view of both camera A and B, the MODIS-derived NDVI are the same in Fig. 5a and b.The PlanetScope-derived NDVI time series show noisier temporal patterns as compared to Sentinel-2 and MODIS-derived NDVI.Nevertheless, the combined use of all PlanetScope acquisitions resulted in a denser NDVI time series as those derived from MODIS and Sentinel-2.Table 1 compares satellite-with camera-derived SOS 50 and EOS 50 ; it shows that MODIS series resulted in SOS 50 retrievals that were more similar to the camera-based SOS 50 retrievals as compared to those from PlanetScope and Sentinel-2 series.The largest discrepancy between PlanetScope-and camera-derived SOS 50 occurred for SR2017 (Fig. 5b).This may be attributed to the unrealistically low NDVI P values for 16 and 24 November 2017, which based on visual examination seem to result from small undetected thin clouds in the PlanetScope imagery for the area that is compared with the camera B field of view.For EOS 50 , retrievals from PlanetScope series were most similar to camera-based retrievals.Nonetheless, a large deviation between EOS 50 from camera, PlanetScope, and Sentinel-2 existed for SR2017 (Fig. 5c); this is because the secondary green-up for that season is contained in the PlanetScope time series, but not captured by Sentinel-2 due to a lack of observations in January.The RMSD for EOS 50 retrievals is approximately twice as large as for SOS 50 retrievals.According to the MSD statistics in Table 1, on average satellite-based SOS 50 and EOS 50 retrievals tend to be later than those from the camera series.

Mapping of phenological metrics
Pixel-level phenology was retrieved over the Kapiti Farm.Table 2 compares the image availability and the fit statistics for PlanetScope, Sentinel-2 and MODIS.For all seasons, the average number of cloudfree observations is almost double for PlanetScope as compared to Sentinel-2 and MODIS.We note however that in case of MODIS the smaller number is due also to the 16-day VI products used (Section 2.5).The maximum gap (in days) between two sequential cloud-free observations in SR2017 is 50% shorter for PlanetScope as compared to Sentinel-2 and MODIS.While maxGap provides information about the maximum lack of observations in a certain period, minGap is not reported here because it would merely reflect the satellites' revisit time.For the dry season, it is relatively easy to find two consecutive cloudfree images, resulting in a minGap of 5 for Sentinel-2 and 1 for Pla-netScope.For each season, Sentinel-2 was the dataset with the largest maximum gap.LR2018 was the season with the largest maxGap (> 30 days for all sensors), which is due to the persistent cloud cover in April and May 2018.In principle, a larger number of valid observations and shorter temporal gaps ensure a better description of the actual vegetation dynamics.In this light, the smaller r value for PlanetScope (Table 2) may be partially explained by more noise, but also by the effective capturing of more natural temporal variation that is not properly modelled by a single season (for example the above-mentioned secondary green-up in SR2017, see Figs. 5 and 6).For Sentinel-2 and MODIS, this second peak was not apparent, resulting in similar r values for all seasons.
Fig. 6 shows the maps of PlanetScope-derived SOS 50 , EOS 50 , maxNDVI, and cumNDVI 50 for the four seasons considered.Spatial patterns as observed on high-resolution WorldView-2 imagery and a DTM (Fig. 1) are reflected in the various maps.The different phenological metrics derived from PlanetScope series are generally spatially smooth (i.e. they do not show unrealistic abrupt spatial variations over homogenous areas).It is worth noting that spatial patterns are different across phenological metrics and seasons.For example, the EOS 50 for LR2017 is relatively early in the southeast, while for SR2017 the EOS 50 is significantly later in the centre.Not all spatial differences depicted in the maps can be attributed to real phenological differences.This is the case for the EOS 50 map for LR2018 (Fig. 6j), where a linear separation is apparent, with later EOS 50 dates west of this line.The main cause was that a cloud-free image acquired on 17 July 2017 only covered the northwestern part of Kapiti.When repeating the phenological analysis without that image, the EOS 50 difference between two neighbouring pixels decreased from 10 to 2 days (Fig. 7ab, points 1 and 2).This demonstrates that despite the many valid observations in NDVI P time series, a single observation can have an important influence on the fitting and resulting phenological metrics.This artefact is not reflected on the maxNDVI and cumNDVI maps for LR2018, suggesting that those metrics are less sensitive to (the lack of) individual observations.Fig. 8 shows the maps of Sentinel-2-and MODIS-derived SOS 50 .Apparent artefacts exist in Sentinel-2-derived SOS 50 maps for LR2018 and SR2018 (Fig. 8b-c).Unlike the straight-line separation in the Pla-netScope-derived map, these artefacts have irregular shapes.An explanation for these artefacts could be the lack of observations from March until the beginning of April in 2018.Fig. 7c-d shows that a single observation in that period, not available for the entire study area due to clouds, strongly influences the fitting.Despite these artefacts, the coarse spatial patterns of the Sentinel-2 and MODIS-derived SOS 50 maps agree with those obtained from PlanetScope.For example, both MODIS and PlanetScope indicate a delayed SOS 50 for LR2017 in the northwestern part of Kapiti.For SR2017, the relatively late SOS 50 in southeastern Kapiti is consistently shown across PlanetScope-, Sentinel-2-and MODIS-derived maps.The MODIS-derived SOS 50 shows less detailed spatial patterns and on average occurs for all seasons earlier than their PlanetScope-and Sentinel-2 equivalents (see reported mean in Fig. 8).Visual assessment of the phenology results for small areas reveals that PlanetScope effectively captures the spatial variation caused by various vegetation covers (e.g.individual trees) better than Sentinel-2, as can be observed from comparison with the 0.5 m resolution WorldView-2 image (Fig. 9).In this example the single tree (Vachellia tortilis) shows an earlier SOS 50 as compared to the grass surrounding it.We can recognize this similar spatial pattern for Sentinel-2, but the tree phenology is less apparent due to the coarser resolution and consequent mixed spectral response; nonetheless the patterns show that despite reported positional errors, PlanetScope and Sentinel-2 retrievals are reasonably aligned.The earlier SOS 50 for this tree is different from what is observed elsewhere, for example in the area of camera C, where both camera (Fig. 4c) and PlanetScope retrievals (Fig. 10) detected a tree cover green-up that is later than for grass.
Fig. 11 compares anomalies of the phenological metrics from Sentinel-2 and MODIS to those from PlanetScope at the 10 × 10 m and 250 × 250 m grid cell level, respectively.For Sentinel-2, and to a lesser extent for MODIS, a stronger correlation is found for maxNDVI and cumNDVI 50 as compared to SOS 50 and EOS 50 .The negative MSD of SOS 50 means that on average SOS 50 from Sentinel-2 was earlier than SOS 50 derived from PlanetScope.This is especially obvious for MODISderived SOS 50 (MSD = −5.8).For EOS 50 , Sentinel-2 shows a positive MSD, which could be partially caused by the large area of artefacts in the PlanetScope-derived EOS 50 map (Fig. 7a).The RMSD values for SOS 50 and EOS 50 indicate that on average Sentinel-2-derived SOS 50 and EOS 50 had 8-and 10-day difference from PlanetScope-derived ones.The difference between PlanetScope and MODIS was 9 days for SOS 50 and 11 days for EOS 50 .The positive MSD values for maxNDVI and cumNDVI 50 indicate that on average NDVI S and NDVI M had larger values than NDVI P series.For all metrics except SOS 50 , PlanetScope showed a stronger correlation with Sentinel-2 as compared to MODIS.This correlation was particularly strong for maxNDVI and cumNDVI 50 , with R 2 above 0.5.

Discussion
Our results showed that fine-resolution satellite systems are effective in retrieving detailed spatial patterns of vegetation phenology for semi-arid rangelands with short vegetation cycles.A curve fitting approach applied to PlanetScope and Sentinel-2 series resulted in similar spatial patterns.The camera series showed that herbaceous vegetation responds rapidly to rainfall and moisture deficits, at times resulting in multiple greenness peaks within one season.Such rapid temporal variability could to some extent be captured by the frequent PlanetScope observations, but not by Sentinel-2.As compared to camera-based in-situ measurements, satellite-derived SOS 50 and EOS 50 were on average within 9 and 18 days of camera-derived SOS 50 and EOS 50 , respectively (Table 1).Other studies also revealed differences between in-situ observations and satellite-retrieved phenological metrics of 10 to 30 days (White et al., 2014).For digital repeat photography, these differences can be partially explained by the use of different vegetation indices (Bolton et al., 2020;Liu et al., 2017;Vrieling et al., 2018).Unlike GCC that only uses the visible spectrum, NDVI uses NIR reflectance, which is sensitive to changes in vegetation structure and leaf area index.In addition, different viewing angles contribute to the inconsistency between satellite-and camera-derived phenological metrics.Vegetation greenness observed from an oblique camera view differs from greenness as observed by the satellites' more nadir view (Bolton et al., 2020) given that non-photosynthetic elements like stems and grass heads are more dominant in the oblique view (Vrieling et al., 2018).In fact, radiative transfer model simulations offered a physicallysound explanation for the earlier EOS 50 observed from oblique versus nadir viewing angles (Vrieling et al., 2018).Despite these inconsistencies of observation characteristics, the satellite retrievals effectively captured the main seasonal changes as observed by the cameras.
Phenology retrieval results over the Kapiti Farm showed that finescale spatial variability in phenological timings exists.Identifying the precise causes for this variability was beyond the scope of this study, but they certainly include terrain characteristics, vegetation composition, water availability, and herbivory.Terrain characteristics (Fig. 1) affect both vegetation composition and water availability, for example due to runoff processes.With regard to vegetation composition we found that the timing of green-up of woody plants versus herbaceous vegetation is not consistent between locations (Figs. 4,9,and 10).This confirms findings by Whitecross et al. (2017) who highlighted that tree and grass phenologies are variable for African savannahs.In addition, we found that spatial patterns were different for the various phenological metrics and between seasons.Because water availability is the main factor influencing vegetation growth in semi-arid areas (Scholes and Walker, 1993), a large part of the spatial and temporal differences in phenology are likely a result of variable plant water availability and intra-and inter-annual variability of precipitation.Given that rainfall variability can be important even at small spatial scales (e.g., Fischer et al., 2016), setting up additional rain gauges across the study area  could help to better understand if precipitation is an important driver of spatial variability of phenology.Moreover, precipitation is not the sole determinant of water availability, but soil properties, topography, and water logging could be further explanatory factors for phenological variability.Lastly, herbivory through grazing by livestock and wildlife could also affect the variability, as it impacts species composition, vegetation structure, and biomass quantity and quality (Altesor et al., 2005;Hiernaux, 1998;Jansen et al., 2016;Metzger et al., 2005;Olsen et al., 2015).However, Olsen et al. (2015) found that grazing-induced differences in vegetation life spans and biomass can hardly be reflected in MODIS-based NDVI time series, possibly due to the coarse resolution.While the simultaneous grazing by livestock and wildlife in Kapiti makes the collection of grazing data cumbersome, it would be interesting to acquire detailed spatio-temporal data on grazing intensity and assess if these may explain part of the phenology variability found in the PlanetScope and Sentinel-2 derived phenology estimates.
The comparison between fine-resolution retrievals of phenology (i.e. from PlanetScope and Sentinel-2) and the more traditional coarseresolution MODIS retrievals confirmed prior findings on spatial scaling of phenology, and illustrated the benefit of phenological retrievals at finer spatial scales.Zhang et al. (2017) found that coarse-resolution phenology is controlled more by the earlier SOS at finer resolution, resulting in a better match between coarse-resolution SOS and the 30th percentile of fine-resolution SOS as compared to its spatial average, particularly in heterogeneous environments.Consequently, we may expect earlier spatial average SOS 50 values from coarse-resolution data, which is confirmed for MODIS by this study (Figs. 8 and 11) and also by Vrieling et al. (2017).Although PlanetScope's spatial resolution is also finer than the one of Sentinel-2, SOS 50 retrievals from Sentinel-2 were not significantly earlier, probably partially because the resolutions are more similar, and additionally because of the effect of spatial artefacts.Further efforts could be made to evaluate scaling effects between Pla-netScope and Sentinel-2 phenology by only selecting pixels/seasons with high-quality retrievals, as was done by Zhang et al. (2017).We nonetheless showed that spatial patterns of the various phenology metrics displayed a reasonable agreement between the retrievals from the various sensors (Figs. 6,8,A3).
A clear advantage of the fine resolution phenology is that it may be linked more directly to specific vegetation communities and in-situ phenological measurements (e.g.visual observations, field cameras, and flux tower data), as was recently also shown for other spatial contexts (Dronova et al., 2020;Granero-Belinchon et al., 2020;Pastick et al., 2020).For example, PlanetScope retrievals (and also those from Sentinel-2, although less clearly) showed that phenological metrics for a single tree differed substantially from its surroundings (Fig. 9).This could help in gaining a more detailed understanding of ecosystem structure and function at the local scale for semi-arid rangelands.For example, the finer-resolution mapping of cumNDVI as a proxy measurement of gross primary productivity may provide a more accurate estimate of rangeland productivity and carbon stock (Schwieder et al., 2018), which in turn may allow a better understanding of its drivers.
Maps of phenological metrics derived from Sentinel-2, and to a lesser extent from PlanetScope, displayed spatial artefacts due to noisy or missing data, which implies that at least for some locations the phenology estimates are inaccurate.These artefacts and the resulting uncertainty in the estimates are larger for SOS 50 and EOS 50 than for maxNDVI and cumNDVI (see also Figs.S1 and S2 in the Supplementary Data for simulation tests on data elimination).This results from the fact that retrievals of timing metrics (SOS 50 and EOS 50 ) are more sensitive to the shape of the fitted curves as compared to maxNDVI and cumNDVI (Myers et al., 2019;Vrieling et al., 2017).Hence, depending on the phenological metric required in a specific application, the presence of artefacts can be more or less critical.For example, despite inaccuracies in SOS and EOS, the spatio-temporal assessment of rangeland productivity based on cumNDVI (Vrieling et al., 2016) could be less problematic, and estimates of its spatial variability or interseasonal variation may notwithstanding be of acceptable quality.
Due to more frequent observations, PlanetScope-derived phenology maps showed fewer artefacts due to clouds and slightly more robust simulation results as compared to Sentinel-2 (Figs.S1 and S2).This higher frequency allows for more frequent observations particularly in periods of more intense cloud cover, and for capturing irregular short vegetation cycles that result from precipitation variability in semi-arid regions.However, PlanetScope surface reflectance series are noisier than those from Sentinel-2 and MODIS.This larger noisiness is a result of two factors: 1) PlanetScope imagery is acquired by many satellites with varying illumination geometry (i.e., sun azimuth and elevation), relatively low radiometric quality, and poor inter-sensor calibration, resulting in less stable surface reflectance retrievals.One option to enhance the quality could be through cross-calibration with other, more stable, satellite sensors like Sentinel-2 and MODIS (Houborg and McCabe, 2018;Wang et al., 2020); 2) our cloud and shadow screening approach does not result in masks of similar quality as available for Sentinel-2 and MODIS, which both have more dedicated spectral bands.A new "usable data mask" (UDM2) was announced by Planet Labs in April 2019 (Planet Labs Inc, 2019), which could improve our present implementation, but was not available for the full study period, nor is the approach used well documented.Alternative deep learning approaches for cloud, shadow, and haze discrimination have recently been developed (Shendryk et al., 2019) and could benefit phenological retrieval (e.g., White et al., 2014).
Nonetheless, even without these potential further improvements, we found that the dense PlanetScope series were effective in retrieving fine-scale phenological patterns.Because currently free access to PlanetScope data is limited to monthly quota-bound academic licenses, and phenological analysis require many images in time, scaling this analysis to larger areas may at present be prohibitive.Therefore, the integration of multiple non-commercial data sources (e.g., the Harmonized Landsat and Sentinel-2 dataset: Claverie et al., 2018) should be further explored.In addition, synergistic use of PlanetScope and Sentinel-2 series may allow to attain more robust fine-resolution estimates of phenology, which requires further research.Precise coregistration through image registration would benefit such synergistic use to ascertain a good spatial fit between data sources across the spatial domain of analysis (Behling et al., 2014;Stumpf et al., 2018).
The upper envelope function fitting approach, used in this paper, was shown to be most sensitive to removing observations in the (short) green-up and senescence phases (Fig. S2), which is in line with the results of similar sensitivity analyses conducted by Vrieling et al. (2018) and Zhang et al. (2009).To mitigate the impact of reduced observations on the estimation of phenology, more efforts can be made to improve the robustness of the curve fitting method used in this study.Recent efforts to achieve this (Bolton et al., 2020;Jönsson et al., 2018) were discussed in the introduction of this paper, and involve the incorporation of data from alternative years to help describe the expected vegetation dynamics, in case the specific year has limited observations.It remains an open question whether such an approach would be effective in systems with short vegetation cycles and large interannual variability in timing and magnitude of greenness as found in Kapiti Farm.A further improvement to our curve fitting approach would be to avoid the manual definition of the seasonal timeframe for which a curve should be fitted.This definition was currently based on visual observation of the time series and prior knowledge.Such a user definition of expected phenology is also part of existing phenology software packages like TIMESAT (Jönsson and Eklundh, 2004), when creating spatial output.More automatic approaches could set breakpoint based on a multi-year VI climatology (Meroni et al., 2014b) or cycle search techniques that first retrieve local maxima and minima from the series (e.g., Bolton et al., 2020;Meroni et al., 2020;Vrieling et al., 2011).

Conclusions
This study demonstrated the potential of using PlanetScope and Sentinel-2 image series to retrieve vegetation phenology at fine (3-10 m) spatial resolution for heterogeneous landscapes with short vegetation seasons.Satellite-derived SOS 50 was on average within nine days and EOS 50 within 18 days of their camera-derived equivalents.Due to its shorter revisit time, PlanetScope had a better density of cloud-free observations as compared to Sentinel-2 NDVI time series.For that reason, phenology retrievals from PlanetScope were less affected by persistent cloud cover during rainy seasons, resulting in more consistent spatial patterns and fewer spatial artefacts as compared to Sentinel-2.Nonetheless, opportunities exist to improve the PlanetScope NDVI stability, and consequently its derived phenological estimates, among others through improved cloud masking and cross-calibration with more stable sensors like Sentinel-2.Our study does not intend to be conclusive about Sentinel-2's potential for studying phenology of short vegetation cycles, but shows that for tropical rangelands its revisit interval cannot ascertain sufficient cloud-free observations across space for every season.While cost may still be prohibitive for scaling phenological analysis with PlanetScope to larger regions, we demonstrated its ability to estimate the phenology of vegetation communities and even individual scattered trees.This may assist in improved understanding of ecosystem functioning at the local scale for heterogeneous landscapes like semi-arid rangelands.

Fig. 1 .
Fig. 1.Overview of the study area.The left map shows the boundary of Kapiti Farm and the location of the three field cameras.The centre map shows a 5 mresolution Digital Terrain Model, which was generated from data acquired with a Leica ALS60 aerial LIDAR survey.The background imagery for both maps is a WorldView-2 scene (natural colour; 0.5 m) acquired on 2 February 2017 provided by DigitalGlobe.The maps on the right show the location of Kapiti Farm in Kenya.
Cameras A and B are in an approximately flat area, and camera C is on a slope.The elevation of camera A, B, and C are 1632 m, 1636 m and 1723 m, respectively.Camera A takes images of an open grassland,

Fig. 2 .
Fig. 2. Precipitation in Kapiti Farm based on manually-collected in-situ daily rain gauge data.The upper panel shows daily and monthly rainfall from March 2017 to February 2018, and bottom panel for March 2018 to February 2019.The average monthly precipitation was calculated from data of the past 18 years (January 2001-February 2019).

Fig. 3 .
Fig. 3. (A)-(C) Sample photos taken by cameras A, B, and C at Kapiti Farm at 12:00 (EAT) on 19 November 2017 (top row), with corresponding locations placed on the same WorldView-2 scene as in Fig. 1 (bottom row).The red polygons in A-C indicate the region of interest (ROI) used for averaging GCC values.A1 refers to ROI1 in camera frame A. ROIs A1, B1, B2, C1, and C2 attempt to isolate different vegetation communities.B0 merges two communities in a landscape-level view to increase comparability with satellite observations.On the bottom row the yellow lines indicate the field of view of the camera, and the 10 m Sentinel-2 grid is shown in black.The blue boxes represent the selected Sentinel-2 pixel locations used for comparison of time series.

Fig. 4 .
Fig. 4. Time series of GCC90 for each vegetation community in the field of view of camera A, B, and C.Tthe lines are fitted curves and the horizontal bars near the xaxis indicate the duration of growing seasons from SOS 50 to EOS 50 , as estimated from the GCC90 series.The blue vertical lines in panel A indicate the daily precipitation from 1 September 2017 to 28 February 2019.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Time series of satellite-derived NDVI and camera-derived GCC90 at three camera locations (A, B, and C in Fig. 1) in SR2017, LR2018 and SR2018.GCC90 was extracted for the entire field of view of each camera location (A1, B0, and C0 in Fig. 2).The lines show the fitted curves and the horizontal bars near the x-axis indicate the duration of growing seasons from SOS 50 to EOS 50 , as estimated from each respective data source.The blue vertical lines in the uppermost panel indicate the daily precipitation.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .Fig. 7 .Fig. 9 .
Fig. 6.SOS 50 , EOS 50 , maxNDVI and cumNDVI 50 retrieved from PlanetScope-derived NDVI time series for Kapiti Farm for four seasons.The mean and standard deviation are reported beside each map.All maps are visualized by showing the difference from the mean.For SOS 50 /EOS 50 , red colours indicate that the date is before the mean, and blue after the mean.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11 .
Fig. 11.Density scatterplots for SOS 50 , EOS 50 , maxNDVI and cumNDVI 50 for the combined SR2017, LR2018, and SR2018 seasons.The x-and y-axis indicate the difference of each phenological metric as compared to the mean SOS 50 , EOS 50 , maxNDVI, and cumNDVI 50 of each considered season calculated from the PlanetScope series.The x-axis represents the PlanetScope-derived phenological metric, resampled to 10 or 250 m to match the spatial resolution of the data source indicated on the y-axis.The y-axis is the phenological metric derived from Sentinel-2 or MODIS.Blue to red colours indicate an increase in frequency.The black line in each plot is the linear regression model, and the dashed line the 1:1 relationship.RMSD, MSD and R 2 are reported in each plot.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. A3 .
Fig. A3.EOS 50 , maxNDVI and cumNDVI 50 retrieved from Sentinel-2-and MODIS-derived NDVI time series for Kapiti Farm for four seasons.The mean and standard deviation are reported beside each map.The maps for SOS 50 and EOS 50 are visualized by showing the difference from the mean.Red colours means that the date is before the mean, and blue colours mean that the date is after the mean.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Comparison of satellite-derived SOS 50 and EOS 50 with camera-derived ones.The RMSD and MSD were calculated by combining the results for three camera locations in three seasons.

Table 2
Image availability and fit statistics for PlanetScope, Sentinel-2 and MODIS; r is the correlation between fitted and original NDVI values.nImages refers to the average number of cloud-free observations in the NDVI time series.maxGap indicates the maximum difference (days) between two consecutive cloud-free observations in the NDVI time series.Mean and standard deviation are reported for each measurement.The statistics were averaged using all pixels within the extent of Kapiti Farm.