Machine learning and Global Vegetation: Random Forests for Downscaling and Gapfilling

. Drought is a devastating natural disaster, where water shortage often manifests itself in the health of vegetation. Unfortunately, it is difficult to obtain high-resolution vegetation drought impact, which is spatially and temporally consistent. While remotely sensed products can provide part of this information, they often suffer from data gaps and limitations in spatial or temporal resolutions. A persistent feature among remote sensing products is tradeoffs between spatial resolution and revisiting times, where high temporal resolution is met by coarse spatial resolution and vice verse . Machine learning methods 5 have been successfully applied in a wide range of remote sensing and hydrological studies. However, global applications to resolve drought impacts on vegetation dynamics still need to be made available, while there is significant potential for such a product to aid improved drought impact monitoring. To this end, this study predicted global vegetation dynamics based on the Enhanced Vegetation Index ( evi ) and the popular Random Forest algorithm (RF) at 0.1 ◦ . We assessed the applicability of RF as a gap filling and downscaling tool to generate spatial and temporal consistent global evi estimates. To do this, we trained 10 an RF regressor with 0.1 ◦ evi data using a host of features indicative of water and energy balances experienced by vegetation and we evaluated the performance of this new product. Next, to test whether the RF is robust in terms of spatial resolution, we downscale global evi , the model trained on 0.1 ◦ data is used to predict evi at 0.01 ◦ resolution. The results show that the RF can capture global evi dynamics at both the 0.1 ◦ (RMSE: 0.02 - 0.4) and at the finer 0.01 ◦ (RMSE: 0.04 - 0.6) resolution. Overall errors were higher in the down-scaled 0.01 ◦ compared to the 0.1 ◦ product. Yet, relative increases remained small, 15 thus demonstrating that RF can be used to create downscaled and temporally consistent evi products. Additional error analysis reveals that errors vary spatiotemporally, with underrepresented landcover types and periods of extreme vegetation conditions having the highest errors. Finally, this model is used to produce global spatially continuous evi products at both the 0.1 ◦ and 0.01 ◦ spatial resolution for 2003-2013 at an 8-day frequency.


