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.

Methods

Image texture measures

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^th 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^th 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).

Ancillary environmental variables

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; available online; 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 (< 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 R² 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 (|r| < 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 (|r| range = 0.1–0.4; Appendix S1: Fig. S1). We evaluated adjusted R² values to compare the predictive performance of models with and without textures. Using the best performing model based on the adjusted R², 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 r = -0.72, between EVI and proportion of shrubland cover, with all other correlations |r| < 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 R² 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 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).