A Random Forest Modelling Procedure for a Multi-Sensor Assessment of Tree Species Diversity

: Earth observation data can provide important information for tree species diversity mapping and monitoring. The relatively recent advances in remote sensing data characteristics and processing systems elevate the potential of satellite imagery for providing accurate, timely, consistent, and robust spatially explicit estimates of tree species diversity over forest ecosystems. This study was conducted in Northern Pindos National Park, the largest terrestrial park in Greece and aimed to assess the potential of four satellite sensors with di ﬀ erent instrumental characteristics, for the estimation of tree diversity. Through ﬁeld measurements, we originally quantiﬁed two diversity indices, namely the Shannon diversity index ( H’ ) and Simpson’s diversity ( D1 ). Random forest regression models were developed for associating remotely sensed spectral signal with tree species diversity within the area. The models generated from the use of the WorldView-2 image were the most accurate with a coe ﬃ cient of determination of up to 0.44 for H’ and 0.37 for D1 . The Sentinel-2 -based models of tree species diversity performed slightly worse, but were better than the Landsat-8 and RapidEye models. The coe ﬃ cient of variation quantifying internal variability of spectral values within each plot provided little or no usage for improving the modelling accuracy. Our results suggest that very-high-spatial-resolution imagery provides the most important information for the assessment of tree species diversity in heterogeneous Mediterranean ecosystems.


Introduction
Forest ecosystems cover nearly one third of the Earth's land surface and contain over 80% of the terrestrial biodiversity [1]. Forest biodiversity expressing the variability among living organisms in forest ecosystems is an important factor for the terrestrial ecosystems' functions and process [2] and several studies have demonstrated the positive linkages between forest biodiversity and ecosystem stability and health [3]. Other studies have also demonstrated the role of forest biodiversity in forest ecosystem services provision, such as biomass production, carbon storage, wild fruits production, nutrient use and retention, and water provision [4,5]. Yet, direct and indirect factors, such as habitat loss, landscape degradation, fire, soil erosion, anthropogenic activities, demographics, and climate change threaten and lead to a decline in the forest biodiversity levels [6]. In particular, Mediterranean forest habitats, characterized as one of the richest biodiversity hotspots worldwide are among the most threatened on Earth [7,8], identified as priority area for conservation [9].
Remote Sens. 2020, 12, 1210 3 of 16 and Landsat-8 OLI satellite sensors, has been recently evaluated by Toressani et al. [2] for tree species diversity estimation in an alpine conifer forest. Forest beta-diversity was evaluated by Khare et al. [28] in a subtropical deciduous forest using images acquired from Pléiades 1A, RapidEye, and Landsat-8 OLI satellite sensors. Ozdemir et al. [29] evaluated diversity correlation with texture measures calculated from RapidEye, SPOT-5, and Aster images. Finally, Chrysafis et al. [30], using a Lasso regression approach, evaluated Sentinel-2 MSI and RapidEye sensors for tree diversity estimation in a heterogeneous forest in northwest Greece.
Since only a limited number of studies have evaluated the effect of image spectral and spatial resolutions in the estimation of tree species diversity, the aim of this study was to assess, at multiple scales, tree species α-diversity, using Landsat-8 OLI, Sentinel-2 MSI, RapidEye, and WorldView-2 optical satellite imagery, in forest habitats of Greece's largest terrestrial National Park.
Since tree species diversity detection from satellite imagery is expected to vary with image resolution [3], we selected a range of satellite sensors with different instrumental characteristics. In addition, Landsat-8 OLI and Sentinel-2 MSI images are accessible cost-free, offering the opportunity to minimize costs for forest ecosystems mapping and monitoring as for any project based on remote sensing [31].
Contrary to the majority of previous studies employing multiple linear regression, we adopted a random forest (RF) regression modelling procedure [32]. RF includes integrated cross-validation for assessing the accuracy of the model and is capable of ranking important variables and addressing data that interacts in a non-linear manner [32][33][34][35][36]. Furthermore, a drawback in the use of multiple linear regression relates to the assumptions of linearity and the absence of collinearity amongst variables. However, these assumptions are rarely observed in forest and remote sensing data, underlining the need for the use of more robust statistical methods [37].
Specifically, this study aims to answer the following questions: (i) Which are the most relevant spectral bands for tree diversity estimation? (ii) What is the strength of the relationship between remotely sensed data and field-measured diversity indices? (iii) Which sensor provides the best remote-sensing-based approach for tree species diversity mapping? And finally, (iv) is random forest regression an efficient approach for modelling tree diversity?