Introduction
Drought is one of the most disruptive natural hazards, causing negative repercussions on the economy, society and the environment, which can affect large areas and populations (Naumann et al., 2014;Vereinte Nationen, 2021).Drought is a complex multivariate phenomenon and its universal identification, and the definition of drought impacts remain a challenge (Blauhut et al., 2016;Vogt et al., 2018;Sutanto et al., 2019).Currently, our understanding of direct drought impacts exceeds those manifested through indirect pathways, and this is primarily due to difficulties in identifying and monitoring indirect drought effects (Vereinte Nationen, 2021).Remotely sensed products that monitor earth system responses during periods of drought are one promising tool that strives towards the alleviation of this issue (AghaKouchak et al., 2015).However, they suffer from tradeoffs between spatial and temporal resolution.The production of high-resolution products will provide a more holistic view of drought responses.
Vegetation is involved in numerous drought-impact pathways (Zhang et al., 2021b).Drought disrupts terrestrial water and carbon cycles, which can reduce the integrity of the ecosystem and associated ecosystem services (Banerjee et al., 2013;Crausbay et al., 2017;Han et al., 2018;Smith et al., 2020).More subtly, vegetation also affects the dynamics of drought propagation; under favourable antecedent conditions, vegetation overshoot may exacerbate and facilitate the onset of rapid and intense droughts (Zhang et al., 2021b).Vegetation is also expected to play a crucial role in shaping drought resistance under future climate change (Vereinte Nationen, 2021).In the absence of such resistance, interventions to alleviate the negative impacts of disrupted ecosystem services can cost up to a billion dollars per drought event (Banerjee et al., 2013;Cammalleri et al., 2020).
In addition to disaster mitigation and management, predicting vegetation-drought dynamics is essential for numerous other applications, such as advances in fundamental ecological theory (Murray et al., 2018;Schwalm et al., 2017;Meza et al., 2020;Chen et al., 2021), studies investigating future change in land use, and climate change interventions (Jiang et al., 2017).
It follows that formulating appropriate responses to drought and alleviating the negative effects of ecosystem disruption during these periods requires accurate predictions.Vegetation indices derived from satellite-based remote sensors, such as the Enhanced Vegetation Index (evi), have proven to be an indispensable tool for monitoring vegetation at multiple scales, from the fine scale, such as crop patches (Moussa Kourouma et al., 2021;Sharifi, 2021) to the global scale (Huang et al., 2021;Vicente-Serrano et al., 2010).In recent decades, numerous satellite-based vegetation indices have been developed (Li et al., 2021a).However, a persistent feature among these products is trade-offs between spatial resolution and revisiting times, where high temporal resolution is met by coarse spatial resolution and vice verse.For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) captures the entire Earth with a high temporal resolution every 1 to 2 days (Zhao and Duan, 2020) with a maximum resolution of 250 m.Landsat and Sentinel-2 data have a higher spatial resolution, 10 and 30 m, but longer revisiting times of approximately 10 and 5 days, respectively (Zhu, 2017;Li et al., 2021a).Revisiting times for Landsat and Sentinel-2 are further prolonged when sensors or retrievals are interrupted by bad weather.Besides temporal frequency, temporal coverage is another important consideration.Coarser scale products are associated with older satellites and have more extended temporal coverage than the newer ones; MODIS products reach as far back as 1999 whereas Sentinel-2 products only go back to 2017.The ideal product for monitoring vegetation dynamics would have global coverage, little to no data gaps, and high spatial and temporal resolution.
Machine learning (ML) methods have been used for downscaling and gap-filling purposes in remote sensing products and can be seen as one tool that may lead to the production of high quality remote sensing products (Zhu et al., 2022;Zeng et al., 2013).Furthermore, ML methods have been successfully applied to a wide range of drought-related (Hauswirth et al., 2021;Shamshirband et al., 2020;Tufaner and Özbeyaz, 2020;Shen et al., 2019;Das et al., 2020;Hauswirth et al., 2022) and remotely sensed vegetation studies (Roy, 2021;Li et al., 2021b;Reichstein et al., 2019).Compared to conventional statistical downscaling techniques, ML is considered the superior alternative; given that no strict statistical assumptions are required, https://doi.org/10.5194/hess-2022-430Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.
were not available at 0.1 • resolution, the nearest-neighbour interpolation scheme from xarray (Hoyer and Hamman, 2017) was used to match the variables to the same spatial resolution.
Vegetation Index -The reference data used in this study is the evi index.evi data provide the observational benchmark for the training and validation of the ML-based products created in this study.The evi can be used as an indicator of overall vegetation status and health, as it is sensitive to chlorophyll content and correlates with primary production, photosynthesis rates, and vegetation physiognomy (Box et al., 1989).Compared to the more widely used Normalized Difference Vegetation Index, evi is considered the superior index, as it is less sensitive to atmospheric conditions and saturation effects in areas of dense vegetation (Gao et al., 2000).These data arise from the Moderate Resolution Imaging Spectroradiometer aboard the Terra and Aqua satellites.Sensors aboard Terra and Aqua are identical, and the 16-day composite images from each sensor are released 8-days apart.In this study, Google Earth Engine's python Application Program Interface (Gorelick et al., 2017) was used to access the terra (MOD13A2.006)and aqua (MYD13A2.006)evi data.These two products were combined to produce a quasi-eight-day time series (Didan, 2015(Didan, , 2021)).
Feature Variables -Global vegetation type patterns are largely driven by terrestrial water and energy balances (Hawkins et al., 2003).Similarly, the responses of vegetation to drought are regulated, in part, by water and energy availability (Xu et al., 2010).Accordingly, a suite of data indicative of terrestrial water and energy balances was selected as potential input variables.
These variables are introduced below, and Table 1 provides an overview.
Meteorology -Hourly data for total precipitation (tp), two-meter temperature (t2m), volumetric soil moisture layer 1 (swvl1), and soil temperature layer 1 (stl1) were retrieved from the hourly ERA5-Land Reanalysis product by the European Centre for Medium-Range Weather Forecasts (Muñoz-Sabater et al., 2021).In addition, potential daily evaporation (pet) was acquired from Singer et al. (2021), pet is calculated following the Penman-Monteith formulation with ERA5-Land as the input data.All meteorological data were resampled to match the 8-day frequency of the evi data.Tp was aggregated by taking the cumulative sum of the previous 8 days, whereas the remainder of the variables were averaged over 8 day windows.
Drought Indices -Aside from short-term changes in water availability, it is also key to understand the long-term dynamics to identify drought legacy effects on the current vegetation states (Schwalm et al., 2017).To this end, the Standardised Precipitation Index (spi) (McKee et al., 1993) and Standardized Precipitation Evapotranspiration Index (spei) (Vicente-Serrano et al., 2010) were used to characterise these legacy effects.The spi and spei were calculated at the 1, 3, 6, 9, 12 and 24-month aggregation lengths.The different lengths of aggregation are related to types of drought: precipitation, soil moisture, and hydrological droughts.Precipitation and soil moisture droughts mostly correlate short-term deficits in soil water (1-3 months), and are important for vegetation with shallow roots; hydrological drought (6-12 months) can be a good proxy for impacts on shrubs, bushes and trees that have deeper roots and are likely to rely on local groundwater for water (12-24 months).In addition, the inclusion of drought indices allows for the characterisation of past climate memory effects on current vegetation growth (Reichstein et al., 2019;Schwalm et al., 2017) associated with past climatic conditions.The equations and steps for calculating spi and spei are detailed in Appendix A2.
In this product, landcover types are classified according to the International Geosphere-Biosphere Programme classification scheme.Barren land, deserts, permanent snow and water bodies were masked in all further analyses.It is important to note that the RF was supplied with the remainder 15 unique landcover types; however, these were collapsed into eight broader classifications for brevity and clarity in the results, discussion and visualisations.To capture the variations in water and energy 130 availability attributable to topographic effects, elevation (elv) and height from the nearest drainage basin (hand) were accessed from MERIT Hydro, a high-resolution global hydrography map (Yamazaki et al., 2019).Lastly, slope and aspect was calculated from elv using the relevant functions in xarray-spatial (Hoyer and Hamman, 2017).While an abundance of ML approaches has been used to predict vegetation status, here the Random Forests Regressor (RF) was selected to link meteorology, land cover, topography, and drought inputs to vegetation health.RF is an ensemble method that fits many decision trees on different subsets of data.RF is advantageous given its relatively straightforwarded implementation, ability to incorporate categorical features, ability to easily identify causal links and limited risk of overfitting.The general pipeline used throughout consisted of five sequential steps (Fig. 1).Here, the RF was implemented in Python 3.9 (Rossum and Drake, 2010) under the scikit_learn framework (Pedregosa et al., 2011).
Feature Selection -In an attempt to include only relevant data in the ML model, the list of potential variables described in Section 2.1.1 and Table 1 was evaluated for their ability to provide meaningful information during model fitting.A pairwise Spearman rank correlation was calculated between all features to ensure that input data correlated with evi.
Those variables that exhibited strong correlations were retained in further analysis, whereas variables that experienced weak correlations were excluded.Aspect did not exhibit strong correlations with evi (Fig. A1).Similarly, spi (at all aggregation times) did not correlate strongly with evi.In addition, spi and were closely correlated with spei, spi was excluded in favour of spei (Fig. A1).spi and aspect were excluded from further analyses; features that were excluded are highlighted in Table 1.
Pre-processing -Given that the RF algorithm accepts 2-dimensional numeric arrays as input, the 3-dimensional data was processed so that each unique latitude and longitude was associated with a time series of each variable.The single categorical feature (lc) was converted to binary numeric.Each unique landcover type is assigned to a new feature, with 1 indicating presence and 0 indicating absence.
Split Strategy and Hyper-parameter optimisation -HalvingRandomSearchCV with a 3-fold cross-validation approach was applied to refine the number of estimators and maximum depth.This hyper-parameter optimisation provides the optimal configuration for the RF so that the critical vegetation dynamics are captured while simultaneously reducing the scores (Fig. 2a).Given that only the risk of overfitting increases with increasing Maximum_depth and number_of_estimators and only marginal increases in test scores are observed past these points, 12 and 15 were identified as the optimal settings.
After determining optimal parameter settings, the data were split into training and validation sets.However, three-dimensional data could conceivably be split along the temporal dimension where the model is trained on all locations with only a subset of the temporal availability (i.e., temporal splitting), or the data can be split according to location where only a subset of the grid pixels are selected for training but over the entire available period (i.e., spatial splitting).Given that previous research has highlighted that RF performance is sensitive to spatial vs temporal splitting, this is especially true for extreme events such as droughts (Hauswirth et al., 2021).We conducted a cursory analysis to determine whether a temporal or spatial splitting approach better balances trade-offs between computational complexity and learning rates.Learning curves for cursory RF models using each splitting approach were quantified and compared.Each model was supplied with increasing training sizes, and test scores were calculated and plotted to visualise learning curves.This cursory analysis revealed that spatial splitting yields faster learning curves than the temporal splitting approach (Fig. 2b); therefore, spatial splitting was identified as the preferred approach.Train -For the final RF model, a spatial split with a (0.06:0.94) (train: predict) ratio was used to train the final model.A 0.06:0.94split was chosen, and there was very little increase in performance past training sixes of 6% (Fig. 2b).Maximum_depth and number_of_estimators were set at 12 and 13, respectively.The parameters that were not subjected to hyperparameter optimisation were set as follows: the squared_error criterion was used to measure the quality of the splits in branches, the maximum number_of_features considered in each split was set at auto, and the minimum and maximum samples_per_leaf_nodes was set at 1 and 2, respectively.
Predict and Evaluate -As a first-pass assessment of overall performance, the model was scored using the validation and the TreeExplainer (Lundberg et al., 2020).

