Large Area Forest Yield Estimation with Pushbroom Digital Aerial Photogrammetry

Low-cost methods to measure forest structure are needed to consistently and repeatedly inventory forest conditions over large areas. In this study we investigate low-cost pushbroom Digital Aerial Photography (DAP) to aid in the estimation of forest volume over large areas in Washington State (USA). We also examine the effects of plot location precision (low versus high) and Digital Terrain Model (DTM) resolution (1 m versus 10 m) on estimation performance. Estimation with DAP and post-stratification with high-precision plot locations and a 1 m DTM was 4 times as efficient (precision per number of plots) as estimation without remote sensing and 3 times as efficient when using low-precision plot locations and a 10 m DTM. These findings can contribute significantly to efforts to consistently estimate and map forest yield across entire states (or equivalent) or even nations. The broad-scale, high-resolution, and high-precision information provided by pushbroom DAP facilitates used by a wide variety of user types such a towns and cities, small private timber owners, fire prevention groups, Non-Governmental Organizations (NGOs), counties, and state and federal organizations.


Introduction
The Forest Inventory and Analysis program (FIA; see also Appendix A, Glossary of terminology and statistical notation) within the US Forest Service (USFS) maintains and manages the national grid of forest inventory plots for the United States [1]. It comprises a grid of fixed area plots designed to cover the USA, with an approximate sampling intensity of one plot per 2430 ha [2]. This rich data resource is designed to support inferences from areas as broad as entire states down to individual counties. There are a variety of applications, however, that would benefit from consistent and precise forest measurements and monitoring across the USA at a finer than county resolution. Applications include (e.g.) mapping and small area estimation of fuels, forest yield, cover, and change. Researchers have mapped forest attributes over broad areas at moderate to fine resolutions with satellite-based imagery, such as MODIS and LANDSAT [3,4]. These products are improvements in the sense of providing higher-resolution spatial maps, but they have limited ability to quantify variation in vertical forest structure relative to what is feasible with remote sensing technologies such as lidar and photogrammetry. Lidar capabilities have been widely demonstrated to model forest height [5], volume [6], biomass [7], carbon [8], canopy fuels [9], and other height and cover related forest attributes [10].
While lidar has proven effective for modeling a variety of vertical structure related forest attributes, a well-known limitation of lidar is its high cost, especially at the scale of entire states. Recent lidar prices in the USA typically range from approximately $0.40/ha to $7/ha. Even at the low end, it would cost approximately 125 million USD to cover all USA forests a single time. While there is currently an effort underway by the USGS [11] to acquire nationwide lidar coverage, it is a long term project, and the return interval may exceed what is useful for monitoring vegetation. Digital Areal Photogrammetry (DAP) is an alternative remote sensing technology which also has the capacity to support large-area forest structure measurements. This technology is commonly referred to by a number of names including Structure from Motion (SfM) and, colloquially, phodar, an inaccurate combination of the term "photo" in photogrammery and Detection and Ranging "DAR" from lidar. The imagery used to support DAP can range in price from less than $0.01/ha as an in an archival project to over $2 USD/ha depending on acquisition specifications. For example, new stereo imagery which is used for the National Agriculture Imagery Program (NAIP) is available every 2-3 years over the entire conterminous USA, and requesting digital surfaces models for the acquisitions costs less than $0.01/ha. DAP-based measurements of forest structure have been shown to be slightly worse than lidar in terms of performance for measurements such as height, basal area, and volume [12][13][14]. However, unlike lidar, a requirement for using DAP effectively is to have an alternative source for a Digital Terrain Model (DTM) because DAP typically measures the first surface, and therefore does not measure the ground surface beneath vegetation. In contrast, lidar can provide both measurements of the vertical forest structure and a detailed representation of the ground surface. Lidar data in the public domain can serve as a reliable source of DTMs for DAP, but much of the USA is not covered with lidar. An alternative source of DTMs for the USA is the US Geological Survey (USGS) National Elevation Dataset (NED) which is derived from diverse sources and "… is a seamless dataset with the best available raster elevation data of the conterminous United States, Alaska, Hawaii, and territorial islands" [15]. The NED rasters for the USA are approximately 10 m horizontal resolution and 1.5 m-2.4 m vertical RMSE [16]. NED elevations are inferior to most modern lidar acquisitions with respect to resolution and vertical accuracy, but they are often the best free source of DTMs for expansive areas in the USA. One question of interest then, is how well the USGS NED DTM can be used with DAP to support forest inventory applications.
For this study, DAP was derived from airborne pushbroom sensor imagery. The motivation to collect data with pushbroom sensors is their wide swath width and multiple look angles, which support stereo (from two look angles) from a single over-flight. A noted limitation with pushbroom imagery relative to frame imagery is increased difficulty to extract a detailed, artifact-free canopy surface model. Artifacts include smoothing or omission of surface features such as lone trees, canopy edges, and canopy gaps. The challenges encountered with pushbroom imagery are discussed in greater detail in the methods section.
The objective of this investigation is to evaluate forest yield estimation with pushbroom DAP and FIA field measurements. The performances of the DAP-based methods were assessed by determining whether estimation with auxiliary height information increased the efficiency (estimated variation relative to sample size) of estimation. Multiple components of the estimation strategy are examined, including the estimator used (regression versus post-stratification), plot positioning accuracy (low versus high), and the source of DTM (lidar versus NED).

