## 3 Methods
### Step 2: Time Series Analysis
#### 3.2.2 Land Surface Phenology for Vegetation Dynamics 2.0  
For the VD 2.0 analysis, additional processing was performed to generate annual time series devoid of seasonal effects that may arise from the static selection. With LGAC data migration and data restoration of imagery lacking payload-correction data, it is generally feasible to extract the timing of specific phenological key points from the full data record on a per-pixel basis, i.e., to derive land surface phenology (LSP) metrics directly [32], or to downscale coarse resolution LSP metrics using data fusion techniques [85]. In a first step, the quality bit product was used to remove all poor-quality observations, i.e., pixels flagged as cloud, cloud shadow, snow, as well as saturated or sub-zero reflectance values. Acknowledging that cloud, and especially cloud shadow masks, are not perfectly accurate, a time series-based outlier detection followed by an inlier restoration was performed based on the noise estimation proposed by [PERSON] et al. [86]. Given a triplet of successive observations, the absolute residual \(r\) between the central observation \(i\) and linear interpolation between its neighbors is given by:  
\[r_{i}=\left|\rho_{i}-\frac{\rho_{i+1}-\rho_{i-1}}{t_{i+1}-t_{i-1}}.(t_{i}-t_{i -1})-\rho_{i-1}\right|, \tag{4}\]  
where \(\rho\) is the reflectance of the shortest wavelength, i.e., blue, which commonly is most severely affected by remaining noise due to atmospheric effects; \(t\) is the time measured in continuous days. The noise \(\varepsilon\) of a time series with \(n\) observations is estimated as:  
\[\varepsilon=\sqrt{\sum_{i=2}^{n-1}r_{i}/n-2}. \tag{5}\]  
We iteratively removed the observation \(i\) with the largest residual \(r_{i,max}\), if \(r_{i,max}/\varepsilon\) exceeded a threshold. Subsequently, we restored all observations that were eliminated via \begin{table}
\begin{tabular}{c c c c c} \hline \hline
**Spectral Band** & **Photosynthetic Active Vegetation** & **Soil** & **Rock** & **Shade** \\ \hline Blue & 3.2 & 7.3 & 26.2 & 0.0 \\ Green & 5.6 & 14.5 & 31.0 & 0.0 \\ Red & 4.5 & 22.4 & 33.4 & 0.0 \\ NIR & 36.7 & 27.5 & 47.0 & 0.0 \\ SWIR1 & 17.0 & 40.2 & 72.4 & 0.0 \\ SWIR2 & 7.1 & 32.2 & 54.9 & 0.0 \\ \hline \hline \end{tabular}
\end{table} 
Table 3: Endmember set derived from [42]; values are in percentage reflectance. the quality bit product, if \(r/\varepsilon\) was below a threshold. We set these thresholds to 3 and 1, respectively, which were found to be reasonable values based on practical experience. Subsequently, the \(f\) time series were interpolated using an ensemble of Radial Basis Function convolution (i.e., smoothing) filters as described by [PERSON] et al. [87], where the observations in the kernel are weighted according to a Gaussian distribution:  
\[w_{t}=\frac{1}{\sigma\cdot\sqrt{2\pi}}\cdot\exp\left(-0.5\cdot\left(\frac{ \Delta t}{\sigma}\right)^{2}\right), \tag{6}\]  
with \(t\) being the time of the observation, \(\Delta t\) the observation time relative to the kernel center, and \(\sigma\) the full-width-at-half-maximum of the Gaussian bell-shaped curve, all values given in days. We used three kernel widths \(\sigma\) with 8, 16, and 32 days, as e.g., employed by [PERSON] et al. [88]. A kernel cutoff value that preserves 95% of the area under the Gaussian curve was used. The estimates from the three kernels were eventually aggregated using a weighted average, wherein the weights correspond to the data available within each kernel. This weighting strategy gives preference to kernels with a higher data density for the final value estimation. Due to the implemented kernel cutoff value, no-data values occasionally remained during the winter season when data gaps are substantially larger than the employed kernel widths, and as such, we subsequently performed a linear interpolation to strictly enforce equidistance and absence of data gaps.