Mapping forest structural heterogeneity of tropical montane forest remnants from airborne laser scanning and Landsat time series

Tropical montane forests are important reservoirs of carbon and biodiversity and have a central role in the hydrological cycle. They are, however, very fragmented and degraded, leaving isolated remnants across the landscape. These montane forest remnants have considerable differences in forest structure, depending on factors such as tree species composition and degree of forest degradation. Our objectives were (1) to analyse the reliability of airborne laser scanning (ALS) in modelling forest structural heterogeneity, as described by the Gini coefficient (GC) of tree size inequality; (2) to determine whether models are improved by including tree speciessensitive spectral-temporal metrics from the Landsat time series (LTS); and (3) to evaluate differences between three forest remnants and different forest types using the resulting maps of predicted GC. The study area was situated in Taita Hills, Kenya, where indigenous montane forests have been partly replaced by single-species plantations. The data included field measurements from 85 sample plots and two ALS data sets with different pulse densities (9.6 and 3.1 pulses m−2). GC was modeled using beta regression. We found that GC was predicted more accurately by the ALS data set with a higher point density (a cross-validated relative root mean squared error (rRMSECV) 13.9%) compared to ALS data set with lower point density (rRMSECV 15.1%). Furthermore, important synergies exist between ALS and LTS metrics. When combining ALS and LTS metrics, rRMSECV was improved to 12.5% and 13.0%, respectively. Therefore, if the LTS metrics are included in models, ALS data with lower pulse density are sufficient to yield similar accuracy to more expensive, higher pulse density data acquired from the lower altitude. In Ngangao and Yale, forest canopy has multiple layers of variable tree sizes, whereas elfin forests in Vuria are of more equal tree size, and the GC value ranges of the indigenous forests are 0.42–0.71, 0.20–0.74, and 0.17–0.76, respectively. The single-species plantations of cypress and pine showed lower values of GC than indigenous forests located in the same remnants in Yale, whereas Eucalyptus plantations showed GC values more similar to the indigenous forests. These results show the usefulness of GC maps for identifying and separating forest types as well as for assessing their distinctive ecologies.


