UAV Remote Sensing for Biodiversity Monitoring: Are Forest Canopy Gaps Good Covariates?

Forest canopy gaps are important to ecosystem dynamics. Depending on tree species, small canopy openings may be associated with intra-crown porosity and with space among crowns. Yet, literature on the relationships between very fine-scaled patterns of canopy openings and biodiversity features is limited. This research explores the possibility of: (1) mapping forest canopy gaps from a very high spatial resolution orthomosaic (10 cm), processed from a versatile unmanned aerial vehicle (UAV) imaging platform, and (2) deriving patch metrics that can be tested as covariates of variables of interest for forest biodiversity monitoring. The orthomosaic was imaged from a test area of 240 ha of temperate deciduous forest types in Central Italy, containing 50 forest inventory plots each of 529 m2 in size. Correlation and linear regression techniques were used to explore relationships between patch metrics and understory (density, development, and species diversity) or forest habitat biodiversity variables (density of micro-habitat bearing trees, vertical species profile, and tree species diversity). The results revealed that small openings in the canopy cover (75% smaller than 7 m2) can be faithfully extracted from UAV red, green, and blue bands (RGB) imagery, using the red band and contrast split segmentation. The strongest correlations were observed in the mixed forests (beech and turkey oak) followed by intermediate correlations in turkey oak forests, followed by the weakest correlations in beech forests. Moderate to strong linear relationships were found between gap metrics and understory variables in mixed forest types, with adjusted R2 from linear regression ranging from 0.52 to 0.87. Equally strong correlations in the same forest types were observed for forest habitat biodiversity variables (with adjusted R2 ranging from 0.52 to 0.79), with highest values found for density of trees with microhabitats and vertical species profile. In conclusion, this research highlights that UAV remote sensing can potentially provide covariate surfaces of variables of interest for forest biodiversity monitoring, conventionally collected in forest inventory plots. By integrating the two sources of data, these variables can be mapped over small forest areas with satisfactory levels of accuracy, at a much higher spatial resolution than would be possible by field-based forest inventory solely.