Study Area
The Northern Pindos National Park in northwest Greece has a total area 1969 sq. km and it is among the largest protected areas of the country, including eleven sites belonging to the Natura 2000 European Network of Protected Areas. The mountainous landscape comprises an intense relief with high topographical diversity. Due to the montane character of the area, the prevailing climate type is transitional from Mediterranean to continental and locally varies depending on elevation and aspect; the amount of precipitation is among the highest in Greece, a fact that is imprinted to the extended plant coverage of the area. Annual rainfall ranges between 1000 and 2250 mm, whereas the annual maximum temperature reaches 34 • C.
The vegetation at lower altitudes (400-1000 m) is composed of thermophilous, deciduous oak species, such as Hungarian oak (Quercus frainetto), Turkish oak (Q. cerris), Macedonian oak (Q trojana), Downy oak (Q. pubescens). Pure or mixed stands of deciduous tree species, such as the hop hornbeam (Ostrya carpinifolia), oriental hornbeam (Carpinus orientalis), common hornbeam (Carpinus betulus), and manna ash (Fraxinus ornus), are also found in this zone. Most important forest ecosystems are found at higher altitudes (900-1600 m) in the National Park. They consist of pure and mixed stands of black pine (Pinus nigra), beech (Fagus sylvatica), and King's Boris fir (Abies borisii-regis). Bosnian pine (Pinus leucodermis) formation occupy the upper mountain slopes and form the tree limit, while sub-alpine grasslands extend above 1800 m to peaks. Riparian vegetation along the rivers and streams within the Park is comprised mainly of oriental plane (Platanus orientialis) galleries.

Remote Sensing Data and Preprocessing
Since the main aim of the study was a comparative multi-resolution assessment of tree species diversity, we used sensor data with different spectral and spatial resolutions. The medium spatial resolution imagery evaluated in this study included Sentinel-2 MSI (S-2) and Landsat-8 OLI (L-8) images acquired on 25th August 2017 and 12th July 2017 respectively. Both sensors have spectral bands covering visible to shortwave infrared (SWIR) spectrum, and a 12-bit radiometric resolution. The spatial resolution of the Landsat-8 OLI bands used is 30m, while Sentinel-2 MSI bands are available at 10 m, 20 m, and 60 m spatial resolution ( Table 1). The images provided by the European Space Agency (Sentinel-2) and United States Geological Survey (Landsat-8) were preprocessed at surface reflectance and orthorectified to correct for relief displacement. Higher spatial, but lower spectral, resolution RapidEye 3A level imagery acquired on 25th August 2017 was also evaluated for tree diversity mapping. RapidEye (RE) is a commercial optical Earth observation mission with a pixel-size of 6.5 m and five optical bands in the visible near-infrared part of the electromagnetic spectrum (VNIR), with a 12-bit dynamic range. The level 3A provides images radiometrically corrected at-sensor, orthorectified to a map projection, and resampled to a 5-m pixel size.
Commercial very-high-spatial-resolution (2 m) WorldView-2 images, available at Ortho Ready Standard Imagery (OR2A) level, acquired on 15th of July 2015, were also used. The OR2A are radiometrically and sensor corrected, mapped to the average base elevation of the terrain covered by each individual scene. For the orthorectification process of the WorldView-2 (WV-2) images, ground control points identified on very large scale orthoimages and a 5-m digital elevation model were used.
Finally, for the conversion of the RE and WV-2 DN values into surface reflectance, the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) algorithm/code was employed [38].