Introduction
The Eastern Arc Mountains (EAM) are a chain of crystalline Precambrian basement mountains, stretching from southern Kenya to eastern Tanzania (Burgess et al., 2007). The tropical montane forests in the EAM are important reservoirs of carbon and biodiversity (Adhikari et al., 2017;Burgess et al., 2007;Lovett and Wasser 1993). Furthermore, the montane forests capture moisture and store precipitation on the hilltops, hence having a central role in the hydrological cycle (Pellikka et al., 2009). However, montane forests of EAM, once dense and continuous stands of indigenous forest, are now threatened and fragmented due to land-use change and conversion to cropland and agroforestry, which leaves a patchy landscape of forest remnants . Depending on such factors as tree species composition and degree of forest degradation, the extant montane forest remnants have considerable differences in forest structure with effects on carbon stocks Omoro et al., 2013), biodiversity (Thijs et al., 2014), and forest function (Pfeifer et al., 2018).
Vertical structural heterogeneity describes the size variability of plant assemblages directly competing for resources at a local scale, thus T reflecting characteristics of the ecosystem's complexity such as conditions of competitive dominance of resource depletion or pre-emption in a given ecosystem (Weiner, 1990), or habitat disturbance history . The structural complexity of ecosystems provides habitat quality for fauna, affecting species diversity (Willson, 1974;Erdelen, 1984). Therefore, changes in morphological traits describing the vertical complexity of the plant ecosystem may be key to determining a threat to the overall biodiversity (MacArthur and MacArthur, 1961;Camargo and Kapos, 2009).
Airborne laser scanning (ALS; a.k.a. airborne LIDAR) provides a method for high-resolution mapping of forest three-dimensional structure via direct measurements of forest canopy height and canopy cover and through predictive modeling of forest structural attributes. When using ALS, some authors express proportions of leaf area density (LAD) along vertical strata, e.g. overstorey versus understorey (Miura and Jones, 2010;Lesak et al., 2011;Hill et al., 2014;Almeida et al., 2019). Understorey density has been found to be valuable in determining the dispersion of key animal species Jung et al., 2012;Zellweger et al., 2017). Indicators to concisely summarize plant height distributions are fundamentally based on either the entropy (structural diversity) or the equitability (size variance) of the plant size distribution (Valbuena et al., 2012).
Many studies have employed ALS derivations of McArthur and McArthur's (1961) foliage height diversity (FHD) as a morphological trait describing the complexity of the plant community (Schneider et al., 2017), showing its relationship to species distributions and biodiversity Bergen et al., 2009). On the other hand, the most widespread alternatives to express size variance are the standard deviation (McRoberts et al., 2008, Coops et al., 2016 and the Gini coefficient (GC) (Valbuena et al., 2013b(Valbuena et al., , 2017cDalponte et al., 2018;Mononen et al., 2018;Zhang et al., 2019;Erfanifard et al., 2019) of ALS heights. Descriptors of tree size variability are also good proxies for biodiversity. Goetz et al. (2007), for instance, found that the total species richness of the community was related to a ALS distribution ratio. GC is calculated based on the basal area of the individual trees, which represent the tree size inequality. Compared to other measures that can be used to evaluate tree size variability, such as the standard deviation, the GC brings the advantage of being independent from the mean tree size, and being bounded by meaningful theoretical values such as such GC = 0 representing total equality and GC = 1 maximum inequality, or GC = 0.5 representing maximum entropy (Valbuena et al., 2012). GC has been used as a proxy to measure forest structural diversity (Pach and Podlaski, 2015) and biodiversity (Nölte et al., 2018). Many studies have examined modeling GC using ALS data for assessing boreal forest structure (Goodwin et al., 2006;Valbuena et al., 2014;Valbuena et al., 2016), but the method has not previously been assessed in tropical montane forest landscapes.
Although ALS can provide highly accurate forest structure estimates, it is possible that the addition of independent phenological data that correlate with forest biophysical properties, e.g. multispectral data, can further improve its ability to predict GC. Incorporating simultaneously acquired ALS and high-resolution multispectral imagery improved species classification at tree level (Ørka et al., 2012). There has been some research on synergistic use of ALS and high spatial resolution spectral metrics from aerial imagery to predict forest structural attributes (e.g. forest canopy cover, canopy height, GC, and density) (Ahmed et al., 2014, Hudak et al., 2002, Manzanera et al., 2016, Valbuena et al., 2017b, but use is limited due to the costs involved in collecting high-resolution data. The open (free) data policy of the Landsat satellite mission has revolutionized the use of Landsat data (Zhu et al., 2019). Multispectral Landsat time-series data have been applied successfully for mapping tree height distributions in sub-Saharan Africa (Hansen et al., 2016) and above-ground biomass in tropical forests (Phua et al., 2017). However, the medium spatial resolution data, e.g. Landsat 8 OLI data, have not been tested for modeling GC to date in the Afromontane forests. Furthermore, there are no studies on combining ALS metrics and spectral-temporal metrics from the Landsat time series (LTS) for GC mapping.
In this study, our objective was to characterize the variation of forest structure in three Afromontane forest remnants in the Taita Hills. More specifically, we aimed (1) to analyse how well tree size inequality, described by GC, can be modeled by two airborne laser scanning data sets acquired from different altitudes using two different sensors; (2) to determine whether models are improved by including tree speciessensitive LTS; (3) to generate GC maps for the forest remnants and evaluate differences between the three forest remnants and indigenous and exotic plantation forests.

Study area
The Taita Hills montane forests in south-eastern Kenya are part of the EAM. The semi-arid plains isolate Taita Hills from the other hills, e.g. Pare and Usambara Mountains in Tanzania and Mount Kasigau situated 50 km to the southeast of Taita Hills in Kenya. An average elevation at the lower plains is 700 m above sea level (a.s.l.) and at the hills 1500 m a.s.l., with the highest hilltop (Vuria) reaching 2208 m a.s.l. .
The landscape is fragmented due to the conversion of forested areas to agricultural lands, while at the same time exotic plantations are established on, for example, erosion-prone areas. The forest is currently limited to small fragments of indigenous and plantation forests between 50 and 200 ha in size. The indigenous cloud forest fragments are restricted to above 1400 m altitude on the southeastern slopes and above 1700 m on the northwestern slopes, in the areas receiving over 900 mm of precipitation annually (Pellikka et al., 2009). In our study, we focused on three forest remnants: Ngangao, Yale, and Vuria ( Fig. 1), located on the mountain ridges with steep slopes at an altitude range of 1700-1952 m, 1750-2104 m, and 1655-2208 m, respectively (Fig. 2).
The study area has a bimodal rainfall annual regime, with long rains from March to June and short rains between October and December. Annual rainfall in Taita Hills is according to Erdogan et al. (2011) between 1100 and 1400 mm, while in the lowlands it is between 400 and 600 mm. The upper montane cloud forests, like Vuria, receive abundant moisture from the low-lying clouds, fog, and moisture-laden southeast trade winds originating from the Indian Ocean (Räsänen et al., 2018). These forests are typically single-layered and covered by epiphytic mosses and lichens. In contrast, the lower montane forests, like Ngangao, are drier and taller with multi-layered tree canopy (Stam et al., 2017). A large part of Yale mountain is barren rock, heathland, or rocks covered by Acacia mearnsii (hereafter wattle tree), while the steep slopes are covered by Eucalyptus spp. (hereafter eucalyptus) and Cupressus lusitanica (hereafter cypress) (Pellikka et al., 2009). Each of the three forest fragments consists of indigenous montane forest parts with several planted exotic trees (Fig. 2).
In Taita Hills, the total forest cover remained virtually unchanged between 1955 and 2004 (Pellikka et al., 2009). However, 50% of the indigenous montane forest was cleared during that period due to climatic and edaphic conditions favorable for agriculture, and because new exotic plantation forests were established (Pellikka et al., 2009). In 1955-2004for Ngangao and 1955-1994 for Yale, the change of nonforested area to the indigenous forest was 12.0 ha and 11.1 ha, and to exotic forest 81.5 ha and 94.2 ha, respectively (Pellikka et al., 2009). The indigenous forests of Taita Hills are rich in species; for example, in Ngangao 73 woody tree species were identified in the study plots by Mbuthia (2003) and 52 tree species in the plots by Schäfer et al. (2016). The most common indigenous tree species in the montane forests include Albizia gummifera, Macaranga capensis, and Maesa lanceolate. The plantation forests established between the 1950s and 1970s include stands of eucalyptus, cypress, and pine (Pinus caribea, Pinus elliottii, and H. Adhikari, et al. Ecological Indicators 108 (2020) 105739 Pinus patula), and between the 1970s and 1980s, Grevillea robusta (hereafter grevillea) and Maesopsis eminii (Adhikari et al., 2017;Pellikka et al., 2009). Exotic trees were planted for soil erosion protection and production of wood and building material (Pellikka et al., 2009). Other important canopy, subcanopy, and understorey tree species are Craibia zimmermannii, Cola greenwayi, Newtonia buchananii, Tabernaemontana stapfiana, Strombosia scheffleri, Syzygium sclerophyllum, Ochna holstii, and Macaranga conglomerate in the moist montane forest (Ngangao) and Prunus africana and Philippia pallidiflora in high altitude and upper montane forests (Yale and Vuria) (Aerts et al., 2011). The high dissimilarity in woody plant communities may be an effect of historical or recent isolation (Aerts et al., 2011), which cannot be confirmed due to the absence of detailed historical forest maps. Selective logging of the most valuable tree species, such as Ocotea usambarensis, has taken place in forests, evidenced in Ngangao by, for instance, hundreds of sawpits (Boström 2010), which have reduced the multi-layer characteristics of the forests.

Field data
The field measurements from 85 circular 0.1 ha sample plots (radius 17.84 m) collected between 2013 and 2015 were used for forest structural diversity modeling (Fig. 1). This plot size secures a robust and unbiased estimation of GC, as described in Adnan et al. (2017). The sampling of forest plot was guided by canopy height model (CHM) based on ALS data set (see Table 1 for more details) and visible to nearinfrared imaging spectroscopy data (AisaEAGLE) acquired in 2013 (Piiroinen et al., 2018). The number of plots in Ngangao, Yale, and Vuria were 37, 27, and 21, respectively. Sample plot centers were positioned using Trimble GeoXH GNSS receiver, and the differential correction were made using a GNSS base station set up at the Taita Research Station, Wundanyi. For all tree stems with a diameter (D) at breast height (1.3 m) > 10 cm, D was measured (Adhikari et al., 2017). Most of the tree species were identified by a local para-taxonomist and a field guide.
GC based on Lorenz curve quantifies tree size inequality among trees within the forest (see also Valbuena et al., 2013a) and was used as a forest structure indicator. Lorenz curves and GCs for each sample plot were calculated based on basal areas (area occupied by a given D) of individual trees . GC is calculated as the area between the line of equality and Lorenz curve divided by the total area below the line of equality. The GC ranges from 0 to 1, where 0 represents perfect equality (all trees of equal size) and 1 represents perfect inequality (few trees have the largest share) (Nölte et al., 2018). According to Valbuena et al. (2012), GC < 0.5, GC = 0.5, and GC > 0.5 represent even-sized forest, irregular forest, and bimodal diameter distributions, respectively.

Airborne laser scanning data
ALS data were collected using two sensors in three different years (Table 1). ALS1 and ALS2 are different in flight altitude, pulse density, and sensor. Both data sets were preprocessed and delivered as georeferenced point clouds in the UTM/WGS84 coordinate system with ellipsoidal heights by the data vendors Topscan Gmbh and Ramani Geosystems . Returns were classified as ground returns and used to produce a digital terrain model (DTM) with a 2 m cell grid using LAStools software (Isenburg, 2014). The ALS point cloud elevations were normalized to heights from ground level using DTM from the corresponding campaign. Noisy returns (e.g. high points) and returns from electric lines were removed manually.
ALS returns corresponding to the field plots were extracted based on plot center coordinates and fixed plot radii. ALS metrics for predictive models of forest structure attributes were calculated for each field plot using FUSION software (McGaughey, 2016). A 3-m height threshold H. Adhikari, et al. Ecological Indicators 108 (2020) 105739 was used to separate understorey and ground returns from canopy returns (Adhikari et al., 2017;Heiskanen et al., 2019). Canopy cover metrics were calculated from all returns (single, first, last), while the height metrics were calculated from the first and last returns (Appendix Table A1) (Gorgens et al., 2017).

Landsat time series
Landsat 8 Table A1). The Japan Aerospace Exploration Agency (JAXA) digital elevation model (DEM) was used to carry out topographic normalization of Landsat images (JAXA, 2015). The LTS included several percentile values (10%, 25%, 50%, 75%, and 90%), trimmed means (10% and 25%), inter-percentile range (10-90), and interquartile range (25-75) for all the bands and five vegetation indices (Adhikari et al., 2016;Potapov et al., 2012). The standard deviation of seasonal features was also calculated for 3 × 3 pixel windows centered at each plot. NDVI and RSR are robust against topographic effects (Adhikari et al., 2016). Therefore, no topographic normalization was performed for these indices.   H. Adhikari, et al. Ecological Indicators 108 (2020) 105739 with the "betareg" package (Cribari-Neto and Zeileis, 2010) and "leaps" package (Lumley, 2017). Beta regression was used for modeling GC because it is a modeling alternative well suited for these types of response variables ranging between 0 and 1 (Valbuena et al., 2013b). Furthermore, beta regression allows regression in conditions when the distributions in the response variables other than normal. A two-step predictor variable selection strategy was used due to a large number of predictor variables within each predictor group. Firstly, out of many highly inter-correlated metrics (Spearman correlation coefficient > 0.9), the one that had higher correlation with the variable of interest (here GC) was retained. Secondly, probable best predictors were identified using "regsubset" function in "leaps" package. Experiments on different link functions indicated the best performance of "loglog" link function and it was used in all the models (Cribari-Neto and Zeileis, 2010). Predictive models were calculated using ALS1 metrics only (ALS1), ALS2 metrics only (ALS2), LTS metrics only (LTS), both ALS1 and LTS metrics combined (ALS1 + LTS), and both ALS2 and LTS metrics combined (ALS2 + LTS). The best model for each combination was identified using the accuracy assessment detailed below.