Introduction
Forest canopy gaps are regarded as hotspots that provide ideal conditions for rapid plant reproduction and growth, maintenance of floristic richness in the understory, and an increase of diversity and structural complexity of the forest habitat [1].Gaps dynamics, associated with the fall of individual or clumped canopy trees, are observed in boreal, subalpine, temperate, and tropical forest types worldwide, especially in late stages of forest development [2][3][4][5].In managed forests, exploiting gaps' dynamic is at the root of gap-based close-to-nature silviculture [6,7].
Small canopy openings may be attributed to tree crown architecture, namely, intra-crown porosity and the space among crowns.Even in fully stocked forest areas, some tree species display so-called 'canopy shyness', characterized by a canopy cover with inter-crown gaps.However, little is known about the relationships between this fine-scale pattern of canopy openings, the sunlight penetration into the canopy, and its impact on the forest understory layer.
Remote sensing from unmanned aerial vehicles (UAV) allows the detection of small size canopy gaps and their spatial patterns [8].Since the optimal spatial resolution of the sensor is suggested to be five times smaller than the monitored object [9], the new versatile imaging platforms offer the potential to explore canopy features at a very high spatial resolution.UAVs have many benefits in forest remote sensing, such as the flexibility of data acquisition with low material and operational costs, and high-intensity data collection [10].However, the use of UAVs is hindered by a set of challenges such as coverage of small areas per flight due to low flight altitude, limited battery capacity, and existing regulatory requirements to keep the UAV in visual line-of sight (VLOS) of the pilot.Nevertheless, lightweight fixed-wing UAVs can cover a larger area with one flight (up to 150 ha) than multicopter drones.Furthermore, the use of multiple batteries on lightweight fixed wing UAVs may enable aerial coverage of up to 1000 ha in one working day in optimal weather conditions [11].
UAV 3D image-based point clouds, normalized by high resolution Digital Terrain Model (DTM), have been explored for the estimation of forest structural variables (e.g., tree density, growing stock, or above-ground biomass) as a cheaper alternative to airborne laser scanning [13,18,19].In contrast, UAV orthomosaic applications cover the detection of canopy openness [20], tree species identification [21], and assessment of infestation and damages at tree level [22][23][24][25].Most studies have focused on wall-to-wall fine scale mapping of forest canopy attributes over small forest areas [13,[26][27][28][29][30][31][32].As an alternative, however, using UAV orthomosaic images to extract covariates has been less explored.Covariates, such as canopy gap-related metrics, are potentially correlated with compositional and structural biodiversity variables.Indeed, the definition of canopy gap remains unclear and inconsistent in the literature [33,34].
Although it is known that gap size and shape considerably influence gap microclimate [1,35], unambiguous indication of thresholds for these gap metrics is lacking [35].Forest canopy gaps depend on the geographical context and the forest composition, in addition to the ability of detecting those gaps.Bonnet et al. [36] defined canopy gaps as 'openings in the canopy with a minimum area of 50 m 2 , a minimum width of 2 m and a maximum vegetation height of 3 m'.Getzin et al [37], adopting a different definition, mapped canopy gaps of 1 m 2 from a true colors UAV product of 7 cm spatial resolution in beech-dominated deciduous and mixed deciduous-coniferous forests in Germany.Furthermore, Getzin et al. [8] extracted forest gaps and their spatial pattern from ten different plots of 1 ha.The authors recommend collecting aerial data in a cloudy condition in order to reduce the effect of the shadow that could lead to the misinterpretation of dark pixels as gaps.
Likewise, gap age is important, since canopy openings fill with time, and biodiversity varies between newly opened gaps versus older gaps [1].In addition, according to Muscolo et al. [6], canopy gap age and size are the primary factors affecting the regeneration of understory besides the suitable substrate.Techniques and packages developed for manned airborne and satellite-borne images hardly suit UAV imagery classification for gap mapping because of different data acquisition parameters [38,39].Although some studies attempted to overcome some of the challenges [8,20,37], they used methods rarely available at low cost or suitable for forest managers.
Since lightweight fixed-wing UAVs equipped with a digital commercial camera with red, green, and blue (RGB) bands are most common in forest remote sensing [12,13], testing their use for exploring correlations between canopy openness and the occurrence of key forest biodiversity features in the understory and overstory is an interesting opportunity.This study explores this potential correlation in a small test site with temperate deciduous forests in central Italy.An integrated survey combining RGB imagery acquisition by UAV with field-based forest inventory was conducted at an experimental site in 2016 in the framework of the LIFE project "FREShLIFE-Demonstrating Remote Sensing integration in sustainable forest management".The subsequent data were processed in the present study.Our objectives were to explore relationships between canopy gaps and certain variables describing the structure and diversity of forest stands and of the understory.Specifically, we examined the following questions: (1) Are there adequate image processing techniques to delineate openings in the canopy cover from UAV imagery?
(2) Are patch metrics of these canopy features correlated with structural and biodiversity-related variables of the forest understory or of the forest stand?Does this correlation vary across stands characterized by different dominant canopy species?

Study Site
The test site extends over 240 hectares on the slopes of Mount Venere (500-800 m a.s.l.) in Caprarola Municipality (Central Italy).In the study, there were three forest types (Figure 1): beech-dominated forest (Fagus sylvatica, 149 ha), turkey oak forest (Quercus cerris, 76 ha), and mixed forest of the two species (14 ha).The vertical profile is, to a large extent, mono-layered in the beech forest and bi-layered in the other two forest types.
The presence of very fertile volcanic soils and the microclimate as a result of the presence of a lake create favorable conditions for the growth of beech, even at a low elevation.In fact, the conformation of the basin, the frequent formation of mists, the high relative air humidity, and the protection from extreme winds ensure a suitable habitat for this species.As a result, the beech forest of Monte Venere is ecologically important, as it grows at much lower altitudes (around 500 m above sea level) than those usually occupied by beech in the Central Apennines (optimum between 1000 and 1700 m a.s.l.).
Beech forest of Monte Venere belongs to the Aquifolio-Fagetum association, which reaches the extreme northern boundary of its ecological range in the Province of Viterbo [40].Turkey oak becomes the dominant species on the south-facing aspects.Turkey oak forests also contain accompanying species such as ash (Fraxinus ornus L.) and hornbeam (Ostrya carpinifolia L.).Turkey oak forests present a well-developed shrub layer including Rosa canina L., Cornus sanguinea L., Crataegus monogyna Jacq., Crataegus oxychanta L., and Ruscus aculeatus L. The test site contains 50 forest inventory plots, each 23 m × 23 m in dimension (area of the plot 529 m 2 ), distributed according to a one-per-stratum systematic sampling design.This sampling scheme partitions the study area into square grids of 529 m 2 plots spatially aggregated into 50 equal size strata.One sample-site location is independently and randomly selected in each stratum, resulting in the spatial distribution displayed in Figure 1.There were 28, 13, and 9 plots in beech, turkey oak, and mixed forests, respectively.

Field Data Processing
The spatial position of the center of inventory plots was recorded with a sub-meter accuracy Trimble GEO7X GNSS receiver and was post-processed with data from local base stations situated approximately 30 km from the study area.At each plot, all plants (trees and shrubs) with a diameter at breast height (DBH) greater than 2.5 cm were inventoried during the 2016 growing season.Key measurements included species identification, diameter at breast height (DBH), tree height, and the presence of microhabitat-bearing living trees according to the reference catalogue provided by Winter and Möller [41].These elements, which are characteristic of natural forests, especially of the old-growth phases, are often absent or rare in managed forests and are, therefore, considered key indicators of biodiversity.
Tree and shrub data were processed in the present study to distinguish the understory species from canopy trees (living trees).The understory layer was identified as the layer including plants with heights up to 50% of the maximum tree height observed in the plot [42].In each plot, we calculated for the understory: (1) density and development, by means of the parameters such as number of plants (N_PLANTS), mean DBH (MEAN_DBH), mean plant height (MEAN_HTOT), total basal area (G_TOT), and total understory volume (V_TOT), which is calculated based on species-specific allometric equations developed using a combination of height and diameter; (2) species diversity, by means of number of species (N_SPECIES), Shannon's index (I_SHANNON, H') [43], and Pielou index (I_PIELOU) [44].
Data collected on living trees were processed to derive structural data (mean DBH, MEAN_DBH; mean total height, MEAN_HTOT) and biodiversity data such as number of habitat trees (HAB), percentage of habitat trees (%HAB) of total number of plants in the plot, vertical species profile by the index of Pretzsch (I_PRETZSCH), Shannon's index (I_SHANNON), and Margalef index (I_MARGALEF) [45].Table 1 displays the summary of biodiversity indices computed for each plot.The test site contains 50 forest inventory plots, each 23 m × 23 m in dimension (area of the plot 529 m 2 ), distributed according to a one-per-stratum systematic sampling design.This sampling scheme partitions the study area into square grids of 529 m 2 plots spatially aggregated into 50 equal size strata.One sample-site location is independently and randomly selected in each stratum, resulting in the spatial distribution displayed in Figure 1.There were 28, 13, and 9 plots in beech, turkey oak, and mixed forests, respectively.