Downscaling globally continuous vegetation status using ML
In this section, the focus shifted toward whether RF can be used to downscale global evi values, that is, whether a model trained on 0.1 • can accurately predict evi at a finer 0.01 • scale.To this end, a 0.01 • data set was compiled.In cases where data were not at the 0.01 • resolution, the nearest neighbour interpolation scheme from xarray (Hoyer and Hamman, 2017) was used to match the variables to the same spatial resolution.This data set was used as new input data to the already trained RF model to predict evi at the 0.01 • scale.The evaluation approach for the downscaled values remained much the same, the overall model accuracy was assessed using (R 2 ) and (RMSE), and Pearson correlation coefficients were calculated for each grid cell.

Applicability of ML informed vegetation status products during periods of drought
One noticeable shortcoming of the RF is its relatively poor ability to predict extreme values depending on the training selection (Hauswirth et al., 2021).To determine to what degree this may influence the generality of the two products mentioned above, we further investigated the accuracy of the predicted evi under low growing conditions by calculating the anomaly correlation index (Eq.1), where eviA i,j denotes evi anomaly for the month j in year i, ēvi ,j denotes the average evi of month j over 2003-2013; σ stands for the standard deviation of evi during the period.

Results
The results here are presented in three parts.First, the results of the model trained on the 0.1 • data are presented; here, the focus is retained on the model's performance and ability to predict the status of the vegetation at the spatial resolution it is trained.
We also touch on which features are most important in predicting the status of the vegetation.Subsequently, we present the model's performance when used to downscale evi and predict 0.01 • data.Lastly, we explore how this module can be used to gain insight into global vegetation dynamics by assessing the accuracy of both products under drought conditions.to be a lower degree of information (Fig. 3).The model was able to reproduce global vegetation patterns by correctly predicting high vegetation density in tropical forests and low vegetation density in arid and urban regions of the world (Fig. 4).When trained on only 6% of the data, the RF was able to predict global evi accurately with a spatial resolution of 0.1 • (R 2 = 0.86; Fig. 4, 5, 6 & 7a).Looking more closely at the distribution of errors, less than 1% of grid cells showed negative correlations and more than 80% showed correlations higher than 0.5 (Fig. 7c) and RMSE ranged between 0.02 and 0.4 (mean: 0.05 ± 0.03; Fig. 7d).However, it is important to note that the accuracy was neither spatially nor temporally uniform.Landcover 220 types were an important feature in determining predictive ability.The predictions of evi in areas dominated by urban and crop landcover types showed the highest degree of error (Fig. 5a & 6a).On the contrary, the most natural types of land cover, such as forests and grasslands, were the most accurately represented by the model (Fig. 5a & 6a).For all types of land cover, the periods of maximum and minimum evi were less accurately predicted than the intermediate periods (Fig. 6a).When the model trained with 0.1 • data was used to predict evi at 0.01 • spatial resolution, there was a slight drop in accuracy.
The predictive capacity was still good but reduced compared to the 0.1 • product, with a median R 2 of 0.75 (Fig. 7b).Again, the RF was able to accurately capture spatial and temporal vegetation dynamics when supplied with 0.01 • data (Fig. 4 & 6b).
The errors also increased, the proportion of grid cells displaying negative correlations now being 5% (Fig. 7c) compared to less than 1% for the 0.1 • product.RMSE ranged between 0.04 and 0.6 (mean: 0.09 ± 0.07; Fig. 7d), with the majority of the grid cells exhibiting RMSE around 0.05.

