Habitat heterogeneity captured by 30-m resolution satellite image texture predicts bird richness across the U.S.
Farwell, Laura et al. (2020), Habitat heterogeneity captured by 30-m resolution satellite image texture predicts bird richness across the U.S., Dryad, Dataset, https://doi.org/10.5061/dryad.1c59zw3s7
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 90th 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 90th 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.
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 R2 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 R2 values to compare the predictive performance of models with and without textures. Using the best performing model based on the adjusted R2, 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 R2 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).
U.S. Geological Survey, Award: 140G0118C0009