Field Data Processing
The spatial position of the center of inventory plots was recorded with a sub-meter accuracy Trimble GEO7X GNSS receiver and was post-processed with data from local base stations situated approximately 30 km from the study area.At each plot, all plants (trees and shrubs) with a diameter at breast height (DBH) greater than 2.5 cm were inventoried during the 2016 growing season.Key measurements included species identification, diameter at breast height (DBH), tree height, and the presence of microhabitat-bearing living trees according to the reference catalogue provided by Winter and Möller [41].These elements, which are characteristic of natural forests, especially of the old-growth phases, are often absent or rare in managed forests and are, therefore, considered key indicators of biodiversity.
Tree and shrub data were processed in the present study to distinguish the understory species from canopy trees (living trees).The understory layer was identified as the layer including plants with heights up to 50% of the maximum tree height observed in the plot [42].In each plot, we calculated for the understory: (1) density and development, by means of the parameters such as number of plants (N_PLANTS), mean DBH (MEAN_DBH), mean plant height (MEAN_HTOT), total basal area (G_TOT), and total understory volume (V_TOT), which is calculated based on species-specific allometric equations developed using a combination of height and diameter; (2) species diversity, by means of number of species (N_SPECIES), Shannon's index (I_SHANNON, H') [43], and Pielou index (I_PIELOU) [44].
Data collected on living trees were processed to derive structural data (mean DBH, MEAN_DBH; mean total height, MEAN_HTOT) and biodiversity data such as number of habitat trees (HAB), percentage of habitat trees (%HAB) of total number of plants in the plot, vertical species profile by the index of Pretzsch (I_PRETZSCH), Shannon's index (I_SHANNON), and Margalef index (I_MARGALEF) [45].Table 1 displays the summary of biodiversity indices computed for each plot.The Shannon index expresses the frequency of the i-th species in a community; its values generally lie between 0 and 3.5; higher values correspond to higher species diversity.Its maximum value (MAX_SHANNON) is given by the natural logarithm of the number of species found in the test area and occurs when all species are equally present.

Pielou index (E)
H / ln(S) [0,1] The Pielou index measures the relative abundance of species groups.The index can take values between 1 (all species are equally abundant) and 0 (there is only one species).
The Pretzsch index summarizes and quantifies species diversity and the vertical distribution of species in a forest.The index is lowest in one-story pure forests, whereas it rises for pure forests with two or more stories.Peak values are reached in mixed woodlands with heterogeneous structures.
It quantifies the presence of a number of species in a community.It also depends on the number of plants found in the sampling area.The index value increases with increasing species diversity.
S = number of species; N = total number of plants; pij = the frequency of species i in the layer j; pi = frequency of the specie i. z = number of layers.