Accuracy under drought conditions
The ACC analysis revealed that the RF was still able to capture evi anomalies (Fig. 8), but to a lesser extent compared to overall performance (Fig. 7c).The majority of grid cells showed positive correlations, with less than 10% displaying negative correlations.At least 50% of grid cells exhibited an ACC of 0.25 for 0.1 • compared to 45% when evi was predicted at 0.01 • (Fig. 8).This indicates that for that 90% of the locations, the RF can reproduce anomalies from the average seasonal cycle and thus can be used to identify periods of negative or positive evi impacts resulting from droughts or more favourable growing conditions.overall performance of the RF is discussed; after that its usefulness as a gap-filling or downscaling tool is critically evaluated.
We then highlight important outcomes from applying the model during periods of drought.Lastly, we discuss the importance of landcover types concerning model accuracy.

Overall performance
The RF successfully predicted evi at 0.1 • scale from meteorological, topography, and landcover types as input data.Of these data, features related to water and energy balances were most important in predicting evi.The RF successfully captured annual vegetation growth cycles and was able to distinguish between the main global biomass with high accuracy.These results add to the current body of research showing that the RF is a powerful technique for predicting temporal and spatial vegetation dynamics from remote sensor data (Roy, 2021;Staben et al., 2018;Wang et al., 2021).Error analysis revealed that prediction accuracy could have been more homogeneous across space and time and varied according to the growing season and land cover type.This behaviour can also be linked to the relative abundance of land cover types, where more dominant land cover types are simulated with higher accuracy.

