Airborne laser altimetry and multispectral imagery for modeling Golden-cheeked Warbler (Setophaga chrysoparia) density

Robust models of wildlife population size, spatial distribution, and habitat relationships are needed to more effectively monitor endangered species and prioritize habitat conservation efforts. Remotely sensed data such as airborne laser altimetry (LiDAR) and digital color infrared (CIR) aerial photography combined with well-designed field studies can help fill these information voids. We used point count-based distance sampling survey data and LiDAR-fused CIR aerial photography to model density of the Golden-cheeked Warbler (Setophaga chrysoparia), an endangered songbird, on the 10 000-ha Balcones Canyonlands National Wildlife Refuge (BCNWR). We developed a novel set of candidate models to explain Golden-cheeked Warbler detection probability and density using habitat covariates characterizing vegetation structure, composition, and complexity as well as habitat fragmentation, topography, and human infrastructure. We had the most model support for covariates calculated using focal means representing a 3.2 ha territory size (100 m radius) vs. 1.8 and 7.0 ha territory sizes. Detection probability decreased with canopy cover and increased with topographic roughness. Golden-cheeked Warbler density increased with canopy cover, was highest at a 7:3 ratio of Ashe juniper (Juniperus ashei) to broadleaf tree canopy cover, and decreased with global solar radiation. Predicted warbler densities using 3 min point counts were similar to six estimates from independently collected warbler territory mapping on BCNWR with a mean difference of 6% and a Root Mean Squared Error of 1.88 males/40 ha. The total population size for BCNWR was estimated at 884 Golden-cheeked Warbler males (95% CI 662, 1206) and predicted densities across the refuge ranged from 0.0 to 0.50 male warblers per ha. On the basis of observed habitat relationships, we defined high quality habitat as having at least 60% canopy cover with Ashe juniper comprising 50-90% of the canopy. We estimated 48% of the area at BCNWR managed for Golden-cheeked Warblers was in high quality habitat conditions and identified patches within the lower habitat quality areas (14% of warbler management areas) that had the greatest potential to become high quality habitat with management. Our approach combined robust wildlife surveys with highly scalable remotely sensed data to examine habitat relationships, estimate population size, and identify existing areas of high quality habitat. This method can be applied to other species of conservation interest and can be used with multiple years of remotely sensed data to assess changes in habitat at local to regional scales.