Accuracy assessment
All of the GC models were assessed and compared by leave-one-out cross-validation (CV). One field plot (i) was taken out each time to eliminate the influence of that plot during model fitting, and the remaining field plots were used to predict a value of GC for that plot (pre i CV ) using beta regression. Hereafter, the abbreviations/subscripts/ superscripts "CV" and "fit" represent measures calculated using CV and non-cross-validated measures (i.e. from model fit residuals), respectively. The model coefficients were averaged from all the iterations to obtain the final model. pre CV was compared against those measured in the field (obs). The relative mean difference (rMD cv ) (Eq. (1)) was used to detect any under-or over-prediction. Prediction precision was evaluated using the cross-validated root mean square error (RMSE CV ) (Eq. (2)) and relative RMSE (rRMSE CV ) (Eq. (3)). The degree of agreement was evaluated using the coefficient of determination (R 2 ) (Eqs. (4) and (5)) . Furthermore, Piñeiro et al., (2008) hypothesis tests were used to evaluate whether the observed and predicted values follow the 1:1 correspondence line, based on the null hypothesis that the intercept ( ) and slope ( ) of the linear regression between the observed and predicted values are H0: = 0 and H0: = 1 (Eq. (6)), respectively (Valbuena et al., 2017a). If the null hypothesis for the slope is rejected, the predictions are not consistent with observed values, while if the slope hypothesis is accepted and the null hypothesis for the intercept is rejected, then the model is biased (Piñeiro et al., 2008). If both (slope and intercept) null hypotheses are rejected due to significant < 1, then we detect an averaging effect caused by over-fitting (Valbuena et al., 2017a). Finally, to avoid the degree of overfitting to the sample employed in the model fitting, the sum of squares ratio (SSR) (Eq. (7)) was used to limit the models to SSR ≤ 1.1 (Valbuena et al., 2017a;Valbuena et al., 2017b). SSR is the ratio between the square root of the residual sums of squares attained in the cross-validation (SS CV ) (Eq. (8)) and that using the whole data set (SS fit ) (Eq. (9)).