Field Data Collection and Tree Diversity Indices
The field data were collected during the 2018 summer period (July -August). A total number of 100 sampling plots, 20×20 m each, were established using the gradsect method of survey [39] (Figure 1). To implement this method, field plots were located randomly along gradsects that were oriented along the strongest gradients of topographic variables. To establish the gradsects, the elevation, slope, and aspect of the area were used along with the road network. Within each plot, the coordinates from Remote Sens. 2020, 12, 1210 5 of 16 both the center and one corner were recorded with high accuracy Global Navigation Satellite System (GNSS) in order to minimize positional inaccuracies. Then, the tree species were identified, and all the individual trees with a diameter at breast height (dbh) larger than 8 cm were recorded. As a result, the recorded trees included all the canopy overstory individuals as well as a number of the midstory and understory individuals whose diameter exceeded 8 cm. Records of the tree species in each plot were used to calculate diversity indices on the basis of the number of individuals (presence-absence).
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 16 along the strongest gradients of topographic variables. To establish the gradsects, the elevation, slope, and aspect of the area were used along with the road network. Within each plot, the coordinates from both the center and one corner were recorded with high accuracy Global Navigation Satellite System (GNSS) in order to minimize positional inaccuracies. Then, the tree species were identified, and all the individual trees with a diameter at breast height (dbh) larger than 8 cm were recorded. As a result, the recorded trees included all the canopy overstory individuals as well as a number of the midstory and understory individuals whose diameter exceeded 8 cm. Records of the tree species in each plot were used to calculate diversity indices on the basis of the number of individuals (presence-absence). In order to compare tree diversity levels among regions, temporal periods, etc., the use of mathematical functions known as diversity indices is well established [40]. We used two indices ( Table 2) to account for richness and evenness components of diversity. Species richness refers to the absolute number of species present in the ecosystem, while evenness considers the relevant abundance of the species. The Shannon Index (H') of diversity [41] takes both species abundance and species richness into account [42]. In addition, we used one other diversity index related to the original Simpson index D giving more emphasis to the evenness component of diversity [43]. The Simpson's diversity (D1) is the complement of Simpson's original formulation and represents the probability that the two individuals represent different species [40]. In order to compare tree diversity levels among regions, temporal periods, etc., the use of mathematical functions known as diversity indices is well established [40]. We used two indices ( Table 2) to account for richness and evenness components of diversity. Species richness refers to the absolute number of species present in the ecosystem, while evenness considers the relevant abundance of the species. The Shannon Index (H') of diversity [41] takes both species abundance and species richness into account [42]. In addition, we used one other diversity index related to the original Simpson index D giving more emphasis to the evenness component of diversity [43]. The Simpson's Remote Sens. 2020, 12, 1210 6 of 16 diversity (D1) is the complement of Simpson's original formulation and represents the probability that the two individuals represent different species [40]. Table 2. Diversity indices used in the study.