RF for Gap-filling
A promising aspect of this study is that the RF can accurately predict evi at unseen geographic locations when trained on relatively few data -only 6% in this case.It follows that this approach can be used to for gap-filling purposes and produce high-resolution vegetation indices from other satellite sources or be used in conjunction with satellite products.For example, Landsat and Sentinel-2 data produce high-resolution data vegetation products; however, retrievals are strongly affected by weather conditions, which results in data gaps.In addition, its relatively low orbiting altitude means that the spatial coverage for each pass over is low.The approach outlined in this study could be applied to Landsat and Sentinel-2, to produce continuous vegetation index data sets at the 30-10m spatial resolution.This approach has been previously used to impute missing values for other remote sensing products like land cover types (Holloway-Brown et al., 2021), leaf traits (Moreno-Martínez et al., 2018) and soil moisture (Nguyen et al., 2022) and can now be the extent to vegetation health or evi.

RF as downscaling tool
The RF accurately predicted evi at finer spatial scales than was trained, successfully predicting evi at a scale of 0.01 • using highresolution auxiliary data.However, it should be noted that this resulted in a reduction in precision compared to the 0.1 • product.
This is an expected result, given that evi at the 0.01 • resolution will exhibit greater variances and more extreme values during periods of high and low growth.Scale-dependent drivers of vegetation dynamics may be another phenomenon that contributes to decreased precision when predicting evi at the 0.01 • using a model trained at a coarser resolution.Meteorology has been shown to be tightly coupled to vegetation at the ecosystem scale but less so at finer scales, where biotic processes, such as competition, herb ivory, disease, and fire, are more important (Franklin et al., 2020).When predicting evi, the relative increases in error remained small.This product can be of particular use in cases where the benefits of having high resolution long-term