Forest structure characterization and analysis
The boundaries of the three forest remnants and different forest types within the forest remnants were mapped using a high-resolution false-colour orthomosaic and CHM based on ALS1. The orthomosaic was based on airborne hyperspectral data (visible to NIR bands, 400-1000 nm) acquired using AisaEAGLE sensor (Specim Ltd., Finland) on 3-8 February 2013 (Piiroinen et al., 2018). AisaEAGLE is a pushbroom scanner with an instantaneous field of view of 0.648 mrad and a field of view of 36.04°. The sensor produced 129 bands with an output pixel resolution of one meter. The bandwidth of each band varies between 4.5 and 5.0 nm (Piiroinen et al., 2018). Only three bands (Green (571 nm), Red (693 nm), and NIR (811 nm)) were used for mapping forest boundaries.
First, AisaEAGLE data were segmented based on the spectral information in ArcGIS 10.3 using Segment Mean Shift tool (ESRI, 2015). Spectral data, tree height (CHM), and field information were then used to identify forest boundaries and to clump the segments manually into stands of the same forest type. Furthermore, orienteering maps of Ngangao (scale 1:10 000) (Boström, 2010), a forest stand map of Vuria (scale 1:10 000) (Boström, 2013), and field plots were used for verifying stand boundaries and dominant species. Each polygon was classified as one of the following types: indigenous, eucalyptus, cypress, pine, bushland, woodland, and rock. For this study, only indigenous, eucalyptus, cypress, and pine were used because these occurred in most of the three forest remnants. Finally, we used the best models to predict GC at 30 m × 30 m grid size over Ngangao, Yale, and Vuria forests.
Based on the maps, we analyzed how GC varies between the moist H. Adhikari, et al. Ecological Indicators 108 (2020) 105739 montane and high altitude montane forests, and between indigenous forests and exotic plantations. Wilcoxon signed-rank test (one sample), Wilcoxon rank-sum tests (also known as Mann-Whitney U test), and Kruskal-Wallis test were used to assess whether the median GC in each forest type is above or below the 0.5 thresholds stated in Valbuena et al. (2012), and whether the GC values were significantly different between the forest types compared (R functions "wilcox.test" and "kruskal.test"). The non-parametric tests were used since normality could not be assumed for each forest. In the case of GC based on the maps, more detailed analyses were done for the main exotic plantation species.