Study Site
This study was conducted for forested areas in the state of Washington (Figure 1), the Northwestern-most state of the conterminous USA. Washington State is a complex mixture of forest and non-forest conditions that spans extreme moisture and elevation gradients, including temperate rainforest, shrub-steppe, desert, alpine forest, five or more mountain ranges, and numerous volcanic peaks, including Mount Rainier, which rises from near sea level to 4352 m.

Field Measurements
Field measurements for this study were collected on Forest Inventory and Analysis (FIA) forest inventory plots measured by the USFS. Measurements in Washington state are taken annually over a 10 year period, but paneled such that inference can be made from FIA measurements taken in any given year (where a panel consists of a regularly spaced subset of the full grid; see Reference ( [2] p. 20) for an expanded explanation of panels). The observations used for this study were measured during the 2015 field season.

Plot Positioning
Starting in 2014 the Pacific Northwest (PNW; including CA, OR, WA, AK, HI and the Pacific islands) FIA unit has implemented a program to collect high-precision locations (HPL) on all PNW field plots. The objective of the HPL program is sub-meter precision for sub-plot positioning. We will henceforth use the term High Precision Locations (HPL) to refer to these updated coordinates. Field personnel collect HPL plot positions on sub-plot centers for a minimum of 15 min while collecting data at 1 Hz with dual-frequency GPS/GLONASS GNSS receivers. The majority of plot positions were obtained with Trimble (The use of trade or firm names in this publication is for reader information and does not imply endorsement by the U.S. Department of Agriculture of any product or service.) Geo 6000 receivers, with a small number of plot positions obtained with Trimble Geo 7X and Javad Triumph 2 receivers. GNSS data for plots were differentially corrected using a single CORS [17] base station. Currently, the majority of positions are obtained with receivers placed directly on the ground, in part due to the fact that the Trimble Geo 6000 and Geo 7x receivers do not have the means to thread onto a tripod directly. We expect the horizontal precision of plot locations surveyed under dense 55 m tall conifer forests with Trimble receivers to be better than 1.85 m on average [18] using the described protocols and receivers. Many plots are measured in the open and under less dense canopy conditions than in Reference [18] which we expect to have better GNSS performance. Given the wide variation in field conditions, we were not able to accurately estimate overall HPL plot precision, except to say that it is expected to exceed 1.85 m.
Prior to implementation of the HPL program, the Standard Precision Locations (SPL) for FIA plots were obtained using a variety of approaches and the precision of SPL coordinates is not known. Previous methods for plot positioning included using the target plot location, digital photo interpretation, use of recreational grade GNSS receivers, and use of Rockwell PLGR GNSS receivers. One of the objectives of this study is to understand the degree of improvement in coordinate precision and in the estimation that is afforded by the HPL program relative to SPL coordinates. For analyses of positioning, we treat the HPL locations as the "true" plot locations. Positional errors are computed as: Since in Equation (1) above HPL coordinates are measured with error, values are slightly upward biased for the "true" positioning errors: , ℎ true coordinates. ( Our discussion of positioning involves only precision and not accuracy because we reasonably assume that, on average, HPL and SPL coordinates are unbiased for the true plot locations, e.g.,

Photogrammetric Point Cloud (Photo Points)
Imagery collected for photogrammetic work was collected as part of the Hexagon Content Program and is used to create orthophotos for the National Agricultural Imagery Program [19]. Photo points were derived from pushbroom imagery collected by a gyro-stabilized Leica ADS100 pushbroom sensor in the summer of 2015 for Washington State. The 4-band sensor has a 62.5 mm focal length, and three look angles including 25.6 degrees forward, nadir, and 19.4 degrees backward. There is a total of 13 lines of pixels in the CCD arrays (each look angle) and they are 20,000 pixels wide. Each pixel is 5 um wide. Only data from the nadir and backward look angles were processed for this analysis. The average flying height was 5000 m which resulted in a 40 cm nominal on-theground pixel resolution. Digital Surface Models (DSMs) were created from the pushbroom imagery using the Auto Spatial Modeler (ASM) module in BAE Socet GXP software. The ASM module produces an irregular TIN surface, which was converted into the LAZ point format [20]. The resulting point cloud was also 40 cm resolution, or approximately 6.25 points per square meter. The horizontal accuracy of pushbroom-derived orthophotos is expected to be better than 2m for the NAIP program [21]: The source of pushbroom imagery for this study. Vertical accuracy is not well documented for point clouds derived from NAIP pushbroom imagery and will vary depending upon a variety of factors including target surface type and illumination.
Although we eventually succeeded in generating a forestry suitable point cloud (e.g., Figure 2), it proved exceedingly difficult to achieve satisfactory performance generating a point cloud from pushbroom imagery with photogrammetric software. In our experience with modern digital framebased cameras, satisfactory performance was achieved with every software we tested, including Trimble Inpho, Erdas Photogrammetry Suite, and Agisoft Photoscan. Not all photogrammetric software natively supports pushbroom imagery, but multiple lines of pushbroom can be combined into pseudo-frames to match the expectations of software that cannot natively import Leica ADS100 imagery.
Amongst the software tested, the BAE Socet GXP ASM module was the only software that created a point cloud which had fine details such as individual tree crowns. Point clouds generated by other software had excessive defects, including heavy smoothing, lack of penetration into canopy gaps, and diagonal interpolation of canopy edges to the ground. BAE Socet performance was initially similar to that of other software until we developed a heuristic approach to search for suitable software parameters in the ASM configuration file. To achieve a satisfactory configuration, we automated the generation of hundreds of thousands of alternative configuration scenarios and analytically compared the resulting point clouds with lidar measurements for the same locations. This enabled us to identify configuration settings optimized for our forest conditions. The same configuration file has also been used satisfactorily for data from acquisitions in other states such as Oregon, California, and Wyoming. Despite having improved upon the original outputs, the final product still exhibits defects, especially the erasure of lone trees in openings. Tree omission errors are also known to occur frequently with frame-based imagery. The authors will make the configuration file containing approximately 100 settings available to readers on request.

Remotely-Sensed Structure Measurements
Processing of remotely sensed information (3-dimensional point clouds) was completed primarily with executables in the FUSION software bundle developed by the US Forest Service PNW Research Station [22]. Height metrics were computed to generate wall-to-wall rasters covering the extent of the remote sensing using the "gridmetrics" tool and height metrics for individual plots were computed using the "cloudmetrics" tool. Point metrics used in this study were computed by differencing canopy elevations from DTM elevations. Metrics evaluated for this study included height percentiles (e.g., 10, 20, …, 100 percentile heights of points above 2 m) and the proportion of returns above 2 m (cover).
The source of canopy surface points used in this study came from stereo pushbroom imagery, however, as is evident in the rightmost side of Figure 2, the ground is often obscured under dense continuous canopies. As a result, we had to rely upon other sources of information for a digital terrain model. We evaluated two sources of DTMs including the National Elevation Dataset (NED) produced by the USGS, and a composite of publicly-available and accessible lidar DTMs for WA State ( Figure  3). The NED DTMs are rasters for the USA with a maximum resolution of 10 m and a vertical accuracy (RMSE) ranging from 1.5m-2.4m for forested areas [16]. The lidar-derived DTM was composited from many projects [23] ranging from the years 2000-2017 with a wide range of acquisition specifications. The composite lidar DTM has a 1 m horizontal resolution. The vertical accuracy depends upon the contributing projects, but typically exceeds 30 cm. With considerable exceptions, lidar acquisitions are concentrated in large part in Western Washington around urban environments and in low-elevation forested areas. Terrain mapping including geologic hazards is the dominant motivation for lidar collections in Washington State.
To facilitate direct comparison between estimation with lidar and NED DTMs, inference was restricted to areas which intersected lidar DTMS and forests ( Figure 3). NED DTMs (extent not shown) are available for the entire study area.

Forest Yield Estimation
Performances were compared for three estimation strategies including regression, poststratification, and for estimation without auxiliary information. Each of the estimators was implemented consistent with data collected as a Simple Random Sample (SRS).
Post-stratification was implemented using a regression estimation approach. Post-strata were treated as categorical variables in a linear regression model. As a result, the same formula for the variance estimator of the regression mean, ̂ , was also used for the variance estimator of the post-stratification mean, ̂ ). Strata mean estimates obtained with the regression approach are identical to strata estimates with a traditional stratified sampling estimator, but variance estimates differ because deviations around the strata means are not weighted by the areas of the strata in the regression approach. The advantage of this approach is that it is easier to compare performances of estimation strategies when estimators are formulated similarly. A convenient result when performing estimation for an SRS design is that formulae for estimators can be simplified because the sampling probability is the same for every element in the population. For example, the formula for an unbiased estimator of the population mean for an SRS does not depend on sample selection probability (shown in Reference [24] for the total): We estimated the variance of the mean estimator (̂ ) with: A regression mean estimator for a probability sample is [24]: When used with Ordinary Least Squares (OLS) regression, the mean estimator ̂ , can be simplified for an SRS design to: The ̂ and ̂ , estimators are equivalent for an SRS design because the second term in ̂ , is equal to zero: We used the following variance estimator for ̂ [25]: There is a clear practical advantage to sampling with equal probability in that the estimators above are simplified. While not explicitly dealt with in this study, it is worth noting that there is an added advantage to regression estimation for SRS in that aggregations of mapped predictions from sub-areas (or sub-domains) can be summed to the total. This is of practical importance because it is common to build raster-based forest yield maps. For a sample collected with an SRS design, aggregates of mapped predictions from an OLS model are compatible with the estimate of the total. For example, to create a synthetic estimate for a county (a synthetic estimator is an estimator which depends upon a relationship for a population to make inferences for a sub-region; e.g., "MLR_A" in Reference [26]), one can simply sum model predictions for a county: If completed for each county in the state, the sum of the county totals will be equal to the modelassisted state total estimate without any adjustments (for an SRS). Although county or other sub-area synthetic estimates are design biased for the sub-area and we cannot use formula (9) to estimate variances for sub-areas, it is convenient that they sum to the model-assisted SRS total ̂ without adjustment, or simply:

Model Fitting and Post-Strata
Post-stratification and model fitting were implemented with simple approaches. For example, we did not attempt to compare all combinations of explanatory variables (e.g., Reference [27]) in an attempt to minimize residual errors. Instead, we used interpretable and defensible explanatory variables including a measure of height and a measure of cover and their interaction to represent the phenomena of increasing forest height and cover, which should describe much of the variation in tree size and presence. Post-stratification was implemented using a raster composed of a single heightbased explanatory variable, the 70th percentile height (ht_p70) for the entire study area, computed after excluding points below 2 m in height. We used 70th percentile height because it has most of the explanatory power of larger height percentiles (e.g., 80th, 90th, 95th), but is more robust to "noise" or blunders that are common in DAP. Bin thresholds are based upon deciles computed for all ht_p70 values in the study area.
In Figure 4 we can see how field measured volume is distributed relative to binned remotely sensed heights. Median volume clearly increases with ht_p70, and the relationship is both non-linear and heteroskedastic (unequal variance). More complex post-strata may use multiple remotely sensed forest attributes [28], or even area-wide volume or biomass predictions when mild gains in estimation are deemed to outweigh the disadvantages of added complexity (References [29,30]).

Efficiency Assessment
We evaluate performances in terms of efficiency for multiple estimation strategies where an "estimation strategy" consists of an estimator paired with a dataset. Efficiency is defined as the sampling variation for a given number of sample observations. An estimation strategy that is more efficient than another will have lower variance than another estimation strategy for the same number of observations, all other things being equal. There is not, however, a single way to measure or report efficiency.
In this study we primarily quantify efficiency with Relative Efficiency ( ). Relative efficiency, as defined here, is the estimated ratio of SRS mean variance without auxiliary information over the mean variance for an alternative estimator: ̂ , is the variance estimator for some alternative design (12) An advantage to using relative efficiency for inference is that it can be interpreted as a function of the sampling effort (e.g., the number of measured field plots). Relative efficiency is the estimated proportion of observations needed with ̂ to achieve the same precision as with ̂ , . For example, if a dataset has = 2.0 for some alternative mean estimator ̂ , then the variance estimated for ̂ with n plots is equivalent to the variance estimated for ̂ with 2 × plots (̂ requires twice as many plots as ̂ , ). This makes it possible to estimate the added value of ̂ in term of the cost required to measure additional plots with ̂ to be as precise as ̂ , . Our inferences rely largely on comparisons between the relative efficiencies of estimation strategies.
A common measure of statistical performance in modeling is the coefficient of determination ( ) which can also be used to measure efficiency. The interpretation of is the (estimated) proportion of variation explained. Using the same notation as for relative efficiency, we define: From a sampling perspective, the main inference that can be drawn directly from values is that the estimation strategy with a higher value is more efficient than another strategy. Advantages of reporting R 2 include that it is commonly used in the literature related to area-based forest modeling with lidar (e.g., Reference [31]), its ease of interpretation, and its stability for high values of equivalent relative efficiencies.
As can be seen in Figure 5, relative efficiency and R 2 are closely related statistics and there is an exponential relationship between relative efficiency and R 2 . While there are advantages to using relative efficiency for inference in a sampling context, the exponential behavior of relative efficiency can be a weakness when interpreting values on the order of 5 or greater. The issue is mitigated somewhat with a sufficiently large sample size (RE becomes more stable), or by simply recognizing that high relative efficiency values may be unstable and building some degree of caution into inferences.