RF for predicting drought effects
Compared to the overall performance, the RF was less capable of capturing extreme values of evi.Increased error among extreme values is a known limitation of the RF (Hauswirth et al., 2021).During RF training, an evaluation metric, in this case squared_error, is used to minimise the error for the model as a whole.In this scenario, optimal fits inevitably result in reduced errors for values close to the mean at the expense of inflated errors for the outliers (Ribeiro and Moniz, 2020).
In the current study, this means that evi during normal growth periods is prioritised over periods of extremely low or high vegetation growth.The production of ML frameworks that accurately reproduce vegetation responses during extreme periods is an essential consideration for future research directives.

Importance of Landcover Types
Varying error according to landcover type in the 0.1 • and 0.01 • is expected for at least two reasons.The first relates to the inherent features of the RF algorithm itself, and the second to the environmental process that affects the dynamics of evi.A limitation of the RF algorithm is that when data is imbalanced, underrepresented groups are less well explained by the algorithm.Accordingly, accuracy varied according to a proportional abundance of landcover types (Jung et al., 2020).
Dominant landcover types, such as forests and grasslands, displayed the least amount of error; in contrast, minority landcover types regions that have undergone human modification (i.e., urban areas and croplands) were associated with the highest error.
Second, the features used in this study may not incorporate processes critical to vegetation status equally among landcover types (Moussa Kourouma et al., 2021).Forests, grasslands, and other natural ecosystems are closely coupled to natural weather processes.However, croplands and urban areas may be less influenced by weather and more influenced by anthropogenic manipulations of water and energy balances (Zhang et al., 2004;Hawkins et al., 2003;Tang et al., 2021).A potential solution to this problem is to rely on Extreme Gradient Boosted Decision Trees, which have been shown to provide more accurate predictions where data are imbalanced (Li et al., 2021b) or include information on human-water management to better represent drought responses (Wanders and Wada, 2015).
Landcover-specific variations in the model's ability to predict vegetation are an important outcome of this study.Apart from the statistical reasons detailed in the previous paragraph as potential mechanisms for this phenomenon, an additional, and most likely compounding explanation is that the data used to predict evi may be more relevant for some landcover types and levels of vegetation growth than others.For instance, vegetation status in urban areas and croplands shows weak correlation or high errors (Fig. 5 & 6).The meteorological data used here to predict evi may not be the only factor driving the vegetation dynamics in human-modified areas.It is possible that irrigation and water management influence vegetation.Indeed, vegetation in urban areas have been shown to grow more rapidly and have a longer growing season than rural counterparts; this is thought to be driven by higher temperatures, high concentration of airborne phosphorous and other aerosol pollutants (Sicard et al., 2018a, b;Pretzsch et al., 2017).In contrast, natural forests and grasslands show high levels of accuracy and correlations, thus suggesting that the data used here is appropriate for the machine learning models to capture vegetation dynamics.Similarly, poor accuracy in wetlands is not unexpected as wetland vegetation is primarily driven by water quality, salinity, and pH (Grieger et al., 2021).
On the contrary, forests and grasslands show high accuracy when using meteorological variables, since these are important drivers of vegetation growth in these areas.Although not directly related to vegetation, Hauswirth et al. (2021) showed that by including water management practices in machine learning models, the predictions of groundwater head and stream flow were more accurately predicted.Therefore, we further iterate the need to provide models with appropriate input data sources.

