Assessing Greater Sage-Grouse Selection of Brood-Rearing Habitat Using Remotely-Sensed Imagery: Can Readily Available High-Resolution Imagery Be Used to Identify Brood-Rearing Habitat Across a Broad Landscape?

Greater sage-grouse populations have decreased steadily since European settlement in western North America. Reduced availability of brood-rearing habitat has been identified as a limiting factor for many populations. We used radio-telemetry to acquire locations of sage-grouse broods from 1998 to 2012 in Strawberry Valley, Utah. Using these locations and remotely-sensed NAIP (National Agricultural Imagery Program) imagery, we 1) determined which characteristics of brood-rearing habitat could be used in widely available, high resolution imagery 2) assessed the spatial extent at which sage-grouse selected brood-rearing habitat, and 3) created a predictive habitat model to identify areas of preferred brood-rearing habitat. We used AIC model selection to evaluate support for a list of variables derived from remotely-sensed imagery. We examined the relationship of these explanatory variables at three spatial extents (45, 200, and 795 meter radii). Our top model included 10 variables (percent shrub, percent grass, percent tree, percent paved road, percent riparian, meters of sage/tree edge, meters of riparian/tree edge, distance to tree, distance to transmission lines, and distance to permanent structures). Variables from each spatial extent were represented in our top model with the majority being associated with the larger (795 meter) spatial extent. When applied to our study area, our top model predicted 75% of naïve brood locations suggesting reasonable success using this method and widely available NAIP imagery. We encourage application of our methodology to other sage-grouse populations and species of conservation concern.


Introduction
Greater sage-grouse (Centrocercus urophasianus; hereafter sage-grouse) were determined to be not warranted for protection under the 1973 Endangered Species Act [1]. Populations have decreased steadily since European settlement in western North America [2], and the overall range of sage-grouse has been reduced to 56% of its pre-settlement distribution [3]. The major reasons for the decline include degradation, fragmentation, and loss of sagebrush (Artemisia spp.) habitats [4][5][6][7]. Sage-grouse are sagebrush obligates and are highly susceptible to changes in sagebrush habitats. Loss or alteration of sagebrush communities has occurred from invasion by native and exotic plants, increased fire frequency and intensity, overgrazing by livestock, energy development, and agricultural or urban development [4,[7][8][9][10][11][12]. As much as 45% of sagebrush communities that originally existed in western North America have been converted to other landcover types [12]. Consequently, identification of preferred habitat characteristics is necessary to inform conservation and management within remaining sagebrush habitat.
This shift in focus to larger spatial extents has been facilitated by the increased availability and functionality of Geographic Information Systems (GIS) analysis and the widespread availability of multispectral satellite and aerial imagery [37][38][39]. Recent studies have utilized satellite imagery acquired by the Landsat Thematic Mapper [18] or Enhanced Thematic Mapper [30] satellites. These sensors acquire data at a minimum spatial resolution of 30 m (with the exception of the 15 m resolution panchromatic band which is of limited use for analyses) and a spectral resolution consisting of 7 unique bands across the electromagnetic spectrum [40]. The large spatial resolution of these sensors allows for analysis of expansive areas; however, the minimum unit size for any analyses conducted is also limited by the 30-m spatial resolution. Although Landsat data is collected at relatively short intervals (i.e., weekly or monthly) and has been somewhat useful in understanding habitat selection patterns [24,[30][31], this coarse (30m resolution) imagery may not provide the resolution necessary to evaluate habitat characteristics at spatial extents relevant to the process of habitat selection [41].
In addition to the widely and freely available Landsat imagery, there is another dataset available through the National Agricultural Imagery Program (NAIP). This free imagery is collected at approximately 3-year intervals by aerial sensors at a spatial resolution of 1 m and a spectral resolution consisting of 4 unique bands. This fine spatial resolution allows researchers to examine habitat relationships undetectable at the coarser Landsat resolution of 30 m. Using NAIP imagery, factors such as edge effects in highly heterogeneous areas, where patch size is often much smaller than 30 m, can be examined while retaining the capability of assessing large landscapes. In the past, ground-derived microhabitat data were used to assess habitat selection, but its utility was limited [42]. Technology now permits the combination of high-resolution aerial imagery and ground verification, high-resolution images combined together [43], or the combination of high-resolution imagery and ground-based imagery [44] to map habitat and assess habitat selection.
Our goal was to evaluate habitat selection of female sage-grouse with broods across a range of spatial extents utilizing 1 m-resolution NAIP imagery. Our specific objectives were to: 1) determine important features of brood-rearing habitat using fine spatial extent NAIP imagery, 2) assess the spatial extent at which sage-grouse selected of brood-rearing habitat in our study area, and 3) create a predictive habitat model that could be applied across our entire study area to identify areas of preferred brood-rearing habitat. We hypothesized that we would be able to identify features selected by brood-rearing sage-grouse using NAIP imagery and that sagegrouse would select habitat characteristics at multiple spatial extents, as demonstrated in other life history stages for this bird [18].