Coordinate Precision
One of the comparisons examined in this study is between estimation with HPL (High Precision Location) versus SPL (Standard Precision Location) coordinates. In Figure 6 we can see that positioning precision is a relevant concern, as 5% of observations appear to exceed 38 m in positioning error. For these plots, there would be no overlap between the remote sensing footprint and the location on the ground where field measurements were obtained. We investigate this issue further in later sections by looking at the impact of coordinate error on estimation. Substantial coordinate error may result in high variability between field and remotely-sensed forest structure measurements.   Figure  7 suggest that the association between volume and DAP height metrics is strongest for HPL with a lidar DTM, with the greatest reduction in performance resulting from using a USGS DTM, followed by using SPL coordinates. Although we can see that DTM source and coordinate quality have demonstrable impacts on the relationship between volume and height, even in the worst case (USGS DTM and SPL coordinates), there is still a strong relationship between the variables. Visually, the relationships and positions of points in the scatterplots appear very similar with only slight differences in spreads and in trends.  Figure 8 demonstrates one of the practical advantages of having high-resolution auxiliary height information to support forest yield estimation. When paired with precise, fine-resolution remote sensing, forest inventory measurements can be used to make precise, fine-resolution maps of forest attributes such as volume. FIA measurements without auxiliary height measurements support county-level estimates. In Figure 8, for example, Mason County (in green) would have a single estimate of total volume. The mapped product, in contrast, provides a very fine level of detail, describing the spatial distribution of volume within the county. While at times useful, it is important to note that aggregates (e.g., mean or sum) of mapped predictions are not design-unbiased for the stand, stratum, county, or any other tessellation of the population. Asymptotically unbiased estimation is feasible at the county level with the appropriate estimator, although biased synthetic estimates may be sufficient for some applications (e.g., Reference [32]).  Table 1 provides the efficiencies estimated for multiple estimation strategies. Standard error (SE) and R 2 values are also provided, although all of our inferences are based on RE. SE and R 2 values in Table 1 are present as a reference to enable comparisons with other studies. The difference in estimation strategies in Table 1 are based upon estimator (regression or post-stratified), type of plot positioning (HPL or SPL), and source of DTM (USGS or lidar). The results with respect to HPL plot locations and use of a lidar DTM are fairly intuitive. Both the lidar DTM and HPL were presumed to be higher quality data sources, which was confirmed when estimation with the lidar DTM and HPL coordinate data resulted in better relative efficiencies than the alternatives. For post-stratification, estimation with SPL caused the greatest singular reduction in efficiency. Estimation with a lidar DTM and HPL is equivalent to having 4.16 times as many plots as an SRS, but estimation with a lidar DTM and SPL is only equivalent to 3.47 times as many plots. Using SPL coordinates had the greatest impact on post-stratification, likely because even minor changes in height due to registration error could move a plot into an entirely different height stratum. The regression estimator was fairly insensitive to HPL versus SPL. The post-stratification approach was more sensitive to changes in estimation strategy: It yielded both the best result (HPL with a lidar DTM) and the worst result (SPL with a USGS DTM).