Variation of Gini coefficient based on field data
The GC values calculated from the field plots clearly showed perstratum distributions other than normal (Fig. 4), and for this reason, we hereby report their medians (henceforth denoted as GC) to describe their location, and their interquartile ranges (IQR) as a descriptor for their dispersion. GC < 0.5 represents even-sized forests, GC = 0.5 represents irregular forests, and GC > 0.5 represents bimodal diameter distributions or negative exponential. Overall, tree size heterogeneity in the study area was above the GC > 0.5 threshold in Ngangao (H 0 : GC ≤ 0.5, H 1 : GC > 0.5, p < 0.001); however, we failed to prove that for Yale (GC = 0.51, IQR = 0.14) or Vuria (GC = 0.50, IQR = 0.18) when using field plot data only. When a plot had > 70% of its basal area covered by the indigenous species, it was considered pure "indigenous", whereas plots were classified as "exotic" when > 70% were exotic species (eucalyptus, cypress, and pine combined). The rest of the plots were considered "mixed". The GC and IQR (in parentheses) for the indigenous, exotic, and mixed plots were 0.58 (0.08), 0.46 (0.14), and 0.52 (0.03) in Ngangao, 0.56 (0.02), 0.44 (0.13), and 0.54 (0.11) in Yale, and 0.54 (0.12), 0.43 (0.10), and 0.41 (0.17) in Vuria, respectively.
Examples of the GC and ALS data are shown for one indigenous forest plot and one exotic plantation forest plot in Fig. 5. The indigenous forest plot has a stratum of smaller trees in DBH range of 10-20 cm along with another stratum of large trees having DBH > 20 cm (Fig. 5a). In these plots, Lorenz curves are shown simultaneously with the diameter distributions. Trees were ordered from the largest to the smallest basal area (left to right, x-axis) and cumulative share of basal area (top to bottom, y-axis) ( Fig. 5a and b). In the case of Fig. 5a showing the indigenous forest, its Lorenz curve illustrates that a smaller quantity of big trees (~30%) have a higher proportion of the total basal area (~76%), whilst a larger quantity of small trees (~70%) have a lower proportion of the total basal area (~24%). On the other hand, the plantation forest plot has a more equal DBH distribution with a single vertical stratum since trees were planted at the same time (Fig. 5b). Fig. 5c-f show the point clouds of ALS1 and ALS2 for the same plots. Due to the lower point density in ALS2, the point clouds have fewer echoes reaching the ground and less backscatter from the understorey (Fig. 5d, f). Table 2 summarizes the GC modeling results based on ALS metrics only, LTS metrics only, and their combination. In total, 23 regression models were obtained with the number of metrics ranging from one to six depending on the sensor combination (see Appendix Table A2 for  selected variables and Table A3 for model parameters). Low values of SSR (SSR < 1.1) indicate no over-fitting problems. For all models, the null hypothesis was rejected (i.e. α ≠ 0 and β ≠ 1), and the p-values were non-significant, which ensures 1:1 correspondence between the observed and predicted values.