IntroductIon
A frequent goal and sometimes legal requirement of wildlife managers is monitoring population size and responses to management for priority species (MacKenzie et al. 2006). Yet, estimates of population size and species-specific habitat relationships are typically unavailable at the spatial grain and extent to sufficiently inform land management decisions (Scott et al. 2002). This is particularly true for the National Wildlife Refuge System (NWRS), a network of >600 000 km 2 of lands specifically managed for conservation of plant and animal species (National Wildlife Refuge System Improvement Act of 1997). Recent US Fish and Wildlife Service (USFWS) initiatives further outline science-based strategies for conservation, emphasizing the role of modeling priority species habitat relationships to formulate land management objectives (Johnson et al. 2009). This presents a significant challenge for refuge and other land managers who must prioritize conservation activities with limited knowledge of existing landscape conditions and areas best suited for sustaining or recovering animal populations (Gregory et al. 2013).
Conversely, the potential for developing detailed information from which to model spatially explicit habitat relationships at a spatial scale and extent desired by land managers has increased substantially in recent years (Craighead and Convis 2013). Remotely sensed data are increasingly available at high spatial resolutions which can be readily combined with wildlife field surveys to examine habitat relationships and generate spatial predictions of density and occurrence (Kerr andOstrovsky 2003, Seavy et al. 2009). A significant advancement in remote sensing technology has been the development of airborne laser altimetry, also known as discrete return light detection and ranging (LiDAR). As an airborne system, LiDAR has achieved operational use for a variety of natural resource applications, much the same as commercial aerial photography during the post-World War II period (Jensen 2007, Hudak et al. 2009). Over one-half of the 47 National Wildlife Refuges in the USFWS Southwest Region (AZ, NM, OK, TX) currently have Li-DAR data available and all are covered by NAIP high spatial resolution (1-m pixels) multispectral aerial photography flown at 3 to 4 yr intervals.
The LiDAR provides an intuitive threedimensional data source (i.e., multiple x, y, and z coordinates in space or a "point cloud") that can be readily converted to vegetation canopy height, cover, density, and other metrics describing habitat conditions and complexity at a desired species-specific spatial scale (Kane et al. 2010. Less common are applications aimed at fusing LiDAR data with high spatial resolution aerial photography for characterizing habitat composition. Inaccurate colocation of aerial image pixels with LiDAR point clouds can hinder cross-sensor applications to estimate parameters such as tree species composition. However, modern digital aerial photography acquired by high geometric quality scanning technology greatly enhances habitat information derived from both sensor types (Hartfield et al. 2011, Chen et al. 2012. Within the United States, the US Department of Agriculture (USDA) Farm Service Agency National Agriculture Imagery Program (NAIP) collects high-resolution (1-m pixels) statewide multispectral and color infrared (CIR) digital aerial photography at 2-and 3-yr intervals that can be paired with LiDAR data to develop wildlife habitat metrics (Chen et al. 2012). A number of recent studies have related woodland composition and structure to Golden-cheeked Warbler (Setophaga chrysoparia) density , Peak and Thompson 2013, Reidy et al. 2015 and occupancy (Collier et al. 2010, Warren et al. 2013a). Thus far, only two studies have used LiDAR data to describe habitat characteristics important to warblers and both found v www.esajournals.org

SESNIE ET AL.
strong support for LiDAR-derived canopy height and cover metrics as predictors of warbler occupancy ) and density (Reidy et al. 2015).
The LiDAR and high spatial resolution multispectral imagery also provide synoptic habitat and natural resource information at a landscape scale that is helpful for prioritizing land management activities (Hudak et al. 2009). Reducing hazardous fuels and susceptibility to droughtinduced tree mortality are increasingly important for maintaining areas of late-successional juniper and oak woodland habitat (Murray et al. 2013, Andruk et al. 2014. Mechanical tree thinning and prescribed burning of understory vegetation are the primary silvicultural practices currently used; however, relatively few treatments have been applied in warbler habitat because of uncertainty regarding post-treatment effects on habitat quality. LiDAR and high spatial resolution multispectral imagery are useful tools for prioritizing land management activities at a landscape scale because fuel and habitat parameters such as canopy cover, base height, and bulk density often collected from extensive field inventories can be accurately estimated from LiDAR (Andersen et al. 2005).
We evaluated NAIP CIR image and LiDARderived habitat metrics as predictors of density for the endangered Golden-cheeked Warbler (hereafter "warbler") on Balcones Canyonlands National Wildlife Refuge (BCNWR) in central Texas. We estimated warbler density using a novel set of covariates including remote sensingbased metrics of vegetation structure, composition, and complexity as well as other spatial data in a geographic information system (GIS) describing fragmentation, topography, and human infrastructure to help determine potential human development and land management impacts on warbler habitat. We developed model covariates based on previously published studies in addition to those related to management activities potentially impacting habitat conditions important to warblers. Our objectives were to: (1) evaluate models comprised of habitat metrics created by integrating NAIP CIR imagery with LiDAR for predicting warbler density, (2) validate the most supported model by comparing predicted warbler density to density estimated from field plots with color-banded birds, and (3) demonstrate how our model can be used to prioritize conservation, monitoring, and management activities to locations with high potential for increasing the amount of warbler habitat. In addition, we assessed the potential for using high spatial resolution and validated remotely sensed habitat data for predicting warbler population size throughout its breeding range in central Texas.

Study species
The warbler is an insectivorous songbird that nests in mature wooded habitat co-dominated by Ashe juniper (Juniperus ashei) (Ladd andGass 1999, DeBoer andDiamond 2006). The warbler's narrow breeding range, confined to the Cross Timbers and Edwards Plateau ecoregions of central Texas (Ladd andGass 1999, Groce et al. 2010), combined with extensive habitat loss and fragmentation has contributed to concern over population status and trends. Central Texas has experienced one of the highest human population growth rates in the United States (4.5% yr −1 ; United States Census Bureau, http://www.census.gov/), resulting in a loss of ~29% of total potential habitat between 2000 and 2010 (Duarte et al. 2013). Other potentially significant threats to warbler populations are increased nest predation by edge-adapted species, which may be more abundant in highly fragmented woodland landscapes (Reidy et al. 2009, Sperry et al. 2009).

Study area
We conducted our study on 47 400 ha in the Edwards Plateau of central Texas, encompassing the 10 000-ha Balcones Canyonlands National Wildlife Refuge (BCNWR) and possible refuge acquisition lands (Fig. 1). This area straddles the division between the Lampasas Cut Plain and Balcones Canyonlands sub-regions (Riskind and Diamond 1988 included a mixture of grasslands, shrublands, and woodlands, with woodlands consisting primarily of mixed age stands of Ashe juniper, Spanish oak (Quercus buckleyi), and plateau live oak (Quercus fusiformis). All woodlands were cleared or thinned prior to refuge establishment in 1992, with the oldest stands in the study area >50 yr old at the time of the study. Since BCNWR establishment, many non-wooded areas have been managed by juniper clearing and prescribed burning to encourage a diversity of plant successional stages and habitat conditions.

Avian surveys
We conducted point-count surveys across the study area using a simple random sample of 250 locations separated by ≥250 m. Survey points intersected all dominant vegetation types on BCNWR, including areas with a low probability of being used by warblers, to characterize warbler habitat relationships across the range of biotic and abiotic conditions in this region.
We rotated six observers among points with each point surveyed for 10 min approximately once every 2 weeks for a total of four surveys per point between 2 April 2012 and 21 May 2012 (with the exception of two points that were only surveyed three times). We conducted surveys from 15 min before sunrise to 5 h after sunrise when air temperature was ≥10 °C, wind was ≤19 km per hr, and there was no rain. We only recorded detections of singing warblers and six other bird species to focus effort on a small number of target species. For warblers, we only recorded males that were sing-ing either of the two most common songs used by this species, described as A and B songs (Spector 1992, Bolsinger 2000, Leonard et al. 2010). These songs are typically loud (mean of 55 dB at 6 m [Warren et al. 2013b], although they may also be "whispered" [Bolsinger 2000]), with the peak intensity at approximately 5000 Hz (Bolsinger 2000). Hearing of all but one surveyor was measured using an audiogram to confirm no significant hearing loss at 4000 and 6000 Hz. All surveyors were trained in distance sampling techniques for 1 week prior to surveying, including several days practicing distance estimation to singing warblers in the field. To assist in determining distances to singing birds, the presumed location of each bird detected was marked in the field on a printed aerial image (1:2000 scale) with point distance bands printed at 25, 50, and 100 m. We also used a laser rangefinder to measure distances to presumed bird locations. This was particularly useful in open and canyon areas where birds could be heard from trees unobstructed from sight. When feasible, we measured distances to singing birds using a global positioning system (GPS) by walking to the bird's location if it was still singing at the end of the 10-min survey. We recorded to the nearest second the time each warbler was first detected to allow the data to be parsed to any survey length and facilitate comparisons with other surveys.

Vegetation surveys
We used georeferenced tree measurements taken in the field to train and validate models used for developing remotely sensed habitat layers such as canopy cover and tree height. We measured vertical canopy cover (sensu Fiala et al. 2006) on 100 × 100-m plots (n = 122) by recording the species of tree canopy intercepted at 10-m intervals (i.e., 100 intercepts/plot) from 1 m above the ground with a densitometer (Geographic Resource Solutions, Arcata, California, USA). We used tree species data recorded on plots to develop model training data sets to classify Ashe juniper and broadleaf canopy from remotely sensed data described in the sections below. We measured mean tree heights on 11.3-m fixed circular radius plots (n = 192) using visual height estimates for trees ≤5 m, accurate to within 0.5 m, to reduce measurement time in excessively dense woodlands and a clinometer for trees >5 m to derive an overall estimate of mean tree height. Summarized field data provided a measure of total, broadleaf, and juniper tree cover and mean tree height at the plot scale. We compared field measurements to remote sensing-based estimates using linear regression and Pearson correlation coefficients for validation of canopy cover and vegetation height data layers. We also used independent image validation pixels to assess classification accuracy for broadleaf and Ashe juniper canopy cover.

Remotely sensed data
We developed spatially explicit habitat metrics to model warbler density from airborne LiDAR altimetry and NAIP CIR aerial photography. We obtained LiDAR data from the Texas Natural Resource Information System (TNRIS) and the Texas Capital Area Council of Governments (CAPCOG). These LiDAR data were collected for the study area during leaf-on periods in 2006, 2007 (Travis and Williamson Counties), and 2011 (Burnet County). LiDAR data collection followed Federal Emergency Management Agency (FEMA) specifications for flood hazard mapping to meet minimum standards for mapping 61-cm (2-foot) contour intervals (FEMA 2010). LiDAR contractors used an Optech ALTM 2050 system for the 2006 and 2007 acquisitions to obtain point cloud data at a 0.7-m (riparian areas) to 1.4-m point spacing and 32° maximum scan angle and an Optech ALTM Gemini system for the 2011 acquisition, that had a lower maximum scan angle (16º) and average point spacing of ≤0.5 m. LiDAR pre-processing steps are described in Appendix A and Fig. A1. We also obtained NAIP CIR orthorectified aerial photography (June/July 2012) from TNRIS. Four band (blue, green, red, and near infrared) NAIP 2012 photographs of high geometric and radiometric quality were acquired using a Leica ADS80 linear array sensor at a 1-m pixel resolution.

Remote sensing-based habitat variables
We hypothesized that variables from remotely sensed data (detailed below) were related to warbler habitat structure and composition based on previous studies (Diamond 2007, Warren et al. 2013a) with the addition of woodland structure and complexity variables potentially altered by management activities (Andruk et al. 2014). We also used LiDAR to map urban development that has been linked to warbler declines (Reidy et al. 2008). We derived all habitat covariate layers (except canopy cover, see below) at a 10-m pixel size. We then calculated focal means of surrounding pixel values at three scales (75, 100, and 150-m radius) corresponding to 1.8, 3.2, and 7 ha, respectively, which spans the range of territory sizes reported for the warbler (Davis et al. 2010, J. Reidy, University of Missouri, unpublished data).
We developed a digital terrain model (DTM) derived from LiDAR (explained in Appendix A) to create topographic variables for global solar radiation (W/m 2 ), slope, and elevation roughness using Spatial Analyst Tools in ArcGIS v. 10.1 (ESRI 2012). Global solar radiation is a measure of the annual amount of radiant energy incident on a site and is important to a number of biophysical processes such as evaporation, transpiration, and photosynthesis affecting woodland composition and structure (Breshears et al. 1997). We calculated elevation roughness as the standard deviation of elevation; higher values indicated steep and irregular terrain.
We used rumple, mean and maximum height, and standard deviation of height as a way to characterize late-successional patterns and woodland habitat heterogeneity potentially associated with warbler densities (Goetz et al. 2007, Vogeler et al. 2013. Rumple is a three dimensional measure of tree canopy heterogeneity and a metric that has shown to be correlated v www.esajournals.org SESNIE ET AL. with vegetation structure and successional stage (Kane et al. 2010). We derived these metrics using FUSION software v. 3.3 (McGaughey 2013) from LiDAR point-cloud data classified as vegetation (Table 1).
To evaluate the potential impacts of woodland management activities such as tree thinning and fuels mitigation, we calculated several variables to examine effects of amount and distribution of over-and understory vegetation on warbler density. We took relative vegetation density measurements at four height strata (1-3 m, 3-5 m, 5-10 m, and 10-15 m) to characterize vegetation understory and overstory conditions. To calculate relative density at the 10-m grid cell scale, we divided the number of LiDAR point returns within each height strata by the total number of returns and multiplied by 100. We also characterized relative density of the overstory as point returns >3 m divided by the total number of returns multiplied by 100. Understory density (point returns ≤3 m) could not reliably be compared between the 2006 and 2007 and the 2011 LiDAR data given differences in the mean number of returns per square meter (<2.0) and scan angles (16º and 32º) between acquisition dates. Instead, we calculated understory distribution at the warbler territory scale by classifying the presence or absence of vegetation with ≤3 m height for each 10-m grid cell; a cell could only be classified as having understory when overstory trees (vegetation ≥3 m) were present.
We used 1-m NAIP 2012 CIR aerial photography to distinguish between broadleaf and juniper canopy cover in combination with LiDAR vegetation height to improve tree canopy classification. NAIP image processing, classification, and validation steps are detailed in Appendix A. We used pixels classified as broadleaf or juniper trees to estimate percent total, broadleaf, and juniper tree cover using the focal sum/total pixels within each focal radius and multiplied by 100. We then resampled canopy cover estimates to a 10-m pixel size to match all other covariate layers. We calculated distance from woodland edge and edge density developed from our remotely sensed canopy cover data to estimate potential habitat fragmentation effects on warbler density. To remove individual or diffuse trees primarily distributed across pasture lands from edge calculations, we removed less continuous areas smaller than 2 ha with less than 20% tree canopy from the edge calculations. This was below the minimum threshold patch size (15 ha) required by warblers for attaining reproductive success (Butcher et al. 2010), but smaller areas may still be utilized by warblers depending on landscape structure and composition (Collier et al. 2012). We converted the remaining patches to a polyline layer to calculate line density and Euclidean distance using Spatial Analyst tools in ArcGIS v. 10.1. We adjusted edge density values to account for the amount of potential edge. For example, areas that were primarily grassland or woodland had much lower amounts of edge habitat than areas with an even mixture of the two habitats. We calculated an adjusted edge density as: ( 1) where, minimum cover was the percentage of woodland or open habitat, whichever was smaller (e.g., minimum cover would be 20 for canopy cover = 20% or canopy cover = 80%). The adjusted edge density represented the level of habitat fragmentation for a given amount of canopy cover (See Appendix B for additional details).
We used building footprints extracted from the LiDAR point cloud data to estimate housing density for areas inside the study area. A point location for every building polygon centroid was used to calculate building density at each of the three spatial scales used above.

Warbler density models
We evaluated candidate models to predict warbler density using variables scaled to three warbler territory sizes (i.e., 1.8, 3.2, and 7 ha) and then compared the models using Akaike's Information Criterion (AIC; Akaike 1974) to determine the scale with the most model support. In addition, a previous study on warblers demonstrated that distance-based models can overestimate the number of warblers present and overestimation increases with point-count duration (Peak 2011). Peak (2011) found that counts limited to 2-min produced the most accurate density estimates. Therefore, we determined the optimal point-count duration for our study by comparing predicted density estimates to independently collected territory mapping data that is described below. We determined the survey duration (between 2 and 6 min) that minimized the overall bias without overestimating density using the Root Mean Squared Error (RMSE).
We developed models for each survey duration, sequentially, by first determining the most supported detection function, and then the most supported covariate model for detection, and lastly the most supported covariate model for bird density. We used AIC to determine model support and only the most supported models from the first and second steps were used in subsequent models. We initially ran our analyses using the "gdistsamp" function in the "unmarked" package v. 0.10-4 (Fiske and Chandler 2011) in program R v. 3.1.0 (R Core Team 2014). The "gdistsamp" function allowed us to model visit-level covariates (e.g., time of day, observer) in the detection function. However, these visitlevel covariates received considerably lower support than site-level habitat covariates such as topography and canopy cover and models with both site-level and visit-level covariates did not converge. Therefore, we simplified analyses by running the models using the "distsamp" function, which only allows for site-level covariates. For "distsamp", we pooled observations from all visits to a point and included the log of the number of visits as an offset term in the density component of the model. Prior to model building, we assessed covariates for normality and, when necessary, transformed covariates using the power transformation recommended by Adjusted edge density = edge density minimum cover + 1 v www.esajournals.org SESNIE ET AL.
the "powerTransform" function from the "car" package v. 2.0-21 (Fox and Weisberg 2011). We scaled all covariates to facilitate model convergence (μ = 0, σ = 1; see Appendix C: Table C1 for covariate summaries). To select the most supported detection function, we first ran the global model (density = vegetation height (ht) + juniper canopy cover (jcc) + broadleaf canopy cover (bcc) + broadleaf canopy cover squared (bcc 2 ) and detection = total canopy cover (cc)) to compare the hazard, uniform, and half-normal functions. We next ran competing models of habitat covariates hypothesized to affect detection. We predicted that vegetation density (e.g., canopy cover, vegetation height) might reduce the distance at which an observer could hear a bird whereas an observer may be able to detect birds from farther away on hillsides vs. plateaus. Therefore, we ran four models containing different metrics related to vegetation density (canopy cover, vegetation height, overstory density, and understory dis-tribution) and two models containing terrain features (topographic roughness and slope). We then combined the most supported vegetation and topographic covariates into one model and, of these seven models, we used the detection model with the most support in all subsequent models.
We developed 23 a priori candidate models representing variation in vegetation structure, tree composition, topography, and woodland fragmentation to evaluate support for habitat factors influencing warbler density (Table 2). We included building density in the candidate set to determine potential impacts of human development on warbler density, although human infrastructure was concentrated in a few minor subdivisions to the east of BCNWR. Explanatory variables were checked for correlations and only variables that were not strongly correlated (<0.7) were included in the same model. To prevent uninformative parameters from being included in Mixed composition and structure model (Structure and composition heterogeneity) cv_horz + jcc + bcc + bcc 2 Succession heterogeneity and composition model ht + jcc + bcc + bcc 2 Succession and composition undst Importance of understory model jcc + bcc + bcc 2 Composition only model cc + jcc_comp + jcc_comp 2 + cc*jcc_comp + cc*jcc_comp 2 Succession and composition edist + edist 2 + bdist Fragmentation and anthropogenic change edist + edist 2 + bdens Fragmentation and anthropogenic change cc + jcc_comp + edens + srad 2 + cc*jcc_comp + cc*edens Structure, composition, fragmentation, physiography (microclimate) cc + jcc_comp + jcc_comp 2+ cc*jcc_comp + cc*jcc_comp 2 + ruf Succession, composition, and terrain. cc + ruf Succession and terrain. cc + bcc + bcc 2 + ruf Succession, composition, and terrain.
the most supported models when using AIC for model selection (Arnold 2010), we sequentially removed the term with the least support for each model (as identified by smallest z-value; absolute value of beta/SE) until removing terms no longer reduced the AIC score (Pagano and Arnold 2009). Although Burnham and Anderson (2002) have criticized this approach because of the potential for model selection bias, we chose this method because our primary objective was to identify the most parsimonious model (Devries et al. 2008, Pagano andArnold 2009) and because of the lack of a perfect solution to addressing the issue of uninformative parameters when using AIC for model selection (Arnold 2010). Therefore, the models in the final set contained simplifications that were a variation in those initially developed in Table 2. The most supported model for each analysis was examined for fit using three bootstrapped goodness-of-fit metrics: sum-of-squares error (SSE), chi-square, and Freeman-Tukey. We developed five spatial predictions of male warbler density using the first 2-6 min of pointcount survey data split into 1 min intervals at the spatial scale with the most model support. We then compared warbler abundance from these five spatial predictions with known abundance from independently collected territory mapping data (see "Territory mapping" below; Robbins 1970, Dobkin andRich 1998) to determine which minute interval had the smallest error.

Territory mapping
As part of a separate project, three surveyors mapped territories and monitored productivity of male warbler populations on six, 40-ha plots from March through June, 2012-2014 (Fig. 1). We captured male warblers using playback of prerecorded warbler vocalizations and banded individuals with a unique color combination. Surveyors visited each plot two to three times per week to locate and record locations of banded and un-banded male and female warblers, locate and monitor nests, re-sight color-banded birds, and determine the number of young fledged per territory. Surveyors updated field maps each week with general locations of territories and traversed each plot at least once a week to identify previously missed males. We used cumulative observations of banded and un-banded males prior to the date each was observed tending fledglings to delineate territories and create minimum convex polygons in ArcMap. For this analysis, we used only data gathered in 2012 and excluded males detected outside the point-count survey dates (e.g., a territory from a male detected during all of March, but not again in April or May would not be included). Surveyors used a combination of band re-sights, contemporaneous singing, plumage or song characteristics, nest locations, breeding activity, and fledgling age to differentiate unique territories. We excluded observations of unidentified males. Each surveyor interpreted their own observations and all data were reviewed and finalized by the lead investigator. We calculated the number of territories in each plot as all of the territories fully in the plot plus half of the partial territories.
We determined broad-scale management implications using density model outputs and habitat layers to identify (1) the extent of existing high quality warbler habitat as defined by modeled habitat relationships and density and, (2) woodland areas with high potential for becoming high quality habitat with the appropriate management action.

Warbler surveys
Overall, we detected 487 singing male warblers during the 998 10-min point counts with warblers detected at least once at 104 of the 250 points. Detections of unique warblers accumulated throughout the 10-min survey with 47%, 58%, 66%, 72%, and 78% occurring during the first 2, 3, 4, 5, and 6 min, respectively. Subsequent analyses were based only on detections that occurred within the first 2-6 min, depending on the analysis (Fig. 2a). The median distance of detected warblers from the observer was 85 m (Fig. 2b). Detections were right truncated to eliminate outliers, with 10% of the farthest observations being discarded following Buckland et al. (2001). Truncation distance differed slightly by point-count duration and ranged from 171 to 180 m.

Remotely sensed habitat variables
We found a high positive correlation between field measurements for LiDAR-derived vegetation height (r = 0.84) and NAIP and LiDAR-derived v www.esajournals.org SESNIE ET AL. canopy cover (r = 0.96). Field and LiDAR heights likely showed lower correlation because observers estimated heights in the field for trees <5 m tall. Combined NAIP imagery and LiDAR-based estimates of broadleaf and evergreen canopy cover were also highly correlated with field measurements (r = 0.91 and r = 0.90, respectively; also see Appendix A: Figs. A2, A3). In addition, validation pixels independent of classification model training showed 89% accuracy for broadleaf and 95% accuracy for Ashe juniper canopy cover. Broadleaf class accuracy improved from 84% to 89% with the addition of LiDAR canopy height, which aided in discriminating broadleaf canopy from bright green cultivated land (Appendix A: Fig. A2).

Territory mapping
We monitored 63 territories (range: 4.5-11.5 territories on a plot) across 6 plots in 2012, of which 30 males were banded (average of 50% per plot; range: 38-60%). We additionally banded 2 females, of which one was paired with an un-banded male and one with a banded male. Territory density averaged 0.22 males per ha across all plots (range: 0.11-0.29 males per ha). The majority of un-banded males was separated spatially or had banded neighbors.

Habitat modeling results
The models developed at the 3.2 ha (100-m radius) scale had the most support for both warbler detection (w i = 0.83) and density (w i = 1.00). We calculated the mean estimates from the 2 to 6 min intervals as 83 to 110% of the territory mapping results, with the smallest RMSE observed using the first 3 min of the point-count survey (estimated mean number of male warblers for the six plots was 94% of the number of mapped territories; RMSE = 1.88 birds/40 ha). Model estimates for four of the 40-ha plots were within one male of the territory mapping results whereas model results from two plots underestimated by about three males/40 ha (Fig. 3). Therefore, we limit reporting to the 3.2 ha (100-m radius) scale and the 3-min pointcount duration. Results obtained at each spatial scale can be found in Appendix C, Table C2. We found the most support for the half-normal detection function (w i = 1.00). For detection covariates, we found more support for terrain variables than vegetation characteristics, but the combined model, comprised of topographic roughness and canopy cover, had the most support (w i = 0.98; Table 3. Detection probability increased with surface roughness and decreased with canopy cover (Fig. 4a,b).
Our most supported density model included canopy cover, proportion of canopy composed of juniper, an interaction between canopy cover and the proportion of juniper cover, and solar radiation (cc + jcc_comp 2 + cc*jcc_comp 2 + srad 2 , Table 3). This model received overwhelming support over our other models (w i = 1.00) and fit the data well based on the three goodness-of-fit tests conducted (all P > 0.19). Warbler density increased with canopy cover, decreased with global solar radiation, and had a non-linear relationship with percent juniper composition (Fig. 5a-c). Although models containing measures of fragmentation were not in the most supported model set, beta estimates for edge density and distance to edge indicated a preference for habitats with lower edge density that were farther from edges (see Appendix C: Table C2).
On the basis of the mapped model-averaged predictions of warbler density across BCNWR, we estimated there were 884 (95% CI 662, 1,206) male warblers on BCNWR during the breeding season of 2012. The highest warbler densities were predicted in a relatively contiguous block of woodland habitat covering the central to southeastern portion of the study area (Fig. 6). Estimates of standard error around these density predictions generally increased with density, but there was more uncertainty on the steeper slopes than on the plateaus for similar density estimates (Fig. 7a,b).

dIscussIon
We were able to estimate warbler density consistent with field observations using our novel set of variables created by integrating aerial imagery and LiDAR. Our results suggest that both CIR aerial photography and LiDAR have high potential for developing large-scale or breeding range-wide warbler density estimates. While LiDAR data were an important component to our study, LiDAR acquisitions are less frequent than CIR aerial imagery and currently cover only about one-third of the breeding range for warblers in Central Texas. However, widely available NAIP CIR photography was, by itself, effective (84% class accuracy for broadleaf trees in this study) for mapping Fig. 3. Comparison of spatial model and territory mapping derived estimates for the number of male Golden-cheeked Warbler territories. Territory mapping data were collected independently from six plots at Balcones Canyonlands National Wildlife Refuges, Texas. Error bars represent 95% CIs, dashed line is 1:1 correlation, and solid line is actual fit. general tree cover categories. In the absence of LiDAR coverage, CIR imagery may be a suitable data source to determine range-wide warbler density and document temporal habitat changes. Misclassification errors we observed such as commissions errors [see Appendix A] can potentially be removed using national or state-wide land cover data sources to mask other vegetation types from woodlands. Spectral and spatial filtering of CIR imagery may also improve classification accuracy of tree species composition (Sesnie et al. 2010) in the absence of LiDAR.
A primary focus of this study was to develop high-resolution digital data layers to better characterize vegetation structure and composition im-portant for maintaining high quality warbler breeding habitat. From validated habitat data, we were able to estimate warbler density and population  numbers that compared well with intensive territory mapping surveys. Warbler locations recorded during the surveys frequently fell within predicted high density areas (Fig. 7a). These results suggest that using the focal means to scale model covariate layers to warbler territory size were robust against potential bias. We found that LiDAR and CIR aerial photography were complementary data  sources, particularly for accurately distinguishing tree canopy composition and thresholds important to warbler density. We were also able to examine habitat complexity not previously explored for this species and found that total canopy cover and the ratio of Ashe juniper to broadleaf tree cover were more important for predicting warbler density than all other variables considered. These findings were consistent with other habitat studies for this species , Warren et al. 2013a, Reidy et al. 2015.
From our models, we further determined that warbler detection increased in rough terrain (Fig. 4). We termed this the "amphitheater effect", where birds may be heard from farther away when the observer is surveying canyon slopes, much the way an amphitheater aids in sound transmission. Conversely, total canopy cover had a negative relationship with warbler detection likely because dense trees were a barrier to song transmission (Fig. 4). Using simulation studies, Pacifici et al. (2008) also found songbird detection rates decreased with increased vegetation density in deciduous forest. Reidy et al. (2015) did not find this type of effect when restricting their evaluation to the amount of area classified as juniper woodland within a 100-m radius of points rather than using a total tree cover estimate. We may have detected warblers at greater distances than other warbler studies such as Peak (2011), because 58% (n = 146) of our points were located in terrain with <10% canopy cover rather than high canopy cover areas typically surveyed for this species. Terrain conditions in combination with vegetation have only been examined in a few studies to determine their influence on warbler detections (Watson et al. 2008, Warren et al. 2013a). Warren et al. (2013a) found that degree slope was positively associated with warbler detection within a 100-m survey radius for occupancy models on Balcones Canyonlands Preserve (BCP). Watson et al. (2008) did not find that slope affected warbler detections and surmised that the effects were small relative to the effects they found due to study areas and season. However, factors affecting detection during occupancy surveys can do so through multiple drivers. For example, if bird density is higher on steeper slopes, bird detection probability may also increase with slope. We found lower support for detection models including percent slope than for topographic roughness (Table 3). Although slope and topographic roughness were highly correlated (Pearson's r = 0.97), models with topographic roughness had two-to-five times greater support than models with slope.
We found the most support for density models that included an interaction between canopy cover and proportion of canopy cover comprised of Ashe juniper. Our characterization of woodland composition using LiDAR and NAIP CIR derived canopy data indicated that peak warbler densities occurred when canopy cover exceeded 80% and supported a ratio of 7:3 of Ashe juniper to broadleaf cover (Fig. 5c). Similarly, Peak and Thompson (2013) and Reidy et al. (2015) reported higher densities of warblers in mixed juniper-oak woodlands than juniper-dominated woodlands and in areas with high canopy cover. Importantly, our results suggest densities fluctuate within woodlands co-dominated by juniper. In general, diverse habitat conditions associated with mixed woodland tree composition have been shown important for meeting warbler life requirements. For example, the interaction between juniper and oak tree species composition has previously been related to greater arthropod abundance and possibly warbler pairing success (Klassen et al. 2012, Marshall et al. 2013. Other model covariates such as tree height and broadleaf canopy cover did not appear in the most supported models because metrics characterizing vegetation structure were often highly correlated with one another (e.g., canopy cover vs. tree height, r = 0.87). As in other woodland habitats, height and canopy cover were also correlated with vegetation density, because vegetation density often increases with the age and successional status, in the absence of major disturbance factors such as wildfire (Miller and Tausch 2001).
Global solar radiation was also included in the most supported density model. The amount of annual solar radiation within our study area varied greatly depending on aspect and adjacent terrain (Appendix C, Table C1). Solar radiation values have been integrated into breeding rangewide models of warbler occurrence (Diamond 2007), but have not previously included in models for predicting warbler density. Lower solar radiation environments likely have higher site moisture content and lower transpiration rates increasing oak regeneration and development (Russell andFowler 2004, Murray et al. 2013) and arthropod abundance (Marshall et al. 2013). No covariates specifically developed to characterize vertical or horizontal woodland complexity, such as canopy height standard deviation, were included in the most supported warbler density models. LiDAR complexity metrics may not be directly applicable to woodland tree canopy architecture in this study area. Kane et al. (2010) noted that LiDAR complexity metrics developed for Pacific Northwest forests required further testing in other vegetation types. While not within the scope of this study, we found low correlation between LiDAR and complexity indicators measured from vegetation plots such as the standard deviation of tree diameters and LiDAR height rumple (r = 0.05) and standard deviation of vegetation height (r = 0.28). Other LiDAR metrics such as understory distribution and coefficients of variation for relative density within height strata were poorly correlated with mean and standard deviation of tree diameters from plots (data not shown). More generally, for this study area mean vegetation height obtained from LiDAR was more positively correlated with the mean (r = 0.48) and standard deviation (r = 0.40) of tree diameters measured on plots. Mixed juniper-oak woodlands with taller, more mature trees that are preferred by the warbler (Pulich 1976, Reidy et al. 2015 likely have greater structural complexity. Candidate models including fragmentation variables such as distance from woodland edge and edge density were not among the most supported models (Table 3). Warbler density may be more influenced by woodland patch composition and height structure than fragmentation in our study area. Reidy et al. (2015) similarly found that woodland composition and canopy height were better predictors of warbler density; however, densities were negatively affected by increasingly open edge.
Overall, warbler density estimates compared well to territorial mapping from a separate field study on BCNWR (Fig. 3). Density estimates were 94% of mean territory map density estimates, with high correlation on four plots and underestimation on two plots. Although some studies have found no statistical difference in abundance estimates between point count and territory mapping (Dobkin and Rich 1998), many studies using 5 or 10 min point counts have found that point counts overestimate true density (Buckland 2006, O'Donnell et al. in press). This overestimation has been documented for warblers in other studies (Peak 2011, Warren et al. 2013b, even when using 2 min point counts (Peak 2011). An essential assumption of distance sampling is that animals are detected at their initial locations (Buckland et al. 2001). Movement during the point count, even random movement, violates this assumption and results in overestimation (Buckland 2006). Conversely, incomplete availability during the point count, that is, not all males present in the point count area sing during the point count, results in underestimation of density. By comparing estimates from different point-count durations to territory mapping results, we attempted to find the duration that balances these opposing biases.

Management implications
We determined woodland composition and structural conditions important for maintaining and managing warbler habitat on BCNWR and identified areas of high quality warbler habitat and woodland areas suitable for developing future habitat. On the basis of observed habitat relationships (Fig. 5a-c), we defined high quality habitat as having at least 60% canopy cover with Ashe juniper comprising 50-90% of the canopy. Within the BCNWR warbler management units, 48% of the area met this definition. As broadleaf cover can be promoted through management activities such as thinning or understory burning (Andruk et al. 2014), we defined potential future high quality habitat on BCNWR as woodlands having close to or greater than 60% canopy cover, but that are currently dominated by juniper (≥90%), with at least some broadleaf species present (≥1%). These conditions occur on 14% of BCNWR warbler management units (1060 ha). These juniper-dominated habitats can potentially maintain successful breeding pairs of warblers (Klassen et al. 2012), but low levels of disturbance favoring broadleaf tree species may increase warbler densities and prevent that habitat from transitioning into pure juniper cover (Murray et al. 2013). Future avian surveys, disturbance, and habitat change information from remotely sensed data will likely help to refine density models and management approaches used to enhance warble habitat conditions.
Mitigating fire hazard and risk is also a concern for Ashe juniper woodlands and warbler habitat areas (Andruk et al. 2014). In many cases, v www.esajournals.org SESNIE ET AL. understory thinning and burning in woodlands and forest are promoted to reduce ladder fuels and the potential for a surface fire to transition to a crown fire (Stephens et al. 2012, Huffman et al. 2013. From our candidate models, we did not find the presence and distribution of understory vegetation within a warbler territory to be important to warbler density. However, reproductive success should be evaluated under varying levels of understory vegetation to determine its relationship with habitat quality. Understory thinning and burning in woodlands can also increase the amount of fine fuels on the ground due to creation of canopy gaps and potential tree mortality. Thus, proposed treatments warrant careful scrutiny to ensure that fire hazard is not exacerbated in these vegetation types. Although this study focused primarily on determining warbler habitat and populations size on BCNWR, we found that high quality habitat also exists on private land proximate to BCNWR boundaries (Fig. 6). Landscape-scale models in conjunction with other field survey data and assessments may help target state and federal conservation programs (e.g., mitigation lands and conservation banks) designed to protect warbler habitat (USFWS 2013).

conclusIons
We developed robust survey methods to estimate warbler density on BCNWR while identifying a repeatable and transferable methodology for monitoring songbirds. Our sampling design using randomly distributed and repeated pointcount surveys produced comparable warbler population estimates when matched to those from intensive territory mapping methods. We successfully combined LiDAR with CIR imagery to develop predictors of warbler density that can be used for managing habitat for this species. Low-cost and high resolution remotely sensed data are increasingly available at the county, state, and national level through a variety of programs (e.g., USGS 3D Elevation Program (3DEP) and USDA Farm Service Agency National Agri culture Imagery Program). Airborne LiDAR will be acquired for all of the conterminous United States, Hawaii, and US Territories within the next 8 yrs through 3DEP. The population monitoring and remote sensing applications developed in this study can help with USWFS-wide efforts to implement strategic habitat conservation planning and meet species recovery goals.