Terrain Models
The results indicate that the best performance is achieved with a lidar DTM, although poststratification on height with at USGS DTM is still 3.57 times as efficient as estimation without DAP heights. This agrees with results from another study which found that for a boreal forest, results using DAP based heights computed with a 10 m DTM compared well with using DAP with a lidar DTM [33]. We did not, however, explore what is the best way to incorporate two separate sources of DTMS into a single large-area remote sensing augmented forest inventory. Two approaches we considered include 1) mosaic lidar and USGS DTMs and give priority to the lidar DTMs, and 2) implement the inventory separately for areas with lidar and USGS DTMs. In the case of the second approach, this means fitting models separately for the two distinct data types, and then mosaicking results and estimates back together at the end. Metrics computed on lidar and USGS DTMs may be sufficiently different that the models may fit better in approach 2, although clearly this warrants further testing.
Another approach worth real consideration based upon the result observed here is implementation using only NED DEMs. This would be the simplest approach and would eliminate any edge artifacts that could arise due to differences between lidar and USGS DTMs. Clearly this would sacrifice a modicum of performance in areas with lidar DTMs and could be especially concerning for fine-resolution map products.

Positional Accuracy
Plot positioning under the canopy has proven to be a difficult endeavor over the history of forest measurements, and even today precise plot locations under canopy requires considerable investments in technology, planning, and GNSS data processing. Even predating the current HPL program, PNW FIA invested heavily in plot positioning, from having field crews attempt accurate pin-prick plot centers on orthophotos to using military-grade GPS equipment. As a result, the distribution of positional errors prior to HPL was actually encouraging. On average, 50% of plots were found to have better than 8 m error and 95% had better than 38 m of error (treating HPL coordinates as "correct"). However, these errors had a demonstrable effect on estimation performance, especially when using post-stratification.
The observed efficiencies for an area-based approach with SPL demonstrate one of the advantages of area-based inference over individual tree analysis approaches in the presence of positioning errors. While PNW has invested in precise plot locations on recent plot measurements, older recordings and other parts of the country do not have precise plot locations. Our results indicate that while performance will decline, plots with poor positioning accuracy can still be used with DAP to improve estimation. However, the observed SPL positing errors are too large to be expected to reliably align with individual trees.
Plot positioning errors (SPL) are sufficiently large that the plots with the worst errors likely contribute a large amount of noise. For example, the worst 5% of plot positions had positioning errors of 38 m or greater. The magnitude of errors observed in plot locations in Figure 7, and the demonstrated loss of estimation efficiency without HPL plot locations demonstrates the importance of HPL positioning to the implementation of remote sensing augmented forest inventory approaches. From this perspective, not having HPL may be equivalent to replacing 5% of observations with random noise, and commonly the worst plot positions fall in areas with the greatest topographic relief the largest trees, where forest types with the largest trees are commonly of the greatest interest.