Modeling results
The rRMSE CV for the best models fitted with ALS1 and ALS2 metrics was 13.9% and 15.1%, respectively. The model using only LTS metrics underperformed relative to those using ALS metrics. The coefficient of determination was R CV 2 = 0.28, and rRMSE CV = 16.5% was the highest uncertainty of all the alternatives tested ( Table 2). The fusion of ALS and LTS metrics improved the models in comparison with using either of them alone. R CV 2 was 0.59 and 0.55 for the models combining LTS with ALS1 and ALS2, respectively. The model combining ALS1 and LTS had the lowest rRMSE CV, 12.5%, and with ALS2 it was as low as 13.0% (Table 2).
We observed differences according to the pulse density with regard to the ALS metrics, which were selected as the best model predictors (see details in Table A2 in the Appendix). In the case of ALS1, the metrics representing canopy cover, the coefficient of variation of the first returns, L-moment coefficient of variation, the 4th L-moment, and the 5th height percentile of the first and last returns were frequently selected in the models when ALS1 was used, either alone or in combination with LTS metrics. However, in the case of ALS2, in addition to those selected in ALS1, additional ALS metrics included 3rd L-moments kurtosis and L kurtosis height of the first returns, height kurtosis, and Lmoment coefficient of variation of the last returns, skewness, and 1th and 99th height percentile for the first return.
The most frequently selected model predictors in case of LTS included, for example, percentile values (10% and 90%), and inter-percentile and interquartile range ('range10_90', and 'range25_75') of the Fig. 4. Distribution of Gini coefficient in each forest remnant and forest type based on the field plots: (a) Ngangao, (b) Yale, (c) Vuria, and (d) all plots pooled together. Indigenous class includes plots dominated by the indigenous species, exotic class includes eucalyptus, cypress, and pine plantations, and mixed class includes both indigenous and exotic species. The white dot represents the median, the thick black bar in the center represents the interquartile range, and the thin black line represents the rest of the distribution (1.5 × interquartile range). Horizontal line represents GC = 0.5.
H. Adhikari, et al. Ecological Indicators 108 (2020) 105739 vegetation indices (NDVI and RSR) and spectral bands (red, NIR, and SWIR). In the case of ALS1 + LTS, high percentile values (90%) for NIR and interquartile range of NDVI and RSR were frequently selected. In the case of ALS2 + LTS, high percentile values (90%) for NIR and NDVI and interquartile range and standard deviation of SWIR and red were frequently selected. In the models combining ALS and LTS metrics, the LTS metrics appeared in the models with more than four predictors. Thus, LTS metrics have less predictive power for the GC, but once they are incorporated with ALS, they explain additional variability. Fig. 6 compares the cross-validated GC predictions obtained when using metrics from the single sensor alone (ALS1, ALS2, LTS) and for the combination of ALS and LTS, against those observed in the field plots. LTS H. Adhikari, et al. Ecological Indicators 108 (2020) 105739 models tend to overestimate small values and underestimate large values (Fig. 6c).

Analysis of Gini coefficient maps
The best beta regression model (ALS1 + LTS) was applied to obtain spatially continuous predictions of GC for Ngangao, Yale, and Vuria (Fig. 7). The GC predictions revealed considerable differences among the three forest remnants and between the different forest types. Forest stands have even-sized, irregular, or bimodal diameter distribution of tree sizes within each forest remnant (Fig. 7). The exotic plantation forests of cypress and pine show low GC (typically < 0.5), meaning stands of even height, while indigenous and eucalyptus stands show high GC (typically > 0.5), meaning more uneven canopy. For example, a pine stand in the middle of Yale shows uniform canopy in Fig. 2b and a GC of less than 0.4 in Fig. 7, while the eucalyptus stands in Fig. 3b shows uneven canopy and a GC of more than 0.6 in Fig. 7. The uneven canopy structure of indigenous forests in Ngangao is shown in Fig. 2c, while in Fig. 7 it is represented as GC between 0.5 and 0.8. In addition, the uneven canopy structure of pine and cypress stands is seen, represented as GC between 0.1 and 0.5 in Fig. 7.
The distribution of GC in the indigenous, eucalyptus, cypress, and pine forests in each forest remnant is shown in Fig. 8. Yale showed relatively high values of GC, indicating emergent trees and several layers. On the other hand, the indigenous forest in Yale and Vuria showed GC prediction even below 0.30, indicating either a seedling stand or even-sized mature trees (Fig. 8b, c). Yale contained more eucalyptus forest than Ngangao or Vuria (Fig. 7). The eucalyptus forest showed higher values of GC than the other plantation forests (Table  A4). This is seen in Fig. 2, in which the right hand top of (a) covered by eucalyptus shows the most uneven canopy.
H. Adhikari, et al. Ecological Indicators 108 (2020)  , for both p < 0.001). This is evident since the indigenous forest is most intact in Ngangao, while in Vuria it is severely degraded due to selective logging, grazing in the forest, and planting of exotic trees within the indigenous stands.

