DOI: 10.5061/dryad.1c59zw3s7
title: Habitat heterogeneity captured by 30-m resolution satellite image texture predicts bird richness across the U.S.

description: Species loss is occurring globally at unprecedented rates, and effective
 conservation planning requires an understanding of landscape
 characteristics that determine biodiversity patterns. Habitat
 heterogeneity is an important determinant of species diversity, but is
 difficult to measure across large areas using field-based methods that are
 costly and logistically challenging. Satellite image texture analysis
 offers a cost-effective alternative for quantifying habitat heterogeneity
 across broad spatial scales. We tested the ability of texture measures
 derived from 30-m resolution Enhanced Vegetation Index (EVI) data to
 capture habitat heterogeneity and predict bird species richness across the
 conterminous U.S. We used Landsat 8 satellite imagery from 2013-2017 to
 derive a suite of texture measures characterizing vegetation heterogeneity
 (available at http://silvis.forest.wisc.edu/webmap/landsat8-evi-textures).
 Individual texture measures explained up to 21% of the variance in bird
 richness patterns in North American Breeding Bird Survey (BBS) data during
 the same time period. Texture measures were positively related to total
 breeding bird richness, but this relationship varied among forest,
 grassland, and shrubland habitat specialists. Multiple texture measures
 combined with mean EVI explained up to 41% of the variance in total bird
 richness, and models including EVI-based texture measures explained up to
 10% more variance than those that included only EVI. Models that also
 incorporated topographic and land cover metrics further improved
 predictive performance, explaining up to 51% of the variance in total bird
 richness. A texture measure contributed predictive power and characterized
 landscape features that EVI and forest cover alone could not, even though
 the latter two were overall more important variables. Our results
 highlight the potential of texture measures for mapping habitat
 heterogeneity and species richness patterns across broad spatial extents,
 especially when used in conjunction with vegetation indices or land cover
 data. By generating 30-m resolution texture maps and modeling bird
 richness at a near-continental scale, we expand on previous applications
 of image texture measures for modeling biodiversity that were either
 limited in spatial extent or based on coarse resolution imagery.
 Incorporating texture measures into broad-scale biodiversity models may
 advance our understanding of mechanisms underlying species richness
 patterns and improve predictions of species responses to rapid global
 change.
description: Methods
 <b><i>Image
 texture measures</i></b> We calculated a suite of texture
 measures based on Landsat 8 EVI composite images available in Google Earth
 Engine (GEE; http://earthengine.google.org). The 30-m resolution,
 eight-day EVI composite product (GEE Image ID LANDSAT/LC8_L1T_8DAY_EVI)
 was generated from Level L1T orthorectified scenes, using
 top-of-atmosphere (TOA) reflectance that accounts for solar angle and
 seasonal variation in Earth-Sun distance (Chander et al. 2009). The EVI
 was calculated based on three bands (near-infrared, red, and blue) of each
 image (Huete et al. 2002). We selected EVI over NDVI because it is less
 sensitive to soil background and atmospheric effects and less prone to
 saturation at high levels of biomass (Huete et al. 2002). To create a
 smooth composite image for texture analysis, we extracted
 90<sup>th</sup> percentile EVI values from available images
 between May–September during 2013–2017, thereby characterizing peak
 greenness of vegetation during the summer growing season while excluding
 spuriously high EVI values (Culbert et al. 2009, Tuanmu and Jetz 2015). We
 masked pixels covered by permanent water bodies using a static water mask
 also derived from Landsat imagery (Hansen et al. 2013). In image texture analysis, central
 pixels within a moving window are assigned a value based on the spectral
 variability of neighboring pixels (Hall-Beyer 2017). First-order textures
 are statistical summaries (e.g., mean, variance) of pixel spectral values
 within the moving window, while second-order textures are based on the
 gray-level co-occurrence matrix (GLCM) and thus take into account the
 spatial arrangement and relationships among neighboring pixels (Haralick
 et al. 1973). Using GEE, we calculated one first-order and six
 second-order texture measures from the 90<sup>th</sup>
 percentile EVI composite (see Table 1 for list of texture measures and
 descriptions), which we selected based on their performance in predicting
 local and regional bird richness patterns (St-Louis et al. 2009, Culbert
 et al. 2012). We used a moving window size of 5x5 pixels (2.25 ha), an
 area large enough to encompass one or more typical breeding bird
 territories (Leonard et al. 2008, Jones 2011) while capturing relatively
 fine-resolution landscape features. We selected a single moving window
 size for our analysis because texture measures across varying window sizes
 tend to be highly correlated and have similar relationships with bird
 richness (St-Louis et al. 2006, Culbert et al. 2012). <b><i>Ancillary
 environmental variables</i></b> To test how well texture measures
 predict bird richness compared to more commonly used abiotic and biotic
 predictors of biodiversity, we calculated two metrics based on topography
 and three metrics based on land cover. We analyzed elevation data from the
 National Elevation Dataset (NED), which provides seamless elevation
 coverage for the conterminous U.S. at 1 arc-second (approximately 30 m)
 resolution (GEE Image ID USGS/NED). We also derived the terrain ruggedness
 index (TRI) from the NED, which quantifies topographic heterogeneity by
 taking the square root of the sum of squared differences between an
 elevation pixel and the eight pixels surrounding it (Riley et al. 1999).
 To characterize land cover composition, we analyzed 30-m resolution 2011
 National Land Cover Data, also derived from Landsat imagery (NLCD; GEE
 Image ID USGS/NLCD). We focused on three dominant land cover classes that
 provide habitat for terrestrial birds: forest (deciduous, evergreen, and
 mixed classes combined), shrubland, and grassland cover. Breeding bird
 data The North American Breeding Bird
 Survey (BBS; <i>available online</i>; Sauer et al. 2017) is a
 long-term, annual survey of ~3,000 routes across the U.S. and Canada.
 Volunteer observers record all birds seen and heard during 50 three-minute
 counts spaced evenly (0.8 km apart) along each 39.4-km route. Surveys are
 conducted during the height of the breeding season between May–July, with
 additional guidelines for time of day and weather conditions intended to
 increase detectability and reduce biases in the data (Robbins et al.
 1986). We removed surveys of non-randomly established routes, surveys
 conducted in inclement weather or outside of established date and time
 ranges, and surveys that did not follow other BBS sampling protocols based
 on the BBS quality codes. We also removed surveys conducted by first-year
 observers on a given route to reduce potential observer effects (Kendall
 et al 1996). We
 restricted our analyses to data collected during 2013–2017, i.e., all BBS
 data collected after the launch of Landsat 8 that were available at the
 time we conducted our analysis. For each BBS survey route within the
 conterminous U.S., we calculated route-level species richness as the
 cumulative number of unique species observed during the 5-year period. We
 excluded rare species with insufficient data (&lt; 30 observations),
 species not adequately sampled by diurnal transect surveys (i.e., raptors
 and crepuscular species), and species associated with habitats
 under-sampled by BBS route locations (i.e., marine, coastal, and
 freshwater species). After pre-processing, 2,749 BBS survey routes
 remained for analysis (Fig. 1), with a total of 535,676 observations of
 332 bird species (Appendix S1: Table S1). We analyzed total bird richness
 and richness of habitat specialists within three broad habitat types for
 birds: forest, grassland, and shrubland (Appendix S1: Table S1). We
 defined habitat specialists as birds known to primarily occur in only one
 habitat type during the breeding season, and determined specialization
 based on habitat designations of major importance from BirdLife
 International (IUCN 2019) and detailed species accounts from Birds of
 North America (Rodewald 2015). Of the 332 species included, we identified
 95 forest specialists (29% of species evaluated), 16 grassland specialists
 (5%), and 38 shrubland specialists (11%). Some species were not affiliated
 with any of these three habitats, or were affiliated with one or more
 habitats but were not identified as specialists. Statistical analysis
 To relate our remotely sensed data
 to bird species richness, we calculated the mean value of each
 environmental predictor within 19.7 km of the centroid of each BBS route
 (sensu Flather and Sauer 1996, Pidgeon et al. 2007). We selected a radius
 half the length of a BBS route to ensure that each sampled landscape would
 contain the entire BBS route, and defined route centroids as the center of
 the minimum bounding rectangle encompassing each route. Because BBS
 richness data conformed closely to a normal distribution (results not
 shown), we used linear regressions to relate environmental predictors to
 bird richness. We included both linear and quadratic terms of predictors
 in regression models to account for potential nonlinear relationships
 between bird richness and environmental variables. For all analyses, we
 fitted models for all species combined and for each species group (forest,
 grassland, and shrubland). To test the performance of texture
 in predicting bird richness, we first evaluated each texture measure
 individually in single-texture models, including both linear and quadratic
 terms. We used coefficient estimates to assess the strength and direction
 of relationships between each texture and bird richness, and adjusted
 <i>R</i><sup>2</sup> values to evaluate the
 explanatory power of models. We also used results of single-texture models
 to inform which textures to include in subsequent multivariate analyses.
 Because different texture measures represent different landscape
 characteristics and might complement each other to predict bird richness,
 we also tested models incorporating multiple textures. However, because
 many texture measures are correlated (Baraldi and Parmiggiani 1995), we
 checked for collinearity by calculating pairwise Spearman’s correlation
 coefficients. We found strong correlations among some textures (Appendix
 S1: Fig. S1), and thus only included uncorrelated textures
 (|<i>r</i>| &lt; 0.7) in multiple-texture
 models. To assess
 whether EVI-based textures of habitat heterogeneity complement EVI as a
 measure of available energy in predicting bird richness, we assessed a
 model including only EVI as a predictor. We then ran a series of multiple
 linear regression models combining EVI with one or more textures as
 predictors. Again, we first checked for collinearity and did not find
 strong correlations between EVI and our suite of textures
 (|<i>r</i>| range = 0.1–0.4; Appendix S1: Fig. S1). We
 evaluated adjusted <i>R</i><sup>2</sup> values to
 compare the predictive performance of models with and without textures.
 Using the best performing model based on the adjusted
 <i>R</i><sup>2</sup>, we produced predictive maps
 of total breeding bird richness and for habitat specialists within each
 habitat group. To
 evaluate the relative importance of texture compared to more commonly used
 measures of environmental heterogeneity in predicting bird richness, we
 fitted a global model including our two topographic metrics (elevation,
 TRI), three land cover metrics (proportion of forest, grassland, and
 shrubland cover), productivity (EVI), and a single texture measure
 (dissimilarity). We opted to use dissimilarity because it was the best
 predictor in single-texture models for total species richness and is an
 intuitive metric of habitat heterogeneity. We centered and standardized
 all predictors to allow for unbiased comparisons of effect sizes. We then
 generated a set of linear regressions including all possible combinations
 of these seven predictor variables and their quadratic terms. Because of
 the large number of variables in the global model, we ranked models using
 the Bayesian Information Criterion (BIC), which penalizes
 over-parameterized models. Prior to fitting models, we assessed
 collinearity among predictors (Appendix S1: Fig. S2). The highest
 correlation was <i>r</i> = -0.72, between EVI and proportion
 of shrubland cover, with all other correlations |<i>r</i>|
 &lt; 0.7. As an additional collinearity check, we calculated variance
 inflation factors (VIFs) for each predictor in top-ranked models and
 removed predictors with VIFs above a cut-off value of 10 (O’Brien 2007).
 We used adjusted <i>R</i><sup>2</sup> values to
 assess the total explanatory power of top-ranked models, and evaluated the
 contribution of each variable in predicting bird richness by plotting
 their effect sizes (standardized regression coefficients) with 95%
 confidence intervals. To further evaluate the relative importance of
 predictors, we used hierarchical partitioning to assess the independent
 and joint contributions of each predictor to overall variance explained
 (Chevan and Sutherland 1991). In hierarchical partitioning analysis, joint
 contributions represent the explanatory power of each predictor that
 cannot be disentangled from other predictors due to multicollinearity,
 while independent contributions represent the variance uniquely explained
 by each predictor (Mac Nally 2000). Lastly, we checked for potential
 biases arising from spatial autocorrelation of BBS route locations by
 calculating Moran’s <i>I</i> and analyzing model residuals in
 correlograms using 500 permutations. All statistical analyses were
 performed in R version 3.5.1 (R Core Team 2018; see Appendix S1: Table S2
 for list of R packages used).