Estimators
The best performance was observed for estimation with post-stratification, although it could not be deemed a clear "winner". The regression estimator appeared to be more stable with respect to positioning errors, which would be advantageous for datasets for which HPL was not available. This makes sense because unlike post-stratification, a well-fit regression model on a continuous explanatory variable does not cause sharp changes in predictions for slight changes in the explanatory variable, which can happen with post-stratification when near the borders of the strata. That said, the magnitude of change between height bins is an artifact of the number of height bins, and for sufficient numbers of bins, regression and post-stratification would essentially converge. However, increasing the number of strata requires a sufficient number of observations to be able to reliably estimate strata means for a large number of height bins, something that will not always be feasible, but the point remains that we may be able to refine our approach to post-stratification to mitigate the sensitivity of the strategy to positioning errors. For example, it may help to post-stratify on predicted volume, or some other linear combination of our explanatory variables, further enmeshing the post-stratification and regression estimation as in References [29] and [34].
From a practical standpoint, post-stratification is much more convenient to implement than regression estimation. FIA already uses post-stratification to deal with forest/non-forest conditions and to adjust for non-response bias. From the FIA perspective, post-stratification could be a drop-in replacement or augmentation for the post-stratification that is already in place.
Regression estimation is more complex to implement and requires making a number of additional decisions. The best performance for a single attribute will be achieved by fitting each attribute separately, i.e., a separate linear regression model with a separate set of predictors. This tends to be unwieldy for an inventory program that estimates hundreds or thousands of attributes. The preferable alternative is to sacrifice slightly on performance and use the regression relationship to estimate a weight for each observation, and then apply the weight identically for all attributes. This approach to estimation is called regression calibration [35]. People have also turned to nearest neighbor algorithms to deal with simultaneously predicting large numbers of attributes from remote sensing [36]. Nearest neighbor approaches are convenient for managing the complexity of forest attributes, but there is currently limited documentation for implementation in a design-based framework. In practice the results (maps and estimates) are likely to be fairly similar between poststratification and regression calibration, suggesting that, initially at least, post-stratification is preferable to incorporate remote sensing into estimation with FIA field measurements.

