Incorporating canopy structure from simulated GEDI lidar into bird species distribution models

The Global Ecosystem Dynamics Investigation (GEDI) lidar began data acquisition from the International Space Station in March 2019 and is expected to make over 10 billion measurements of canopy structure and topography over two years. Previously, airborne lidar data with limited spatial coverage have been used to examine relationships between forest canopy structure and faunal diversity, most commonly bird species. GEDI’s latitudinal coverage will permit these types of analyses at larger spatial extents, over the majority of the Earth’s forests, and most importantly in areas where canopy structure is complex and/or poorly understood. In this regional study, we examined the impact that GEDI-derived Canopy Structure variables have on the performance of bird species distribution models (SDMs) in Sonoma County, California. We simulated GEDI waveforms for a two-year period and then interpolated derived Canopy Structure variables to three grid sizes of analysis. In addition to these variables, we also included Phenology, Climate, and other Auxiliary variables to predict the probability of occurrence of 25 common bird species. We used a weighted average ensemble of seven individual machine learning models to make predictions for each species and calculated variable importance. We found that Canopy Structure variables were, on average at our finest resolution of 250 m, the second most important group (32.5%) of predictor variables after Climate variables (35.3%). Canopy Structure variables were most important for predicting probability of occurrence of birds associated with Conifer forest habitat. Regarding spatial analysis scale, we found that finer-scale models more frequently performed better than coarser-scale models, and the importance of Canopy Structure variables was greater at finer spatial resolutions. Overall, GEDI Canopy Structure variables improved SDM performance for at least one spatial resolution for 19 of 25 species and thus show promise for improving models of bird species occurrence and mapping potential habitat.