Conclusions
The results from this study reveal that the RF is an appropriate method for predicting evi at the global scale, at the 0.1 • and downscaled 0.01 • resolution.In general, RF was capable of predicting evi dynamics with high accuracy; global patterns of vegetation and temporal dynamics were well captured.However, it is essential to note that higher errors were associated with underrepresented landcover types and periods of extreme vegetation growth, such a drought periods.Lower accuracy for underrepresented classes in unbalanced data sets and a hampered ability to predict extreme values is a common criticism of the RF.In accordance with this study, landcover types that account for a smaller fractional cover of the earth's surface, and periods of extreme vegetation growth, were associated with the highest error.Predicting evi at a finer resolution resulted in increased errors.This is attributable to higher variances in the 0.01 • product compared to 0.1 • and it is important to note that the relative increases remained small.
The results here also highlight the use of RF for efficiently and accurately predicting missing data and downscaling, ultimately allowing for the production of spatially continuous evi data at very high spatial and temporal resolutions.Toward this end, this study produces spatially continuous evi product at 0.1 • and 0.01 • resolution, which could be used to fill existing gaps in satellite observations or in conjunction with satellite data to have improved monitoring of drought impacts on vegetation.
This study adds to previous research efforts that have successfully applied the RF in predicting vegetation status.Here the RF was used to produce a global spatial and temporally continuous evi product at 0.1 • and 0.01 • , with a median R 2 of 0.86 & 0.75, respectively.The approach outlined in this study could be applied to Landsat and Sentinel-2, to produce continuous vegetation index data sets at the 30-10m spatial resolution.The RF algorithm is a powerful technique for predicting temporal and spatial vegetation dynamics from remote sensor data, as well as those using RF for gap filling purposes on remote sensing products.The novelty of this product, compared to previous studies, is that it has global coverage, high spatial and high temporal resolution.For the calculation of spi: where i is the month in question and m = i − scale.
For the calculation of spei: where: D i = tp i − pet i and

Figure 1 .
Figure 1.The five sequential steps followed during the RF fitting and evaluation.

Figure 2 .
Figure 2. (a) Evolution of RF performance during HalvingRandomSearchCV hyper-parameter optimization of: maximum_depth (blue) and number_of_estimators (red).(b) RF performance following the incremental increase of train set size using a location (red) based split approach compared to a temporal (blue) based split approach.
default coefficient of determination (R 2 ) scorer in the scikit RF implimentation.The model predictions were further evaluated https://doi.org/10.5194/hess-2022-430Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.bycalculating the root mean squared error (RMSE) and Pearson correlation coefficients.These were calculated independently for each grid cell to provide information on the spatial variation of errors.Last, to gain insight into which features were the most essential for predicting evi, global feature importance was calculated using Shapley Additive exPlanations' (SHAP)

3. 1
Random Forest Regressor at the 0.1 • resolution SHapley Additive exPlanations values provided an understanding of the relative importance of each feature in predicting evi.The most important features were those associated with meteorology, landcover type and topography; drought indices proved https://doi.org/10.5194/hess-2022-430Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.

Figure 3 .
Figure 3. Feature importance for the RF predicting evi at 0.1 • .The features are ordered by level of importance, with higher mean SHAP values indicating higher importance.

Figure 4 .
Figure 4. Mean (2003 -2013) observed (top) and mean predicted evi values by the model at the 0.1 • (left) and 0.01 • (right).Barren land, deserts, permanent snow, and water bodies were masked and represented by black.

Figure 5 .
Figure 5. RMSE calculated for the Amazon Basin, Great Lakes and Western Europe for predicted evi values by the model at the (a) 0.1 • and (b) 0.01 • .

Figure 7 .Figure 8 .
Figure 7. : (a) Scatter plot of observed and predicted evi at 0.1 • and (b) 0.01 • ; Cumulative distribution function for (c) Pearson Correlation Coefficients overall grid pointsat 0.1 • (blue) and 0.01 • (orange), (d) violin plot of RMSE for all grid points at 0.1 • and 0.01 • https://doi.org/10.5194/hess-2022-430Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.evi products outweigh the limitations associated with error increases.The product presented here can and should be used in further studies investigating global vegetation dynamics.
Figure A1.Correlation Matrix of pairwise Spearman rank correlation coefficients between all potential variables tp i−j,l + j l=1 tp i,l , if j < k j l=j−k+1 tp i,l , if j ≥ k (A3)This time series is then fitted to a gamma distribution taken the following steps: First α and β fitting parameters as calculated as: https://doi.org/10.5194/hess-2022-430Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.