Estimating Ladder Fuels: A New Approach Combining Field Photography with LiDAR

Forests historically associated with frequent fire have changed dramatically due to fire suppression and past harvesting over the last century. The buildup of ladder fuels, which carry fire from the surface of the forest floor to tree crowns, is one of the critical changes, and it has contributed to uncharacteristically large and severe fires. The abundance of ladder fuels makes it difficult to return these forests to their natural fire regime or to meet management objectives. Despite the importance of ladder fuels, methods for quantifying them are limited and imprecise. LiDAR (Light Detection and Ranging), a form of active remote sensing, is able to estimate many aspects of forest structure across a landscape. This study investigates a new method for quantifying ladder fuel in the field (using photographs with a calibration banner) and remotely (using LiDAR data). We apply these new techniques in the Klamath Mountains of Northern California to predict ladder fuel levels across the study area. Our results demonstrate a new utility of LiDAR data to identify fire hazard and areas in need of fuels reduction.


Introduction
Remote sensing has been key to many areas of environmental science, and wildland fire is no exception [1][2][3][4]. Wildland fires consume homes, endanger lives, and consume a large number of taxpayer dollars every year [5,6]. During the past century, many forests in the Western United States have undergone unprecedented change in forest structure and fire behavior due to forest and fire management practices, including resource extraction and fire suppression [7]. One dramatic change in Western forests that historically experienced frequent, low-to moderate-severity fire has been the infilling of ladder fuels that help facilitate the movement of fire from the forest floor to tree canopies [7,8].
Imagery acquired from the Landsat program has informed the study of wildland fire severity [9,10]. Analyses of these Landsat-based fire severity assessments have shown that wildfires over the last 30 years have produced an uncharacteristically large proportions of high severity effects (where over and herbaceous cover, but also small trees and lower tree limbs. We use the term ladder fuel cover to convey our intention to estimate the connectivity from the surface fuel layer to the canopy fuel layer. Although there is overlap in methodology with some previous efforts to estimate understory cover, our use of the term ladder fuel cover draws the distinction between our estimates and the more commonly used approach for estimating understory cover based on areal extent within a given plot.
Aerial LiDAR is an active remote sensing technology that is able to penetrate the upper canopy, allowing accurate estimations of forest attributes [46,47] including canopy structure [48,49], canopy bulk density [50], and leaf area index [51][52][53]. Airborne LiDAR presents a distinct advantage over terrestrial and space borne platforms for characterizing landscapes because it is able to cover a wide area, while still maintaining high pulse coverage. LiDAR has been applied to wildlife habitat [54,55] and wildland fire behavior modeling, providing fine-scale input across the landscape for terrain, canopy cover, tree height, crown bulk density, and fuel model [26,[56][57][58] with varying degrees of accuracy. However, even though we can now model many inputs to fire behavior models, ladder fuels remain underrepresented. Because LiDAR is able to penetrate the canopy, it presents an excellent data source from which to derive additional information about these ladder fuels. LiDAR has been used to calculate metrics similar to ladder fuels, such as understory vegetation cover [59] and shrub biomass [60], although both of these studies necessitated intense sampling efforts, and neither captured the density and vertical continuity of ladder fuels as a whole.
Here we developed a ground-based methodology for predicting ladder fuel cover that provides more precise, quantitative estimates than the categorical system derived by Menning and Stephens, yet does not require the extensive field data collection efforts required by others [13,59,60]. We then investigate whether LiDAR can predict ladder fuel cover. Specifically, our main research questions for this project are: (1) Does the photo-banner methodology adequately sample ladder fuel cover within each plot? (2) How robustly can aerial LiDAR predict ladder fuel cover across a range of conditions?

Study Area
The Klamath River Basin is located in Northern California, and includes diverse land ownership and management. The 2622 ha (6480 ac) study area is located between Orleans, CA and Happy Camp, CA, centered at 41.5 • North, 123.5 • West ( Figure 1). While LiDAR was collected over a 14,323 ha (35,394 ac) area, four focal areas are utilized by the overarching Western Klamath Restoration Partnership project to concentrate sampling effort in areas of maximum project interest for fuels assessment, and limiting the study area to 2622 ha. The four focal areas represent local forest diversity around geographically distant communities of homes. Topography is steep and complex, with elevation in our study area ranging between 160 and 520 m. Canopy cover ranged between 0 and 100%, but was often quite dense, averaging 90% across the study area (derived from the proportion of LiDAR first returns over breast height, or 1.37 m).
The Klamath has a Mediterranean climate: hot, dry summers, and cool, wet winters, although there is substantial variability in climate across the region [8,61]. Normal daily maximum and minimum temperatures range between 11 • C and 35 • C, and between 1 • C and 11 • C, respectively [8]. The majority of precipitation falls during the winter months, with an average of 143 cm every year [8]. The Klamath includes the highest diversity of conifer forests in North America [62]. Our study area contains areas of forest, shrub, and meadow, with the dominant tree species being tanoak (Notholithocarpus densiflorus) and Douglas-fir (Psedotsuga menziesii). It is also common to find pacific madrone (Arbutus menziesii) and canyon live oak (Quercus chrysolepis). The study area also includes rarer abundance of golden chinquapin (Chrysolepis chrysophylla), California black oak (Quercus kelloggi), California bay (Umbellularia californica), pacific dogwood (Cornus nuttallii), big-leaf maple (Acer macrophyllum), red alder (Alnus rubra), ponderosa pine (Pinus ponderosa), and sugar pine (Pinus lambertiana). Many of these tree species, such as California bay and multiple species of oak, are important sources for tribal gatherers, both historically and in the present [63]. Coniferous species, as well as a high diversity of oaks and other hardwoods present in this system, provide tribal community uses and important habitat for many species [8].
Remote Sens. 2016, 8,766 4 of 23 as well as a high diversity of oaks and other hardwoods present in this system, provide tribal community uses and important habitat for many species [8]. Fire has long been an important part of the Klamath, with many ignitions coming from lightning strikes [8] as well as the indigenous people [64]. Before the 1900s, fires were ignited and managed by indigenous tribes, ranchers, and foresters, maintaining a frequent, mixed severity fire regime [8,64,65]. In 1911, fires were discouraged by the Weeks Act, and by the 1920s, fires were actively suppressed in accessible areas, with remote fires remaining unchecked until 1945, with the advent of enhanced firefighting equipment [8]. Since the implementation of fire suppression, fire frequency decreased dramatically, resulting in a sharp decrease in annual area burned [8,61,66].
The change in fire frequency also led to a change in forest structure, with oaks becoming overtopped, and meadows getting encroached upon by conifers. This resulted in a decrease in forest complexity, with denser, more homogenous forests, and smaller forest openings [61,67]. Fire has long been an important part of the Klamath, with many ignitions coming from lightning strikes [8] as well as the indigenous people [64]. Before the 1900s, fires were ignited and managed by indigenous tribes, ranchers, and foresters, maintaining a frequent, mixed severity fire regime [8,64,65]. In 1911, fires were discouraged by the Weeks Act, and by the 1920s, fires were actively suppressed in accessible areas, with remote fires remaining unchecked until 1945, with the advent of enhanced firefighting equipment [8]. Since the implementation of fire suppression, fire frequency decreased dramatically, resulting in a sharp decrease in annual area burned [8,61,66].
The change in fire frequency also led to a change in forest structure, with oaks becoming overtopped, and meadows getting encroached upon by conifers. This resulted in a decrease in forest complexity, with denser, more homogenous forests, and smaller forest openings [61,67].

Field Data
Field plots were established and measured during the summer of 2015. Plot locations were chosen based on a stratified random sample, which sorted stands in the combined four focal areas by: (1) management history-based on the CalVeg database [68]; (2) insolation; and (3) quadratic mean diameter of the stand. This yielded 10 different strata combinations, or stand types, which are shown in Figure 2.

Field Data
Field plots were established and measured during the summer of 2015. Plot locations were chosen based on a stratified random sample, which sorted stands in the combined four focal areas by: (1) management history-based on the CalVeg database [68]; (2) insolation; and (3) quadratic mean diameter of the stand. This yielded 10 different strata combinations, or stand types, which are shown in Figure 2. More common stand types were assigned more plots, but rare types still contained at least three plots. Plot locations were limited to areas with <45% slope, between 25 m and 200 m from the nearest road, at least 25 m from all stand boundaries (where forest type changes), and with at least 100 m between plots. These limits were established to ensure the safety of the field crew and reduce travel time to plots to maximize plots sampled. Sixty plots were established and measured during the 2015 field season.
Plot centers were recorded with a Trimble Geo 7X GPS receiver with a Tempest antenna. One thousand points per plot were collected, followed by differential correction with a PDOP filter of 4. This resulted in plot center coordinates with an average horizontal precision of 87 cm. While the positional accuracies of the differential GPS and of the LiDAR data are high (sub-meter), their combined error creates an inconsistent misalignment between the two datasets. Because the accuracy of our predictive models relies heavily on the ground data lining up with LiDAR, we tested and adjusted plot center coordinates using mapped tall tree crowns. The distance and azimuth from plot center were measured to at least five tall, dominant conifers whose crowns reached above the main canopy and therefore would be apparent in the LiDAR canopy height model (CHM). Many of these identification trees fall outside the plot boundary. The distances to trees were measured using a TruPulse 360 laser rangefinder with an average error of ±30 cm, thus the local accuracy of the stem map's geometry, as well as the accuracy of the GPS and LiDAR measurements, all fall below 1 m, and should coincide reasonably well. The stem maps-plot centers and locations of non-leaning trees-were incorporated into a GIS database based on the plot centers' GPS positions, and then compared to the CHM generated by the LiDAR point cloud. All tree characteristics recorded in the More common stand types were assigned more plots, but rare types still contained at least three plots. Plot locations were limited to areas with <45% slope, between 25 m and 200 m from the nearest road, at least 25 m from all stand boundaries (where forest type changes), and with at least 100 m between plots. These limits were established to ensure the safety of the field crew and reduce travel time to plots to maximize plots sampled. Sixty plots were established and measured during the 2015 field season.
Plot centers were recorded with a Trimble Geo 7X GPS receiver with a Tempest antenna. One thousand points per plot were collected, followed by differential correction with a PDOP filter of 4. This resulted in plot center coordinates with an average horizontal precision of 87 cm. While the positional accuracies of the differential GPS and of the LiDAR data are high (sub-meter), their combined error creates an inconsistent misalignment between the two datasets. Because the accuracy of our predictive models relies heavily on the ground data lining up with LiDAR, we tested and adjusted plot center coordinates using mapped tall tree crowns. The distance and azimuth from plot center were measured to at least five tall, dominant conifers whose crowns reached above the main canopy and therefore would be apparent in the LiDAR canopy height model (CHM). Many of these identification trees fall outside the plot boundary. The distances to trees were measured using a TruPulse 360 laser rangefinder with an average error of ±30 cm, thus the local accuracy of the stem map's geometry, as well as the accuracy of the GPS and LiDAR measurements, all fall below 1 m, and should coincide reasonably well. The stem maps-plot centers and locations of non-leaning trees-were incorporated into a GIS database based on the plot centers' GPS positions, and then compared to the CHM generated by the LiDAR point cloud. All tree characteristics recorded in the field, including a unique ID number, tree height, species, DBH, and lean, were added to the GIS to help with the CHM comparison. The CHM was generated at 25 cm pixel resolution to correct for possible GPS-LiDAR misalignment.
In general, the tree positions in the CHM were close to the tree positions derived from the GPS; however, in a number of instances, the tree configuration apparent in the CHM was clearly different than the ground data. Plot locations were only shifted if an improved location fit was apparent. The combination of individual plot center and trees associated with that plot were moved manually and as a unit. Only x-y shift was applied (no stretching or warping) until the new fit between GPS trees and the CHM was optimal. We cross-checked the alignment of the shifted plot using tree heights as measured on the ground vs. the local maxima of the LiDAR-derived CHM. Table A1 shows the horizontal precision reported for each plot, as well as the distance of plot shift. Of the 60 plots, 45 were shifted an average of 3.3 m. The average shift distance for all plots, including those that were unshifted, was 2.5 m. Figure 3 shows an example case of plot shifting. In this example, some trees and part of the unshifted plot fall on the road, whereas trees in the shifted plot line up with the LiDAR-derived canopy height model. Each plot center was intentionally located at least 25 m from roads, showing that the plot center as measured by the GPS does not line up with the LiDAR data. Had the unshifted plot been used, model accuracy would have been reduced due to the inclusion of road in the LiDAR point cloud. field, including a unique ID number, tree height, species, DBH, and lean, were added to the GIS to help with the CHM comparison. The CHM was generated at 25 cm pixel resolution to correct for possible GPS-LiDAR misalignment. In general, the tree positions in the CHM were close to the tree positions derived from the GPS; however, in a number of instances, the tree configuration apparent in the CHM was clearly different than the ground data. Plot locations were only shifted if an improved location fit was apparent. The combination of individual plot center and trees associated with that plot were moved manually and as a unit. Only x-y shift was applied (no stretching or warping) until the new fit between GPS trees and the CHM was optimal. We cross-checked the alignment of the shifted plot using tree heights as measured on the ground vs. the local maxima of the LiDAR-derived CHM. Table A1 shows the horizontal precision reported for each plot, as well as the distance of plot shift. Of the 60 plots, 45 were shifted an average of 3.3 m. The average shift distance for all plots, including those that were unshifted, was 2.5 m. Figure 3 shows an example case of plot shifting. In this example, some trees and part of the unshifted plot fall on the road, whereas trees in the shifted plot line up with the LiDAR-derived canopy height model. Each plot center was intentionally located at least 25 m from roads, showing that the plot center as measured by the GPS does not line up with the LiDAR data. Had the unshifted plot been used, model accuracy would have been reduced due to the inclusion of road in the LiDAR point cloud. Plots used a 16.93 m radius to produce a plot area of 900 m 2 . While a suite of measurements was taken at each plot, the only measures used by this study were a set of four photographs. A banner, measuring 4 m tall × 0.5 m wide with 1 m vertical gradations, was placed at plot center and a photograph was taken from plot edge at each cardinal direction towards the center with a T70h Android Quadcore Rugged Tablet. The banner was similar to a vegetation cover board used in many wildlife studies [39,40]. The photos were taken from eye-level, so camera height varied between 1.5 and 2 m. Figure 4A shows an example of one of these photos in an area with low ladder fuel cover. Plots used a 16.93 m radius to produce a plot area of 900 m 2 . While a suite of measurements was taken at each plot, the only measures used by this study were a set of four photographs. A banner, measuring 4 m tall × 0.5 m wide with 1 m vertical gradations, was placed at plot center and a photograph was taken from plot edge at each cardinal direction towards the center with a T70h Android Quadcore Rugged Tablet. The banner was similar to a vegetation cover board used in many wildlife studies [39,40]. The photos were taken from eye-level, so camera height varied between 1.5 and 2 m. Figure 4A shows an example of one of these photos in an area with low ladder fuel cover.  Figure 4A,B shows plots with low ladder fuel cover. Figure 4A's banner is unobstructed, while the banner in Figure 4B is partially obscured by tanoak leaves close to the camera lens, leading to only part of the banner being used for analysis, shown in Figure 4D.
Ladder fuel cover in each 1 m vertical segment of each photo was separately analyzed by four different technicians. Although we attempted photo processing automation similar to that used by Bennett, Judd, and Adams [69], we found that we were not able to accurately assess photos where large foreground objects occluded the camera lens, such as that shown in Figure 4B. We therefore utilized a standard ocular estimation technique for estimating canopy cover as a proxy for estimating ladder fuel cover from each photograph. Cover estimates were binned into five classes between 0 and 100%, shown in Figure 5. One cover class was assigned to each 1 m vertical segment of banner that best approximated ladder fuel cover as shown by the photo for that 1 m tall × 0.5 m wide section of banner.
To maintain consistency between estimates, photo-analysis periods were limited to 45 minutes. At the start of each session, technicians had to pass a cover estimation calibration test-scoring five correct estimates in a row. Correct classifications for analysts' tests were determined by an accurate assignment of a sample cover image, (one of the 30 shown in Figure 5), to a class designation (the six classes shown in Figure 5). To assist technicians with their estimates, they were provided with a printout of Figure 5, which served as a reference for their analysis.  Figure 4A,B shows plots with low ladder fuel cover. Figure 4A's banner is unobstructed, while the banner in Figure 4B is partially obscured by tanoak leaves close to the camera lens, leading to only part of the banner being used for analysis, shown in Figure 4D.
Ladder fuel cover in each 1 m vertical segment of each photo was separately analyzed by four different technicians. Although we attempted photo processing automation similar to that used by Bennett, Judd, and Adams [69], we found that we were not able to accurately assess photos where large foreground objects occluded the camera lens, such as that shown in Figure 4B. We therefore utilized a standard ocular estimation technique for estimating canopy cover as a proxy for estimating ladder fuel cover from each photograph. Cover estimates were binned into five classes between 0 and 100%, shown in Figure 5. One cover class was assigned to each 1 m vertical segment of banner that best approximated ladder fuel cover as shown by the photo for that 1 m tall × 0.5 m wide section of banner.
To maintain consistency between estimates, photo-analysis periods were limited to 45 min. At the start of each session, technicians had to pass a cover estimation calibration test-scoring five correct estimates in a row. Correct classifications for analysts' tests were determined by an accurate assignment of a sample cover image, (one of the 30 shown in Figure 5), to a class designation (the six classes shown in Figure 5). To assist technicians with their estimates, they were provided with a printout of Figure 5, which served as a reference for their analysis. Binned cover estimates, used to approximate ladder fuel cover. A canopy cover image series was adapted from Forest Inventory and Analysis (FIA) protocol [70], and binned into six cover classes. The cover estimates for each image set is shown at the bottom, and the average cover in each bin is shown at the top of each class. This was also the diagram used to assist technicians in assigning cover estimates.
To maintain a focus on ladder fuel fuels and eliminate the complications of perspective (a leaf near the camera lens sometimes covered an entire 1 m section of the banner), a 10% rule was established, whereby, if the banner area covered by a single object was more than 10% of a 1 m banner section, that area was removed from analysis. This addressed potential error introduced by tree boles in front of the banner, as well as leaves and smaller branches that were close to the camera lens and thus took up a larger than appropriate area on the banner. One example of this rule is demonstrated by Figure 4B (the raw photo) and Figure 4D (with unanalyzed area greyed out), where the banner is partially obscured by tanoak leaves close to the camera lens and the 10% rule was used. In this case, all large tanoak leaves were eliminated from analysis (shown in Figure 4D) because each leaf obscured over 10% of a banner 1 m segment. Ultimately, Figure 4A was analyzed as having 5% cover (the lowest class) in all height bins; Figure 4B was analyzed as having 20% cover from 0-1 m, 33% cover from 1-2 m, too little banner to determine cover from 2-3 m, and 27% cover from 3-4 m.
Because four independent technicians estimated the percent cover in each 1 m section of each photo, estimates often were not in perfect agreement. In cases where estimates differed by more than two consecutive bins (15% of samples), photos were re-assessed by a pair of technicians, who added a fifth measure to the set of estimates. Estimates for each 1 m segment of each photo were then averaged and used for further analysis.

LiDAR Data and Processing
Quantum Spatial collected aerial LiDAR between 23 May and 26 May 2015. A Leica ALS70 laser system mounted on a Cessna was used to collect LiDAR points utilizing a scan angle of ±15° from nadir. A pulse rate of 172 kHz, with unlimited returns per pulse (though often no more than 5) was used, averaging eight pulses per m 2 . IPAS TC v.3.1, ALS Post Processing Software v.2.75, Waypoint Inertial Explorer v.8.5, Leica Cloudpro v.1.2.1, and the TerraSolid software suite v.14 and v.15 were used to calculate point positions, classify points, and test spatial accuracy of the point cloud by Quantum Spatial [71]. The vendor completed real time kinematic (RTK) and post-processed kinematic (PPK) surveys and reported that average vertical and horizontal accuracy were 5.7 cm and 0.9 cm, respectively, based on the mean divergence of points from ground survey point coordinates (60 ground survey points were collected and compared to measured LiDAR ground points across the study area).
LiDAR metrics were extracted for each 900 m 2 plot using LasTools [72]. Metrics included basic point statistics (six metrics), percentile heights (eight metrics), and relative percent cover values (40 metrics), yielding a total of 54 LiDAR-derived metrics (see Appendix Table A2 for a detailed breakdown). Relative percent cover values were calculated for a range of 20 strata each for first Figure 5. Binned cover estimates, used to approximate ladder fuel cover. A canopy cover image series was adapted from Forest Inventory and Analysis (FIA) protocol [70], and binned into six cover classes. The cover estimates for each image set is shown at the bottom, and the average cover in each bin is shown at the top of each class. This was also the diagram used to assist technicians in assigning cover estimates.
To maintain a focus on ladder fuel fuels and eliminate the complications of perspective (a leaf near the camera lens sometimes covered an entire 1 m section of the banner), a 10% rule was established, whereby, if the banner area covered by a single object was more than 10% of a 1 m banner section, that area was removed from analysis. This addressed potential error introduced by tree boles in front of the banner, as well as leaves and smaller branches that were close to the camera lens and thus took up a larger than appropriate area on the banner. One example of this rule is demonstrated by Figure 4B (the raw photo) and Figure 4D (with unanalyzed area greyed out), where the banner is partially obscured by tanoak leaves close to the camera lens and the 10% rule was used. In this case, all large tanoak leaves were eliminated from analysis (shown in Figure 4D) because each leaf obscured over 10% of a banner 1 m segment. Ultimately, Figure 4A was analyzed as having 5% cover (the lowest class) in all height bins; Figure 4B was analyzed as having 20% cover from 0-1 m, 33% cover from 1-2 m, too little banner to determine cover from 2-3 m, and 27% cover from 3-4 m.
Because four independent technicians estimated the percent cover in each 1 m section of each photo, estimates often were not in perfect agreement. In cases where estimates differed by more than two consecutive bins (15% of samples), photos were re-assessed by a pair of technicians, who added a fifth measure to the set of estimates. Estimates for each 1 m segment of each photo were then averaged and used for further analysis.

LiDAR Data and Processing
Quantum Spatial collected aerial LiDAR between 23 May and 26 May 2015. A Leica ALS70 laser system mounted on a Cessna was used to collect LiDAR points utilizing a scan angle of ±15 • from nadir. A pulse rate of 172 kHz, with unlimited returns per pulse (though often no more than 5) was used, averaging eight pulses per m 2 . IPAS TC v.3.1, ALS Post Processing Software v.2.75, Waypoint Inertial Explorer v.8.5, Leica Cloudpro v.1.2.1, and the TerraSolid software suite v.14 and v.15 were used to calculate point positions, classify points, and test spatial accuracy of the point cloud by Quantum Spatial [71]. The vendor completed real time kinematic (RTK) and post-processed kinematic (PPK) surveys and reported that average vertical and horizontal accuracy were 5.7 cm and 0.9 cm, respectively, based on the mean divergence of points from ground survey point coordinates (60 ground survey points were collected and compared to measured LiDAR ground points across the study area).
LiDAR metrics were extracted for each 900 m 2 plot using LasTools [72]. Metrics included basic point statistics (six metrics), percentile heights (eight metrics), and relative percent cover values (40 metrics), yielding a total of 54 LiDAR-derived metrics (see Table A2 for a detailed breakdown). Relative percent cover values were calculated for a range of 20 strata each for first returns and all returns.
A visual representation of this suite of values is shown in Figure 6. Rasters of each LiDAR-derived metric were also generated at a spatial resolution of 30 m to be used to make predictions across the landscape. returns and all returns. A visual representation of this suite of values is shown in Figure 6. Rasters of each LiDAR-derived metric were also generated at a spatial resolution of 30 m to be used to make predictions across the landscape.

Statistical Analysis
Our workflow through plot and LiDAR processing is summarized in Figure 7. Average ladder fuel cover was calculated for each 1 m height band in each photo by averaging the estimates of the four photo-interpretation technicians. To determine whether plots were sufficiently sampled, the R software program [73] was used to calculate the standard error of ladder fuel cover among the four photos for a given plot, for each 1 m height stratum. This was done for a random sample of two-, three-, and four-photo samples to examine whether repeated sampling per plot decreased standard error and increased the reliability of this measure.
For all subsequent analyses, the four estimates for each 1 m height band (one from each cardinal direction) were averaged to produce a single estimate of ladder cover for each 1 m height band in each plot. These data were used to examine the distribution of ladder fuel cover in each 1 m height band across the study area. Ladder fuel cover estimates between 1 and 4 m were then averaged for each plot to produce a dependent variable, upon which to build a regression model using LiDAR. Cover below 1 m was left out of this measure since fuels in this range are commonly defined as surface fuels [74]. LiDAR is also less precise at differentiating ground points from low points within 1 m of the ground, and these points are often removed from analysis, though many different thresholds have been used [75][76][77][78].
Because LiDAR-derived independent variables were highly correlated, an iterative model building approach was taken. (1) The best regression models using one to four independent variables were chosen using the leaps package in R [79]. (2) Each model was evaluated for variable collinearity (correlation >0.6 in either Pearson or Spearman correlations was used to indicate collinear variables) and variable significance (p < 0.05). Only models where all independent variables were significant and none were collinear with one another were kept. (3) The leaps package was

Statistical Analysis
Our workflow through plot and LiDAR processing is summarized in Figure 7. Average ladder fuel cover was calculated for each 1 m height band in each photo by averaging the estimates of the four photo-interpretation technicians. To determine whether plots were sufficiently sampled, the R software program [73] was used to calculate the standard error of ladder fuel cover among the four photos for a given plot, for each 1 m height stratum. This was done for a random sample of two-, three-, and four-photo samples to examine whether repeated sampling per plot decreased standard error and increased the reliability of this measure.
For all subsequent analyses, the four estimates for each 1 m height band (one from each cardinal direction) were averaged to produce a single estimate of ladder cover for each 1 m height band in each plot. These data were used to examine the distribution of ladder fuel cover in each 1 m height band across the study area. Ladder fuel cover estimates between 1 and 4 m were then averaged for each plot to produce a dependent variable, upon which to build a regression model using LiDAR. Cover below 1 m was left out of this measure since fuels in this range are commonly defined as surface fuels [74]. LiDAR is also less precise at differentiating ground points from low points within 1 m of the ground, and these points are often removed from analysis, though many different thresholds have been used [75][76][77][78].
Because LiDAR-derived independent variables were highly correlated, an iterative model building approach was taken. (1) The best regression models using one to four independent variables were chosen using the leaps package in R [79]; (2) Each model was evaluated for variable collinearity (correlation >0.6 in either Pearson or Spearman correlations was used to indicate collinear variables) and variable significance (p < 0.05). Only models where all independent variables were significant and none were collinear with one another were kept; (3) The leaps package was used to find the best non-collinear independent variable to add to each model; (4) The process was repeated for all potential models until independent variables were no longer all significant; (5) The model with the best adjusted R-squared value was chosen from these models; (6) This model was tested for heteroscedasticity (where the variability across the range of a regression is not consistent across the range of values) using the Breusch-Pagan test [80]. R-squared, root mean squared error (RMSE), and a 10-fold cross-validation error using the CVTools package in R [81] were calculated to evaluate the fit of the final model. We repeated this model building for shifted, as well as unshifted plots to determine the difference in model fit between the two methods for this study.
used to find the best non-collinear independent variable to add to each model. (4) The process was repeated for all potential models until independent variables were no longer all significant. (5) The model with the best adjusted R-squared value was chosen from these models. (6) This model was tested for heteroscedasticity (where the variability across the range of a regression is not consistent across the range of values) using the Breusch-Pagan test [80]. R-squared, root mean squared error (RMSE), and a 10-fold cross-validation error using the CVTools package in R [81] were calculated to evaluate the fit of the final model. We repeated this model building for shifted, as well as unshifted plots to determine the difference in model fit between the two methods for this study. After a model was created, LiDAR-derived rasters of each independent variable used by the model were generated at 30 m resolution. This spatial resolution was chosen to match the plot size of 900 m 2 . A new 30 m raster of predicted ladder fuel cover was then generated for the entire study area.

Assessment of Photo-Banner
To determine whether each plot had been sufficiently sampled, we calculated the standard error between two-, three-, and four-photo estimates of ladder fuel cover at the same height within a single plot. Table 1 shows the distribution of these standard error measures of increasing sampling density across all plots for each height band. Table 1. Average standard error between measures of percent ladder fuel cover within a given plot and height bin are shown. We report the average standard error for two-, three-, and four-photo samples for each 1 m band in each plot. The variability of samples in the lowest height class (0-1 m) did not change much with increased sampling. With the lowest sampling density of only two photos per plot, standard error After a model was created, LiDAR-derived rasters of each independent variable used by the model were generated at 30 m resolution. This spatial resolution was chosen to match the plot size of 900 m 2 . A new 30 m raster of predicted ladder fuel cover was then generated for the entire study area.

Assessment of Photo-Banner
To determine whether each plot had been sufficiently sampled, we calculated the standard error between two-, three-, and four-photo estimates of ladder fuel cover at the same height within a single plot. Table 1 shows the distribution of these standard error measures of increasing sampling density across all plots for each height band. Table 1. Average standard error between measures of percent ladder fuel cover within a given plot and height bin are shown. We report the average standard error for two-, three-, and four-photo samples for each 1 m band in each plot.

Photos 3 Photos 4 Photos
Height The variability of samples in the lowest height class (0-1 m) did not change much with increased sampling. With the lowest sampling density of only two photos per plot, standard error was highest for the highest height class, at 9.94. However, for all but the height class within 1 m of the ground, standard error decreased with every subsequent photo sample added. Since we do not have more than four photos per plot, it is difficult to speculate where standard errors level off. However, based on the decrease in standard error of less than 1% from three to four photos in all height strata above 1 m, it appears that four photos may be sufficient. Average standard error for plots over 1 m, with four photo samples per plot, was only 6.98%. Considering that the estimates of ladder fuel cover were derived from six classes (shown in Figure 5), with most classes representing a range of 20%, the low standard error was encouraging. These low values suggest reasonable consistency between measurements, even between only two samples per plot.
Average ladder fuel cover is summarized in Figure 8 for each height class, as well as a combination of height classes from 1-4 m. Average 1-4 m ladder fuel cover was used in subsequent analyses to build a predictive model from LiDAR-derived metrics. The majority of plots had 1-4 m ladder fuel cover between 21% and 45%, with a maximum value of 91%, and minimum value of 7%, showing a good representation of sparse ladder fuels, and a weak, though still present, representation of dense ladder fuels.
was highest for the highest height class, at 9.94. However, for all but the height class within 1 m of the ground, standard error decreased with every subsequent photo sample added. Since we do not have more than four photos per plot, it is difficult to speculate where standard errors level off. However, based on the decrease in standard error of less than 1% from three to four photos in all height strata above 1 m, it appears that four photos may be sufficient. Average standard error for plots over 1 m, with four photo samples per plot, was only 6.98%. Considering that the estimates of ladder fuel cover were derived from six classes (shown in Figure 5), with most classes representing a range of 20%, the low standard error was encouraging. These low values suggest reasonable consistency between measurements, even between only two samples per plot.
Average ladder fuel cover is summarized in Figure 8 for each height class, as well as a combination of height classes from 1-4 m. Average 1-4 m ladder fuel cover was used in subsequent analyses to build a predictive model from LiDAR-derived metrics. The majority of plots had 1-4 m ladder fuel cover between 21% and 45%, with a maximum value of 91%, and minimum value of 7%, showing a good representation of sparse ladder fuels, and a weak, though still present, representation of dense ladder fuels. Figure 8. A box-and-whisker plot that displays the distribution of ladder fuel cover within each 1 m class across all plots is shown. We also include the average ladder fuel cover between 1 and 4 m in each plot.

Linking LiDAR to Ground-Based Measures
Using the iterative model-building approach described, we built a multiple regression using the shifted and unshifted plot centers. Shifted plot centers produced the model shown in Equation 1, where independent variables included the percentage of all LiDAR returns between 1 and 8 m (COV1_8), the standard deviation of LiDAR point heights above two meters (STD), and the percentage of first return LiDAR points between 8 and 16 m (FCOV8_16). Table 2  (1) Figure 8. A box-and-whisker plot that displays the distribution of ladder fuel cover within each 1 m class across all plots is shown. We also include the average ladder fuel cover between 1 and 4 m in each plot.

Linking LiDAR to Ground-Based Measures
Using the iterative model-building approach described, we built a multiple regression using the shifted and unshifted plot centers. Shifted plot centers produced the model shown in Equation (1), where independent variables included the percentage of all LiDAR returns between 1 and 8 m (COV1_8), the standard deviation of LiDAR point heights above two meters (STD), and the percentage of first return LiDAR points between 8 and 16 m (FCOV8_16). Table 2   For the above model, all independent variables, as well as the model as a whole, were significant at p < 0.05, with an R-squared value of 0.73, root mean squared error of 9.92, and 10-fold cross-validation error of 10.57. Figure 9 shows a scatterplot of the predicted versus measured values of ladder fuel cover, using Equation (1). For the above model, all independent variables, as well as the model as a whole, were significant at p < 0.05, with an R-squared value of 0.73, root mean squared error of 9.92, and 10-fold cross-validation error of 10.57. Figure 9 shows a scatterplot of the predicted versus measured values of ladder fuel cover, using Equation (1). Figure 9. Ladder fuel cover derived from the LiDAR point cloud is compared to ground-based measurements. This model is significant at the p < 0.05 level, with an R-squared value of 0.73, root mean squared error of roughly 10%, and 10-fold cross-validation error of roughly 11%. The ground-based measurements range between 7% and 91% ladder fuel cover.
A model was also created using data from the unshifted plot locations. Model fit was nearly identical, with significance at p < 0.05, an R-squared value of 0.73, and 10-fold cross-validation error of 11.02. While percentage of all LiDAR returns between 1 and 8 m still emerged as the primary model driver, other variables composed the rest of the model: height of the 90th percentile of LiDAR returns above 2 m, percentage of LiDAR first returns over 8 m, and percentage of LiDAR returns over breast height (1.37 m).
Equation (1) was used to predict ladder fuel cover across the entire study area, and is shown in Figure 10. Predicted ladder fuel cover ranged from −8% to 99% across the four focal areas, with an average of 36.6% and standard deviation of 15.9%. Although some areas had predicted ladder fuel cover under 0%, these only represented 0.12% of the total study area. While most areas have moderate to low ladder fuel cover, a few areas have pockets with very high ladder fuel cover. The top 10% of ladder fuel cover in the study area encompassed predicted densities between 58% and 100%, and were concentrated in focal areas C and D. Table 3 displays the average ladder fuel cover for each focal area, and the proportion of each focal area with exceptionally dense ladder fuel cover (in the top 10% of predicted values for the study area). Average ladder fuel cover ranged between 31% and 40% across the four focal areas, but area D had, by far, the greatest proportion of dense ladder fuel cover, having over four times as much as area A where ladder fuel cover was over the Figure 9. Ladder fuel cover derived from the LiDAR point cloud is compared to ground-based measurements. This model is significant at the p < 0.05 level, with an R-squared value of 0.73, root mean squared error of roughly 10%, and 10-fold cross-validation error of roughly 11%. The ground-based measurements range between 7% and 91% ladder fuel cover.
A model was also created using data from the unshifted plot locations. Model fit was nearly identical, with significance at p < 0.05, an R-squared value of 0.73, and 10-fold cross-validation error of 11.02. While percentage of all LiDAR returns between 1 and 8 m still emerged as the primary model driver, other variables composed the rest of the model: height of the 90th percentile of LiDAR returns above 2 m, percentage of LiDAR first returns over 8 m, and percentage of LiDAR returns over breast height (1.37 m).
Equation (1) was used to predict ladder fuel cover across the entire study area, and is shown in Figure 10. Predicted ladder fuel cover ranged from −8% to 99% across the four focal areas, with an average of 36.6% and standard deviation of 15.9%. Although some areas had predicted ladder fuel cover under 0%, these only represented 0.12% of the total study area. While most areas have moderate to low ladder fuel cover, a few areas have pockets with very high ladder fuel cover. The top 10% of ladder fuel cover in the study area encompassed predicted densities between 58% and 100%, and were concentrated in focal areas C and D. Table 3 displays the average ladder fuel cover for each focal area, and the proportion of each focal area with exceptionally dense ladder fuel cover (in the top 10% of predicted values for the study area). Average ladder fuel cover ranged between 31% and 40% across the four focal areas, but area D had, by far, the greatest proportion of dense ladder fuel cover, having over four times as much as area A where ladder fuel cover was over the 90th percentile. Since area D has the highest human population density, it is especially critical to identify these pockets of dense ladder fuel cover for fire planning and mitigation.
Remote Sens. 2016, 8, 766 13 of 23 90th percentile. Since area D has the highest human population density, it is especially critical to identify these pockets of dense ladder fuel cover for fire planning and mitigation. Figure 10. Ladder fuel cover, as predicted from the LiDAR point cloud using Equation (1), is displayed across the study area. Table 3. Average ladder fuel cover in the four focal areas, as well as the proportion of each focal area with ladder fuel cover in the top 10% (58%-100% ladder fuel cover) for the entire study area, are shown.

Discussion
We developed an improved ground-based method for quantifying ladder fuels using an ocular estimate of ladder fuel cover as a surrogate. LiDAR-derived variables included in our model were Figure 10. Ladder fuel cover, as predicted from the LiDAR point cloud using Equation (1), is displayed across the study area. Table 3. Average ladder fuel cover in the four focal areas, as well as the proportion of each focal area with ladder fuel cover in the top 10% (58%-100% ladder fuel cover) for the entire study area, are shown.

Discussion
We developed an improved ground-based method for quantifying ladder fuels using an ocular estimate of ladder fuel cover as a surrogate. LiDAR-derived variables included in our model were the percentage of all points between 1 and 8 m, the standard deviation of point heights above two meters, and the percentage of first return points between 8 and 16 m. The primary model driver was the percentage of all points between 1 and 8 m, a good representation of understory and midstory vegetation cover. This variable was the single best predictor of ladder fuel cover in models derived for both shifted and unshifted plots. It is curious that the most logical indicator of our field-derived ladder fuel cover, percentage of points between 1 and 4 m, was not a driver of the model. While the correlation between the LiDAR-derived 1-4 m and 4-8 m relative cover was not strong (see Table 2), the correlation of both of these variables with 1-8 m relative cover is high, indicating that this could better capture both the understory and midstory. This midstory may influence vegetative ladder fuels through light availability and the amount of competition with other plants. Light penetration can also influence the number of LiDAR pulses that reach the lower stratum with enough energy to register a return from understory vegetation. Less important independent variables are more contentious, since many LiDAR-derived variables are highly correlated (and one can become a better predictor than another with just a small change in plot data). However, standard deviation of point heights is a good indication of the overall forest structure, and has been used to predict forest structure [82]. Percentage of first return points is commonly used to estimate cover [60,75], so a stratum between 8 and 16 m is a good representation of upper canopy cover, since all but a few trees were under 16 m in height. While this is not a direct representation of ladder fuel cover, upper canopy cover could be important for light availability and LiDAR pulse penetration for the same reasons described above. Shifting plot centers to better match the LiDAR point cloud did not improve model accuracy, likely due to the large size of the plots and raster cells, and the small shift in most plots. However, we emphasize the importance plot shifting could have for smaller plots and still believe it to be a good check for larger plot shifts or those that could negatively influence the model. This method is a more quantitative approach to ladder fuels measurement than the five-class system developed by Menning et al. [13], and is also a promising method for change detection. Because our estimation relied on calibrated field photographs, our methods could be repeated in plots over time or after treatment or fire. The resulting series of photos taken in field plots could be utilized to record a quantitative change in ladder fuel cover, but could also serve as a visual indication of change, with a consistent and apparent measure of scale. In our study, area ladder fuel cover was fairly consistent within most plots, with low standard error between estimates. However, to minimize standard error in all plots in our study area, our results suggested collecting at least three samples per plot to bring the average standard error in plots (for 1 m bands over 1 m) below 8% and bring the maximum standard error near 20% (the size of most of our ladder fuel cover bins).
Using the relative ladder fuel cover derived from this ground-based measure, we built a model that predicted this value from the LiDAR point cloud with reasonable precision. This landscape-scale prediction also identified areas where there are pockets of uncharacteristically dense ladder fuels, indicating dense ladder fuels. Although it is out of the scope of this study to further investigate these areas, they indicate locations where efforts to mitigate ladder fuel hazard could be prioritized. While our measurement is not a direct input for current fire models, we believe that it holds much potential for landscape monitoring, fuels management, cultural resource identification and protection, and the evaluation of wildland fire behavior. This method is robust, quantitative, repeatable, fairly inexpensive (10 person-minutes in the field plus 16 in the office per plot (assuming four independent assessments per photo)), and is based on direct observation (i.e., does not rely on allometry).

Immediate Implications for Managers
Managers can use our methods for ground-based ladder fuel characterization to establish a baseline surrogate for ladder fuels across their management area and return to these plots to remeasure ladder fuel cover after change, similar to photographic monitoring used for evaluating fuels treatment and wildland fire effects [83]. Although our photographic methods require some analysis in the office, they are efficient to implement in the field, and provide a visual reference for each plot for further analysis or evaluation of change over time. Mechanical and/or prescribed fire fuel reduction treatments are planned throughout this study area, likely including several areas that overlap our plots. Plot remeasurement is planned post-treatment to assist in evaluating treatment effectiveness. Not only is this a quantitative measure of change, but also a visual indication, since plot photos, with consistent scale, are byproducts of plot measurement. This will allow for robust quantification of change in ladder fuel cover at the plot scale, and if LiDAR data collection is repeated, at the landscape scale as well.
Knowledge of the locations of areas with especially dense ladder fuels (such as those with ladder fuel cover in the 90th percentile for the study area) can be invaluable for prioritizing fuel reduction to facilitate community safety and fire preparedness. If areas of dense ladder fuels are near communities or evacuation routes (especially along routes that represent the only evacuation corridor for residents), these areas of dense ladder fuel cover could be targeted for fuel reduction or flagged as potentially hazardous for crown fire if a wildland fire were to occur in the area. Making these high-risk areas known to residents and fire managers could improve safety on large fires that threaten these communities in the wildland urban interface. Information about dense pockets of ladder fuels can also be useful to firefighters on the ground, for risk assessment of crown fire, ease of cutting fireline, and accessibility to and mobility through these areas, especially if they are to be traversed to get to a safety zone. Our results directly inform land management, including treatment prioritization at the landscape scale within the Western Klamath Restoration Partnership's Geographic Zones of Agreement [84].
Our work can augment work that addresses other aspects of wildland fire risk and modeling that can be predicted by LiDAR. Aerial LiDAR is an excellent predictor of topography (thus elevation, slope, and aspect), tree height, and canopy cover, all used extensively by wildland fire models [1]. Canopy bulk density has also been predicted by LiDAR [50,58,60,85]. Surface fuel model has been predicted in some areas [56,57,86], but still presents a challenge in denser forest types on steep, topographically complex terrain [26]. CBH has also been estimated from LiDAR [25,87], completing the suite of landscape metrics required by many wildland fire behavior models. However, because neither fuel model nor CBH is a direct measure of ladder fuels, we encourage managers and modelers to consider incorporating a metric that is better able to account for ladder fuels, both for land management and fire modeling. While ladder fuel cover is also not a perfect surrogate, it represents an alternative method of measuring ladder fuels that focuses on the fuels, rather than tree crowns, to estimate ladder fuels.
For wildlife that require understory structural complexity such as fishers, these data could also be valuable [88]. From a cultural perspective, which is especially important in the tribal communities surrounding the study area, these results will be used to mitigate potential impacts from fire near valuable cultural resources, both historic and present [64,89]. These resources could include archaeological, heritage, and tribal sacred sites, as well as legacy oaks and pines that provide a traditional food source, but are also critical habitat elements [90,91]. Thinning and prescribed fire treatments to enhance these cultural resources that are at elevated risk of intense wildland fire will be facilitated by our predictive map of ladder fuel cover across the landscape [92]. Because areas with lower densities of surface and ladder fuels facilitate enhanced accessibility, these areas are preferred by traditional gatherers, as suggested by Hummel and Lake [89]. Our work to predict ladder fuel cover could be used to identify especially open understories that may be valuable for identifying gathering locations for (non-timber) forest products and associated ecological services for the community, tribes and public [93,94].
Beyond the immediate implications of this work, this photographic technique, as well as the relationship we built between ground measures and the LiDAR point cloud, represents a new method to quantify an important fuel source that has been so elusive to previous fire behavior models that an allometric proxy was implemented in the place of a quantitative measure based on fuels actually present. Although Kane et al. [95] were not able to find a significant relationship between standard forest structure metrics and fire severity, we encourage further investigation into whether this estimation of ladder fuel cover has any predictive power over fire behavior and effects. It will be particularly useful to test whether these methods for estimating ladder fuel cover can be used to explain observed fire severity in areas where pre-fire LiDAR data are available.

Study Limitations
While photo interpretation methods were robust, we did not show that the estimate of ladder fuel cover relates to any specific measure on the ground, such as ladder fuel bulk density. Therefore, this model was only useful for predicting the relative ladder fuel cover, and should not be used to predict specific quantities of fuel. This method did not account for vertical continuity of fuels (such as a tree crown reaching all the way to the ground and providing a continuous "ladder" for fire to climb). Vertical continuity was addressed by Menning and Stephens [13], but remains a difficult parameter to incorporate, both for sampling methods and sampling density, and was out of the scope of this study. Photo interpretation methods also contained numerous inherent sources of error, including: (1) slope distorting the apparent height of vegetation; (2) objects in the foreground appearing disproportionately large; (3) using standard focal length settings, which could alter relative object sizes in the photo [96]; and (4) photo interpreter error. We also noticed a moderate degree of inaccuracy, averaging 2.5 m (up to 15 m for some plots), between the LiDAR and the ground-based measures of tree locations. While we attempted to mitigate this error by shifting plot center GPS coordinates based on the mapped height of tall trees in relation to plot center, there are likely still inaccuracies between the plot measured on the ground and the LiDAR point cloud extracted from each plot perimeter. Furthermore, only a single study area was used, introducing the potential for locational bias in our analysis. Because our study area was a single sample of the landscape, results from this study should not be applied to other areas without first collecting ground reference data, testing the relationship, and adjusting where necessary to ensure accurate model estimates.
This research also emphasizes the need for better accuracy assessment of plot center coordinates for LiDAR-based studies. LiDAR vertical accuracy was reported to be 5.7 cm. Despite our use of a high-accuracy GPS with antenna, collecting at least 1000 estimates per plot, using differential correction, and receiving an average horizontal precision of 87 cm, plot centers were still mis-aligned with the LiDAR data by an average of over 2.5 m, and up to at least 15 m. We urge those that work with LiDAR to assess how well recorded plot centers line up with the LiDAR point cloud and take measures to better match the plot to the LiDAR before undertaking further analyses. This is especially important when plots are located in heterogeneous landscapes or near roads, and when plot sizes are small (so a small shift in plot center equates to a dramatically different area of the LiDAR point cloud). While some LiDAR-derived metrics are not as sensitive to shifts in plot center location [97], many likely are, and plot accuracy should be better accounted for in all cases.

Future Research
This research demonstrated the importance of matching plot data to the LiDAR point cloud. We urge others to pursue further research into the best practices of measuring and addressing error when working with the LiDAR point cloud.
We are eager to further explore the ability of LiDAR to predict ladder fuel cover, and plan to work more extensively with land managers and community members in this area to discuss ways to best utilize the products from these analyses. We urge others to experiment with our methods, not only as a way to estimate ladder fuels, which are important to the field of wildland fire and tribal cultural resources [92], but also as a monitoring method for plot vegetation and fuel loading, as well as for the enhancement of access to and the quality of cultural resources.
We encourage others to test this metric against fire behavior to determine its influence on the likelihood of crown fire, and hope that fire behavior models will incorporate some version of this metric to further refine future wildland fire behavior and effects models.

Conclusions
Ladder fuels are an important piece of fire behavior and effects modeling. However, because they are difficult to quantify in the field, ladder fuels are indirectly measured through allometric equations that estimate canopy base height, leaving out the shrubs and small trees that are a large component of ladder fuels in many forests. Our work presents a new ground-based method for estimating ladder fuels and establishes a strong predictive link with aerial LiDAR. We were able to predict our ground-based estimates of ladder fuel from aerial LiDAR with an R-squared value of 0.73. We encourage others to explore this relationship in other areas and use our estimates of landscape scale ladder fuel cover to reduce wildfire risk to communities in our study area through targeted treatment. We also demonstrate that plot center accuracy must be independently tested, and urge LiDAR users to critically compare their plot locations to the LiDAR point cloud.

Abbreviations
The following abbreviations are used in this manuscript:

Variable Group Individual Variable Description
Basic point statistics