## 2 Materials and Methods
### Variable Selection  
Variable selection has an important influence over land cover classification exercises, particularly as a method to reduce data dimensionality and redundance and improve classification accuracy and computational efficiency and interpretability [48, 87, 88]. As expected, we found high (Pearson's) correlations between several of the predictors (Figure S3), supporting the need for a variable selection procedure.  
Variable selection in RFs often depends on variable importance measures. Among the several measures available (see [87, 89]), the permuted variable importance measure (hereafter referred to as the PVIM) is the most popular for RFs [90]. It measures the mean decrease in out-of-bag accuracy caused by permuting each variable over all the trees; the variables causing larger mean decreases are more important and can thus be ranked based on the PVIM. Performance-based variable selection algorithms implemented by R libraries, such as VSURF and varSelRF, automate this process to optimize RFs to achieve the highest accuracies [89, 91, 92]. However, these algorithms select the best variable set according to the predictors that together minimize the out-of-bag error. As a performance validation measure, out-of-bag error is less reliable than totally independent validations, given that each new sample is only validated by the trees that did not include them in the training phase. Moreover, out-of-bag error measures may perform inconsistently with our dataset because the spatial independence within the training dataset is greater among pixels of different polygons than among pixels of the same polygon. For these reasons, we applied a modified variable selection procedure inspired by the backward recursive feature elimination process implemented in varSelRF, but with a stricter and more demanding independent validation and variable elimination method throughout each step of the process, as described below.  
The overall framework of our step-by-step variable selection process is the following: (1) start with a full model, including all the predictors, and compute an RF; (2) perform the validation of the model (using the totally independent validation dataset) and record the Overall Accuracy and Kappa Coefficient (see Section 2.12 for details on these validation statistics); (3) compute the PVIM and select the least important variable to be eliminated; and (4) repeat the process until only two predictors are left. The PVIM has, however, been shown to be potentially biased in the presence of a large number of correlated variables [48, 87]. To account for this, we computed the PVIM at each step of the recursive elimination, instead of only at the beginning. Moreover, to ensure robustness, we repeated each of the RF models with the same set of predictors 25 times and excluded the variable that featured more times as the least important at each step. The best predictor set was then chosen based on a plot of the mean Overall Accuracies and Kappa Coefficients achieved in each step, averaged over the 25 RF replicates. The Overall Accuracy and Kappa Coefficient threshold were set considering the balance of keeping the prediction accuracy high while using the minimum number of predictors, reducing redundancy, and increasing interpretability.