UAV Data
Aerial images of the study area were collected near peak greenness during two consecutive days of the last week of May 2016 by a fixed-wing eBee UAV (SenseFly, Cheseaux-Losanne, Switzerland).The eBee was equipped with a commercial SONY WX 18.MP RGB camera, recording in blue (450 nm), green (520 nm), and red (660 nm) wavelengths.The eBee, a hand-launched fixed-wing UAV with an electric motor-driven pusher propeller, has a 96 cm wingspan and a weight of about 700 g including camera, inertial measuring unit, GPS, and battery payloads.The maximum area covered in a single flight is about 140 ha at 120 m altitude or 9 km 2 at 2,000 m altitude, with a maximum flight time of 45 min.
Before UAV image acquisition, 15 ground control points (GCPs) were placed in open areas on the ground using 50 × 50 cm targets with a black and white chalkboard pattern to ensure high visibility in the images.The coordinates of the center of the GCPs were recorded with a Trimble Geo 7X receiver, with a 2-s logging rate for approximately 15 min each.The recorded coordinates were post-processed with data from the nearest ground station, located at 30 km from the study area, using the Pathfinder software.The post-processed GCP coordinates resulted in an average standard deviation from north and east and heights of 0.9 m, 0.7 m, and 1.9 m, respectively.
The flight altitude was 145 m above ground level, and the flight lines were set to obtain a longitudinal overlap of 85% and lateral of 75%.The total flight time was 2 h and 43 min among eight flights, launched from two different areas.Flight line spacing was 42 m, and the distance between two adjacent photos was 34.6 m.The focal length of the camera was set to 4 mm, and the ISO sensibility was ISO-100 with a shutter speed of 1/2000 s.A total of 433 images was acquired with a field of view of 200 m × 150 m.
Visual inspection of the acquired images revealed no problems related to light and atmospheric conditions or saturation.Orthomosaic images were derived from the UAV images using the SfM photogrammetry software Agisoft PhotoScan Pro [46].This software allows one to fully automate the photogrammetric workflow to process aerial images and produce 3D (e.g., image-based point cloud) and 2D (e.g., Digital Surface Model and Ortomosaic) models, thanks to the multi-view stereo-matching and SfM algorithm implemented.This software has been already used for forest analysis [11,13,[47][48][49][50], and the UAV images were processed using the following steps: (a) image alignment, (b) mesh building, (c) guided marker positioning and optimization of camera alignment (georeferencing of created scene), (d) dense cloud building, (e) raster grid DSM generation with a resolution of 0.2 m × 0.2 m, and (f) building orthomosaic.From the SfM photogrammetric workflow, we obtained a raster grid orthomosaic with a spatial resolution of 10 cm.

Image Processing and Variable Selection
We tested different settings of the Contrast Split segmentation algorithm to map forest canopy openings from the RGB orthomosaic images.The red band was the only one that significantly changed with openings in the canopy cover and the crowns.The Contrast split algorithm is considered effective in separating dark objects from bright objects [51].The purpose of testing different settings was to extract forest canopy openings as accurately as possible, which are expected to appear darker, compared to surrounding canopy pixels.This process was accomplished using eCognition Developer (http: //www.ecognition.com)software.We validated the gap delineation from contrast split segmentation by visual interpretation of the orthomosaic image to identify omission or commission errors.The visual validation techniques have been often used in forest canopy gap mapping [33].We later tested different thresholds to filter small openings that might not affect ecological phenomena under investigation and, therefore, constitute a noise masking of relevant data.We only considered mapped features falling within the plot boundary and gaps within the plot boundary for at least 50% of their area.For each feature, we calculated the shape and size indices and other gap patch metrics commonly used in landscape-scale analyses.In total, 19 patch metrics (see Supplements; more details in [52]) were processed for each gap in each plot.This suite of metrics reflects different nuances of two-dimensional gap properties that may be potentially important for linking image-detected higher-level structures to dependent lower-level processes of the biota [37].
The size or extent gap patch metrics are the area (A), length, border length (or the perimeter), ratio of length/width, and the width.The shape metrics include border index (B.Index), asymmetry (Asym), roundness (Round), compactness (Comp), shape index, density, rectangular fit (Rect.fit),radius of the largest enclosed ellipse (RLE), radius of the smallest enclosed ellipse (RSE), elliptic fit, and other composite indices such as the shape complexity index (GSCI).The gap shape complexity index is an important measure of forest gaps [53].It is the ratio of a gap's perimeter to the perimeter of a circular gap of the same area.A value of 1 describes a perfect circle, while increasing values indicate increasing shape complexity.For example, values of 1.20 and 2.70 have 20% and 170% complexity, respectively.The last three gap metrics are the patch fractal dimension (PFD) [54], the fractal dimension (FD) [55], and the fractal dimension index (FDI) [56].A description of each metric and its range of variation is given in Table A1 (Appendix A).
For each patch metric, we calculated the median (mdn), mean (avg), standard-deviation (SD), sum (SUM), and the coefficient of variation (cv) per plot.We tested gap metrics with two different gap size thresholds, namely, greater than 1 m 2 and greater than 2 m 2 .This leads to two times 95 variables for each plot.These two thresholds were set, because some ecological phenomena, such as understory structure dependencies, are only quantifiable if small gaps are taken into account, yet very small gaps may not affect at all the lower dependency phenomenon and therefore constitute noise [57].

Statistical Analysis
First, we performed exploratory analysis with the field data using ANOVA (Analysis of Variance) tests when the distribution was normal, and KRUSKAL-WALLIS analysis otherwise with the categorical variable being the forest type characterizing each single plot.Statistical analyses were performed with R Foundation 3.2.3(https://www.r-project.org/foundation/).The code used in this paper is available at https://github.com/MartinBagaram/forest_gaps.
We used all gap variables for Pearson's product-moment correlation analysis and Spearman's rank correlation with the field data, using post-stratification of sampling units by forest types where appropriate.Usually, the two correlations give approximately the same results (see [58]).To avoid the fact that the observed correlation, although significant according to the p-value lower than 5%, could have occurred by chance, its significance was further tested using the 9999 permutation method proposed by Legendre and Legendre [59].This method has been successfully used in investigating the relationship between forest canopy gap metrics and biodiversity [37].
For gap variables correlated with field data, we tested the correlation among significant potential predictors.This third step allowed a decrease in the redundant information brought by potentially correlated predictors and eliminated the collinearity or singularity.To determine the most significant variables to be used as possible predictors, we used a forward stepwise analysis, since this approach selects only the minimum significant number of predictors.There exist three methods for model selection known as best, forward, and backward subset selection.In general, the three methods lead to similar though not identical results [60].Owing to the huge number of predictors (190 in this study), the best subset selection requires for each field parameter to compare up to two of 190 models, which deterred us from using the best set selection.
When the forward stepwise yielded only one predictor variable as sufficient to model the field parameter (which was often the case), we compared the correlation coefficients of Spearman and Pearson.If these two coefficients were significant and close to one another, then the patch metric was selected as a predictor if it passed the 9999 permutations test of Legendre.Otherwise, the second variable with the highest values of both Spearman and Pearson coefficients was selected, and then we repeated the 9999 permutations method.This approach has the advantage of, on one hand, excluding variables depicting high Pearson's correlation due to extreme values, and on the other hand, not solely relying on the Spearman's correlation, which is a measure of monotony [61] and not of linearity.Furthermore, using 9999 permutation assures that the p-value observed does not occur by chance.
Finally, we performed linear regression tests using the selected significant variables.Following Møller and Jennions [62] and Getzin et al. [37], who considered a coefficient of determination R 2 > 0.25 as meaningful, since the predictor value leads to a great change if it explains over 25% of variance, we used a threshold of R 2 > 0.5.We then validated the regression by checking regression quality assumptions such as the normality of residuals, the homoscedasticity of residuals, and that the mean of the residuals is zero.For forest parameters that can be predicted with R 2 greater than 50%, we generalized the model to the whole forest type extent using the grid of 529 m 2 plots.We performed the validation with cross-validation (Leave-One-Out Cross-Validation), which produced root mean squared errors (RMSE).
For inference, it is important to provide the variance, confidence interval, and bias of R 2 .Bootstrapping is a method commonly used for that purpose [63][64][65], especially in the case of linear regression [66].Bootstrap is a resampling method developed by Efron and Tibshirani [67].We computed those statistics for variables with R 2 greater 0.50 using formulas given by Mura et al. [63].Although it is suggested that the bootstrap sample size should be big enough (greater than 200), there is no consensus yet on what should be the actual size.Following Mura et al. [63], we set the bootstrap sample size to 500.The mapping of biodiversity attributes was achieved with ArcGIS 10.2.A graphical abstract of the methodology is presented in Figure 2.

Canopy Gaps Mapping
The contrast split algorithm based on the red band accurately differentiated dark objects (shaded canopy gaps) from bright objects, which, in most cases, corresponded to forest canopy.The mapping faithfully delineated shaded canopy gaps but performed poorly in illuminated gaps when the bare soil was apparent (Figure 3).Detected gaps from plots ranged from 100 cm 2 (1 pixel) to 122 m 2 .This range encompassed small openings or inter-crown cracks in the canopy cover and larger gaps generated by the fall of one or more canopy trees.The gap size distribution is roughly the same in the three forest types.When considering gaps greater than 1 m 2 , over 75% of gaps in the three forest types were less than 5 m 2 .When considering gaps greater than 2 m 2 , over 75% of gaps in the three forest types were less than 7 m 2 (Figure 3).
The collinearity analysis of gap patch metrics indicated that patch metrics were strongly correlated to each other.A sample result of collinearity analysis of predictor variables is reported in Figure 4.The lowest correlations were obtained with the coefficient of variation, whereas the highest correlations were associated with the plot-level sum of patch metrics.These findings indicate that gap patch metrics are not suitable for a multiple linear regression.

Canopy Gaps Mapping
The contrast split algorithm based on the red band accurately differentiated dark objects (shaded canopy gaps) from bright objects, which, in most cases, corresponded to forest canopy.The mapping faithfully delineated shaded canopy gaps but performed poorly in illuminated gaps when the bare soil was apparent (Figure 3).Detected gaps from plots ranged from 100 cm 2 (1 pixel) to 122 m 2 .This range encompassed small openings or inter-crown cracks in the canopy cover and larger gaps generated by the fall of one or more canopy trees.The gap size distribution is roughly the same in the three forest types.When considering gaps greater than 1 m 2 , over 75% of gaps in the three forest types were less than 5 m 2 .When considering gaps greater than 2 m 2 , over 75% of gaps in the three forest types were less than 7 m 2 (Figure 3).
The collinearity analysis of gap patch metrics indicated that patch metrics were strongly correlated to each other.A sample result of collinearity analysis of predictor variables is reported in Figure 4.The lowest correlations were obtained with the coefficient of variation, whereas the highest correlations were associated with the plot-level sum of patch metrics.These findings indicate that gap patch metrics are not suitable for a multiple linear regression.

Correlation between Gap Metrics and Understory Variables
The understory layer exhibits density, development, and diversity features that seem to vary with forest types cover.Namely, variables related to density, development, and diversity of the understory had significant differences amongst the forest types, except for Pielou index, MEAN_HTOT, and MEAN_DBH (Table 2).
All eight field parameters presented significant correlations (p < 0.05), with gap patch metrics in at least one of the three forest types (Table 3).Tables B1-B3 of Appendix B present the threshold, Pearson's, and Spearman's correlations of gap patch metrics and understory.The regression coefficient had the lowest standard error in mixed and Fagus forests compared to Quercus forest; MEAN_DBH in mixed forest had a very low variability as shown by the coefficient of the regression predicting MEAN_DBH in mixed forest is merely equal to zero.
As shown in Table 3, the highest coefficient of determination was in mixed forest (MEAN_HTOT, R 2 = 0.87, p < 0.000) and the lowest was in Fagus forest (N_SPECIES, R 2 = 0.15, p < 0.05).The best results were found in Quercus and mixed forests, whereas the worst ones were in Fagus forest.I_SHANNON and MEAN_DBH in Quercus forest, and MEAN_HTOT and N_PLANTS in mixed forest were

Correlation between Gap Metrics and Understory Variables
The understory layer exhibits density, development, and diversity features that seem to vary with forest types cover.Namely, variables related to density, development, and diversity of the understory had significant differences amongst the forest types, except for Pielou index, MEAN_HTOT, and MEAN_DBH (Table 2).
All eight field parameters presented significant correlations (p < 0.05), with gap patch metrics in at least one of the three forest types (Table 3).Tables B1-B3 of Appendix B present the threshold, Pearson's, and Spearman's correlations of gap patch metrics and understory.The regression coefficient had the lowest standard error in mixed and Fagus forests compared to Quercus forest; MEAN_DBH in mixed forest had a very low variability as shown by the coefficient of the regression predicting MEAN_DBH in mixed forest is merely equal to zero.
As shown in Table 3, the highest coefficient of determination was in mixed forest (MEAN_HTOT, R 2 = 0.87, p < 0.000) and the lowest was in Fagus forest (N_SPECIES, R 2 = 0.15, p < 0.05).The best results were found in Quercus and mixed forests, whereas the worst ones were in Fagus forest.I_SHANNON and MEAN_DBH in Quercus forest, and MEAN_HTOT and N_PLANTS in mixed forest were

Correlation between Gap Metrics and Understory Variables
The understory layer exhibits density, development, and diversity features that seem to vary with forest types cover.Namely, variables related to density, development, and diversity of the understory had significant differences amongst the forest types, except for Pielou index, MEAN_HTOT, and MEAN_DBH (Table 2).
All eight field parameters presented significant correlations (p < 0.05), with gap patch metrics in at least one of the three forest types (Table 3).Tables A2-A4 of Appendix B present the threshold, Pearson's, and Spearman's correlations of gap patch metrics and understory.The regression coefficient had the lowest standard error in mixed and Fagus forests compared to Quercus forest; MEAN_DBH in mixed forest had a very low variability as shown by the coefficient of the regression predicting MEAN_DBH in mixed forest is merely equal to zero.
As shown in Table 3, the highest coefficient of determination was in mixed forest (MEAN_HTOT, R 2 = 0.87, p < 0.000) and the lowest was in Fagus forest (N_SPECIES, R 2 = 0.15, p < 0.05).The best results were found in Quercus and mixed forests, whereas the worst ones were in Fagus forest.I_SHANNON and MEAN_DBH in Quercus forest, and MEAN_HTOT and N_PLANTS in mixed forest were predicted with R 2 > 0.50.All patch metrics, except one (cv_Length in Fagus forest), correlated with the field parameters were related to gap shape.The spatial patterns of the understory models to the extent of their respective forest types underscored some key features (Figure 5).In Quercus forest, I_SHANNON ranged from just above 0 to 1.8 with RMSE of 0.21.These values are in the range of variation of data collected on plots.The variability in MEAN_DBH in Quercus forest was from 2 to 17 cm with an RMSE = 1.3 cm.
In mixed forest, MEAN_HTOT went up to 11 m, with the majority of grids falling between 5 and 9 m, and the RMSE was 0.43.The number of plants per grid ranged from less than 5 to 33 with an RMSE of 3.84 (Figure 5).We did not perform the spatialization of MEAN_DBH in mixed forest, even though the R 2 was greater than 0.5, because the linear model slope was close to zero.The empty cells in the figure correspond to cells in which the prediction gives unreasonable values.
The spatial patterns of the understory models to the extent of their respective forest types underscored some key features (Figure 5).In Quercus forest, I_SHANNON ranged from just above 0 to 1.8 with RMSE of 0.21.These values are in the range of variation of data collected on plots.The variability in MEAN_DBH in Quercus forest was from 2 to 17 cm with an RMSE = 1.3 cm.
In mixed forest, MEAN_HTOT went up to 11 m, with the majority of grids falling between 5 and 9 m, and the RMSE was 0.43.The number of plants per grid ranged from less than 5 to 33 with an RMSE of 3.84 (Figure 5).We did not perform the spatialization of MEAN_DBH in mixed forest, even though the R 2 was greater than 0.5, because the linear model slope was close to zero.The empty cells in the figure correspond to cells in which the prediction gives unreasonable values.

Correlations between Gap Metrics and Living Tree Variables
All living tree diversity and structural variables exhibited significant differences among the three forest types, except total basal area (G_TOT) and total volume (V_TOT) (Table 4).As with the understory, some gap patch metrics appear to be significantly correlated with these field parameters (p < 0.05).Tables C1-C3 in Appendix C summarize the threshold, Spearman's, and Pearson's correlations associated with gap patch metrics and living trees parameters.In Quercus forest, eight field parameters yielded good results, compared to mixed forest which had four.For Fagus forest, habitat trees and percentage of habitat trees only denoted significant correlations.Coefficients of the linear regression had the lowest standard errors in Quercus forest.

Correlations between Gap Metrics and Living Tree Variables
All living tree diversity and structural variables exhibited significant differences among the three forest types, except total basal area (G_TOT) and total volume (V_TOT) (Table 4).As with the understory, some gap patch metrics appear to be significantly correlated with these field parameters (p < 0.05).Tables A5-A7 in Appendix C summarize the threshold, Spearman's, and Pearson's correlations associated with gap patch metrics and living trees parameters.In Quercus forest, eight field parameters yielded good results, compared to mixed forest which had four.For Fagus forest, habitat trees and percentage of habitat trees only denoted significant correlations.Coefficients of the linear regression had the lowest standard errors in Quercus forest.
The highest coefficient of determination was in mixed forest (HAB, R 2 = 0.79, p < 0.000) and the lowest was in Fagus forest (%HAB, R 2 = 0.11, p < 0.05).Best results were in Quercus forest (I_PRETZSCH, HAB, and %HAB) and mixed forest (HAB, %HAB, MEAN_DBH, and MEAN_HTOT).HAB, %HAB, MEAN_HTOT in mixed forest and HAB and I_PRETZSCH in Quercus forest had R 2 exceeding 0.50 (Table 5).Among the patch metrics correlated with living tree parameters, all but the length (cv_Length in mixed forest), which is gap size related patch metric, were gap shape-related patch metrics.The spatial patterns of regression models for living trees parameters are presented in Figure 6.In Quercus forest, I_PRETZSCH varied from 1 to 2.4 with RMSE of 0.40.A majority of cell grids had indexes higher than 1.7, indicating a complex vertical structure.In the same forest type, the number of habitat trees reached 20 with RMSE of 3.89.
In mixed forest, however, most of cell grids had the number of habitat trees higher 9.In this forest, the prediction error is relatively smaller with an RMSE of 1.6.The percentage of habitat trees per grid is more variable with values ranging from 0 to 100%.The MEAN_HTOT varied from 6 to 28 m with an RMSE of 1.86.

Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 27
The spatial patterns of regression models for living trees parameters are presented in Figure 6.In Quercus forest, I_PRETZSCH varied from 1 to 2.4 with RMSE of 0.40.A majority of cell grids had indexes higher than 1.7, indicating a complex vertical structure.In the same forest type, the number of habitat trees reached 20 with RMSE of 3.89.
In mixed forest, however, most of cell grids had the number of habitat trees higher than 9.In this forest, the prediction error is relatively smaller with an RMSE of 1.6.The percentage of habitat trees per grid is more variable with values ranging from 0 to 100%.The MEAN_HTOT varied from 6 to 28 m with an RMSE of 1.86.

Quality Assessment
The bootstrap R 2 was relatively close to the experimental R 2 (computed from the actual data) for both the understory and the living trees.The standard error and the bias are very low.The highest negative bias was related to I_SHANNON with a bias of −0.05.This suggests that the actual R 2 of I_SHANNON is lower than the one observed.Similarly, the highest positive bias was recorded with HAB in mixed forest (living trees with a value of 0.03)suggesting that the actual R 2 associated with habitat trees in mixed forest is higher than 0.79 (Table 6).

Mapping Forest Canopy Gaps
Gap mapping from UAV RGB 10 cm orthomosaic images produced very promising results.In this study, we mapped openings in canopy as small as 100 cm 2 , corresponding either to intra-crown openings or to inter-crown cracks in the canopy cover.This fine-scale mapping of canopy openings is important, because some ecological phenomena, such as understory structure dependencies, are only quantifiable if small gaps are taken into account.Though it is known that different forests, depending on ecological and geographical contexts, exhibit different canopy structure, there is not yet a rule for a minimum gap size.As a consequence, different authors set, arguably, minimum gap size of 1 m 2 [8,37,57,[68][69][70], 5 m 2 [71], 10 m 2 [33,72], and 50 m 2 [36].To circumvent this problem, we did not set any gap size limit but let the ecological phenomenon under investigation dictate the appropriate gap size limit.Hobi et al. [73] used the same approach by not setting any gap size limit.The gap mapping used in this paper does not take into account the vegetation height.Therefore, our mapping, although consistent with the definition given by Brokaw [74], focuses only on shaded gaps or dark objects.Zielewska-Büttner et al. [33] reported as well that in their attempt to map forest canopy gaps, the shadow occurrence and forest height affected the mapping accuracy.
Furthermore, the forest canopy mapping was highly affected by the quality of the orthomosaic and image acquisition conditions.In this study, the orthomosaic co-registration presented some flaws in certain areas.Those areas displayed a low quality in gap delineation and even in the visual interpretation of what constitutes a gap.The second limitation came from the difference in reflectance between the datasets acquired in the two consecutive days.Although the images were acquired at noon to reduce the bidirectional reflectance effect [22], the sun illumination between the two days was not the same, leading to two sets of images with slightly different brightness.The best practice would be to collect data covering the site in the same light conditions and, if possible, in the same day.
Finally, the contrast split algorithm failed to detect illuminated canopy gaps where the bare soil is visible, because the bare soil reflectance is as high as the one of the vegetation.Therefore, the algorithm missed big gaps with no vegetation but with apparent bare soil.As an alternative, multiresolution segmentation would allow one to detect both types of canopy gaps, though it requires all image objects to be subsequently classified by the user (more time consuming than contrast split).The second alternative would be to use image-based 3D point clouds.This dataset, possessing a third dimension, is sufficient for detecting local minima in the digital crown model, e.g., by means of hotspot analysis [82].

Modeling Understory Variables through Canopy Gaps Covariates
The very high-resolution mapping of gaps in the canopy cover provided a set of metrics that proved to be useful as covariates in regression estimation of the density (number of plants), the diversity (Shannon index) and the development (mean DBH, mean HTOT) of the understory in Quercus and mixed forest types.The same cannot be said for Fagus forest.One explanation for Fagus' poorer performance lies in the different light conditions under Quercus and Fagus canopy cover.The crown of Turkey oak, when grown in fully stocked high forest, is oval-shaped, with high canopy base height and medium-textured.Due to this structure, oak does not cast a dense shadow, and therefore allows other trees or shrubs to grow in lower layer.By contrast, beech casts a dense shadow under which relatively few plants can grow (e.g., shade-tolerant seedlings of beech).As a result, the presence of relatively small openings in the canopy cover, as those mapped here (75% smaller than 7 m 2 ), might have an influence on the growth and diversification of the understory in Quercus and mixed forest types.The same level of canopy openness is not sufficient to modify light availability under beech canopy cover.This is confirmed by the relatively low-density values of beech understory.
The patch metrics correlated to the understory parameters, in mixed and Quercus forest types, are all shape metrics (except for the mean DBH in mixed forest).The density and development of the understory in mixed forest tend to increase as gap shape fits to an ellipse.A similar trend is observed in Quercus forest, where the mean DBH positively correlates with gap shape, when it fits to a rectangular geometry.On the other hand, the negative correlation between the diversity of the understory in Quercus forest and its covariate (standard deviation of Rectangular fit) does not have a straightforward ecological interpretation.
It seems reasonable that the diversity and development of the understory can be affected by the micro-porosity pattern of canopy cover, which light can shine through.The cumulative effect of differently shaped small gaps affects understory development under forest canopies dominated by Turkey oak and mixtures between Turkey oak and beech.In fact, two different gaps with the same extent can have different shape metrics.Our findings suggest that regular shapes (rectangular or elliptic) are influential for understory density and diameter growth.This is a relatively new finding, as most previous studies focused on relationships between gap size and understory features.Popma and Bongers [83], for instance, demonstrated that the growth and morphology of seedlings in a tropical rain forest depended on the size of gaps; likewise, the dominance of understory vegetation by shade tolerant species or shade intolerant species [84].

Modeling Living Trees Biodiversity through Canopy Gaps Covariates
As for the understory, strong correlations were found between canopy shape gap metrics and density of micro-habitat bearing living trees and vertical species profile in Quercus and mixed forest types.The patch metric correlated with the density of micro-habitat bearing trees in canopies dominated by Turkey oak are negatively related to the presence of gaps fitting into a rectangular shape (Sum_rect fit), whereas vertical diversification in species profile is negatively correlated with the variability of gap density (Sd_Density).This reflects variability from filament shape (low value) to high dense shape (square).In mixed forest stands, the density of micro-habitat bearing trees is negatively correlated to the variability of gap length (CV_Length), while the mean height of the stand is positively correlated with the variability in gap shape asymmetry (Sd_Asym).
There might be different interpretations for these findings.Firstly, one could assume that the presence of canopy openings affects the growth in diameter and height of living trees.Schliemann and Bockheim [72] stipulated that after the opening of the canopy, the increased solar radiation allows plants to grow faster.Accordingly, the observed micro-porosity of canopy cover might favor sunlight penetration into the canopy layer(s), enabling conditions for species distribution along the vertical profile.In addition, the increased exposure to weather conditions and biotic agents might influence micro-habitat differentiation along the trunk and branches (e.g., cavities and holes, injuries and wounds, trunk and crown breakage, epixilic species colonization).This idea was primarily developed by Runkle [85] who contended that, in a tree fall gap, the microclimate changes associated with the gap formation extend to the bases of trees surrounding the canopy gap.Secondly, one might argue that micro-habitat bearing trees, which are frequently old large trees, create a structural diversity in the forest canopy because of their multiple dead tops, bole, and top decays.This specific structure influences the canopy arrangement and, therefore, creates small canopy gaps with specific patch metrics that can be used as a footprint of habitat trees occurrence.

Comparison with Other Studies and Implications
In the Mediterranean region characterized by complex environmental conditions and high seasonal variability of vegetation properties, UAV remote sensing offers advantages for capturing canopy parameters correlated with biodiversity metrics, in terms of costs and precision compared to alternative remote sensing technologies.A potential niche UAV fills-by providing covariates combined with field-based forest inventory data-can be used in regression models for spatially-explicit estimation of forest biodiversity variables over small forest areas.This study showed that canopy openness metrics can be used to predict, by linear regression in certain types of broadleaved deciduous forests, understory features (density, development, and diversity), vertical species profile, and functional forest ecosystem attributes such as density of micro-habitat bearing trees, which are considered difficult to measure directly [86].This latter result is particularly interesting considering the recent attention given to the inventory of tree microhabitats [87].
The results obtained in this study are far from being random or happening by chance.The correlation was controlled using permutation tests [59], and all linear regression models satisfied linear regression assumptions.The confidence of the results was further bolstered by bootstrap resampling technique.To the best of our knowledge, this study is the first investigating the relationship between living tree parameters such as habitat trees and gap patch metrics from UAV RGB images.
Gap size distribution was roughly the same in the three forest types, with over 75% of gaps being smaller than 7 m 2 , when gaps smaller than 2 m 2 are not considered.Strong relationships between this micro-porosity of the canopy layer and ecological phenomena beneath and inside the forest canopy were found only in two out of the three examined forest types.This finding suggests that variability in canopy cover composition and the size of field inventory plots are to be taken into account when dealing with phenomena associated with spatial patterns of sunlight penetration into the canopy.
In this regard, the poor results of our study in beech-dominated forest types likely stemmed partially from the field inventory plot size, which might have been too small to capture significant relationships.Fagus forest is constituted with few big trees per plot.These horizontal and vertical structures make it impossible to observe a great variety of canopy gaps in a plot of only 529 m 2 .In their study, Getzin et al. [37] found relevant relationships in forest areas with beech as canopy-dominant species using plots almost twenty times larger than the plot size adopted in this study.
Furthermore, gap age is another important characteristic that could be taken into account to improve the quality of correlations.A newly opened gap has different biodiversity characteristics [1], and its effect on understory and living trees would be minimal compared to an old gap of the same size and shape.Therefore, gap age can be used as another variable to be tested in future studies in which forest canopy gaps are used as covariates of biodiversity-related variables of interest.

Conclusions
This research highlights that UAV remote sensing can potentially provide covariate surfaces of variables of interest for forest biodiversity monitoring, which are conventionally collected in forest inventory plots.By integrating the UAV images and field-based data, these variables can be mapped over small forest areas with satisfactory levels of accuracy, at a much higher spatial resolution than would be possible by field-based forest inventory solely.From a statistical perspective, it matters little if the covariate that best correlates with the variable of interest is gap shape or size metrics.Nonetheless, knowledge of the nature and signs of these correlations can help inform gap-based silvicultural approaches aimed at the development and temporal continuity of specific biodiversity components.This knowledge is crucial in the new silvicultural paradigm, in which managers adopt operations mimicking nature.The Radius of Largest Enclosed Ellipse feature describes how similar an image object is to an ellipse.The calculation uses an ellipse with the same area as the object and is based on the covariance matrix.This ellipse is scaled down until it is totally enclosed by the image object.The ratio of the radius of this largest enclosed ellipse to the radius of the original ellipse is returned as a feature value

Radius of Smallest Enclosing Ellipse (ε min
The Radius of Smallest Enclosing Ellipse describes how much the shape of an image object is similar to an ellipse.The calculation is based on an ellipse with the same area as the image object and based on the covariance matrix.This ellipse is enlarged until it encloses the image object in total.The ratio of the radius of this smallest enclosing ellipse to the radius of the original ellipse is returned as a feature value Rectangular Fit -[0,1] ; where 1 is a perfect rectangle.
The Rectangular Fit feature describes how well an image object fits into a rectangle of similar size and proportions.While 0 indicates no fit, 1 indicates for a complete fitting image object.The calculation is based on a rectangle with the same area as the image object.The proportions of the rectangle are equal to the proportions of the length to width of the image object.The area of the image object outside the rectangle is compared with the area inside the rectangle Beech is usually accompanied by other deciduous trees like Quercus cerris L. (turkey oak) in the upper layer and only at lower altitudes, and by Carpinus betulus L. (hornbeam), Castanea sativa Mill.(Chestnut), Acer opalus L. (maple), Corylus avellana L. (hazel), and Ilex aquifolium L. (holly) in the dominated layer.

Figure 1 .
Figure 1.Study area map with three main forest types and 50 sample plots.

Figure 1 .
Figure 1.Study area map with three main forest types and 50 sample plots.

Figure 2 .
Figure 2. The overall workflow and methods implemented.

Figure 2 .
Figure 2. The overall workflow and methods implemented.

27 Figure 3 .Figure 4 .
Figure 3. Example of gap delineation using the Contrast Split algorithm (A) and the size distribution of gaps limited to plots with minimum threshold of 1 m 2 (B) and minimum threshold of 2 m 2 (C).

Figure 3 . 27 Figure 3 .Figure 4 .
Figure 3. Example of gap delineation using the Contrast Split algorithm (A) and the size distribution of gaps limited to plots with minimum threshold of 1 m 2 (B) and minimum threshold of 2 m 2 (C).

Figure 4 .
Figure 4. Example of the correlation matrix of the most correlated and the least correlated predictor variables: (A) Pearson's Correlation of patch metrics variables associated with the sum and (B) Pearson's correlation from the patch metrics associated with coefficient of variation.

Figure 5 .
Figure 5. Spatial patterns of understory parameters with R 2 greater than 50% in mixed and Quercus forests.The predictor variables were (A) standard deviation of rectangular fit, (B) average of rectangular fit, (C) average of asymmetry, and (D) average of roundness.

Figure 5 .
Figure 5. Spatial patterns of understory parameters with R 2 greater than 50% in mixed and Quercus forests.The predictor variables were (A) standard deviation of rectangular fit, (B) average of rectangular fit, (C) average of asymmetry, and (D) average of roundness.

Figure 6 .
Figure 6.Spatial representation of living tree parameters with R 2 greater than 50% in mixed and Quercus forests.The predictor variables were (A) sum of rectangular fit, (B) standard deviation of density, (C) coefficient of variation of length, (D) standard deviation of asymmetry, and (E) average of radius of the smallest enclosed ellipse.

Figure 6 .
Figure 6.Spatial representation of living tree parameters with R 2 greater than 50% in mixed and Quercus forests.The predictor variables were (A) sum of rectangular fit, (B) standard deviation of density, (C) coefficient of variation of length, (D) standard deviation of asymmetry, and (E) average of radius of the smallest enclosed ellipse.
∞); 0 = ideal The Roundness feature describes how similar an image object is to an ellipse.It is calculated by the difference of the enclosing ellipse and the enclosed ellipse Shape Index b v /4 √ A [1,∞) ; 1 = ideal The Shape index describes the smoothness of an image object border.The smoother the border of an image object is, the lower its shape index ∞) ; 1 = perfect circle It is the ratio of a gap's perimeter to the perimeter of a circular gap of the same area Patch fractal dimension (PFD) 2 • ln(b v )/ln(A)

Table 1 .
Summary of diversity indices.

Table 2 .
Summary of exploratory statistics on understory data.

Table 3 .
Results of linear regression with understory data (B = slope of linear regression; SE = standard error, N = sample size).Bold values refer to adjusted R 2 higher than 50%.

Table 4 .
Summary of exploratory statistics on living trees.

Table 6 .
Results of bootstrapping for parameters with R 2 greater than 0.50.

Table A6 .
Coefficient of correlations of Pearson and Spearman for some selected explicative and living trees dependent variables in mixed forest.

Table A7 .
Coefficient of correlations of Pearson and Spearman for some selected explicative and living trees dependent variables in Fagus forest.