Pushbroom Photogrammetric Heights Enhance State-Level Forest Attribute Mapping with Landsat and Environmental Gradients

Abstract

We demonstrate the potential for pushbroom Digital Aerial Photogrammetry (DAP) to enhance forest modeling (and mapping) over large areas, especially when combined with multitemporal Landsat derivatives. As part of the National Agricultural Imagery Program (NAIP), high resolution (30–60 cm) photogrammetric forest structure measurements can be acquired at low cost (as low as $0.23/km² when acquired for entire states), repeatedly (2–3 years), over the entire conterminous USA. Our three objectives for this study are to: (1) characterize agreement between DAP measurements with Landsat and biophysical variables, (2) quantify the separate and combined explanatory power of the three auxiliary data sources for 19 separate forest attributes (e.g., age, biomass, trees per hectare, and down dead woody from 2015 USFS Forest Inventory and Analysis plot measurements in Washington state, USA) and (3) assess local biases in mapped predictions. DAP showed the greatest explanatory power for the widest range of forest attributes, but performance was appreciably improved with the addition of Landsat predictors. Biophysical variables contribute little explanatory power to our models with DAP or Landsat variables present. There is need for further investigation, however, as we observed spatial correlation in the coarse single-year grid (≈1 plot/25,000 ha), which suggests local biases at typical scales of mapped inferences (e.g., county, watershed or stand). DAP, in combination with Landsat, provides an unparalleled opportunity for high-to-medium resolution forest structure measurements and mapping, which makes this auxiliary data source immediately viable to enhance large-scale forest mapping projects.

1. Introduction

The complexity and scope of forest planning and monitoring across large ownerships (e.g., USDA Forest Service National Forest Systems) is motivated by various management issues, such as climate change, wildfire risk, wood production, and threatened and endangered species, and is made more robust by regionally consistent, multifaceted forest attribute datasets. Spatially explicit information for detailed forest attributes (e.g., structure and composition) can support landscape analyses across a variety of ownerships and biophysical conditions. Thus, wall-to-wall mapping based on time-series of Landsat imagery, environmental gradients, and forest inventory data have become increasingly common as a basis for planning and monitoring [1,2,3,4,5]. While such modeling and mapping efforts provide valuable information from landscape to regional levels, they may not be sufficiently accurate at finer scales (e.g., stands or blocks) needed for tactical management decisions, as indicated, e.g., by coefficients of determination increasing from 0.37 at plot-scales (approximately 1 ha) to >0.90 at landscape-scales (216,500 ha hexagons) in Riemann et al. [6].

Remote sensing of vertical forest structure can improve estimation and mapping in support of planning and monitoring [7]. For example, the addition of detailed forest height information from airborne lidar (henceforth just “lidar”) has been shown to dramatically improve accuracy for forest structural, but not compositional, attribute modeling compared to Landsat-based approaches [8]. Underestimation of biomass and related forest attributes in high biomass ecosystems associated with spectral saturation is common for models based on multispectral imagery [9,10]. This issue is often minimized when using lidar height-based measurements as predictor variables in modeling [11]. However, lidar coverages are not spatially complete, nor do they have consistent repeat acquisitions. As a result, available lidar data do not have the spatial and temporal coverage to fully support regional estimation and monitoring efforts.

The development of canopy surface height measurements from National Agricultural Imagery Program (NAIP) pushbroom-sensor image acquisitions, known as Digital Aerial Photogrammetry point clouds (or simply DAP, for short), may address some of the shortcomings of lidar acquisitions related to cost and frequency in the continental United States (CONUS), albeit with fine-scale artifacts (e.g., missing trees in openings and edge effects) and less precision than either lidar or frame camera DAP [12,13]. Statewide acquisitions occurring at 2–3-year intervals—the current cadence of NAIP imagery acquisition—could form the basis for broad-scale monitoring if DAP height data provide adequate information about forest canopy structure. Integration of lidar structure with other auxiliary data sources has been investigated in the peer-reviewed literature [14,15,16,17,18], but it is not known how the incorporation of DAP (which has less explanatory power than lidar) with Landsat and environmental gradients affects forest attribute modeling and whether the fusion of multiple datasets synergistically improves model performance. Examples of applying DAP based on pushbroom imagery, as is used for the NAIP program, to model forest attributes are relatively uncommon and generally focus on small study areas <20,000 km² [19]. Therefore, even though DAP is far less expensive than lidar ($0.23/per km² vs. $200 per km²) while providing some of the same information, understanding the capacity for DAP to improve mapping and estimation of forest attributes will be required before state or federal agencies will consistently support acquisition of such data.

The objective of this study is to assess whether the vertical structure information provided by NAIP DAP improves predictive power for forest attributes over what is provided by Landsat and environmental gradients. The explanatory power of NAIP DAP is evaluated when modeling a suite of 19 field-derived forest attributes over a large area (Washington State, USA; Figure 1) for a single year, 2015. We also assess whether there is increased explanatory power provided when multiple auxiliary data sources are used together, by how much, and for which forest attributes. Our hypothesis is that NAIP DAP can appreciably improve the explanatory power of forest attribute models relative to using only Landsat and environmental gradients, especially forest attributes closely related to height and cover such as height, volume, biomass, and carbon. These models can then be used for precise large-area forest mapping.

2. Methods

2.1. Study Area

The study was conducted for forest areas in the state of Washington (WA) in the USA (Figure 1). The study area encompasses 184,000 km² and includes diverse biophysical conditions, such as marine forests, montane forests of up to 4392 m elevation (Mount Rainier), rainforests with more than 2.5 m of rainfall, and arid desert conditions with less than 10 cm of rainfall per year. The state is 53 percent composed of forests (90,000 km²), of which 89% are dominated by conifer trees [20].

The most common tree species in WA is Douglas-fir (Pseudotsuga menziessi (Mirb.) Franco). Other common conifers include western hemlock (Tsuga heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata Donn ex D. Don), and various species of true fir (Abies), pine (Pinus), and spruce (Picea). Common broadleaf species include western red alder (Alnus rubra Bong.), bigleaf maple (Acer macrophyllum Pursh), black cottonwood (Populus balsamifera Torr. & Gray), oak (Quercus garryana Douglas ex Hook.), and aspen (Populus tremuloides Michx.).

2.2. Field Sample Data

Field measurements were collected on USDA Forest Service Forest Inventory and Analysis (FIA) plots, the national grid of Continuous Forest Inventory (CFI) plots for the USA [21]. FIA field plots use a clustered 0.067-hectare plot design composed of 4 0.017-hectare circular subplots over an approximately 1.0-hectare area (Figure 2).

We summarized 573 FIA field plots collected in 2015 across Washington. To be included in modeling, plots had to have at least half of their sampled area in forested conditions as determined by field crews. On each of four subplots, crews tallied individual live trees ≥2.5 cm and standing dead (snag) trees ≥12.7 cm using three different fixed-radius plot sizes: microplot ≈1/750 ha for trees 2.5–12.5 cm diameter at breast height (DBH), subplot ≈1/60 ha for trees and snags ≥12.7 cm DBH, and macroplot ≈1/10 ha for trees and snags ≥61.0 cm DBH (western Washington) and ≥76.2 cm DBH (eastern Washington). From individual tree records and their associated expansion factors, we calculated the suite of plot attributes in Table 1 which are summarized in Table 2. We adjusted tree expansion factors to include only sampled areas that were in forested condition classes, which includes recently disturbed (e.g., harvest or fire) areas.