Study Area
Our study area was an 817 km 2 area surrounding Strawberry Reservoir in north-central Utah (Fig 1). We delineated this area by running a fixed-kernel density estimator (using leastsquares cross validation (LSCVh) to select the smoothing parameter (h)) on 3,865 locations (nest, brood, and non-brood) of female sage-grouse collected from 1998 to 2008. We then used Home Range Tools (http://www.blueskytelemetry.com) for ArcGIS version 9.3 1 (ESRI, Inc., Redlands, CA) to create a 95% polygon surrounding these locations. This polygon (Fig 1) contained 824 of 836 (98.5%) brood locations in our dataset. We removed the remaining 12 brood locations from our analysis, considering them to be outliers.
The study area defined by the LSCVh fixed-kernel density estimate was a high mountain valley that transitioned to lower elevations moving eastward. Elevation ranged from 1,946 to 3,150 m. Average annual precipitation varied widely from 43 cm in the lower elevations to 84 cm at the highest elevations (www.ncdc.noaa.gov). Vegetation consisted of shrublands dominated by big sagebrush (A. tridentata). Silver sagebrush (A. cana) occurred in the more mesic areas and black greasewood (Sarcobatus vermiculatus) was found in some areas in the eastern part of the study area at lower elevations. On slopes at higher elevations, tree communities consisted of quaking aspen (Populus tremuloides), Gambel's oak (Quercus gambelii), and various conifers (e.g. Abies spp., Picea spp., and Pseudotsuga spp.). The tree community at lower elevations was dominated by juniper (Juniperus spp.) with scattered pinyon pine (Pinus edulis). Common forbs found in the study area included longspur lupine (Lupinus arbustus), silky lupine (Lupinus sericeus), sticky purple geranium (Geranium viscosissimum), and sulphurflower buckwheat (Eriogonum umbellatum). Common grasses included Kentucky bluegrass (Poa pratensis), and smooth brome (Bromus inermis). In the lower elevations, cheatgrass (Bromus tectorum) occurred in the understory, but this exotic species was largely absent from the study area. Riparian areas were dominated by willow species (Salix spp). No fires occurred within our study area during our study years and grazing by domestic livestock was absent. Similarly, habitat enhancement via mechanical removal of sagebrush or pinyon-juniper did not occur in the study area until 2009.

Data Collection
We captured female sage-grouse annually by netting them with the aid of all-terrain vehicles on and around leks during the months of March through May using a modified spotlighting method [45]. Once captured, we fitted sage-grouse with necklace style radio transmitters (Advanced Telemetry Systems, Inc., Isanti, MN) and tracked them using a 4-element Yagi antenna and either a Telonics TR 2 (Telonics, Inc., Mesa, AZ) or Communication Specialists R-1000 (Communication Specialists, Inc., Orange, CA) digital telemetry receiver. We monitored females approximately twice weekly from April through August. Nest initiation occurred in late April or May for females in our study area. Females that hatched at least one egg and were observed with at least one chick, were considered brooding. We monitored broods at least twice weekly during daytime hours through 56 days post-hatch. We also opportunistically encountered brooding and non-brooding females during the study. If we were unable to visually detect chicks with a female, we left the immediate area and observed the location for 20 minutes or until the female returned. We classified females that flew long distances and did not return to the area, or were located twice consecutively without chicks as non-brooding and did not include them in our brood sample. For more information on trapping and collection of telemetry data see Baxter et al. 2008 [46] or Peck 2011 [47]. Our sample of broods also included those associated with sage-grouse translocated from four neighboring populations. Because these sage-grouse flocked with resident grouse and demonstrated little to no difference in movements, reproduction, or survival, we lumped their locations with those of resident grouse [46,48].

Ethics Statement
Trapping and handling of sage-grouse was permitted and approved by the Utah Division of Wildlife Resources under a Certificate of Registration (#1COLL6817) and by Brigham Young University's Institutional Animal Care and Use Committee (IACUC approval #08-0402).

Imagery Classification
To characterize vegetation in our study area, we performed a supervised classification on 1-m resolution NAIP imagery collected in 2006. This year represented the first that statewide coverage of 1-m resolution imagery was available for Utah. Thereafter, NAIP imagery was collected every 3 years. Although we could have benefited from imagery collected more frequently-particularly during the early years of our study-Strawberry Valley experienced very little of the habitat change that has impacted sagebrush systems in much of western North America. Consequently, we viewed any potential bias associated with collection of sage-grouse locations in years before or after the 2006 image as unlikely to influence our results. We used ENVI EX Feature Extraction 1 (Exelis Visual Information Solutions, Inc. McLean, VA) to classify NAIP imagery. Using this classification, as well as digitization to manually identify road classes in ArcGIS version 10 1 (ESRI, Inc., Redlands, CA), we generated a landcover layer that divided the landscape into the following 10 classes: paved roads, high-use or major dirt roads (graveled/wide enough for two-way traffic), low-use or minor dirt roads (two tracks), bare soil, shrubs, trees, grass, water, riparian areas, and agricultural areas. Our shrub landcover class consisted of almost entirely sagebrush species; however, due to the limited spectral bands available in NAIP imagery we were unable to differentiate between species. In order to ensure the accuracy of our landcover layer and prior to assessment of sage-grouse selection, we performed an on-the-ground accuracy assessment. Using ArcGIS 10, we randomly distributed 502 points across the study area. In the summer of 2011, we visited 202 of these points and recorded which of the 10 landcover classes best described each location. Using this information and our aerial imagery, we visually interpolated the landcover classes for the remaining 300 locations that we were unable to access for a variety of reasons (e.g. private property). We then used these data (S1 File) to calculate accuracy statistics for our landcover classification [49].

Statistical Analysis
Following accuracy assessment, we developed a list of 86 explanatory variables ( Table 1) that may have influenced selection of brood-rearing habitat by sage-grouse in our area based on previous literature [13,18,[21][22][23][24][25][26][27][28][29][30][31] and our own experience in the study area since 1998. We then divided the variables into two groups: those that would be best examined at multiple spatial extents, and those for which a single spatial extent was adequate. Variables that we evaluated at a single spatial extent (n = 44) included distances to various features and variables for which only the values at the actual use or random site were relevant ( Table 1). The remaining 42 variables (Table 2) were dependent on spatial extent. For variables dependent on spatial extent, we calculated values at three spatial extents by generating circular buffers with radii of 45,200, and 795 m surrounding each site. The 200 and 795-m spatial extents were selected based on the lower and upper end of daily brood movements [22]. The 45-m spatial extent was representative of common patch sizes available to broods in our study area. Prior to modeling,  we tested for multicollinearity between explanatory variables and did not combine in a single model any variables with a correlation coefficient > |0.6|.
To determine the variables that best differentiated use from random sites, we used a multistaged information theoretic approach [50] within a mixed-effects logistic regression [51], using a random intercept to account for individual heterogeneity. We scaled all variables to have a mean of zero and a standard deviation of one prior to analysis. We then used ArcGIS 10 to calculate values for all of our explanatory variables at each spatial extent for 675 brood locations collected from radio-marked females between 1998 and 2008 (remaining locations collected between 1998 and 2008, n = 149, were from unmarked females and we withheld them for accuracy assessment along with locations collected between 2009 and 2012). We then generated an equal number of random (i.e. available) locations from within the study area after masking out Strawberry Reservoir. Because our random locations were cast within the boundary of the study area and not associated with individual home ranges, our modeling of resource selection generally corresponded to Johnson's 2 nd order of selection [52] To ensure that 675 random locations adequately characterized our study area, we calculated the true mean (i.e. mean of all pixels/resource units) for continuous variables and compared our sample means with 95% CIs to these values [53]. In every case, the confidence intervals and even standard errors of our sample overlapped the true mean values suggesting that 675 random locations was adequate to characterize our study area.
Next, we developed 35 a priori, univariable and multivariable models (Table 3) and used model selection within each of our three spatial extents based on previous literature [13,18,[21][22][23][24][25][26][27][28][29][30][31] and our own experience (>15 years) to determine which variables best differentiated use from random locations [51]. To evaluate relative model support, we judged models based on minimization of Akaike's Information Criterion (AIC) [54]. We followed this same procedure for the spatial extent invariant variables with a set of 34 a priori models (Table 4). For each of these four groups (45,200, 795-m spatial extents and spatial extent invariant variables), we advanced the top model and any competitive models ( 2 ΔAIC) to a second stage of analysis.
In our second stage of analysis, we combined the models that were advanced from stage 1 into 7 new models (Table 5). We created these models by combining the top models from each spatial extent with the variables in the top model from the spatial extent invariant group. In these 7 models, for spatial extent dependent variables, we used the spatial extent at which the univariable model for that variable had the lowest AIC in the first stage of analysis.
We evaluated beta coefficients based on their standard errors and 85% confidence intervals [52]. To evaluate effect sizes for variables in our top models, we calculated a resource selection function (RSF) by holding all other variables constant at their mean. We used variance inflation factors (VIF) to test for multicollinearity among variables in our final models. We considered VIF > 5 to indicate multicollinearity. To assess predictive ability of our final models, we performed a k-folds cross validation [52] where k = 5. We sorted the data into 5 partitions, with

TreeCover
The proportion of the "tree" landcover class in a circular buffer

SoilCover
The proportion of the "bare soil" landcover class in a circular buffer

ShrubCover
The proportion of the "shrub" landcover class in a circular buffer GrassCover The proportion of the "grass" landcover class in a circular buffer

RiparCover
The proportion of the "riparian" landcover class in a circular buffer

MajDirtRoadCover
The proportion of the "major dirt road" landcover class in a circular buffer

MinDirtRoadCover
The proportion of the "minor dirt road" landcover class in a circular buffer an approximately equal number of locations in each partition. In each iteration of our procedure, four partitions (80% of the data) were used as the training set, while the remaining partition (20% of the data) was used as the test set. We repeated this procedure until all data were used both as the test set and as part of the training set. We regressed the number of locations from the test (used) dataset in each bin against the median RSF value of the random locations in each bin. We report mean coefficient of determination and slope, and we considered the combination of a high coefficient of determination and a positive slope to be indicative of a model that differentiated between use and available locations well [55].

Results
Our NAIP classification showed composition of landcover in our study to be 45.5% shrubs, 28.7% trees, 11.9% grass, 7.6% water, 3.0% bare soil, 1.7% riparian, 0.7% agriculture, 0.7% lowuse dirt roads, 0.1% paved roads, and 0.1% high-use dirt roads. Accuracy assessment of our classification [49] yielded an overall accuracy of 78.5% and a Kappa value of 0.707. We utilized 675 relocations of sage-grouse broods from 120 females over 15 years. The number of locations per female ranged between 1-28, with a mean of 5.6 locations per female. Unique females were relocated from 1-3 years depending on battery life of transmitters and length of time each bird survived.
The top-ranked model at the 45-m spatial extent included the combination of percent shrub, percent grass, percent riparian, percent paved road, meters of shrub/tree edge, and meters of riparian/tree edge (Table 5). There were no other competitive models. Results from the 200-m spatial extent were the same as the 45-m spatial extent (Table 5). At the 795-m spatial extent, we had three competitive models. The top model included the combination of  (Table 5) and only accounted for 14% of model weight. The spatial extent-invariant group also had two competitive models. The top model consisted of distance to trees, distance to transmission lines, and distance to permanent structures. The next best competitive model, with a ΔAIC of 2.2 and 25% of model weight, was the same as the top model with the exclusion of distance to permanent structures (Table 5).
In the second stage of analysis our top model consisted of 5 variables from the 795-m spatial extent (percent shrub, percent grass, percent tree, percent paved road, and meters of sage/tree edge), 1 from the 200-m spatial extent (meters of riparian/tree edge), 1 from the 45-m spatial extent (percent riparian), and 3 from the spatial extent invariant group (distance to tree, distance to transmission lines, and distance to permanent structures) ( Table 6). This model suggested that 7 variables were negatively correlated with selection of brood-rearing habitat based on their coefficients and RSFs: percent grass, percent tree, percent paved road, meters of shrub/ tree edge, meters of riparian/tree edge, and increased distance from transmission lines (Fig 2;  Table 7). The remaining 3 variables were positively correlated with selection: percent shrub, percent riparian, distance to tree, and distance to permanent structure (Fig 2). From the k-folds cross validation, our mean adjusted r-squared was 0.96, mean slope was 17.67, and Pearson's rank correlation coefficient was 0.95. The top model successfully predicted 75% (n = 84) of the 2009 to 2012 brood locations naïve to development of the models (Fig 3). Variance inflation factors of our top models were < 5.

Discussion
Due to the fine spatial extent (1 m) of our input image, we were able to examine influences of habitat edge on habitat selection by female sage-grouse with broods in a way that has not been done before across such a large area. Our top model included two edge-associated variables that sage-grouse appeared to avoid when selecting brood-rearing habitat: shrub/tree edge, and riparian/tree edge. As we could not differentiate between shrub types, the shrub/tree edge may have been avoided because the shrubs nearest the trees may not have been sagebrush but rather a mountain shrub community. These edge-associated variables provide more information on the relationship of selection for areas with high percent shrub and riparian and low percent tree landcover. Sage-grouse not only avoided areas with a high percentage of trees but also areas that consisted of a patchy mosaic of trees and sage-grouse habitat types at the 795 m spatial extent. In addition, sage-grouse broods were found farther from trees, representing an inverse relationship to distance to trees (Fig 2). The variables in our top model were similar to Atamain et al. 2010 [30]. They identified "xeric mixed sagebrush" as a vegetation type that was selected during early brood-rearing and "moist sites with riparian shrubs" and "montane sagebrush" as areas that were selected during late brood-rearing. While we did not make the distinction between early and late brood-rearing habitat due to the mesic nature of our study site, we did identify percent shrub at the 795 m spatial extent and percent riparian at the 45 m spatial extent as factors selected by brood-rearing sage-grouse. We also identified a negative relationship with percent tree at the 795 m spatial extent where Atamain et al. 2010 [30] reported avoidance of pinyon/juniper woodlands. Dzialak et al. 2011 [31] showed a positive relationship between brood-rearing habitat and percent shrub at the 90 m spatial extent. Our results indicated that sage-grouse selected areas with a higher percentage of shrubs, albeit at a much larger spatial extent. Dzialak et al. 2011 [31] and Dinkins et al. 2014 [32] both showed mixed effects of distance to mesic habitat with sagegrouse showing an aversion to mesic areas during early brood-rearing and a selection for mesic areas during mid and late brood-rearing periods. Our top model did not contain distance to mesic areas. Nonetheless, we did show selection for areas with higher proportions of riparian habitat at the 45-m spatial extent.
Anthropogenic structures such as well pads have been identified as influencing habitat selection by brooding sage-grouse, but negatively influencing survival rates [18]. Well pads did not occur in our study area; however, numerous permanent structures (largely cabins) were located in otherwise suitable brood-rearing habitat. Sage-grouse with broods avoided areas close to permanent structures. The difference between our findings and previous research [18] could be due to higher human activity at the permanent structures in our area compared to well pads or some other difference between these structures and how sage-grouse perceived them.
Distance to transmission lines was another anthropogenic structure included in our top model. Sage-grouse with broods in our study area were found closer to transmission lines than random locations. One possible explanation for this is that the right-of-way, cleared for the transmission lines in our study area, may have created desirable microsite conditions for brood rearing sage-grouse. Another possible explanation is that transmission lines in our study area happened to be located in quality brood-rearing habitat and brood rearing sage-grouse did not actively avoid them. Nonetheless, transmission lines are thought to have a negative influence on sage-grouse for a variety of reasons including provision of raptor perches which has the potential to negatively influence survival rates [16]. We did not measure survival in our analyses and it is possible that there could be decreased fitness of broods that selected areas near transmission lines. A similar phenomenon has been demonstrated with other anthropogenic disturbance features. In two separate studies, sage-grouse selected areas closer to anthropogenic disturbance, but exhibited decreased fitness in these areas [18,33]. Sage-grouse with broods selected habitat characteristics at a variety of spatial extents with at least one variable included from each of the three spatial extents we examined. The majority of these variables in the top model were best at the largest spatial extent, which reemphasizes the need to examine sage-grouse habitat selection at larger extents [16,[35][36]. However, inclusion of percent riparian at the 45-m spatial extent in the top model illustrates the importance of examining small spatial extent habitat characteristics as well. While the spectral resolution of NAIP imagery limits the specificity of the classes to broad categories (i.e. all shrubs vs. sagebrush only), the accuracy for these broad classes was sufficient to create a model that successfully predicted 75% of the 2009 to 2012 brood locations. With the success of classified NAIP imagery in this study, we suggest it may be applicable to other sage-grouse populations, species of conservation concern or any other species in habitat conducive to these methods. The accuracy of our NAIP classification was of critical importance to the validity of our model selection results. While overall accuracy of 78.5% and a kappa of 0.707 could be improved, it is within the range of other widely used vegetation classifications in the western United States (e.g., LANDFIRE = 69.3% to 87.1% overall accuracy across all western United States super zones, no Kappa reported; SWReGAP = 0.60 Kappa; SAGEMAP = 0.537 to 0.878 Kappa across all mapping zones).
The benefits or tradeoffs of using our methods with NAIP imagery are dependent upon the question of interest. The rapid and complex changes to modeling techniques as well as the availability, cost, and quality of imagery influence project applicability. Early efforts by Homer et al. in 1993 [56] demonstrated the immediate utility of remotely sensed data in mapping current-day sagebrush habitats as well as future applications in multispatial-extent modeling. Since then many others have assessed and classified habitat on a national [57], state-wide [58], or project-level [59] using medium or lower resolution Landsat imagery. One advantage to our approach was that we were able to determine the importance of edge at the patch scale on brooding female habitat selection. This may not be possible with more course resolution. Newer, more sophisticated methods combine the use of imagery, often at different spatial extents, with ancillary geospatially explicit data sets [43,[60][61], and at times, ground based verification to assess habitat change. A pitfall to our approach (using higher resolution imagery-1m), is that currently, it cannot be used at larger scales [57][58][60][61] due to processing time and extensive field verification. Other studies [62][63] used a similar approach to assess rangeland tree cover characteristics with relatively accurate classification and similar results. Though its efficacy may be limited to project level analyses, the classification of NAIP imagery as a base layer for landscape spatial extent analyses was a cost-effective method for examining habitat selection at varying spatial extents.
Supporting Information S1 File. This table contains the data used to inform the habitat modeling and resource selection function processes. (CSV)