Species Diversity Index Formula Reference
Shannon Index (H') Simpson's diversity (D1) where p i is the proportion of the ith species in the sampling plot

Spectral Information Extraction
Spectral information inherent in the original image bands of each sensor was considered for tree-diversity modelling. For the Landsat-8 OLI bands and Sentinel-2 MSI 20 and 60 m bands, pixel values corresponding to each plot center were extracted and used for the analysis. For all other images, the mean spectral values of the pixels falling within each plot were computed. Furthermore, for these images, having a finer pixel size than plot size, the Coefficient of Variation (CV) was computed as a measure for the quantification of the amount of the variability in forest vegetation within each plot [2,22].

Statistical Modelling
According to previous studies, linear regression models have been widely used for tree diversity patterns prediction [14]. Therefore, we originally examined the data in terms of the assumptions required for implementing a linear regression approach. In a preliminary step using the Shapiro-Wilk test and QQ normal probability plots, we identified that the requirements of linear regression for normality and homoscedasticity could not be met for our data. Therefore, in order to identify the correlation between spectral information and diversity indices, we used the Spearman rank coefficient, a non-parametric measure of correlation [45]. Subsequently, image-based predictive models for the tree-diversity indices were developed using a Random Forest (RF) regression algorithm [32], as implemented in the randomForest package [46].
Originally, using a Spearman's rank analysis, we eliminated highly correlated variables in order to minimize instability during variable importance computation. Then, we applied a variable selection procedure based on backward elimination process for determining the variables of the highest predictive value [47]. We selected this approach since several studies have demonstrated that the identification and inclusion of the most important variables in a reduced parsimonious model improves the efficiency of prediction, minimizes the influence of noisy predictors and redundant data, and increases the interpretability of the final model [48,49].
More specifically, the variables of each image were compared one by one through the Spearman rank matrices. We then grouped these variables according to the area of the electromagnetic spectrum (i.e., visible near infrared and shortwave infrared) recorded by the respective sensor band, and we removed from further analysis the highly correlated variables (r higher or equal to 0.7 with a p-value smaller than 5%) with the highest ranking predictor in the respective part of the spectrum [47,50].
The remaining variables were used for the initial training of the RF models for each image, and the variable importance rankings were generated considering the percentage increase in the mean squared prediction error (%IncMSE) [32,51]. Finally, a backward elimination selection procedure was followed so as to produce the most accurate model with the smallest subset of variables retained for each sensor [52].
For all models, the number of regression trees grown (ntree) was optimized based on an out-of-bag estimate error, the number of random variables used in each tree (mtry) and the minimal size of the Remote Sens. 2020, 12, 1210 7 of 16 terminal nodes of the trees (nodesize) was set to default values [53]. All models were evaluated based on their performance to predict the field measured indices. Thus, the coefficient of determination (R 2 ), the Root Mean Square Error (RMSE), and the Relative Squared Error (RSE) have been calculated following a 10-fold cross-validation.

Results
The identification of the variables with the most predictive power was carried out for each image using the percentage increase in the mean squared prediction error (%IncMSE) of RF models (Figures 2-5).

Results
The identification of the variables with the most predictive power was carried out for each image using the percentage increase in the mean squared prediction error (%IncMSE) of RF models (

Results
The identification of the variables with the most predictive power was carried out for each image using the percentage increase in the mean squared prediction error (%IncMSE) of RF models ( Figures  2-5).        Figure 2 presents the hierarchical ranking of the relative importance for the predictor bands of Landsat-8 dataset using the percentage increase in mean squared prediction error (%IncMSE). The near-infrared band (B05) was the most important variable for both diversity indices (H', D1). The coastal band (B01) band was ranked second for the prediction of H' and D1 indices.
In the case of the Sentinel-2 MSI dataset (Figure 3), the coastal band (B01) band was the most important variable for both diversity indices, followed by the green band (B03). The results also indicate that mean spectral values are in general more important than the coefficient of variation quantifying spectral variation within each plot. The coefficient of variation in the blue band (B02) was the most important among these spectral variation measures. The mean blue band (B02) was also ranked higher than most of the variables in both diversity indices.
Results on the importance of RapidEye bands (Figure 4) indicate that the blue (B01) and red (B03) bands were the most significant variables for both diversity indices. Moreover, in RapidEye models, spectral variation within each plot was less important than the average spectral values recorded.
Finally, in the case of WorldView-2 dataset, the blue (B02) and the second NIR band (B08) were ranked higher from all other WorldView-2 bands. Even though none of the CV features are ranked first, the plots indicate an increased importance of spectral variation compared to lower spatial resolution Sentinel-2 MSI and RapidEye sensors. Figure 5 also illustrates the non-linear nature of the relationship between the second NIR band (B08_m) of WorldView-2 and (c) Shannon's (H') and (d) Simpson's diversity (D1) indices, confirming the use of random forest approach for developing remote-sensing-based predictive models for tree diversity. Based on the results of variable selection using the backward elimination method for the different images (Figure 6), the higher relationship between field measured diversity indices and the remotely sensed data is identified for Shannon's H' index (Table 3).
Landsat-8 dataset using the percentage increase in mean squared prediction error (%IncMSE). The near-infrared band (B05) was the most important variable for both diversity indices (H', D1). The coastal band (B01) band was ranked second for the prediction of H' and D1 indices.
In the case of the Sentinel-2 MSI dataset (Figure 3), the coastal band (B01) band was the most important variable for both diversity indices, followed by the green band (B03). The results also indicate that mean spectral values are in general more important than the coefficient of variation quantifying spectral variation within each plot. The coefficient of variation in the blue band (B02) was the most important among these spectral variation measures. The mean blue band (B02) was also ranked higher than most of the variables in both diversity indices.
Results on the importance of RapidEye bands (Figure 4) indicate that the blue (B01) and red (B03) bands were the most significant variables for both diversity indices. Moreover, in RapidEye models, spectral variation within each plot was less important than the average spectral values recorded.
Finally, in the case of WorldView-2 dataset, the blue (B02) and the second NIR band (B08) were ranked higher from all other WorldView-2 bands. Even though none of the CV features are ranked first, the plots indicate an increased importance of spectral variation compared to lower spatial resolution Sentinel-2 MSI and RapidEye sensors. Figure 5 also illustrates the non-linear nature of the relationship between the second NIR band (B08_m) of WorldView-2 and (c) Shannon's (H') and (d) Simpson's diversity (D1) indices, confirming the use of random forest approach for developing remote-sensing-based predictive models for tree diversity. Based on the results of variable selection using the backward elimination method for the different images (Figure 6), the higher relationship between field measured diversity indices and the remotely sensed data is identified for Shannon's H' index (Table 3). The most accurate model for H' is derived from the WorldView-2 data (Table 3) considering two variables (R 2 = 0.44, RMSE=0.29). The sentinel-2 MSI model included five variables (R 2 = 0.29, RMSE=0.35), slightly higher than the models of Landsat-8 (R 2 = 0.27, RMSE=0.34) and RapidEye models (R 2 = 0.21, RMSE=0.35).

Discussion
Remote sensing has become increasingly popular for mapping and monitoring the status of tree diversity [3,54]. Nowadays, with the increased number of remote sensing systems available, their improved instrumental characteristics, easier access, and increased affordability, the challenge relates to the efficient exploitation of the information content and better integration of the remotely sensed data in the diversity monitoring and conservation process [16]. Research is needed for assessing the scientific applicability of available sensors and the appropriate scale of measurement over different forest ecosystems [13] for complementing and sustaining, or even replacing [55] the field based measurements of tree diversity.
The comparative analysis of the models developed from WorldView-2, RapidEye, Sentinel-2 MSI and Landsat-8 OLI imagery, indicated that the WorldView-2 imagery was able to better approximate the tree species diversity indices in this study site in Northwest Greece. Sentinel-2 MSI models outperformed Landsat-8 OLI image models, while RapidEye models were ranked fourth in terms of productive accuracy. The advantage of the WorldView-2 dataset can be attributed to the very high spatial resolution of the images, which allows individual tree crowns to be visually identified [56]. Previous studies have also indicated that tree diversity prediction through spectral responses varies with spatial scale [3]. Khare et al. [28] calculated various diversity indices for different spatial resolutions and found that tree diversity was better explained at the finest (i.e., 0.5 m) pixel size evaluated. However, in another study in a dry tropical forest environment, it was found that Landsat imagery presented better correlation with tree diversity compared to high-spatial-resolution IKONOS data [27]. The authors related this finding to the very fine scale of IKONOS data for the specific landscape, decomposing tree crowns to sunlit and shadow-covered pixels, suggesting that larger pixel sizes than a single crown reduce the spectral variability as the reflectance signals is the average of larger ranges of surface conditions [27].
Large-scale simulation studies through field spectroscopy can be found in other ecosystem settings (i.e., grasslands), where studies suggested that the detection of α-diversity declined with lower spatial resolution and the optimal pixel size for α-diversity prediction approximates the size of an individual plant leaf or crown [57]. Yet again, in such ecosystem types, other experimental studies present findings indicating a controversy regarding the optimal pixel size [58]. In forest areas, the challenge is aggravated by the fact that optimal spatial resolution is not constant, being primarily affected by the spatial and structural parameters of the forest stands [59]. Therefore, there is a need for experimental studies across biomes, ecosystems, and sites for developing, if possible, operational, generic, applicable approaches for the remote sensing of forest biodiversity [3].
Sentinel-2 MSI and Landsat-8 OLI images, having a spatial resolution not allowing for the direct delineation of individual tree crowns, presented relatively lower levels of accuracy than the WorldView-2 images. The Landsat-8 OLI models developed in our study, presented similar levels of accuracy with the models developed by Madonsela et al. [22] in a savanna woodland area and by Nagendra et al., [27] in a dry tropical forest. Sentinel-2 MSI models, despite the lower accuracy in comparison to WorldView-2 models, presented satisfactory results (R 2 =0.31 for Simpson's diversity). This could be explained by the "surrogacy" concept [3], which suggest that spectral diversity is linked with various types of field tree diversity through physical and ecological rules. This results in a satisfactory correlation between spectral measurements and field-measured diversity, even when pixel size is larger than the size of individual trees.
Medium spatial resolution remotely sensed data, such as the Sentinel-2 MSI imagery available free to the research community, can have practical application as a screening tool to identify critical tree diversity areas [22]. The same data can be fused with information on habitat types, soil, climate, etc., for developing complex predictive models of tree diversity [27]. The RapidEye image was found less suited from all images for assessing tree species diversity despite having a finer spatial scale than Landsat-8 OLI and Sentinel-2 MSI images. A possible reason for this might be the lower spectral resolution of the RapidEye image. Previous studies suggested that increased spectral resolution resulted in an increase in accuracy of diversity estimation [60] and that a large number of spectral bands is preferable to few bands [61]. But this might be a weak hypothesis since the higher spatial resolution models of the WorldView-2 images relied on a limited number of variables (two to four). Potential radiometric distortions caused by resampling the original nominal 6.5-m spatial resolution to a 5-m pixel size (Level 3A) might be another reason driving the lower accuracy of RapidEye diversity models.
The lower accuracy of the RapidEye image might also be related to the absence of information on the SWIR part of the spectrum. This band has been incorporated in almost all Sentinel-2 MSI and Landsat-8 OLI diversity models. SWIR has essential spectral information useful for characterization of vegetation [22] and has been strongly correlated with species richness in a tropical forest ecosystem [62]. SWIR spectral bands have been also reported to provide critical information for enhancing the predictive accuracies regarding forest attributes along with the lower the visible part of the spectrum (notably blue band) and forest structure [63,64]. The blue band was also a very important variable for local diversity estimation in all sensors. The linkages the blue part of the electromagnetic spectrum and forest structure were also noted in Mediterranean [63] and boreal forest environments [65]. In general spectral features in the visible region (400-700 nm) are reflective of leaf chemistry and pigment content, which results in strong absorption features [66]. Specifically, the relative high importance of the CV of the blue band for the Sentinel-2 MSI and WorldView-2 data, might be related to the absorption in the blue by carotenoids [67].
Blue band, adjacent to the coastal one, was also a very important variable for the final Sentinel-2 MSI, RapidEye and Worldview-2 models, similar to the findings of Nagedra et al. [27]. Interestingly, the blue band presented very high correlation with the coastal one, only in the case of the Landsat-8 OLI image (Spearman rank correlation = 0.966) and therefore it was removed from the subsequent RF modelling procedure.
The NIR band, which is widely considered as a key band for diversity estimation [60], has been also included in all diversity models. Information from the coastal part of the electromagnetic spectrum was also a very important for the Sentinel-2 MSI and Landsat-8 OLI models. The reflectance measured in the coastal band is related to the chlorophyll content of plants and often used for forest biomass and carbon estimations [68] and its' strength for species richness estimation has been demonstrated in a savanna woodland [22].
In a previous study, it was identified that the inclusion of the red-edge spectral information was an essential parameter for tree diversity estimation, equal or even higher than the spatial resolution [28]. In our study, the red-edge bands of the Sentinel-2 MSI, RapidEye, and WorldView-2 images presented a weak relationship with all diversity indices and high correlation with the NIR bands, as indicated by the Spearman measure, and therefore were not included in the RF modelling procedure.
In order to quantify the spectral heterogeneity within each plot, the coefficient of variation was calculated for the 10-m Sentinel-2 MSI, RapidEye, and Worldview-2 bands. However, both the variable importance plots and the final variable selected within each model confirmed the weak relationship of this measure of spectral variation with tree diversity. Madonsela et al. [22] also noted that the relationship between spectral indices and measures of tree species diversity declined significantly when surrogate measures of variability were used as predictors of tree species diversity in a tropical savanna environment. Nagendra et al. [27] also reported that spectral/spatial variation of reflectance values within field sampling plots provides either non-significant or weakly significant correlations with species richness and diversity. Meng et al. [37], using multispectral (10-m) and panchromatic (2.5-m) Spot-6 data for forest diversity estimation, also noted that no measures related to the internal variability of forest stands (i.e., texture) were included in the final models. They suggested that this is related to the fact that diversity indices measure only a quantitative aspect of diversity, not being able to correlate with the qualitative spectral variability resulting from different types of mixture, i.e., the pronounced spectral differentiation between an evergreen conifer and a deciduous species in contrast with the similar spectral values between two conifer species. Nevertheless, as Lu and Weng underline [69], the importance of adding information related to the variability of spectral values, increases as the spatial resolution increases.
Unlike to the majority of the previous studies, we did not use a linear regression approach for assessing the relationship between remote sensing data and field the measurements of tree species diversity. RF regression does not require any data distribution assumptions and can detect interactions and higher order relationships between independent variables without a priori specification of these terms [50]. Furthermore, RF are particularly appealing due to their ability to generalize even under a limited training samples regime, as often is the case in remote sensing applications [48]. In our study, the relatively modest or low accuracy achieved by the RF models, might also be related to the field data collection strategy for the tree diversity indices calculation. Since we considered all trees with dbh larger than 8 cm, a small number of these individuals included in the sample may not be directly observed from optical sensor systems. Nevertheless, increased presence of understory species is closely related to canopy gaps, resulting also in an increased spectral contribution of lower canopy stratums and understory to forest reflectance [70].
We adopted a two-step process, removing highly intercorrelated variables for reducing the computational complexity of the algorithm and the dimension of the input data. Several related publications have shown that the predictive power of RFs may benefit from variable selection [50,71]. Finally, in terms of the most appropriate diversity index to be used with remote sensing data, the results indicate that, in our site, a remotely sensed spectral signal correlates better with Shannon's Index than Simpson's diversity, since, across all sensors, the Sentinel-2 MSI model provided the highest coefficient of determination. However, similar to the insights from previous in-depth comparative studies [44], the identification of one ideal diversity index might not be feasible.

Conclusions
In this study, we compared remotely sensed data from four different sensors, namely Landsat-8 OLI, Sentinel-2 MSI, RapidEye, and WorldView-2, for tree species local diversity estimation over the forest habitats of the Northern Pindos National Park in northwest Greece using a random forest modelling procedure.
The models relying on the very-high-spatial-resolution WorldView-2 image provided the highest accuracy for the prediction of the richness and evenness components of diversity. The use of the Sentinel-2 MSI imagery also provides an attractive option for tree species diversity estimation considering the accuracy results, when compared to the relative higher spatial resolution RapidEye imagery.
The variable selection procedure improved the predictive power and robustness of RF models compared to the original full dataset models and allowed for the identification of the important spectral bands for tree diversity estimation. The Coefficient of Variation (CV), quantifying internal variability within each plot, provided little or no usage for improving the modelling accuracy. This research confirmed that the α-diversity is scale and sensor dependent. Freely available Sentinel-2 MSI images characterized by a good compromise between spatial and spectral resolutions, seems promising for providing accurate α-diversity index estimates through remote sensing techniques. This study, as well as previous ones, emphasized the tradeoffs existing between spatial and spectral resolutions for local tree diversity estimation. Further research is required to determine the optimal combination of spectral and spatial resolutions over different biomes and forest stand characteristics.