Locations of FIA subplots were obtained as part of the high-precision positioning program used area-wide (CA, OR, WA, AK, HI, pacific islands) by the Pacific Northwest FIA program PNW FIA [25]. PNW FIA has invested in dual frequency (L1/L5) GPS and GLONASS GNSS receivers with positioning accuracy estimated to be better than 2 m on average under dense canopy [26,27,28]. This level of positioning accuracy may not be needed for intersection with Landsat and environmental gradients, but Strunk et al. [13] showed that high performance positioning provides measurable improvements in models fit with auxiliary NAIP DAP variables.

2.3. Auxiliary Data Sources

Three sources of auxiliary data were used in this investigation, including environmental gradients (Table 3), Landsat derivatives (Table 4), and NAIP Digital Aerial Photogrammetry (DAP; Table 5). Climate attributes were derived from the Parameter-elevation Relationships on Independent Slopes Model (PRISM) using 30-year monthly normals (1981–2010) (http://www.prism.oregonstate.edu/normals/ accessed on 28 September 2016). Ecologically significant climate indices were first calculated at 800 m from monthly temperature, precipitation, and vapor pressure deficit rasters, then resampled to 30 m using bilinear interpolation. Location rasters were generated using pixel center’s geographic coordinates. Location rasters include a coastal proximity index for air temperature that measures the optimal path of air parcels to each pixel [29]. Continuous soil characteristics were downloaded from Soil Information for Environmental Modeling and Ecosystem Management using USDA-NRCS STATSGO soil data [30]. We combined separate soil horizons into composite measures using depth-weighted averaging.

We derived topographic covariates by first mosaicking 10 m (1/3 arc-second) digital elevation models (DEMs) downloaded from the National Elevation Dataset [34], then resampling to 30 m using bilinear transformation. All elevation derivatives were created from the 30 m DEM to avoid inconsistencies stemming from calculating indices from 10 m and then resampling to other distances.

For Landsat derivatives, we used variables derived from a multispectral ensemble modeling approach to the LandTrendr algorithm [35,36]. Methods are described in detail in Bell et al. [37], but we briefly describe the key processing steps. Yearly summer medoid composites of Landsat 5 (TM), 7 (ETM+), and 8 (OLI) were created for 1985–2017 and run through LandTrendr using a variety of spectral bands and indices to delineate segment endpoints. The output of each LandTrendr run is a set of temporal vertices and segments for each pixel that describe change over all years. The sets of LandTrendr outputs then are used as inputs in a random forest model [38] using TimeSync pixel trajectories as training data [39]. As with the inputs, the output of the ensemble modeling is a revised set of vertices and segments that are used to create fitted spectral bands and indices. We applied fitted tasseled-cap values estimated by Crist and Cicone [40] and computed normalized burn ratio [41] indices in this analysis (Table 4).

Stereo imagery used to generate photogrammetric point clouds was collected in the summer of 2015 by Northwest Geomatics Ltd. (Calgary, Alberta, Canada) using a Leica ADS100 sensor. Rear and nadir view angles were used for point cloud creation (forward angle not used). Data were collected from a nominal 5000 m flying height for a target 40 cm GSD. Photogrammetric point clouds were generated from NAIP pushbroom imagery for the state of WA by WA Department of Natural Resources using BAE Socet software [42] at a nominal 40 cm resolution. The point clouds are stored in the laz format [43] where vertical information are elevations. While lidar is typically also processed to create a Digital Terrain Model (DTM), NAIP DAP point clouds were not processed to create DTMS since they generally only represent the canopy surface structure. The DTM used in this study was instead taken from the National Elevation Data (NED) [44]. This causes a reduction in accuracy over what can be achieved with a lidar DTM e.g., [13].

While the DAP point clouds generated from the pushbroom imagery collected for the NAIP program offer many opportunities due to their cost–value relationship (<$0.23/km²) and the relatively high temporal frequency of NAIP (2–3 years), and the data bear a superficial resemblance to lidar, pushbroom DAP point clouds should not be considered equal to airborne lidar with respect to their ability to characterize forest structure. Pushbroom DAP point clouds have many limitations in terms of structure and artifacts relative to lidar data. Pushbroom DAP resembles a coarse, smoothed lidar-derived canopy surface model (elevations) with artifacts resulting from shadows, edges, and single trees in openings (e.g., Figure 3). DAP point clouds can be an excellent source of general canopy structure information, but they do not have the detail, reliability, intra-canopy structure, and (typically) ground returns below vegetation that are provided by lidar [12,13]. Pushbroom DAP point clouds appear to be especially poorly suited for individual tree detection (ITD) for the example datasets that we have worked with thus far: CT, NH, OR, TN, WV, WA. Hopefully, with camera and software updates this will change. Frame-camera-derived DAP point clouds, in contrast, contain more fine canopy details that may better support ITD (not shown). Despite the limitations relative to lidar, the DAP point clouds derived from NAIP stereo provide a unique opportunity to measure forest structure over the entire conterminous USA at regular intervals (2–3 years) at fine spatial resolutions (30–60 cm).

The DAP point cloud was processed to create plot metrics using the U.S. Forest Service FUSION software [45]. FUSION was used to extract DAP points for FIA plot footprints and then FUSION was used to calculate a suite of DAP height-derived metrics for individual plot locations including mean, standard deviation, quantiles, and cover (Table 5). See the FUSION documentation for a complete list of metrics and tools available in the FUSION suite.

2.4. Orthogonal Transformations

One of the objectives of the study is to assess similarities and differences between DAP height, environmental gradients, and Landsat imagery in their ability to explain variation in forest attributes. However, each of the auxiliary data sources is itself a multivariate dataset with varying amounts of explanatory information. To assess the agreement between auxiliary data sources using, e.g., scatterplots, would require hundreds of scatterplots, and even that could be misleading, as some relationships are multivariate in nature. Instead, to simplify comparisons, each of the auxiliary data sources and the response variables were separately orthogonally transformed. In the case of DAP height, environmental gradients, and the response variables, principal component transformations were used. In the case of Landsat spectral data, a tasseled cap transformation was used, which is very similar to principal components [46]. These transformations consolidate unique linear information into the first few columns of data. Visual assessment of the presence of linear trends observed between datasets can then be inferred from plotting the first few bands against each other from each transformed dataset. In the case of auxiliary data sources, the strength of the trends represent the amount of redundant information, and for the response, the strength of the trends represent the amount of explanatory power. We recognize that using ordination to generate orthogonal predictor variables that explain a large portion of the variance in each auxiliary data source does not necessarily provide the optimal set of predictors for any specific model. Reducing the complexity of these datasets to just a few dimensions is a convenient method to assess the potential contribution of DAP, Landsat, and/or environment to the modeling of forest attributes possible.

2.5. Model Exploration

We used Ordinary Least Squares (OLS) to explore the linear agreement between auxiliary variables and response variables. Linear methods were convenient in the context of this study as they facilitate interpretation of the linear agreement between multiple variables. Thus, we could explore the relative explanatory power of each of the three auxiliary data sources. Linear models were fit using the R statistical programming environment [47] with the RStudio GUI [48]. Linear regression was used to relate individual response variables with up to three predictor variables from each of the predictor sets. A best subsets regression approach [49] was used to select up to three predictor variables that achieve the highest value of Bayesian Information Criterion BIC [50], a common model scoring criterion.

Relative model performances and the explanatory information in the three auxiliary data sources were evaluated using the coefficient of determination (R²). While the same data used for fitting were used for evaluation, in preliminary exploration we observed that OLS with large numbers of training observations (e.g., 500+) and few numbers of covariates is fairly robust to overfit and optimistic R² values. Models with higher R² values are assumed to have higher performances, and the contributing auxiliary data sources are assumed to provide more explanatory information than alternative groups of auxiliary data sources. This approach is clearly limited to monotonic linear relationships, where more complex and non-monotonic relationships are not considered. This type of investigation is exploratory in nature and provides indication of the relative contributions of the three auxiliary data sources, but this study is not designed to identify “best” prediction strategies.

Spatial correlation amongst residuals was evaluated both visually by using mapped residuals (clusters of similar residuals) and quantitatively by using the Moran’s I statistic [51], a measure of global spatial autocorrelation. We used Moran’s I with resampling simulations to test the probability that the observed level of spatial autocorrelation could arise through random chance—probability is computed as the proportion of resamples (taken without consideration for spatial arrangement) with the same or larger absolute values of Moran’s I (for a two-sided test).

3. Results

3.1. Linear Agreement between Variable Types

To assess the linear agreement among different sets of predictor variables, we examined scatterplots between DAP HT, ENV (first principal components), and LS (Tasseled Cap 1) (Figure 4). The principal components were fit using all the variables from each type, although the results were insensitive to subsets of variables as the many variables contain redundant information (not shown). In the first panel of Figure 3, we see little or no evidence of an association between DAP HT and ENV—only 8% of variation in ENV.PC1 is explained by HT.PC1 (R² = 0.08). In the second panel, there is evidence of a weak association between LS and DAP HT (R² = 0.26). In the third and final panel of Figure 4, there is no evidence of an association between ENV and LS (R² = 0.01). These results suggest there is little or no linear agreement between the three types of auxiliary data sources, and they may contribute additively to models for the response variables. Scatterplots for additional principal components (e.g., HT.PC2, not shown) demonstrated no visible trends for any combinations of auxiliary variables (data not shown).

We assess the relationship between the three auxiliary data sources and the response variables using the scatterplots in Figure 5 which indicate that the first principal component from forest attributes were weakly to strongly related to each of the auxiliary datasets following orthogonal transformations. The scatterplots suggest that DAP height predictors have the greatest amount of explanatory power for forest attributes (R² = 0.54), followed by Landsat (R² = 0.29), with environment variables having the least apparent explanatory power (R² = 0.12). The first principal component for DAP height (HT.PC1) has a tight grouping with the response vector (RESP.PC1). Landsat (LS.TC1), in contrast, has a curvilinear relationship with response PC1, which appears to be grouped less tightly than the response PC1 vs. DAP height metrics. The strength of the relationship for Landsat as indicated by R² for the simple linear regression fit is negatively affected by the curvilinear relationship. The scatter for ENV.PC1 has a weak linear relationship with the response data. Other principal component auxiliary vectors (e.g., PC2, PC3, TC2, TC3, etc.) demonstrated no linear relationships with response principal components.

3.2. Linear Modeling Performance

The explanatory power of each of the three auxiliary data sources is inferred from the coefficient of determination values (R²) reported for the 19 forest attributes in Figure 6. As was suggested previously in scatterplots of principle components, the OLS results indicate that DAP height derivatives have the greatest explanatory power (i.e., greatest R² for a model with only one auxiliary data source) for the greatest number of variables, followed by Landsat, with the least explanatory power coming from environmental variables. There are response variables for which Landsat exceeds the explanatory power of DAP height, especially SDI and CANCOV. We should note that CANCOV is not a field-measured attribute, it is an allometric prediction from field-measured attributes [52] and the performance may not be representative of the predictive performances for field-measure canopy cover or canopy closure. There were no response variables for which environmental variables alone exceeded the explanatory power of either DAP height or Landsat alone.

The performances for models fit with variables from multiple types of predictors suggest that the components of variation explained by each type of auxiliary variable were at least partially independent, as indicated by their additive contributions to performances when used together. Performance inproves when more than one type of predictor is available, and the best performance is generally achieved using all three predictor types, although the improvement is generally small for environmental gradients. Response variables associated with canopy trees (e.g., BA_GE_3, CANCOV, BPH_GE_3_CRM, DDI, STNDHGHT, VPHC_GE_3) generally exhibit the greatest R² when either a DAP or Landsat remote sensing auxiliary data source is included. Response variables associated with species composition (BAH_PROP), standing dead trees (SVPH_GE_25), and tree densities (TPH_GE_3, TPH_GE_75) had lower R² values across all models.

3.3. Spatial Analysis

Given the large extent of our project area (184,000 km²), a linear model with a fixed set of covariates for the entire area may perform poorly for local predictions. The primary concern is that spatially grouped observations might have correlated residuals. This would result in bias in the aggregation of pixels to polygons such as stands, watersheds or counties. To assess the need for model localization, we evaluate spatial patterns in the residuals for volume (VPHC_GE_3; Figure 7)—a variable that is commonly of interest for land managers in the Pacific Northwest (USA). The spatial patterns evident in the residuals in Figure 7 suggest that a generic model fit for the entire state of WA is not appropriate for mapping purposes and that a localization process of some kind is needed. For example, it is evident that there are clusters of blue, orange, and red residuals at various locations throughout the state. A formal test of spatial correlation between the residuals using Moran’s I indicates that it would be highly unlikely to see the observed residual pattern from spatially independent observations (p-value << 0.001).

4. Discussion

Our investigation demonstrates that, while DAP provides the greatest explanatory power overall for forest attributes (but not for all attributes) of the three auxiliary data sources, there is room for dramatic improvement when DAP and Landsat are used in combination (up to 20 percentage point increase in R² and a mean increase compared to DAP alone of 9.5 percentage points). Environmental gradients, in contrast, provide little or no explanatory over what is provided by DAP or Landsat. Pairwise scatterplots of their principal components suggest that some of the information in the three auxiliary data sources is redundant, but Landsat and DAP especially provide independent explanatory information, as indicated by their weak linear agreement with each other (DAP PC1 and Landsat TC1) and stronger agreement with PC1 from the response data. These findings are encouraging given that Landsat imagery is freely available at high temporal resolution (>1 overflight/month) and NAIP imagery is collected nationwide (conterminous USA) every 2–3 years. These data can readily be combined in large area mapping efforts for considerably less cost than a similar mapping effort with lidar—although with reduced precision relative to lidar.

The combination of NAIP DAP and Landsat provide an exciting opportunity for high fidelity forest mapping and estimation, but there is evidence that careful consideration should be given to spatial correlation between prediction errors. These data can nominally be used in combination for fine-scale predictions (e.g., 30 m or finer), but our exploratory findings suggest that due to the spatial correlation between residuals, aggregates to, e.g., stand or watershed polygons could be highly biased. Observed spatial correlations are consistent with the use of longitude and latitude as predictors as well as modeling and mapping at sub-regional scales (e.g., partitioning the data into 3–4 ecoregions for the state of Washington, USA) in regional forest attribute mapping and monitoring frameworks utilized in OR and WA for the past two decades [1,37,53]. Our results are also consistent with evidence that biomass modeling and mapping based on spectral (Landsat) and height (lidar) that does not account for spatial error can lead to underprediction of estimate error in small-area applications [54]. Additional investigation would be needed to determine the degree to which spatial modeling could reduce local bias and the amount of bias present for various scales of aggregation. Alternative modeling approaches such as errors in variables models [55], spatial models [56], and non-parametric approaches [38,57] should be investigated for predictive performance with respect to small-area inferences.

As has been observed in other studies [8,17], we observed that height-derived attributes contribute the greatest explanatory power in models of forest attributes, followed by spectral data and last, environmental gradients. Forest attribute mapping for the Deschutes National Forest in central Oregon indicated that the addition of lidar to k nearest neighbor (kNN) imputation models using Landsat and environmental predictors increased R² for many structural variables, including QMD (0.1 to 0.55), HCB (0.1 to 0.5), and BA (0.3 to 0.7) [8]. Those improvements are more dramatic than the improvements in QMD (0.17 to 0.29), HCB (0.28 to 0.41), and (0.55 to 0.63) in our study. Differences may be attributed to their use of more accurate height information (lidar with a lidar DTM vs. pushbroom DAP with the USGS NED DTM), their smaller study area (a single National Forest vs. all of Washington state), and possibly due to their use of a more flexible modeling method (kNN vs. OLS). The importance of inclusion of additional auxiliary data sources appears to be greater for pushbroom-DAP-derived height metrics than those derived from lidar, as pushbroom DAP generally explains less variation in forest attributes than would be expected with lidar i.e., there is more variance to explain, e.g., [12].

The R² values observed for modeling forest attributes from pushbroom DAP in this study are lower than has been observed for other studies with pushbroom DAP in the NW USA. The study by Strunk et al. [12] used circular plot footprints over a much smaller 500 km² extent, which may be indicative that localization and a consolidated plot footprint could result in better modeling performances. Another study that had better model performances than we observed also used FIA plots; however, the study used plots that were predominantly located in lowlands in western WA forests and only for areas with a lidar DTM, supportive of the hypotheses that a better DTM and localization may both play a role in predictive performance [13].

The performances of forest attribute modeling studies with Landsat vary dramatically, depending upon modeling approach, forest type, plot type, and data processing steps. For example, in a single study comparing kNN, reduced major axis regression, and random forest for Arizona and Minnesota for the prediction of biomass on FIA plots from Landsat, they observed R² values (converted from variance ratio or VR) ranging from 0.25 (VR = 1.25) to 0.7 (VR = 0.58) [58]. Similar levels of variation were observed in forest modeling results for the Pacific Northwest (USA) as part of a USFS led Landsat and biophysical variable kNN mapping product [1,59,60]. For example, biomass models in OR and WA by the LEMMA Team achieved R² values in the range of 0.41–0.54, which is like the value observed here (0.42). However, it is difficult to provide direct comparison between the Landsat and biophysical variable results observed here and other studies, given the variety of study differences and our emphasis on exploring linear variable agreement instead of optimal modeling approaches. One consistent theme in the prediction of forest attributes from Landsat is that there is saturation at higher levels of biomass, which is less evident from remote sensing technologies that measure (or infer) height and cover such as lidar and DAP (e.g., Figure 4).

One limitation of this study is that the pushbroom DAP data were processed shortly after the 2015 collection, and our familiarity with processing ADS100 pushbroom stereo with BAE Socet software has improved between 2015 and today (2022). There are numerous errors and artifacts in the 2015 dataset that are not present in later pushbroom DAP point cloud datasets for WA (e.g., 2017 and 2019 data). In addition, due to the much greater lidar extent today than in 2015 (lidar DTMs yield improved tree heights relative to heights compute from NED DTMs) and improvements in photogrammetric point cloud extraction, it is likely that new pushbroom DAP datasets (e.g., 2021) will perform better than those reported here, which will also affect the relative contributions of Landsat and environmental gradients. The results are also specific to the FIA clustered subplot footprint, where the performances and relative contributions of the three auxiliary data sources to modeling and prediction may differ for single-square or circle-plot footprints (rather than clustered, as with the FIA plot). The results found here are also not conclusive of the potential for DAP in general, as frame camera DAP has been observed to perform much better than Pushbroom DAP and even similar to lidar [12,61].

Another limitation of this study is that results for cover provided here are not indicative of performances that may be achieved with other approaches to quantify tree canopy cover, such as with a line transect, spherical densiometer, moosehorn, or hemispherical photos e.g., [62]. For this analysis we implemented the allometric cover prediction approach described by Hann [23], which may be more or less correlated with remote sensing than other methods. For example, in Strunk et al. [12], DAP explained greater than 90% of variation in lidar-derived cover, and in unpublished analyses associated with the same investigation, DAP explained 76% or 86% of variation in field-measured cover when measured with a moosehorn or spherical densiometer, respectively. However, the cause of lower cover performance observed here relative to that study is not clear; the lower DAP performance achieved for cover in this study may be due to the use of allometric cover, our use of a different plot footprint composed of widely spaced (≈30 m) subplots or the much larger geographic extent of this study relative to Strunk et al. [12].

Lastly, although we found comparatively less contribution to models from environmental gradients, inferences from our results are limited to structural forest attributes and linear relationships in a global (statewide) modeling approach. Environmental gradients may aid in the prediction and localization of models using other modeling approaches, responses or response types [1,63]

5. Conclusions

The pushbroom imagery collected in support of NAIP, the USA’s national imagery program, enables frequent (2–3 years) high-resolution (30–60 cm) photogrammetrically derived forest structure measurement over broad areas, including the entire conterminous USA. This study demonstrates that DAP is a powerful auxiliary variable in the prediction of many forest attributes, especially in combination with Landsat-derived auxiliary variables. While the explanatory power afforded by pushbroom DAP is substantially lower than lidar or frame camera DAP, modeling with both DAP and Landsat in combination resulted in up to 20 percentage point improvements (R²) over either auxiliary dataset alone. Environmental gradients, in contrast, provided limited additional explanatory contribution in our exploratory investigation—more robust analyses involving small area and spatial modeling methods may yield different results with respect to the importance of environmental gradients.

Since Landsat and DAP are readily available over large areas, their combined integration in mapping efforts could be used to facilitate creation of frequent, precise, spatially consistent forest maps for, e.g., entire states, but informative at fine scales such as counties, watersheds or stands. The resulting maps are not free of local biases at these scales, however, as there was evidence of spatial correlation in mapped residuals. Future research will look at leveraging more sophisticated modeling strategies which can accommodate the complex nonlinear and spatial interdependence among forest attributes and covariates, and leverage other sources of auxiliary data such as satellite DAP and Sentinel imagery.

Author Contributions

Conceptualization, J.L.S.; methodology, J.L.S.; validation, J.L.S.; formal analysis, J.L.S.; investigation, J.L.S.; resources, J.L.S.; data curation, J.L.S. and M.J.G.; writing—original draft preparation, J.L.S., D.M.B. and M.J.G.; writing—review and editing, J.L.S., D.M.B. and M.J.G.; visualization, J.L.S.; supervision, J.L.S.; project administration, J.L.S.; funding acquisition, J.L.S., D.M.B., and M.J.G. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

We would like to express our appreciation to Caleb Maki at the WA DNR for processing NAIP stereo to generate the DAP point clouds used in this study—this study would not have been possible without him and the encouragement from WA DNR to use these data. We would also like to thank Yang Zhiqiang for contributing auxiliary LCMS variables, and Hans-Erik Andersen and Derek Churchill for their suggestions on drafts. We thank the USDA Forest Service Northwest Region and the Pacific Northwest Research Station for funding that supported this work. We also thank the USDA Forest Service Forest Inventory and Analysis program for collecting the field observations used in this study.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The location of the state of Washington, our study area, in the northwestern conterminous USA.

Figure 2. Schematic of FIA field plot layout [21].

Figure 3. Comparison of pushbroom NAIP DAP point cloud (red) and lidar point cloud (blue) heights above lidar DTM for a 1 m wide 170 m long strip, where the tallest tree is approximately 30 m according to the lidar.

Figure 4. Scatterplots and coefficients of determination (R2) of first principal components for Landsat (LS), DAP Height (HT), and environmental gradients (ENV). The black is an OLS fit between variables on x and y axes.

Figure 5. Scatterplots and coefficients of determination (R2) of the response first principal component (RESP.PC1) vs. first principle components for Landsat (LS), DAP Height metrics (HT), and environmental gradients (ENV). The black line is an OLS fit between the RESP.PC1 and respective auxilary variables on the x axes.

Figure 6. Coefficients of determination (R2) in percent (%) for linear models fit relating DAP heights (Ht), environmental gradients (Env), and landsat (LS) to field-measurement-derived forest attributes. The color gradient represents low (blue) to high (red) R2 values.

Figure 7. Volume (VPHC_GE_3) residual pattern (m3/ha in percent of the mean) for OLS model fit with best performing set of height, environment, and Landsat predictors (p-value < 0.0001 for Moran’s I test).

Table 1

Inventory plot attributes for accuracy assessment. Plot attributes were compiled from tree-level records and summarized to plot level (DBH = diameter at breast height).

Attribute	Units	Description
AGE_DOM	years	Basal area weighted stand age based on field recorded or modeled ages of open grown, dominant, and co-dominant trees (FIA crown classes 1–3)
BA_GE_3	m2/ha	Basal area of live trees ≥2.5 cm DBH
BAH_PROP	proportion	Proportion of total basal area that is hardwood
BPH_GE_3_CRM	kg/ha	Component Ratio Method biomass of all live trees ≥2.5 cm DBH [22]
CANCOV	percent	Canopy cover of all live trees [23]
DDI	none	Diameter diversity index [24]
HCB	m	Average height to crown base
QMD_GE_3	cm	Quadratic mean diameter of trees ≥2.5 cm DBH
SC	none	A composite measure of stand size class and stand cover class.
SDBA	m2/ha	Standard deviation of basal area of all live trees
SDDBH	cm	Standard deviation of diameter of all live trees
SDDBHC	cm	Standard deviation of diameter of all live conifers
SDI	none	Stand density index, defined as sqrt (trees per hectare × basal area per hectare)
SIZECL	none	Ordinal size class based on quadratic mean diameter of dominant and co-dominant trees (QMD).
STNDHGT	m	Average tree height of dominant and codominant trees
SVPH_GE_25	m3/ha	Volume of snags ≥25.0 cm DBH and ≥ 2.0 m tall
TPH_GE_3	trees/ha	Density of live trees ≥ 2.5 cm DBH
TPH_GE_75	trees/ha	Density of live trees ≥ 75.0 cm DBH
VPHC_GE_3	m3/ha	Volume of live conifers ≥ 2.5 cm DBH

Table 2

Summary of field measurement derived forest attributes.

	Min	Max	Mean	Median	sd	cv%
VPHC_GE_3	0	1747.45	348.87	232.82	351.82	100.84
TPH_GE_75	0	113.57	8.13	0	17.37	213.53
TPH_GE_3	0	7466.23	1009.72	700.67	1052.42	104.23
SVPH_GE_25	0	1119.22	85.02	26.69	143.92	169.29
STNDHGT	0	64.68	20.12	19.54	11.5	57.17
SIZECL	1	6	3.23	3	1.23	38.13
SDI	0	693.67	169.37	152.06	120.9	71.38
SDDBHC	0	44.24	12.28	11.41	7.65	62.3
SDDBH	0	42.25	12.23	11.45	7.39	60.44
SDBA	0	0.61	0.07	0.05	0.08	103.33
SC	1	7	4.15	4	1.71	41.26
QMD_GE_3	0	100.57	24.31	22.86	13.34	54.88
HCB	0	30.42	7.1	5.89	5.62	79.25
DDI	0	10	3.99	3.58	2.25	56.45
CANCOV	0	99.57	63.67	70.79	27.95	43.9
BPH_GE_3_CRM	0	996,580.22	199,977.61	141,155.54	198,912.17	99.47
BAH_PROP	0	1	0.05	0	0.16	313.39
BA_GE_3	0	128.54	35.67	30.15	26.41	74.05
AGE_DOM	0	526	111.36	94	82.03	73.67

Table 3

Environmental gradients used in modeling.

Attribute	Units	Description
Climate (based on PRISM 30-year normals—1981–2010)
ANNPRE	ln mm	Total annual precipitation
ANNTMP	°C	Mean annual temperature
AUGMAXT	°C	Mean August maximum temperature
CONTPRE	%	Percentage of annual precipitation falling in June–August
CVPRE	none	Coefficient of variation of mean monthly precipitation of December and July
DECMINT	°C	Mean December minimum temperature
DIFTMP	°C	Difference between AUGMAXT and DECMINT
SMRMAXVPD	hPa	Maximum summer vapor pressure deficit
SMRMNVPD	hPa	Mean summer vapor pressure deficit
SMRPRE	ln mm	Mean precipitation from May–September
SMRTMP	°C	Mean temperature from May–September
Location
LAT	none	Geographic latitude
LON	none	Geographic longitude
COASTPROX	none	Coastal proximity for temperature
Soils
AWCL1	none	Available water capacity up to one meter
ROCKDEPTH	cm	Soil rock depth
BD	none	Soil bulk density
SAND	%	Soil percent sand
SILT	%	Soil percent silt
CLAY	%	Soil percent clay
PERM	none	Soil permeability
PH	none	Soil pH
POROS	none	Soil porosity
RVOL	none	Soil rock volume
Topography
DEM	meters	Elevation
HILL	none	Elevation hillshade from azimuth = 315, altitude = 45
MLI	none	McComb’s Landform Index [31]
SLPPCT	%	Slope (percent) (rounded to nearest integer)
PRR	none	Potential relative radiation [32]
TPI150	none	Topographic position index, calculated as difference between cell’s elevation and mean elevation of cells within a 150-m-radius window [33]
TPI300	none	Topographic position index, calculated as difference between cell’s elevation and mean elevation of cells within a 150–300-m-radius annulus
TPI450	none	Topographic position index, calculated as difference between cell’s elevation and mean elevation of cells within a 300–450-m-radius annulus

Table 4

Landsat spectral derivatives using ensemble LandTrendr.

Attribute	Units	Description
TC1	none	Fitted tasseled-cap axis 1 (brightness) based on ensemble Landtrendr segmentation
TC2	none	Fitted tasseled-cap axis 2 (greenness) based on ensemble Landtrendr segmentation
TC3	none	Fitted tasseled-cap axis 3 (wetness) based on ensemble Landtrendr segmentation
NBR	none	Fitted normalized burn ratio based on ensemble Landtrendr segmentation

Table 5

Example DAP height metrics computed using FUSION.

Attribute	Units	Description
Min, max, median, mean, sd, skewness, kurtosis	meters	Various statistics computed on vertical distribution of point heights above 2 m. Points below 2 m were not considered.
P05, P50, P90	meters	Height quantiles computed on the vertical distribution of point heights above 2 m. Points below 2 m were not considered.
Cover	%	Percent of points above 2 m relative to all points

References

1. Ohmann, J.L.; Gregory, M.J. Predictive Mapping of Forest Composition and Structure with Direct Gradient Analysis and Nearest- Neighbor Imputation in Coastal Oregon, U.S.A. Can. J. For. Res.; 2002; 32, pp. 725-741. [DOI: https://dx.doi.org/10.1139/x02-011]

2. Tomppo, E.; Olsson, H.; Ståhl, G.; Nilsson, M.; Hagner, O.; Katila, M. Combining National Forest Inventory Field Plots and Remote Sensing Data for Forest Databases. Remote Sens. Environ.; 2008; 112, pp. 1982-1999. [DOI: https://dx.doi.org/10.1016/j.rse.2007.03.032]

3. Wilson, B.T.; Lister, A.J.; Riemann, R.I. A Nearest-Neighbor Imputation Approach to Mapping Tree Species over Large Areas Using Forest Inventory Plots and Moderate Resolution Raster Data. For. Ecol. Manag.; 2012; 271, pp. 182-198. [DOI: https://dx.doi.org/10.1016/j.foreco.2012.02.002]

4. White, J.C.; Coops, N.C.; Wulder, M.A.; Vastaranta, M.; Hilker, T.; Tompalski, P. Remote Sensing Technologies for Enhancing Forest Inventories: A Review. Can. J. Remote Sens.; 2016; 42, pp. 619-641. [DOI: https://dx.doi.org/10.1080/07038992.2016.1207484]

5. Lister, A.J.; Andersen, H.; Frescino, T.; Gatziolis, D.; Healey, S.; Heath, L.S.; Liknes, G.C.; McRoberts, R.; Moisen, G.G.; Nelson, M. et al. Use of Remote Sensing Data to Improve the Efficiency of National Forest Inventories: A Case Study from the United States National Forest Inventory. Forests; 2020; 11, 1364. [DOI: https://dx.doi.org/10.3390/f11121364]

6. Riemann, R.; Wilson, B.T.; Lister, A.; Parks, S. An Effective Assessment Protocol for Continuous Geospatial Datasets of Forest Characteristics Using USFS Forest Inventory and Analysis (FIA) Data. Remote Sens. Environ.; 2010; 114, pp. 2337-2352. [DOI: https://dx.doi.org/10.1016/j.rse.2010.05.010]

7. Huang, W.; Swatantran, A.; Johnson, K.; Duncanson, L.; Tang, H.; Dunne, J.O.; Hurtt, G.; Dubayah, R. Local Discrepancies in Continental Scale Biomass Maps: A Case Study over Forested and Non-Forested Landscapes in Maryland, USA. Carbon Balance Manag.; 2015; 10, 19. [DOI: https://dx.doi.org/10.1186/s13021-015-0030-9]

8. Zald, H.S.J.; Ohmann, J.L.; Roberts, H.M.; Gregory, M.J.; Henderson, E.B.; McGaughey, R.J.; Braaten, J. Influence of Lidar, Landsat Imagery, Disturbance History, Plot Location Accuracy, and Plot Size on Accuracy of Imputation Maps of Forest Composition and Structure. Remote Sens. Environ.; 2014; 143, pp. 26-38. [DOI: https://dx.doi.org/10.1016/j.rse.2013.12.013]

9. Steininger, M.K. Satellite Estimation of Tropical Secondary Forest Above-Ground Biomass: Data from Brazil and Bolivia. Int. J. Remote Sens.; 2000; 21, pp. 1139-1157. [DOI: https://dx.doi.org/10.1080/014311600210119]

10. Zhao, P.; Lu, D.; Wang, G.; Wu, C.; Huang, Y.; Yu, S. Examining Spectral Reflectance Saturation in Landsat Imagery and Corresponding Solutions to Improve Forest Aboveground Biomass Estimation. Remote Sens.; 2016; 8, 469. [DOI: https://dx.doi.org/10.3390/rs8060469]

11. Sheridan, R.D.; Popescu, S.C.; Gatziolis, D.; Morgan, C.L.S.; Ku, N.-W. Modeling Forest Aboveground Biomass and Volume Using Airborne LiDAR Metrics and Forest Inventory and Analysis Data in the Pacific Northwest. Remote Sens.; 2015; 7, pp. 229-255. [DOI: https://dx.doi.org/10.3390/rs70100229]

12. Strunk, J.L.; Gould, P.J.; Packalen, P.; Gatziolis, D.; Greblowska, D.; Maki, C.; McGaughey, R.J. Evaluation of Pushbroom DAP Relative to Frame Camera DAP and Lidar for Forest Modeling. Remote Sens. Environ.; 2020; 237, 111535. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111535]

13. Strunk, J.; Packalen, P.; Gould, P.; Gatziolis, D.; Maki, C.; Andersen, H.-E.; McGaughey, R.J. Large Area Forest Yield Estimation with Pushbroom Digital Aerial Photogrammetry. Forests; 2019; 10, 397. [DOI: https://dx.doi.org/10.3390/f10050397]

14. Hudak, A.T.; Lefsky, M.A.; Cohen, W.B.; Berterretche, M. Integration of Lidar and Landsat ETM+ Data for Estimating and Mapping Forest Canopy Height. Remote Sens. Environ.; 2002; 82, pp. 397-416. [DOI: https://dx.doi.org/10.1016/S0034-4257(02)00056-1]

15. Andersen, H.E.; Barrett, T.; Winterberger, K.; Strunk, J.; Temesgen, H. Estimating Forest Biomass on the Western Lowlands of the Kenai Peninsula of Alaska Using Airborne Lidar and Field Plot Data in a Model-Assisted Sampling Design. Proc. IUFRO Div.; 2009; 4, pp. 19-22.

16. Ahmed, O.S.; Franklin, S.E.; Wulder, M.A. Integration of Lidar and Landsat Data to Estimate Forest Canopy Cover in Coastal British Columbia. Photogramm. Eng. Remote Sens.; 2014; 80, pp. 953-961. [DOI: https://dx.doi.org/10.14358/PERS.80.10.953]

17. Erdody, T.L.; Moskal, L.M. Fusion of LiDAR and Imagery for Estimating Forest Canopy Fuels. Remote Sens. Environ.; 2010; 114, pp. 725-737. [DOI: https://dx.doi.org/10.1016/j.rse.2009.11.002]

18. Singh, K.K.; Vogler, J.B.; Shoemaker, D.A.; Meentemeyer, R.K. LiDAR-Landsat Data Fusion for Large-Area Assessment of Urban Land Cover: Balancing Spatial Resolution, Data Volume and Mapping Accuracy. ISPRS J. Photogramm. Remote Sens.; 2012; 74, pp. 110-121. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2012.09.009]

19. Goodbody, T.R.H.; Coops, N.C.; White, J.C. Digital Aerial Photogrammetry for Updating Area-Based Forest Inventories: A Review of Opportunities, Challenges, and Future Directions. Curr. For. Rep.; 2019; 5, pp. 55-75. [DOI: https://dx.doi.org/10.1007/s40725-019-00087-2]

20. Holgerson, J.; Stanton, S.; Waddell, K.; Palmer, M.; Kuegler, O.; Christensen, G. Washington’s Forest Resources: Forest Inventory and Analysis, 2002–2011; Gen. Tech. Rep. PNW-GTR-962 US Department of Agriculture, Forest Service, Pacific Northwest Research Station: Portland, OR, USA, 2018; 112. [DOI: https://dx.doi.org/10.2737/PNW-GTR-962]

21. Bechtold, W.A.; Patterson, P.L. The Enhanced Forest Inventory and Analysis Program: National Sampling Design and Estimation Procedures; US Department of Agriculture Forest Service, Southern Research Station: Asheville, NC, USA, 2005.

22. Heath, L.S.; Hansen, M.; Smith, J.E.; Miles, P.D. Investigation into Calculating Tree Biomass and Carbon in the FIADB Using a Biomass Expansion Factor Approach. Proceedings of the Forest Inventory and Analysis (FIA) Symposium 2008; US Department of Agriculture, Forest Service, Rocky Mountain Research Station: Park City, UT, USA, 2009; Volume 56, 26.

23. Hann, D.W. Equations for Predicting the Largest Crown Width of Stand-Grown Trees in Western Oregon; Oregon State University: Corvallis, OR, USA, 1997.

24. Franklin, J.F.; Spies, T.A.; Van Pelt, R. Definition and Inventory of Old Growth Forests on DNR-Managed State Lands (Section One). Report to the Washington Department of Natural Resources: Olympia, WA, USA, 2005; 22.

25. Andersen, H.-E.; Strunk, J.L.; McGaughey, R.J. Using High-Performance Global Navigation Satellite System Technology to Improve Forest Inventory and Analysis Plot Coordinates in the Pacific Region. Gen. Tech. Rep.; 2022; 1000, 444.

26. McGaughey, R.J.; Ahmed, K.; Andersen, H.-E.; Reutebuch, S.E. Effect of Occupation Time on the Horizontal Accuracy of a Mapping-Grade GNSS Receiver under Dense Forest Canopy. Photogramm. Eng. Remote Sens.; 2017; 83, pp. 861-868. [DOI: https://dx.doi.org/10.14358/PERS.83.12.861]

27. Andersen, H.-E.; Clarkin, T.; Winterberger, K.; Strunk, J.L. An Accuracy Assessment of Positions Obtained Using Survey-and Recreational-Grade Global Positioning System Receivers across a Range of Forest Conditions within the Tanana Valley of Interior Alaska. West. J. Appl. For.; 2009; 24, pp. 128-136. [DOI: https://dx.doi.org/10.1093/wjaf/24.3.128]

28. Clarkin, T. Modeling Global Navigation Satellite System Positional Error under Forest Canopy Based on LIDAR-Derived Canopy Densities. Master’s Thesis; University of Washington: Seattle, WA, USA, 2007.

29. Daly, C.; Halbleib, M.; Smith, J.I.; Gibson, W.P.; Doggett, M.K.; Taylor, G.H.; Curtis, J.; Pasteris, P.P. Physiographically Sensitive Mapping of Climatological Temperature and Precipitation across the Conterminous United States. Int. J. Climatol. J. R. Meteorol. Soc.; 2008; 28, pp. 2031-2064. [DOI: https://dx.doi.org/10.1002/joc.1688]

30. Miller, D.A.; White, R.A. A Conterminous United States Multilayer Soil Characteristics Dataset for Regional Climate and Hydrology Modeling. Earth Interact.; 1998; 2, pp. 1-26. [DOI: https://dx.doi.org/10.1175/1087-3562(1998)002<0001:ACUSMS>2.3.CO;2]

31. McCombs, J.W.I. Geographic Information System Topographic Factor Maps for Wildlife Management. Ph.D. Thesis; Virginia Tech: Blacksburg, VA, USA, 1997.

32. Pierce, K.B.; Lookingbill, T.; Urban, D. A Simple Method for Estimating Potential Relative Radiation (PRR) for Landscape-Scale Vegetation Analysis. Landsc. Ecol.; 2005; 20, pp. 137-147. [DOI: https://dx.doi.org/10.1007/s10980-004-1296-6]

33. Weiss, A. Topographic position and landforms analysis. Proceedings of the Poster Presentation, ESRI User Conference; San Diego, CA, USA, 9–13 July 2001; Volume 200.

34. National Elevation Dataset (NED)|The Long Term Archive. Available online: https://lta.cr.usgs.gov/NED (accessed on 7 August 2018).

35. Kennedy, R.E.; Yang, Z.; Cohen, W.B. Detecting Trends in Forest Disturbance and Recovery Using Yearly Landsat Time Series: 1. LandTrendr — Temporal Segmentation Algorithms. Remote Sens. Environ.; 2010; 114, pp. 2897-2910. [DOI: https://dx.doi.org/10.1016/j.rse.2010.07.008]

36. Cohen, W.B.; Yang, Z.; Healey, S.P.; Kennedy, R.E.; Gorelick, N. A LandTrendr Multispectral Ensemble for Forest Disturbance Detection. Remote Sens. Environ.; 2018; 205, pp. 131-140. [DOI: https://dx.doi.org/10.1016/j.rse.2017.11.015]

37. Bell, D.M.; Acker, S.A.; Gregory, M.J.; Davis, R.J.; Garcia, B.A. Quantifying Regional Trends in Large Live Tree and Snag Availability in Support of Forest Management. For. Ecol. Manag.; 2021; 479, 118554. [DOI: https://dx.doi.org/10.1016/j.foreco.2020.118554]

38. Breiman, L. Random Forests. Mach. Learn.; 2001; 45, pp. 5-32. [DOI: https://dx.doi.org/10.1023/A:1010933404324]

39. Cohen, W.B.; Yang, Z.; Kennedy, R. Detecting Trends in Forest Disturbance and Recovery Using Yearly Landsat Time Series: 2. TimeSync–Tools for Calibration and Validation. Remote Sens. Environ.; 2010; 114, pp. 2911-2924. [DOI: https://dx.doi.org/10.1016/j.rse.2010.07.010]

40. Crist, E.P.; Cicone, R.C. A Physically-Based Transformation of Thematic Mapper Data—The TM Tasseled Cap. IEEE Trans. Geosci. Remote Sens.; 1984; 3, pp. 256-263. [DOI: https://dx.doi.org/10.1109/TGRS.1984.350619]

41. Lutes, D.C.; Keane, R.E.; Caratti, J.F.; Key, C.H.; Benson, N.C.; Sutherland, S.; Gangi, L.J. FIREMON: Fire Effects Monitoring and Inventory System; US Department of Agriculture, Forest Service, Rocky Mountain Research Station: Fort Collins, CO, USA, 2006; 164.

42. Walker, S.; Pietrzak, A. Remote Measurement Methods for 3-D Modeling Purposes Using BAE Systems’ Software. Geod. Cartogr.; 2015; 64, pp. 113-124. [DOI: https://dx.doi.org/10.1515/geocart-2015-0009]

43. Isenburg, M. LASzip: Lossless Compression of LiDAR Data. Photogramm. Eng. Remote Sens.; 2013; 79, pp. 209-217. [DOI: https://dx.doi.org/10.14358/PERS.79.2.209]

44. Gesch, D.B.; Evans, G.A.; Oimoen, M.J.; Arundel, S. The National Elevation Dataset. 2018; pp. 83-110. Available online: https://pubs.er.usgs.gov/publication/70201572 (accessed on 7 August 2018).

45. McGaughey, R.J. FUSION/LDV: Software for LiDAR Data Analysis and Visualization [Computer Program]; USDA, Forest Service Pacific Northwest Research Station: Washington, DC, USA, 2014.

46. Kauth, R.J.; Thomas, G.S. The Tasselled Cap—A Graphic Description of the Spectral-Temporal Development of Agricultural Crops as Seen by Landsat. Proceedings of the LARS Symposia; Purdue University: West Lafayette, Indiana, 1976; 159.

47. R Core Team. R: A Language and Environment for Statistical Computing; R Core Team: Vienna, Austria, 2020.

48. RStudio Team. RStudio: Integrated Development Environment for r; RStudio Team: Boston, MA, USA, 2020.

49. Furnival, G.M.; Wilson, R.W., Jr. Regressions by Leaps and Bounds. Technometrics; 1974; 16, pp. 499-511. [DOI: https://dx.doi.org/10.1080/00401706.1974.10489231]

50. Schwarz, G. Estimating the Dimension of a Model. Ann. Stat.; 1978; 6, pp. 461-464. [DOI: https://dx.doi.org/10.1214/aos/1176344136]

51. Moran, P.A.P. Notes on Continuous Stochastic Phenomena. Biometrika; 1950; 37, pp. 17-23. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15420245][DOI: https://dx.doi.org/10.1093/biomet/37.1-2.17]

52. Riemann, R.; Liknes, G.; O’Neil-Dunne, J.; Toney, C.; Lister, T. Comparative Assessment of Methods for Estimating Tree Canopy Cover across a Rural-to-Urban Gradient in the Mid-Atlantic Region of the USA. Environ. Monit. Assess.; 2016; 188, 297. [DOI: https://dx.doi.org/10.1007/s10661-016-5281-8]

53. Davis, R.J.; Ohmann, J.L.; Kennedy, R.E.; Cohen, W.B.; Gregory, M.J.; Yang, Z.; Roberts, H.M.; Gray, A.N.; Spies, T.A. Northwest Forest Plan–the First 20 Years (1994-2013): Status and Trends of Late-Successional and Old-Growth Forests. Gen Tech Rep PNW-GTR-911 Portland US Dep. Agric. For. Serv. Pac. Northwest Res. Stn. 112 P; 2015; 911.

54. Babcock, C.; Finley, A.O.; Andersen, H.-E.; Pattison, R.; Cook, B.D.; Morton, D.C.; Alonzo, M.; Nelson, R.; Gregoire, T.; Ene, L. et al. Geostatistical Estimation of Forest Biomass in Interior Alaska Combining Landsat-Derived Tree Cover, Sampled Airborne Lidar and Field Observations. Remote Sens. Environ.; 2018; 212, pp. 212-230. [DOI: https://dx.doi.org/10.1016/j.rse.2018.04.044]

55. Fuller, W.A. Measurement Error Models; John Wiley and Sons: Hoboken, NJ, USA, 1987; ISBN 978-0-471-86187-4

56. Hengl, T.; Heuvelink, G.B.M.; Rossiter, D.G. About Regression-Kriging: From Equations to Case Studies. Comput. Geosci.; 2007; 33, pp. 1301-1315. [DOI: https://dx.doi.org/10.1016/j.cageo.2007.05.001]

57. Packalén, P.; Temesgen, H.; Maltamo, M. Variable Selection Strategies for Nearest Neighbor Imputation Methods Used in Remote Sensing Based Forest Inventory. Can. J. Remote Sens.; 2012; 38, pp. 557-569. [DOI: https://dx.doi.org/10.5589/m12-046]

58. Powell, S.L.; Cohen, W.B.; Healey, S.P.; Kennedy, R.E.; Moisen, G.G.; Pierce, K.B.; Ohmann, J.L. Quantification of Live Aboveground Forest Biomass Dynamics with Landsat Time-Series and Field Inventory Data: A Comparison of Empirical Modeling Approaches. Remote Sens. Environ.; 2010; 114, pp. 1053-1068. [DOI: https://dx.doi.org/10.1016/j.rse.2009.12.018]

59. Bell, D.M.; Gregory, M.J.; Kane, V.; Kane, J.; Kennedy, R.E.; Roberts, H.M.; Yang, Z. Multiscale Divergence between Landsat-and Lidar-Based Biomass Mapping Is Related to Regional Variation in Canopy Cover and Composition. Carbon Balance Manag.; 2018; 13, pp. 1-14. [DOI: https://dx.doi.org/10.1186/s13021-018-0104-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30218413]

60. Kennedy, R.E.; Ohmann, J.; Gregory, M.; Roberts, H.; Yang, Z.; Bell, D.M.; Kane, V.; Hughes, M.J.; Cohen, W.B.; Powell, S. An Empirical, Integrated Forest Biomass Monitoring System. Environ. Res. Lett.; 2018; 13, 025004. [DOI: https://dx.doi.org/10.1088/1748-9326/aa9d9e]

61. Noordermeer, L.; Bollandsås, O.M.; Ørka, H.O.; Næsset, E.; Gobakken, T. Comparing the Accuracies of Forest Attributes Predicted from Airborne Laser Scanning and Digital Aerial Photogrammetry in Operational Forest Inventories. Remote Sens. Environ.; 2019; 226, pp. 26-37. [DOI: https://dx.doi.org/10.1016/j.rse.2019.03.027]

62. Fiala, A.C.S.; Garman, S.L.; Gray, A.N. Comparison of Five Canopy Cover Estimation Techniques in the Western Oregon Cascades. For. Ecol. Manag.; 2006; 232, pp. 188-197. [DOI: https://dx.doi.org/10.1016/j.foreco.2006.05.069]

63. Henderson, E.B.; Ohmann, J.L.; Gregory, M.J.; Roberts, H.M.; Zald, H. Species Distribution Modelling for Plant Communities: Stacked Single Species or Multivariate Modelling Approaches?. Appl. Veg. Sci.; 2014; 17, pp. 516-527. [DOI: https://dx.doi.org/10.1111/avsc.12085]