Comparison with Other Studies
Research studies which have evaluated lidar and DAP for forest inventory have tended to focus their inferences about performance on R 2 and RMSE instead of relative efficiency. The emphasis on residual variation is likely due in part to the prevalence of linear regression as a tool for inferential modeling and mapping, and the relative obscurity of sampling methods. In an exception that discussed efficiency directly [10], lidar was used to estimate volume in addition to biomass, basal area, and number of stems, and achieve a relative efficiency of 4.5 for volume (reported in the study as design effect, which is 1/ RE), which is similar to the efficiencies observed in this study. It is fairly common for studies using lidar to exceed the R 2 values observed here. In this study we observed at best R 2 values of 0.76, where it is not unusual to see R 2 values which exceed 0.8 or even 0.9 for volume in the literature when using lidar-derived height attributes(e.g., References [37,38]), and values greater than 0.75 for frame-based DAP (e.g., References [13,14,39]).
Our findings with respect to the effects of using SPL versus HPL on efficiency are similar to what was observed in another study [40]. The related study examined the effects of up to 6m in positional error on lidar based regression predictions of aboveground biomass for multiple plot sizes. For the plot size most similar to the area of an FIA plot (707 m 2 in the related study versus 675 m 2 for an FIA plot [40] the effect of 6m positional errors on R 2 values was on the order of a 2-3 percent, or roughly the same as the difference in R 2 values that we observed between HPL and SPL coordinates with a regression estimator. The relative efficiencies observed in this study for post-stratification are comparable to those found for post-stratification of volume using auxiliary lidar heights in Reference [34]. This study differs from Reference [34] in that we used more strata (10 versus 6), and we stratified on a single auxiliary variable. Our finding that the relative efficiencies of post-stratification and model-assisted estimation agrees with another study [41], although our studies used different sources of auxiliary information (DAP versus NLCD [42]). Two other studies, References [29] and [28], disagreed with our findings and found that model-assisted estimation was more efficient than post-stratification. Key distinctions between this study and References [29] and [28] include that the method we used only measured data and post-stratified on a single auxiliary variable where they used simulated data and looked a post-stratification using combinations of multiple auxiliary variables. We are uncertain as to the root cause of our different findings, but they may be attributable to these differences in data and methodology.

Feasibility
A driving motivation to pursue alternatives to lidar to explain vertical forest structure is that lidar is very expensive. We envision DAP as a lower-cost technology which can also provide forest height. Point clouds derived from DAP are typically much less expensive than lidar, even if it is necessary to collect new imagery to support the effort. For example, new imagery and an associated point cloud may cost on the order of $0.3 USD per hectare, and less than $0.01 per hectare if the imagery has already been flown and the vendor has historical imagery in their archive. Area-wide DAP for Washington State (186k km 2 ) can be obtained for under $200k USD when derived from existing imagery. At this kind of pricing, frequent, high-precision, high-resolution forest inventory becomes more readily feasible.

Monetary Value of Efficiency
A practical implication of increased relative efficiency with DAP is that it is feasible to achieve higher estimation precision for the same number of plots. For a large area, e.g., all of Washington State, the effect may not be very impactful, but for sub-regions such as a single county, it can mean the difference between usable and unusable forest yield estimates. For sub-region estimates without auxiliary information, a small number of plots within the region may result in unreliable estimates. Integration of remote sensing information with field measured plot data for estimation is equivalent to having additional plots (for attributes related to height and cover). Here we briefly quantify the added monetary value of DAP when performing an estimation with FIA field plots using some simplifying assumptions. A major caveat to this exercise is that DAP provides little or no explanatory power for many field-measured forest attributes. For these attributes, DAP provides no added value, and it can prove difficult to untangle added value when some estimates are improved by DAP and others are not.
The use of a permanent grid of field plots is fairly standard practice and individual plot costs can be fairly high (100's to 1000's of dollars), often because of the high cost of travel between plots (and often because lots of kinds of measurements are taken on permanent plots). For this exercise, we will assume that the majority of the value of the plot is associated with forest attributes which can easily be related to remotely measured height, attributes such as height [43], basal area [44], volume [45], biomass [43], and carbon [8]. For the sake of simplicity we will assume that the hypothetical cost per plot is $1000 for some inventory (actual FIA plot costs in WA vary greatly but are expected to be ≈ $2000/plot on average), that our area of interest has 100 field measurement plots, and that the values in Table 1 are applicable. Based on our results in Table 1, our estimation strategies with auxiliary DAP ranged from 3-4 times as efficient as estimation without DAP (assuming SRS). In the case of our area with 100 field plots, this is equivalent to having 3-4 times as many plots or an additional 200-300 plots. With our assumption of $1000/plot, this is an added value of $200,000-$300,000, a substantial increase.
This same approach can be used to quantify the value associated with differences amongst the estimation strategies that we presented in Table 1. In the case of post-stratification with HPL versus SPL, there is a 0.7 unit difference in relative efficiency. For our same target area with 100 plots, this is equivalent to having 0.7 × 100 = 70 additional plots, or an added value of $70,000 at $1000/plot. If we divide this added value amongst the existing 100 plots in the county, HPL provides an added value of $70,000 ÷ 100 = $700 per plot. Based on the provided assumptions, as long as the added cost of HPL does not exceed $700/plot, the HPL program adds value to the inventory.

Scope of Inference
Estimates in this study use data that are restricted to the extent of the lidar DTM. As a result, our inferences are only about the relative performances of estimation strategies, and do not represent point estimates for all of Washington State. The tables of results provided in this study are only appropriate to guide the selection of an estimation strategy, not as stand-alone values.

Study Limitations and Future Direction
This study focuses on attributes and measurements which are associated with forest structure. Forest inventory crews typically measure a wide variety of forest attributes, many of which cannot be easily modeled from remotely-sensed forest structure measurements. This oversight (not explicitly addressing other forest attributes) is especially apparent when we monetize the added value of upgrading our inventory with remote sensing. The added value is only applicable to attributes for which remote sensing can aid in estimation, and in our toy example we assign the full value of an additional plot ($1000) to the added efficiency, whereas some of that value is in fact lost due to the fact that remote sensing was not used to aid in (e.g.) estimation of species proportions. In practice, consideration must be given to the fact that many of the attributes which are not easily captured with height-based measurements are also important for a wide variety of important inferences about our nation's forest. Most forest inventory applications are dependent on measurements such as species and damage which are not easily accounted for by the approach used here.
In some instances, estimates for difficult to measure attributes (with remote sensing) can still see some improvement when combined with height-related attributes. Biomass by species (e.g) will not have improvement in the species proportion component of the estimate, but the biomass component of the estimate is improved, which can result in a lower average error for biomass by species.
Future research efforts will include investigation of a greater variety of forest attributes, and integration with additional sources of information such as weather, elevation, and spectral information (especially Landsat time series), which may aid in the modeling and estimation of a wider variety of forest attributes.

Conclusions
This study demonstrates that pushbroom imagery (DAP) can be used to increase both the efficiency and resolution of forest yield estimates. While high-resolution area-wide height measurements have been available from lidar for some time, low-cost pushbroom DAP makes for a compelling case for wide-scale and repeated acquisitions to support forest yield mapping and estimation. The maturation of this technology provides incredible new opportunities to examine vegetative structure such as height, biomass, and volume for entire states at a level of precision and resolution that typically are only available for the extents of relatively small and infrequent lidar projects. One could imagine, for example, that state-wide fine-resolution maps of canopy fuel model parameters such as crown bulk density could prove advantageous for planning and risk assessment in the current era of increasing fire size and severity. For the few exceptions where states have complete lidar coverage, these states are ideally situated to leverage lidar DTMs with frequent lowcost statewide DAP acquisitions.

Acknowledgments:
We would like to acknowledge the extensive efforts of FIA field crews and staff for collecting, curating, and distribute the field measurement data used in this study. We would also like to thank WA DNR and other private and public organizations which shared their lidar data and made it possible to create a lidar DTM over an extensive portion of Washington State. Finally, we would like to thank 3 anonymous reviewers whose detailed comments and suggestions greatly enhanced the quality of this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Glossary of Terminology and Statistical Notation
DAP -Digital Aerial Photogrammetry, see photogrammetry DSM -Digital Surface Model, often stored as a raster DTM -Digital Terrain Model, an approximate representation of the ground surface, often stored as a raster Efficiency -The ratio of sample variation to sample size FIA -Forest Inventory and Analysis, a USFS inventory program based on a nationwide grid of plots covering the conterminous USA, parts of Alaska, and the pacific islands.
Frame camera DAP -Photogrammetric measurements made using a digital frame camera, where the frame camera uses a rectangular sensor array, and stereo or image overlap is achieve by flying overlapping flight lines HPL and SPL -High Precision Locations and Standard Precision Locations Landsat -The longest running earth observing satellite program (1973 -present) with sensor resolutions ranging from 15 -100 m depending on satellite and band Lidar -Light detection and ranging (airborne scanning), an active measurement system which uses a scanning laser and sensors in addition to navigational positioning technologies to precisely map the earth's surface SPL and HPL -Standard Precision Locations and High Precision Locations MODIS -Moderate Resolution Imaging Spectroradiometer, two moderate resolution (250m -1km, depending on band) space-borne sensors which cover the earth every 1-2 days NED -National Elevation Dataset, a national ground model for the USA derived from various sources Phodar -A colloquial term for the point cloud generated by DAP, see photogrammetry Point cloud -a series of XYZ coordinates, or the 3D representation of the physical structures, where the coordinates are commonly acquired with a remote sensing technology Photogrammetry -The processes and technologies used to make measurements from photography [46], most notably here, the ability to make horizontal and vertical measurements from stereo imagery. Photogrammetry is used in this study to generate a 3D point cloud.
Pushbroom DAP -Photogrammetric measurements made using a pushbroom camera, where the pushbroom camera uses a wide linear sensor array, and stereo or image overlap is achieved by having multiple sensor arrays simultaneously collecting data while pointing e.g. forward, down and backwards SfM -Structure from Motion (DAP), see photogrammetry Stereo -The overlap between two or more images taken from multiple camera positions, and in the context of this study is typically achieved by taking a series of images from an aerial platform such as an airplane or helicopter USDA -United States Department of Agriculture, the branch of the federal government which includes the USFS USFS -United States Forest Service, a branch of the federal government which is responsible for the FIA program in addition to numerous other responsibilities such as research and managing some federal lands USGS -United States Geological Survey, a federal agency focused on the measurement and monitoring of environmental conditions

Statistical Notation
DE -Design effect, or the ratio of variance of some design relative the variance estimator for a simple random sample RE -Relative efficiency or ratio of variances for two sampling designs -A single measurement of the response attribute, e.g. the volume measured on a single FIA plot.
̂ -A single deviation between an observation and the fitted regression line.
-A prediction from the regression line for observation i.
-The sample variance of the response attribute, in our case volume.
-The sample variance of residuals around the regression line. ̂ -Unbiased estimator for a simple random sample, which is simply the unweighted sample mean.
̂ -A regression estimator of the mean of some attribute for some domain, e.g. mean forest volume for WA state.
-Coefficient of determination, or the proportion of variation explained. ̂ -A regression estimator of the total of some attribute for some domain, e.g. total forest volume for WA state.
̂ ) -A variance estimator for ̂ , or the expected variance amongst mean estimates when using ̂ for repeated samples from the same population. ̂ -A variance estimator for ̂ , or the expected variance amongst mean estimates when using ̂ for repeated samples from the same population.