Discussion
The results show that the GC of tree size heterogeneity can be predicted with reasonable accuracy, using ALS data and beta regression, in tropical montane forests. The resulting maps can be used to demonstrate the structural difference between forest remnants and forest types better than using field plots alone since they provide an idea of GC variation at the forest fragment level. The accuracy of the GC estimates improved significantly when using a combination of airborne LIDAR and satellite imagery (ALS + LTS), compared with predicting GC using either of these alone. Significant differences were found between forest stands of indigenous and plantation species within and between forest remnants, some of which could not be detected using the field information, which demonstrates the practical advantages of employing these data to study the ecology of the tree assemblages involved.
This was the first study carried out in the tropics showing how ALS can provide reliable and accurate predictions of GC. The accuracy of GC prediction depends on the pulse density (Adnan et al., 2017). ALS2 provided lower R CV 2 than ALS1, which can partially be explained by the lower point density as a result of the higher platform height. Furthermore, the lower frequency of returns may lead to a lower frequency of pulses detecting the understorey, underestimating the overall forest structure heterogeneity (Goodwin et al., 2006;Valbuena et al., 2017c). R 2 for GC prediction based on ALS alone in Valbuena et al., (2017b) Goodwin et al. (2006), where higher platform altitudes underestimated forest structural properties. The difference in ALS sensor and platform height might have had an impact on the final model, warranting additional investigations. Wall-to-wall prediction of GC at grid sizes other than 30 m could have provided different conclusions but was consistent with 0.1 ha field plot size. Further research should, however, implement this method in predicting tree size inequality in other forest environments to test whether the ALS and Landsat combination yields similar synergies for other forest structure attributes. The metrics derived from Landsat imagery, the LTS, alone showed little explanatory capacity for predicting forest heterogeneity. Valbuena et al. (2017b) predicted GC using a multispectral sensor (MS) alone with R 2 of 0.06, which is lower than 0.30 in our case. LTS is limited to spectral information from the top of the closed tree canopy and lacks vertical structural information. However, LTS has the potential to complement the information from the ALS in describing the structural complexity of the forest. Incorporation of spectral-temporal information from Landsat with ALS metrics improved prediction of tree size inequality. For example, with LTS metrics the inter-percentile range of NDVI and interquartile range of RSR showed significant explanatory potential for GC prediction. The key for these can be in the relationships between these remote sensing metrics and the forest density (Adhikari et al., 2016) since Adnan et al. (2017) showed that stand density is a confounding factor in the relationships between ALS metrics and the GC of tree size heterogeneity. This is possibly the reason why Manzanera et al. (2016) detected good potential in NDVI metrics for predicting GC. Our study is the first to exploit such potential with satellite imagery of coarser resolution. Also, by including LTS metrics sensitive to tree species with ALS2 metrics, GC prediction was improved relative to prediction without LTS. This shows that even with lower investments on ALS measurements and including LTS metrics, we can achieve higherR CV 2 . Valbuena et al. (2017b) achieved R CV 2 of 0.45 when combining metrics from ALS and MS (NDVI), which is lower than what we achieved (0.59) by including inter-percentile range NDVI and interquartile range RSR along with ALS metrics (Model 18 in Table A2). Our results suggest that if heterogeneity of tree sizes must be determined using ALS alone, a high point density is needed. In other cases, low-density ALS data sets could be supported by satellite data providing information about variation in forest density to obtain reliable results. If LTS metrics are included along with ALS sensor metrics, an ALS campaign with low point density may suffice to get higher R CV 2 and lower rRMSE CV . LTS data can thus compensate for lower point density. Combining ALS1 with LTS reduced rRMSE CV to 12.5%. Synergistic use of ALS and LTS has been documented for mapping biomass (Phua et al., 2017), canopy height (Hudak et al., 2002), and canopy cover (Ahmed et al., 2014), but this is the first study to map forest structure H. Adhikari, et al. Ecological Indicators 108 (2020) 105739 heterogeneity based on GC using ALS and LTS in the tropics. ALS and LTS metrics-based wall-to-wall prediction of GC provides an opportunity to visualize tree size inequality among different forest types and stands. Based on a range of GC predication, different forest stands within Ngangao, Yale, and Vuria have even-sized, irregular, and bimodal diameter distributions for indigenous and exotic plantation species. As a result, we were able to make comparisons and analyse the ecology of the species in these forests that we would have been unable to do from field information alone.
During 1955 and 2004, a total of 260 ha (50%) of indigenous tropical forest were deforested and degraded by conversion into cropland and bushland. In the same period, plantation and indigenous forest were established on the barren lands, resulting in a balanced total forest area (Pellikka et al., 2009). In 1955 indigenous forest developed in nonforested areas, and some areas were covered by the fast growth of succession species, e.g. Phoenix reclinata in abandoned open areas (Pellikka et al., 2009). Furthermore, indigenous species, e.g. Ocotea usambarensis and Podocarpus latifolius, were extracted from Ngangao for timber in the past, leaving small diameter classes (Aerts et al., 2011). For these reasons, the indigenous forest is a mixture of climax forest (undisturbed) with bimodal diameter distribution and primary successions with even-sized diameter distributions.
Our results were similar to those obtained by Valbuena et al. (2016), where old protected areas similar to the indigenous forest had the highest GC values. Eucalyptus, cypress, and pine were planted between the 1950 s and 1970 s in non-forested areas or degraded areas due to deforestation. As now mature stands were planted, it is evident that their canopy structure is of even height relative to indigenous forests, whose structures were more governed by disturbance ecology (Mbuthia, 2003). These are fast-growing species and can be an alternative to reforestation projects.
It is important to emphasize that many of these conclusions that we were able to draw from the remote sensing prediction map of GC could not be concluded using the field plot information alone. For example, at the level of plot data, the indigenous forest has bimodal diameter distribution in all forest remnants. However, when using the remote sensing prediction, we were able to prove that the indigenous forest in Ngangao is bimodal and in Vuria even-sized. These led us to the inference about their distinctive ecology and historical land use, details that would not have come to light without the assistance of remote sensing. Moreover, the exotic forest (eucalyptus) in Vuria has bimodal diameter distribution at the remote sensing level, while at plot level the exotic forest is even-sized. Furthermore, the forest of Yale has a bimodal distribution, while Vuria is even-sized, which was not possible to statistically prove based on field plot information only due to a much smaller sample size that can be employed from discrete plot samples. We thus emphasize the potential of RS for increasing the statistical power in pursuing ecological hypotheses and showing patterns on landscape scales that cannot be detected otherwise.

