## 3 Methods
### Step 1: Analysis-Ready Data (ARD)
#### 3.1.2 Dataset Reduction for Vegetation Dynamics 1.0  
Before the opening of the Landsat archive, VD 1.0 studies were limited to few Landsat images that needed to be handpicked carefully. [PERSON] et al. [49] manually selected 14 Landsat TM/ETM+ images from 1984-2005 based on cloud cover and seasonal constraints with a maximum gap of one year between acquisitions (Table 2). They aimed to pick images at peak vegetation vigor, which was assumed to approximately correspond to mid-May to mid-June in the mountainous rangelands of Crete. The seasonal validity of each selected image was compared with a 10-day AVHRR NDVI time series [80] for selected plots in these rangelands. To mimic this hand-picking of images, but for the whole  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Symbol** & **Description** & **Symbol** & **Description** \\ \hline \(\Delta D\) & Difference between selected & [x,y] & Cartesian coordinates \\  & observation and target date & [x,y] & Long-term average vector \\ \(\rho\) & Reflectance & [\(\langle r\rangle\),\(\langle f\rangle\)] & \\ \(\Delta\rho\) & Residual between observed and interpolated reflectance & [\(\langle x\rangle\),\(\langle y\rangle\)] & Long-term average of Cartesian coordinates \\ b & Spectral band & \(\langle\theta\rangle\) & Typical (long-term) start of the phenological year \\ F & Fractional cover & s, ns & Season, number of seasons \\ f & Shade-normalized fractional cover of & [\(\overline{r}\), \(\overline{f}\)] & Seasonal average vector \\ E & Photosynthetically active vegetation & [\(\overline{r}\), \(\overline{y}\)] & Seasonal average of Cartesian coordinates \\ t & Time & \(\overline{b}\), & Start of the phenological year in season \\ At & Time difference & y, my & Calendar year, number of years \\ n & Length of the time series & \(\mathrm{f}_{\mathrm{VFS,y}}\) & Annual peak of season fractional cover \\ \(\varepsilon\) & Noise of the time series & \(\mathrm{f}_{\mathrm{VFS,y}}\) & Annual fractional cover seasonal base level \\ w & Weight & \(\mathrm{f}_{\mathrm{VFS,y}}\) & Annual fractional cover seasonal amplitude \\ \(\sigma\) & Sigma of Gaussian bell & a, b & Regression intercept and slope \\ DOY & Day of the Year & CC & Relative cover change \\ [z,f] & Polar coordinates (Day of the Year in radians, fraction) & CC\({}_{\mathrm{act}}\) & Absolute net cover change \\ \hline \hline \end{tabular}
\end{table}
Table 1: Mathematical symbols.  
of Crete, we used a selection algorithm on the full ARD dataset. We generated pixel-based composites with a parametric weighting scheme (15; 22) with the FORCE Higher Level Processing System (v. 3.5) using the seasonal proxy as the only determinant, i.e., the closest observation to a target date was used to populate each pixel in the composite. The selected acquisition dates (Table 2) were used as temporal targets, and a narrow compositing period of \(\pm\)10 days was used to restrain the selection of observations to dates that are as close as possible to the target dates. For most pixels, the exact date was selected, with secondary peaks at +7, \(-\)1, and \(-\)9 days, which are related to Landsat's orbital paths. Note that these composites still contain some clouds as only the seasonal suitability was used as the selection criterion. This was done to ensure that the compositing algorithm is biased towards selecting full images as this matches the manual selection process from the reference studies. Still, clouds were masked out for the analysis using the automated cloud masks, whereas clouds were manually identified in the reference studies.