Introduction
Due to the widespread impacts of land cover change, pollution, wildlife exploitation, and climate change on species and ecosystems across the world, biodiversity assessment and monitoring is imperative. The United Nations' latest global assessment report on biodiversity and ecosystem services estimates that the number of species threatened from extinction now reaches an unprecedented 1 million, rates of extinction are accelerating, and current actions are insufficient to reverse these trends (IPBES 2019). Monitoring biodiversity across large areas through rapid, scalable, and accurate remote sensing methods will greatly improve our understanding of species distributions, loss and resilience, drivers of changes, and impacts at multiple ecosystem levels (Pereira et al 2013, Corbane et al 2015, He et al 2015, Rocchini et al 2015. Proenca and colleagues (2017) reviewed scientific approaches to global biodiversity monitoring, with a focus on GEO BON's Essential Biodiversity Variables, and highlight the use of remote sensing for characterizing ecosystem structure and function. In terrestrial applications, remote sensors provide measurements of vegetation properties, including chemistry, structure and phenology, which are related to species habitat requirements, specifically shelter and food resources. In this regard, many unexplored opportunities with new technologies remain open to investigation. For example, NASA's Global Ecosystem Dynamics Investigation (GEDI; Dubayah et al 2020) is a light detection and ranging (lidar) sensor that was installed on the International Space Station (ISS) in December 2018. This sensor makes fine-scale measurements of 3D vegetation structure, such as overall canopy height and vertical distributions of plant material. The acquisition of these structural measurements over the majority of the land surface is expected to lead to significant advances in biodiversity monitoring and modeling (Bergen et al 2009, Goetz et al 2014, Bakx et al 2019. In this study, our primary goal is to examine the importance and added value of 3D canopy structure data derived from simulated GEDI lidar waveforms in regional-extent bird species distribution models (SDM).
Bird species have been well-monitored relative to other taxa (Troudet et al 2017), especially in North America, and recent work incorporating long-term surveys and radar remote sensing suggests a drastic decline of 2.9 billion birds since 1970 (Rosenberg et al 2019). Similar to other taxa, habitat loss is a major driver of reduction in bird species abundance and range . Suitable bird habitat includes necessities like food, water, shelter, and nesting sites. Habitats are characterized at various scales, from broad associations to more descriptive microhabitats. Habitat associations are typically related to land-cover classes, frequently broad vegetation types (e.g. coniferous forest, grassland), which can be extracted from existing national or global land-cover maps (e.g. MODIS C6; Sulla-Menashe et al 2019). Habitat associations can be further subdivided into microhabitats (e.g. shaded ground beneath mature canopy) which are characterized by fine-scale topography, climate, and vegetation structure. Measurements of vegetation structure are not widely available at fine scales, but previous work suggests that structural complexity is associated with a greater variety of microhabitats which can lead to higher diversity and abundance due to more opportunities for foraging, shelter, and nesting (Cody 1985, Hunter 1999, Whittaker et al 2001, Hill et al 2004.
Previous SDM applications have mainly used remote sensing to derive spatially-explicit environmental predictors, such as topography and bioclimatic conditions (Wilson and Jetz 2016), or land-cover classes (He et al 2015, Bradie and Leung 2017). When considering bird diversity, continentalextent distributions are largely modeled with bioclimatic variables (e.g. WorldClim; Fick and Hijmans 2017), which can predict overall physiological constraints, and correlate with broad vegetation patterns that influence habitat (Stralberg et al 2009, Lawler et al 2009. However, at regional extents and finer spatial resolution (grain), information on vegetation structure is useful for characterizing microhabitats, particularly in areas with heterogeneous vegetation patterns, for example from topographic variation (Stralberg et al 2009) or land-use/landcover change (Sohl 2014). Bird species diversity is known to respond to both vegetation structure and floristic diversity (Macarthur and Macarthur 1961, Rotenberry 1985, Adams and Matthews 2019, both factors that have been missing in models operating beyond the plot scale. Recent developments in airborne remote sensing have provided missing links at regional to local extents, implying possibilities for broader-extent applications when these technologies are elevated to space. Airborne lidar surveys of topography and vegetation structure are becoming more frequent, but are still limited in areal extent due to operational costs. While most airborne lidar sensors typically capture discrete returns from small footprints (<1 m diameter), airborne and spaceborne waveform sensors capture the full vertical profile of vegetation structure, typically using relatively larger footprints (>20 m diameter) collected over a wider area (Lim et al 2003;Anderson et al 2016;Hancock et al 2019). Lidar data are already improving our understanding of how vegetation structure influences bird species distributions at multiple spatial scales (Tattoni et al 2012, Rechsteiner et al 2017, Carrasco et al 2019 and several review papers have examined lidar-based studies that consider avian distributions (Vierling et al 2008, Davies and Asner 2014, Bakx et al 2019. In general, current research tends to use airborne lidar and derived metrics related to canopy height, vertical and horizontal variability, and cover, as well as understory density. Most studies report a positive effect of lidar metrics in explaining bird distributions and/or patterns of richness, however the extent and scale of analysis are important factors. Species sense and respond to a range of temporal and spatial scales (Wiens 1989, Levin 1992, and this has important implications for parameterizing and interpreting SDMs. (Mayor 2009) andMcGarigal et al (2016) review how hundreds of studies address the issue of scale dependence, finding no consensus. At a minimum, they recommend testing SDMs at a variety of spatial scales, ideally optimizing the scale of each variable (multi-scale optimization).
In this study, we focus on implementing simulated GEDI lidar-the first spaceborne lidar mission designed specifically to measure vegetation structure. The instrument is scheduled to collect data for a minimum of two years, from 2019 to 2020. Initial GEDI data were released on January 22, 2020, and so our analysis uses 2 years of GEDI data simulated from airborne lidar data. Previous research using NASA's airborne full waveform Land Vegetation and Ice Sensor (LVIS), a precursor to GEDI, showed that lidar metrics were useful for modeling bird species richness, abundance, and habitat quality at local extents (Goetz et al 2007(Goetz et al , 2010. More recent research has used space-based ICESat-1 waveform lidar at broader scales (Goetz et al 2014), corroborating other studies showing vegetation properties derived from multispectral satellites (e.g. MODIS), such as percent tree cover and life form type, along with bioclimatic variables were generally sufficient predictors of breeding bird species richness.
We created SDMs for 25 species of birds in Sonoma County (figure 1) using canopy structure predictor variables as well as additional climate, phenology, and other auxiliary geospatial predictors. Our three objectives are to: (1) generate bird SDMs and determine the importance of individual predictor variables and variable groups for 25 bird species; (2) determine whether including Canopy Structure variables improve SDM performance, while also exploring the impacts of spatial scale; and, (3) map probability of bird species occurrence and associated uncertainty of these predictions.

Spatial and temporal domain of SDMs
We focus on Sonoma County because there are a relatively large number of geolocated bird species observations and the entire county is covered by airborne lidar data (figure 1) which we used for the GEDI simulation. A previous effort, the Sonoma County Breeding Bird Atlas (U.S. Geological Survey 2019), was based entirely on field observations and mapped occurrence of bird species in Sonoma County at 5-km spatial resolution. We sought to create finer spatial resolution maps showing probability of occurrence by generating three regularly-spaced square grids covering the county at 250-, 500-and 1000 m spatial resolution.
The airborne lidar collection occurred from September to November 2013 and is the primary temporal constraint on our analysis. We chose water years (starting in October and ending in September) 2013 to 2015 to bracket the lidar data collection. Additional remote sensing datasets used for SDMs were temporally-aggregated to the entire three-year period. Due to the limited spatial coverage of species observations from this three-year period, we extended the species observation window to ten years (2006 to 2015).

Bird species data 2.2.1. Selecting species of interest
We gathered field observations from point count surveys collected by Point Blue Conservation Science (figure 1; appendix A) from the California Avian Data Center (http://data.prbo.org/cadc2/) between 2006 and 2015, along with data from the citizen science observation network eBird (Sullivan et al 2009), and the Breeding Bird Survey (Sauer et al 2017). eBird data comprised the majority of records. We used data from all detailed survey events that recorded all species present, and thus from which we could discern detection or non-detection. If presence detections occurred in a given grid cell at a particular resolution then that cell response was classified as presence (value = 1). Grid cell surveys in which the observer did not encounter the species of interest were classified as absence (value = 0). We selected the top 25 species based on number of grid cell presences at 250 m spatial resolution (see table  1). Main species habitat associations (Conifer, Oak, Shrub, Riparian, Grass, Urban, and Variable; table 1) were determined from the Sonoma County Breeding Bird Atlas, Audubon Guide to North American Birds (https://www.audubon.org/bird-guide) and The Cornell Lab of Ornithology Birds of North America guide (https://birdsna.org; Rodewald 2015). A variable habitat association means that the species does not spend a majority of time in a given habitat.

GEDI lidar and canopy structure predictor variables 2.3.1. GEDI Instrument Description
The GEDI lidar sensor measures vertical canopy structure of forests between ±51.6 degrees latitude. The instrument emits pulses of energy (1064-nm) which reflect off of various canopy layers and/or the Earth's surface, and the returned energy of each pulse is recorded as a short duration time series. Using the speed of light, ISS orbital geometry, and precise positioning, this time series is geolocated on the Earth's surface and transformed into a waveform-a profile of returned energy discretized into 15-cm vertical height bins for each ∼22-to 25 m diameter footprint (Dubayah et al 2020). Footprints are separated by approximately 60 m along each of the 8 ground tracks and by 600 m between tracks. GEDI, like other lidar systems, is not able to make measurements through clouds.

Canopy Structure predictor variables from simulated GEDI
The primary GEDI data product (L1B) is a geolocated lidar waveform. Variations of the main waveform and waveform-derived metrics (L2; Tang et al 2012) are output for each footprint.
We simulated two years of GEDI lidar observations from airborne lidar data by using the methods   (TMn), average maximum temperature (TMx), potential evapotranspiration (PET), actual evapotranspiration (AET), and climatic water deficit (CWD). The monthly data for water years 2013, 2014 and 2015 were resampled from a native spatial resolution of 270 m to our three analytical spatial scales using an averaging method. Variables for each water year were then grouped and averaged into four three-month periods: October-December (Q1), January-March (Q2), April-June (Q3), and July-September (Q4). Final SDM climate variables were averages of the three water years (2013-2015) for each three-month BCM variable, resulting in 24 total 'Climate' variables for each spatial analysis scale.

Auxiliary variables
We used four 'Auxiliary' predictor variables as inputs to SDMs: elevation, distance to ocean coast, distance to streets, and distance to streams. Sonoma County has a diverse topography, ranging from Pacific coastlines to mountains greater than 1000 m in elevation. Elevation can be a proximate determinant of a species' presence (Hof et al 2012). Elevation rasters were derived from the USGS National Elevation Dataset (NED; 1/3 arcsecond native spatial resolution). In Sonoma County, distance to coast is correlated with a longitudinal pattern of habitats. Closer to the coast are mixed hardwood and conifer forests, such as Coastal redwoods and Douglas fir, with interspersed patches of chaparral, and further inland occur conifer patches, oak woodlands, chaparral and grasslands. Distance to coast was also calculated initially as a 30 m raster from a coastal vector layer using Arc-GIS. Distance to streets was incorporated as a proxy gradient for human influence. We used the Sonoma County streets vector layer to calculate distances in a 30 m raster. Distances to streams were incorporated as an additional predictor related to moisture and water availability. Numerous species of birds prefer habitats along or near streams (Rottenborn 1999, Mcclure et al 2015. We used the lidar-derived streams vector layer from the Sonoma County Vegetation Mapping and Lidar program (https://sonomavegmap.org) to calculate distances in a 30 m raster. All auxiliary variables were resampled to the spatial analysis scales (250-, 500-, and 1000 m) using an averaging method.

Phenology
Multispectral satellites like the Moderate Resolution Imaging Spectroradiometer (MODIS) are particularly useful for monitoring vegetation presence, vitality, and phenology. We calculated Dynamic Habitat Indices (DHI) Phenological predictor variables from MODIS, most recently described by Hobi et al (2017) and Radeloff et al (2019). These indices are measures of vegetation productivity over the course of a year. We were limited to using MOD13 NDVI since it is the only index measured at (approximately) our finest analysis scale. The DHIs we calculated include: (1) DHI sum-the area under the phenological curve of a year, (2) DHI min-the minimum value of the phenological curve of a year, (3) DHI var-the coefficient of variation of the phenological curve of a year, (4) DHI median-the median value of the composite time series, (5) DHI 95p-the 95th percentile value of the composite time series, and (6) DHI seasonal difference-the difference between mean June and mean November NDVI, a potential proxy for deciduousness in this area. The six DHI metrics were calculated at their nominal spatial resolution (232-, 463-, 927 m) and then resampled to our analysis scales. See appendix C.c. for additional details.

Species distribution model approach 2.5.1. Individual SDMs
We combined Climate (24), Canopy Structure (20), Phenology (7), and Auxiliary (4) variables for a total of 55 predictors covering Sonoma County at each spatial resolution. We rescaled each predictor variable by subtracting the mean and dividing by the standard deviation. To account for potential multicollinearity we used a variance inflation factor (VIF) analysis to remove the most correlated variables. We used our highest spatial scale scenario (250 m) variable values and set a VIF threshold value greater than or equal to 10-equivalent to excluding variables that have an R-square of 0.9 or higher when regressed against all other variables. We selected the same variables for the other two model scales that use Canopy Structure predictors. This model scenario is referred to as 'All with Canopy Structure' . For comparison we ran a scenario referred to as 'All without Canopy Structure' which excluded Canopy Structure variables. We ran seven individual SDMs (Random Forests, Support Vector Machine, Boosting, Extreme Gradient Boosting, Neural Network, Net Regularized Generalized Linear Models, and K-Nearest Neighbors) 500 times (bootstraps) for each species, resolution, and model scenario combination. Each bootstrap corresponded to a random sample (80% training, 20% testing) of spatially-thinned species observations (Aiello-Lammens 2015) and associated predictor variables. Cross-validation was used to tune models and prevent overfitting. Additional details related to individual model tuning and bootstrapping can be found in appendix D. We focused on the area under the receiver operator characteristic curve (AUC; Fielding and Bell 1997) for individual model evaluation. AUC values equal to 0.5 mean that the predicted values are equivalent to random guesses. AUC values greater than 0.5 and up to 1 (perfect) indicate that the model has some level of predictive power.

Variable importance
We used a sensitivity analysis method from the R package rminer (Cortez andEmbrechts 2013, Cortez 2016) to assess variable importance. We selected the data-based sensitivity analysis (DSA) method and associated average absolute deviation from the median (AAD) for measuring importance. We focused on the model scenario 'All with Canopy Structure' and assessed variable importance for species both in terms of individual variables and variable groups (Climate, Canopy Structure, Phenology, and Auxiliary). We also assessed variable importance aggregated by habitat association.

Ensemble SDMs
We generated two different types of ensembles, an individual iteration ensemble (IIE) and combined iteration ensemble (CIE). The IIE for each species calculates a weighted average prediction (Marmion et al 2009) value from up to seven individual model predictions for a single bootstrap (figure 2). To create this ensemble we only selected individual models with AUC greater than 0.5 and then calculated a weighted average prediction for each grid cell, where the weights were based on adjusted AUC (AUC score minus 0.5). All prediction averaging was done following a logit transform of the model predictions. Averaged values were then transformed back to probability (0 to 1). Our second ensemble method, CIE, aggregates species model predictions from individual models across all 500 bootstraps, again weighting the final prediction value by adjusted AUC. The maximum number of individual prediction values was 3500 (7 models * 500 iterations), assuming all individual models had AUC greater than 0.5. We used the weighted average prediction value and associated weighted standard error to display probability of occurrence and associated uncertainty for each species (figure 2). All SDM analyses were performed in R (R Core Team 2019). Appendix H lists specific R packages used and associated references.
2.6. Impact of Canopy Structure variables on model performance 2.6.1. Ensemble SDMs Differences in predictive accuracy when contrasting model constructs can be highly informative with respect to the importance of a particular predictor for understanding the niche of a species, and for management and decision-making (Austin and Van Niel 2011Mod et al 2016, Bell and Schlaepfer 2016. However, consistency in identification of important environmental variables varies with the algorithm used and other factors, so caution must be taken to ensure models do not over-fit when tested with and without a variable, as they may include variable interactions that increase model complexity and may not correctly reflect the ecology of the species (Aguirre-Gutiérrez et al 2013, Bell and Schlaepfer 2016). This was noted by (Araujo and Guisan 2006), who recommend the use of carefully chosen predictors (see also Araújo et al 2019). We postulated that GEDI variables provide ecologically-meaningful information about vegetation structural complexity that partly determines the niche of several of our study species. We visually compared the performance of SDMs that included and omitted Canopy Structure variables. We expected the algorithms used would adjust model complexity to optimally predict with and without the Canopy Structure variables. However, we also expected Canopy Structure variable importance would still be evident as an increase in performance when they are included in the model.
Though information from the full range of the species has been noted as very important for SDM performance (Kadmon et al 2003, Syphard and Franklin 2010, Martínez-Freiría et al 2016see review in Engler et al 2017), the geographic coverage of the airborne lidar dataset precluded us from running models outside Sonoma County. Relatively less important than geographic bias and sample size (Thibaud et al 2014), spatial autocorrelation will also affect SDM performance and was not included in our models. Since our goal was to evaluate the relative importance of the Canopy Structure variables in SDMs, we focused on building relatively accurate models and comparing among constructs that vary only with respect to the variables used.

Logistic regression analyses
In order to provide additional evidence of the impact of Canopy Structure variables, we used all the variables remaining after VIF to construct an additive logistic regression model (the 'All with Canopy Structure' model scenario), and also fitted the 'All without Canopy Structure' model scenario. For each bootstrap sample, we then calculated the likelihood ratio test statistic of the full model vs the restricted model (i.e. with vs without Canopy Structure), and report the range of values of the statistic in relation to the value at which it is expected to occur fewer than 5% of the times (the 'statistically significant' departure value). Results of this logistic regression analysis are provided for comparison in appendix E.c. since our main focus is the use of Canopy Structure variables in ensemble SDMs.

Predictor variables and importance
The VIF method used to remove highly-correlated predictor variables at 250 m spatial resolution Figure 2. Overview of SDM framework for the 'All with Canopy Structure' modelling scenario. VIF reduced the total number of predictor variables to 23. Then for each of the 25 species we used seven different machine learning algorithms at three different spatial scales. We ran 500 bootstraps for each species and scale. We calculated the aggregated variable importance for each species and also produced weighted prediction and uncertainty maps. reduced the number of predictor variables from 55 to 23. This reduction left 38% of Climate, 35% of Canopy Structure, 100% of Auxiliary and 43% of Phenology predictor variables. The following Climate variables remained: Actual Evapotranspiration (AET) from Q1, Q2, Q3, and Q4, Potential Evapotranspiration (PET) from Q1 and Q2, Precipitation from Q1 (PptQ1), Maximum Temperature from Q1 (TMxQ1), and Minimum Temperature from Q1 (TMnQ1). The following Canopy Structure variables remained: Biomass (BM; referred to as BM_2_1 in appendix Table C1), Leaf Area Index 0-to 10 m, 10-to 20 m, 20-to 30 m, and 30-to 40 m, Vertical Distribution Ratio Middle (VDRM) and Vertical Distribution Ratio Bottom (VDRB). All Auxiliary variables remained after VIF. The following Phenology variables remained: NDVI Annual 95th Percentile (NDVI_Ann_95p), NDVI Seasonal Difference (NDVI_Seas_Diff), and NDVI Variance (NDVI_Var). The correlation matrix for the remaining variables is shown in appendix figure C3 .
For all 500 bootstraps we summarized variable importance by predictor variable group for habitat associations (figure 3) and individual species (appendix E). We ranked DSA variable importance for each individual bootstrap and then selected the top-5 most important variables from each bootstrap. Figure 3 shows the relative distribution of all top-5 variables by variable group and habitat association. From this perspective at 250 m, Canopy Structure variables most frequently ranked in the top-5 most important variables for Conifer, Shrub, and Urban specialists. The importance of Canopy Structure diminishes at coarser spatial resolutions, such that Climate variables most frequently ranked in the top-5 when considering Shrub and Urban specialists at 1000 m. Across all spatial resolutions and habitat associations, Climate variables most frequently ranked in the top-5, but the relative frequencies generally matched the proportion of Climate variables available relative to the total number of variables. However, Oak, Riparian, and Grass specialists had a higher frequency of Climate variables ranked in the top-5. For all habitat associations, Phenology variables ranked in the top-5 more frequently than would be expected based on proportion of Phenology variables available relative to the total number of variables. Across all habitat associations at 250 m, Auxiliary variables comprised 12.9% of the top-5 (vs. 17.4% expected), Phenology 19.3% (vs. 13.0% expected), Canopy Structure 32.5% (vs. 30.4% expected), and Climate 35.3% (vs. 39.1% expected). These importance patterns are corroborated by another perspective in which we summed total importance for each resolution and predictor variable group (supplementary figure E2 (stacks.iop.org/ERL/15/095002/mmedia)).
The DSA importance method also provided insight into which individual predictor variables were most important (appendix figure E3). Individual Phenology variables were relatively more important than median importance of variables from all models, regardless of habitat association. Canopy Structure variables were usually about as important as the median variable importance. The two vertical distribution ratios (VDRB and VDRM) were more important than other canopy structure variables, with VDRB showing noticeably higher relative importance for all habitat associations. Canopy Structure variables BM, LAI 0-to 10 m, and LAI 10-to 20 m were relatively more important for Conifer specialists, particularly at finer spatial resolution. Climate variables were relatively homogenous in their importance, but TMxQ1 showed noticeable deviation from the median importance value. Precipitation from the first quarter (PptQ1) was another relatively important Climate variable. Importance of Climate variables generally increased going from fine-to coarse-scale.

Model performance 3.2.1. Comparing performance of different model scenarios
Ensemble models were most frequently the best performing models over all bootstraps (appendix F.a.). Figure 4 shows notched boxplots corresponding to Ensemble Weighted Average (EWA) AUC values from 500 bootstraps for each species and resolution. These boxplots are useful for visualizing differences among bootstrap population medians. The 'All with Canopy Structure' scenario showed median AUC improvements (i.e. non-overlapping notches) for at least one Every Conifer specialist showed a median AUC improvement for at least one spatial resolution. Five of the six Oak specialists showed an improvement in performance for at least one spatial resolution. Few of the Grass or Riparian specialists showed improvement after incorporating Canopy Structure variables. All three Shrub specialists showed a median AUC improvement for at least one spatial resolution. Lastly, two of four Urban species showed a median AUC improvement for at least one spatial resolution. Thresholded differences in AUC are shown in Supplementary Table F1. For example, at 250 m, the median AUC value increased by at least 0.02 AUC for 10 of 25 species, while the 95th percentile AUC value increased for 14 of 25 species. The maximum increase in median AUC after incorporating Canopy Structure variables was 0.04 for CBCH.

Influence of spatial resolution
We summarize 'All with Canopy Structure' median model performance by spatial resolution in Supplementary Table F2. Considering all 25 SDMs, most (11) have the highest median AUC value at 250 m. Grass and Riparian habitat specialists are notable exceptions, with all of their highest median AUC values at 500 m and 1000 m, respectively. Conifer, Oak, and Shrub specialists generally have the highest median AUC at 250 m spatial resolution. Supplementary Table F2 also shows the difference between the maximum and minimum AUC medians across all spatial resolutions, providing insight into the scale variability of each SDM. The average scale variability is 0.03 AUC, while the maximum variability is 0.067 AUC. Therefore, the effect of scale on median model performance is similar to the effect of incorporating Canopy Structure variables.

Ensemble maps
We calculated CIE weighted average predictions from individual models which had AUC values greater than 0.5. Figure 5 shows species distribution maps for birds from a range of habitat associations as predicted by the ensemble weighted average as well as the corresponding uncertainty. The uncertainty estimates include variance due to the model used, the bootstrap sampling, and unexplained variance by the model. Because of the small sample sizes of the bird observation data, the total predicted error is relatively large. All species maps are presented in appendix G.

Discussion
Our ability to describe, parameterize, and model species' habitat continues to improve with the addition and refinement of relevant geospatial datasets. Our central objective was to determine the added benefit of Canopy Structure variables from simulated GEDI lidar in SDMs that also include relevant Climate, Phenology, and Auxiliary variables. We quantified the added benefit of these Canopy Structure variables by calculating variable importance of the 'All with Canopy Structure' model scenario and by comparing the performance of this model scenario against the 'All without Canopy Structure' model scenario using 500 bootstraps. The spatial scale at which SDMs are generated is a very important consideration and impacts our estimates of variable importance and comparative assessment of model per-formance. For this analysis we selected three practical SDM spatial scales which coincided with readily available and ecologically-meaningful Climate, Phenology, and Auxiliary predictor variables. Below we discuss the variables that remained after the VIF process, their importance, the change in SDM performance after incorporating Canopy Structure variables, current constraints on the analysis, and ideas for building on this work in the future.

Model variables and importance
Methods to reduce redundant information across predictor variables, such as VIF, are useful for creating models that are parsimonious and interpretable in terms of variable importance. Furthermore, the elimination of certain variables can provide insight regarding unique information content relative to the entire stack of environmental predictors. Considering all variable groups, we found Climate variables to be the most important overall. Maximum temperature was the most important Climate variable for Grass and Oak species, which agrees with previous studies (Howard et al 2015, Bradie andLeung 2017). Six of nine Climate variables were associated with either potential or actual evapotranspiration (ET), which is related to a combination of temperature, moisture availability, and plant productivity. Hawkins et al (2003) and Coops et al (2018) found ET variables to be important for predicting species richness at broad spatial scales. Barbet-Massin and Jetz (2014) showed that annual PET was among the most relevant variables for predicting bird species distributions in the conterminous U.S. The ET variables do not appear to be relatively as important when grouped by habitat association, but still contain some unique information. For example, PETQ2 is one of the most important variables for Bushtit and Mourning Dove.
Seven Canopy Structure variables remained following VIF reduction, and as a group were usually the second most important variable group depending on habitat association. Canopy structure is an important determinant of avian diversity, including the partitioning of niches for species coexistence, changes in microclimate (Zellweger et al 2019) and refuge from predation (Davies and Asner 2014). Müller et al (2010) found that 3D structure measured with airborne lidar was the main statistical determinant of bird assemblages in a German temperate mixed forest. Macarthur and Macarthur (1961) first reported the link between bird species diversity and foliage height diversity (FHD). However, neither version of the FHD metrics remained following VIF reduction. This is not necessarily because the variable is without utility, but rather some of its information content varied linearly with other variables that also include information about the vertical distribution of plants, like vertical distribution ratios (VDRB and VDRM), LAI profiles, and NDVI. VDRB was a particularly important variable for all habitat associations. This variable was initially included to characterize how evenly plant material is distributed, as it compares the overall canopy height to the height of 50% of the cumulative waveform signal. In theory, values close to 0.5 indicate an even distribution, while lower values should indicate a higher concentration of plant material in the understory, and vice versa. In actuality, when examining the kriged raster we observe that the variable shows spatial correspondence with tall vegetation or built-up structures, both of which could be useful for perching or nesting. Some species SDMs may find this variable useful for discriminating between canopy structure characteristics. For instance, this was one of the most important variables for California Towhee, which commonly uses dense shrub/scrub areas that are low to the ground.
Canopy Structure variables were most important for Conifer specialists (figure 3), suggesting that these species seem to prefer locations with certain Canopy Structure characteristics (i.e. based on importance metrics, more biomass at the canopy level and denser foliage, or LAI, in the 0-to 10 m and 10-to 20 m strata). This agrees with our observations of these conifer-associated species in Sonoma County. For example, Dark-eyed Junco is often found in areas with complex vegetation structure less than 5 m from the ground, but also forages in midcanopy. During the breeding season, males often sing from exposed perches near the tops of conifers and snags. Chestnut-backed Chickadee prefers structurally-complex mature conifer forests, as corroborated by the importance of the variable LAI 20to 30 m, and is often found gleaning small insects and other arthropods from bark and twigs.
Phenology variables mainly capture temporal cycling of top-of-canopy foliage status, but in some cases (low to moderate forest cover) may contain information about foliage density associated with the near-infrared component of NDVI which is related to the phenomenon of leaf additive reflectance (Jensen 1996). The Phenology variables had the highest top-5 importance difference relative to their expected importance (+6.3% overall at 250 m) suggesting that these variables contain particularly useful information for bird SDMs. Visually, the three Phenology variables show some correspondence with broad land cover and vegetation classes (figure 2), which may provide sufficient predictive information for some SDMs. Canopy Structure variables from GEDI also had a positive top-5 importance difference (+2.1% overall at 250 m) for all habitat associations. Canopy Structure variables likely complement Climate and Phenology variables by adding 3D information related to microhabitat structure that is not captured by the other variables we used. This 3D perspective has the greatest benefit for Conifer specialists (+9.5% overall at 250 m).

SDM performance 4.2.1. Ensemble model performance perspective
Fewer than half of the ensemble SDMs had median AUC greater than 0.7, a threshold for assessing SDM quality (Pearce and Ferrier 2000, Duan et al 2014. The overall performance of these 25 SDMs was likely related to the bird observation data and the large distribution range of the individual species (Mcpherson and Jetz 2007). Although Sonoma County is relatively diverse in terms of climate, topography, and vegetation type, more bird observations spanning a wider range of the environmental predictor variables could lead to better presence/absence separation and thus improved SDM performance. Furthermore, the best performing SDMs were generally associated with Conifer and Oak species, which are relatively more selective and less likely to occur across multiple habitats in Sonoma County. Fairly common species SDMs, such as Mourning Doves, Song Sparrows, and House Finches did not benefit much from the incorporation of Canopy Structure information due to their high prevalence in multiple habitats (Kadmon et al 2003, Brotons et al 2004 see Gavish et al 2017 for a counter-example).
Larger observation datasets and more certainty in presence/absence records both help models to more accurately depict the niche of common species (Kadmon et al 2003). While there are a large number of survey events in Sonoma County (n = 10606), these occurred largely in the same locations-when aggregated within cells, only a few hundred cells were surveyed (see table 1). However, because most cells had many survey events, our determination of presence/absence within each cell is probably very accurate, the more so for these common and easy to identify generalist species. It is therefore likely that the relatively poor performance of most models is due to low sample sizes and associated lack of high-density observations from different land-cover classes. These are common limitations associated with citizen science observations, such as those used in this study (Boakes et al 2010, Beck et al 2014. Although Canopy Structure variables frequently improved median model performance, especially at finer spatial resolution, the improvement in performance was relatively modest considering the number of new variables (and potential amount of new information) added to the models. The various machine learning models used in the ensemble performed well at discerning patterns and optimizing predictions with relatively few predictor variables. Another reason for modest improvements in performance was related to correlation and information-overlap between variables remaining following VIF. Even after using VIF, there was still correlation among the predictor variables (see appendix figure C3) suggesting that some canopy structure information could be explained by other variables. For examples we found Pearson correlation coefficients greater than 0.6 when regressing NDVI95p against three Canopy Structure variables (Biomass, LAI 0-to 10 m, and LAI 10-to 20 m).
Lastly, spatial scale dependence cannot be ignored in SDM analyses (Wiens 1989, Jackson andFahrig 2015). In this study, we found an average scale variability of 0.03 AUC (Supplementary Table F2) which is similar to the effect of incorporating Canopy Structure variables in some SDMs. When going from fine to coarse spatial scale, we found Canopy Structure variable importance decreased and incorporation of Canopy Structure variables less frequently improved model performance. Mertes and Jetz (2018) found resampling environmental predictor variables with fine spatial structure decreased the performance of SDMs. Curry et al (2018) found that spatially-explicit ensemble models with finerscale predictor variables consistently performed better at predicting occurrence for 3 of 11 grassland bird species. Our results, and those from others, are related to the spatial structure and autocorrelation of Canopy Structure variables (Mertes and Jetz 2018). We did not explicitly factor in spatial autocorrelation of environmental predictor variables into SDMs, but semivariograms generated for each Canopy Structure variable confirmed the spatial structure (semivariogram range) of each variable was generally within 250to 1000 m (Supplementary Table C2). The majority of semivariograms generated from GEDI simulation subsets within Conifer and Oak forests have ranges below 500 m. Additional SDM performance improvements may be possible through exploration of multi-scale optimization (e.g.

Constraints on the current analysis
The primary limitations of this study are the quantity and geographic distribution of bird observations, as well as our ability to represent canopy structure at optimal scales using simulated spaceborne lidar data. The biggest limitation associated with the bird species observation data is a lack of high-density, evenlydistributed spatial observations. As shown in figure 1, most observations occurred in urban/suburban areas. Since we are primarily interested in how forest Canopy Structure variables improve SDMs, this was not the ideal network of observations because relatively few observations were distant from roads and inside of intact forests. The impact of these Canopy Structure variables is likely more difficult to resolve when urban features are mixed with tree canopies at the grid cell resolutions used in this analysis. We would expect to receive a more consistent signal from a 'natural' forest canopy compared with a suburban or urbanwildland interface canopy interspersed with built-up structures.
Other limitations are associated with our ability to represent canopy structure using simulated spaceborne lidar data. First, the simulated lidar data are hypothetical geospatial observations-actual mission footprint locations will vary as a function of orbital geometry, which cannot be predicted exactly, and local cloud cover patterns. Another major limitation concerns interpolation from footprints (L2 products) to grids (L3 products). There are numerous options for interpolating point data to continuous grids depending on the nature of the dataset. As a starting point, we used a kriging methodology similar to the one outlined in the GEDI L3 Algorithm Theoretical Basis Document (Dubayah et al 2020). Over the course of two years, GEDI is expected to provide sufficient coverage for generating L3 products at 1000 m spatial resolution. However, coverage gaps will still be present in 1000 m spatial resolution grids and would be even more prevalent at the finer spatial resolutions used in this study (see appendix figure C2). For natural and continuous forest types, we expect this kriging method to be a good option for interpolation as neighboring grid cells tend to have similar geostatistical characteristics. For fragmented or urban/suburban forests, this method may not be as effective. More advanced interpolation methods such as co-kriging (Tsui et al 2013) or fusion  with auxiliary data with continuous spatial coverage, like Sentinel-1, Sentinel-2 or Landsat satellite imagery, may prove to be better options for interpolating large areas or filling in sparsely-sampled grid cells.

Future directions
We expect canopy structure measurements from GEDI will be informative for modelling habitat for a variety of forest-dependent taxa. Our modeling approach is applicable to larger spatial extents and taxa, but one important consideration would be the source of Climate data, as BCM is a Californiawide climate product with an unusually fine spatial resolution (270 m). Global climate variables from WorldClim (Fick and Hijmans 2017), which have approximately 1000 m spatial resolution, are more commonly used for larger SDM extents outside of the conterminous U.S. We note that this resolution matches the planned GEDI L3 products. It will be necessary to test the utility of these global climate variables and GEDI L3 metrics in SDMs at their native resolution. However, for some species the tradeoff in model performance as a function of spatial scale may preclude the use of relatively coarse vegetation structure and/or climate data (Manzoor et al 2018). Alternatively, shifting to more of a microhabitat focus, recent studies have proposed frameworks for dynamically modeling microclimate using fine-scale topography and vegetation structure datasets derived from remote sensing (Kearney et al 2019, Lembrechts et al 2019. Using these finer spatial and temporal scale climate datasets in conjunction with vegetation structure measurements will improve our understanding of the mechanistic interplay between the two, as well as how species and ecosystems respond to changes in climate (Zellweger et al 2019).
It will also be beneficial to explore data and methods for increasing the quantity and coverage of canopy structure measurements. GEDI is scheduled for two years of data collection on the ISS but a different type of spaceborne lidar data, being collected simultaneously from ICESat-2, may have some utility in filling GEDI canopy structure coverage gaps. The Advanced Topographic Laser Altimeter System (ATLAS) instrument on-board ICESat-2 uses photon counting technology which has limitations in dense vegetation, especially during daytime observations (Popescu et al 2018). Although vegetation measurement is not the primary focus of the ICESat-2 mission, we still expect some useful canopy structure observations in forests of sparse to moderate canopy cover (Neuenschwander and Pitts 2019).

Conclusions
We incorporated Canopy Structure variables derived from simulated GEDI lidar into 25 bird SDMs in Sonoma County, CA. We found Canopy Structure variables were generally the second most important group of variables across a variety of habitats. Canopy Structure was most important for Conifer and Shrub habitat associations. The incorporation of Canopy Structure variables into ensemble SDMs resulted in performance improvements for the majority of species across all spatial analysis resolutions. More SDMs showed higher, albeit modest, performance improvements at the finest spatial analysis scale suggesting that Canopy Structure variables are more beneficial at smaller scales which are more in line with microhabitats. Although this study covers a regional spatial extent and limited number of species, it demonstrates that Canopy Structure variables derived from spaceborne lidar have the potential to improve our ability to map species distributions.
Our understanding of global forest canopy structure will grow by orders of magnitude over the next decade as a result of the GEDI and ICESat-2 spaceborne lidar missions. We encourage future SDM efforts to explore the utility of operational lidar products from these sensors across natural and anthropogenic landscapes. Future SDM studies may also benefit from fusion of spaceborne lidar with other remote sensing data as this combination can provide more continuous spatial coverage across a broader range of spatial scales. SDMs generated at even finer spatial scales may benefit more from lidarderived Canopy Structure variables, resulting in more detailed map products for land management and conservation of species habitat.
Sun conducted the link margin analysis we used to set the noise levels in the GEDI simulation. Scott Luthcke produced orbital simulations and gave expert suggestions regarding expected data quality/quantity.
Bird observations were largely collected and organized by the invaluable contribution of a number of citizen scientists, through eBird (Sullivan et al 2009) and the BBS programs.
We used Northern Arizona University's high performance computing system (Monsoon) to run all species distribution models. Maxwell Gilbert, Emma Forester, and Mariana Palacio Salazar assisted with pre-processing of predictor variable datasets. Sean Mahoney provided contextual knowledge about bird species of interest. Wendy Schackwitz provided a helpful review.
We acknowledge funding support from NASA's Citizen Science for Earth Systems Program (Cooperative Agreement 80NSSC18M0107 to MC) and the GEDI Science Definition Team (NNL15AA03C to SJG).
Lastly, the manuscript benefited substantially from two anonymous reviews.

Code and Data Availability
The data that support the findings of this study are available upon request from the authors.