Conclusion
We explored the reliability of predicting the GC of tree size heterogeneity in indigenous tropical montane forests and compared the results with those obtained in plantations of eucalyptus, cypress, and pine in the same area. For this purpose, we employed two ALS data sets acquired from different altitudes using two different sensors, and also explored the potential for improving the accuracy of the GC models by incorporating a set of LTS predictors based on Landsat time series. The results showed the inclusion of LTS spectral-temporal metrics sensitive to tree species and forest density improved the modeling accuracy significantly. The inclusion of NDVI inter-percentile range and RSR interquartile range along with ALS metrics had the highest R CV 2 among the models. Finally, we predicted maps of GC for the three forest remnants using the best model and compared forest structural heterogeneity among and within the forest types. Results showed that these maps can be employed to demonstrate ecological hypotheses that would not be plausible using field data alone. Using these GC maps, we observed a high degree of structural heterogeneity in these montane forest remnants. Based on the ranges of GC values observed in these maps, we were able to identify different areas within these forest remnants having even-sized, irregular, and bimodal diameter distributions for indigenous and exotic species. For modeling GC, high point density ALS is required if ALS metrics is used alone. Nonetheless, low scanning densities may suffice when combining the ALS data with LTS metrics. The method showed the potential to discriminate between indigenous and plantation forests.
Modeling results Table A.1 explains the predictor metrics used in modeling GC. Fig. A.1 shows R CV 2 for GC prediction using a different number of variables. The highest R CV 2 obtained using six metrics together from ALS1 and LTS metrics was 0.59. Maximum predictors in ALS1 were four and in ALS2 five; the R CV 2 for four metrics in ALS1 are still higher than five metrics in ALS2 (Fig. A.1, Table A.2). Table A.2 shows all of the models with different numbers of metrics possible in each sensor combination or alone. The estimates of the model coefficients (beta distribution parameters, and ) are given in Table A.3. In the final choice of models, we selected those with lower uncertainty while securing their significance (both in model fit - Table A.2and observed versus predicted -Eq. (6)) and avoiding overfitting to the sample (SSR, Eq. (7)). To avoid overfitting of the models, models with SSR greater than 1.1 were removed. Table A.4 presents the summary statistics of different forest pixels (indigenous, eucalyptus, cypress, and pine) in the three forest remnants.