Evaluating Statewide NAIP Photogrammetric Point Clouds for Operational Improvement of National Forest Inventory Estimates in Mixed Hardwood Forests of the Southeastern U.S.

Abstract

The U.S. Forest Service, Forest Inventory and Analysis (FIA) program is tasked with making and reporting estimates of various forest attributes using a design-based network of permanent sampling plots. To make its estimates more precise, FIA uses a technique known as post-stratification to group plots into more homogenous classes, which helps lower variance when deriving population means. Currently FIA uses a nationally available map of tree canopy cover for post-stratification, which tends to work well for forest area estimates but less so for structural attributes like volume. Here we explore the use of new statewide digital aerial photogrammetric (DAP) point clouds developed from stereo imagery collected by the National Agricultural Imagery Program (NAIP) to improve these estimates in the southeastern mixed hardwood forests of Tennessee and Virginia, United States (U.S.). Our objectives are to 1. evaluate the relative quality of NAIP DAP point clouds using airborne LiDAR and FIA tree height measurements, and 2. assess the ability of NAIP digital height models (DHMs) to improve operational forest inventory estimates above the gains already achieved from FIA’s current post-stratification approach. Our results show the NAIP point clouds were moderately to strongly correlated with FIA field measured maximum tree heights (average Pearson’s r = 0.74) with a slight negative bias (−1.56 m) and an RMSE error of ~4.0 m. The NAIP point cloud heights were also more accurate for softwoods (R²s = 0.60–0.79) than hardwoods (R²s = 0.33–0.50) with an error structure that was consistent across multiple years of FIA measurements. Several factors served to degrade the relationship between the NAIP point clouds and FIA data, including a lack of 3D points in areas of advanced hardwood senescence, spurious height values in deep shadows and imprecision of FIA plot locations (which were estimated to be off the true locations by +/− 8 m). Using NAIP strata maps for post-stratification yielded forest volume estimates that were 31% more precise on average than estimates stratified with tree canopy cover data. Combining NAIP DHMs with forest type information from national map products helped improve stratification performance, especially for softwoods. The monetary value of using NAIP height maps to post-stratify FIA survey unit total volume estimates was USD 1.8 million vs. the costs of installing more field plots to achieve similar precision gains. Overall, our results show the benefit and growing feasibility of using NAIP point clouds to improve FIA’s operational forest inventory estimates.

1. Introduction

1.1. National Forest Inventory

Around the world, National Forest Inventory (NFI) programs collect information on the location, composition, and distribution of forest resources using a combination of field measurements and remote sensing data. This information, which is used to estimate forest characteristics at various points in time, is one of the key sources of data for forecasting the effects of climate change and anthropogenic greenhouse gas emissions [1]. To help guide sustainable management, NFI programs often use design-based sampling to obtain unbiased estimates of timber resources in their respective countries [2,3,4]. Traditionally, estimates have focused on basic attributes, such as forest volume and area. However, in recent years, NFI programs have evolved to provide information on a wider range of forest characteristics, from damage and disease to ecosystem services such as biodiversity and carbon pools [5,6]. To meet the growing needs of local land managers and to improve inputs into international reporting (e.g., United Nations Intergovernmental Panel on Climate Change or IPCC), NFI programs employ techniques such as sample intensification [7], small area methods [8], model-assisted and hybrid-estimation [9], and post-stratification [10] to increase the scale and precision of statistical estimates derived from forest inventory data.

1.2. FIA and Nationwide Forest Attribute Estimation

In the U.S., the Forest Inventory and Analysis (FIA) program, managed by United States Forest Service (USFS), has been designated as the design-based NFI to determine the extent, condition, volume, and growth of forests, and monitor changes in forest health. FIA data are collected on a permanent network of fixed sampling plots set approximately every 2400 ha across the U.S. using a random, spatially balanced design [11]. The tree (e.g., DBH, height, species) and stand-level attributes (e.g., basal area, trees per acre) collected by FIA have a wide range of applications [12], such as assessment of carbon pools [13], analysis of timber product outputs [14], calculation of forest growth trends [15,16], and study of species’ range shifts in response to climate change [17,18,19]. Although FIA is the most inclusive and longest running forest monitoring program in the U.S., its statistical estimates still contain uncertainty due to its relatively sparse sampling frame (which is designed to produce strategic-level estimates), number of unsampled plot locations resulting from denied access and hazardous conditions [20], and natural variability in forest composition and structure among sampled plot locations [21].

To increase the precision of forest inventory estimates without growing the number of sampling plots, a statistical technique known as post-stratification can be used. This involves aggregating measurements of forest attributes (response variables) into strata that are more homogeneous than the full dataset, typically through the use of wall-to-wall land cover maps and other sources of classified remote sensing data (e.g., see McRoberts et al., 2002). Reductions in estimate variance are obtained by deriving population means as the weighted mean of each stratum [22,23]. The FIA program has long used post-stratification to increase the precision of its forest inventory estimates starting in the 1980s with double sampling with aerial photography [24]. Over the years much research has been done on the topic [25,26,27], leading to the use of various auxiliary remote sensing datasets to form strata [28]. Currently, all four FIA units use some form of the Landsat-based, National Land Cover Database (NLCD) Tree Canopy Cover (TCC) map [29] for post-stratification.

1.3. Use of Point Clouds in Forest Inventory Applications

Although the spectral information from optical satellite imagery provides greater precision of strata for canopy-level attributes such as the proportion of forest area, it is less effective for attributes related to below-canopy structure and quantity, such as aboveground biomass, because the spectral signal tends to saturate, making it ineffective at distinguishing between mature and old-growth forests [30,31,32]. In contrast to optical remote sensing, estimating canopy-level forest attributes with airborne laser scanning (ALS) (or airborne LiDAR) data has been more successful [33,34]. Recently, many countries have acquired nationwide ALS LiDAR to produce high-resolution topographic elevation data. As a result of this wall-to-wall availability, point cloud data from ALS has become a viable data source for statewide estimation of forest attributes [2,5,35]. Despite the growing availability of ALS data in the public domain, there are still significant challenges to using it in operational reporting workflows within FIA. First, since most statewide LiDAR acquisitions are focused on deriving topographic elevation, most are acquired during leaf-off conditions, which limits the data’s ability to accurately resolve canopy cover, height, and other structural attributes in hardwood-dominated forests. Second, since the topography is relatively stable and ALS is still expensive to acquire, most states and countries do not have plans for repeat collections, further limiting temporal availability of ALS data for ongoing, broadscale forest monitoring programs.

An affordable alternative for wall-to-wall forest structure data at the state level are point clouds derived from Digital Aerial Photogrammetry (DAP; [36]). Using photogrammetric processing techniques (also known as Structure from Motion, [37]), DAP can be used to produce 3D point clouds from overlapping stereo images. Unlike ALS, which uses an active sensor to penetrate the canopy to obtain information about understory vegetation structure and ground surface topography, DAP point clouds only capture the top surface of 3D objects, such as trees and buildings. Because DAP does not fully resolve the ground surface, a high-resolution Digital Terrain Model (DTM) is needed to obtain normalized surface heights, referred to here as a Digital Height Model (DHM). Here we differentiate a DHM, which represents a rasterized surface of normalized heights for all 3D objects on the landscape including trees, vegetation, and non-forest structures such as houses, buildings, and, in their absence, ground and pavement, from a Digital Surface Model (DSM), a commonly used term for a rasterized surface of non-normalized elevations for all 3D objects; and a Canopy Height Model (CHM), a raster of normalized heights for trees and forests, where heights for non-forest objects have been removed, typically using a mask based on land cover data derived from classified satellite imagery [38].

Although ALS captures more within canopy detail in forests and is the primary source for high-resolution DTMs, DAP offers certain advantages to operational forest inventory programs. For example, DAP point clouds are approximately 1/2 to 1/3 the cost of ALS acquisitions [39,40] and in many cases, the stereo images required to produce them are becoming increasingly available from continuous monitoring programs. For example, in the U.S., the National Agriculture Imagery Program (NAIP) collects new, high-resolution stereo imagery over the conterminous U.S. every 2–3 years, and the acquisition cost of its DSM product is less than $0.01 USD/ha [27]. Thus, point clouds derived from NAIP DAP (NAIP point cloud hereafter) are an affordable and repeatable source of high-resolution forest structure data for all states in the conterminous U.S. While shorter update intervals are preferred for tracking changes in areas where forests are disturbed more frequently, the potential 2–3-year return interval for NAIP DAP is less than half the 10-year shelf-life of most ALS datasets used for model-assisted estimation [41]. In addition, recent studies performed using small and heterogeneous datasets showed that the precision of height estimates from DAP point clouds is comparable to, albeit slightly lower than, the height estimates from ALS [42]. For example, ref. [43] reported that DAP and ALS metrics showed similar levels of correlation in terms of forest structural attributes such as height, volume, and biomass. Slope, canopy cover, and image acquisition date can generate significant differences between the height metrics derived from DAP and ALS point clouds [44]. When considering combinations of several biophysical forest characteristics, metrics derived from ALS tend to perform better than those from DAP [45]. Other factors which limit the accuracy of DAP DHMs include the selected DTM cell resolution (e.g., 1 m vs. 10 m DEM), the geo-precision of plot coordinate locations (e.g., use of recreational grade vs. survey-grade, HPGNSS receivers, see [27]), forest type, stand density, and topographic complexity (e.g., slope, aspect, and elevation) [44]. Despite being slightly less accurate, DAP point clouds have been successfully used in several nationwide applications to estimate forest canopy height [46,47], and are thus a viable and affordable option to consider for use in operational forest inventory programs.

1.4. Using NAIP Digital Height Surfaces to Improve FIA Forest Inventory Estimates

Given its ability to enable large-scale, frequent, wall-to-wall assessments of forest structure, there is growing interest in learning more about the quality of NAIP DAP DHMs (NAIP DHM hereafter) and their potential for improving operational forest inventory estimates derived by the FIA program. In the Pacific Northwest region of the U.S., ref. [27] found that using NAIP DHMs for post-stratification reduced the variance of mean forest volume estimates by as much as four times compared to estimates derived solely from simple random sampling. This success, along with other studies indicating the potential for using DAP-derived DHMs to improve estimates of height, density, and volume [48,49], suggest they may also be useful for improving traditional FIA estimates developed for statewide reporting. While the study by [27] directly points to the potential benefits of using NAIP DHMs in the context of FIA estimation, it does so with a limited number of sample plots (only 1/10 of the FIA plots in WA state were used) in a predominantly conifer-based forest system. To date, the effectiveness of NAIP DHMs has not been evaluated in southeastern mixed hardwood forests, where seasonality (e.g., fall leaf senescence), compositional heterogeneity [50], and disturbance rates are higher than in other regions of the country [51,52].

1.5. Objectives

Here, our study aims to evaluate NAIP DHMs in an operational context, using thousands of FIA plots to quantify their accuracy in resolving tree heights, as well as to assess their ability to improve FIA estimates in diverse, southern mixed hardwood forests (Figure 1). To achieve these objectives, we use FIA field-measured trees and ALS LiDAR data to validate the accuracy of height estimates derived from NAIP DHMs collected in 2018 for the states of Tennessee (TN) and Virginia (VA). Using a series of scatter plots and linear regression models, we evaluate the effect of forest type (hardwood vs. softwood), resolution of rasterized DHM (1 m vs. 10 m) and the timing of FIA field measurement on the level of variability and bias found in the three sources of height information (FIA field measurements, ALS LiDAR, and NAIP point clouds). Additionally, we also assess the utility of NAIP DHMs when used for post-stratification of FIA forest area, total volume (in m³), and volume per area (in m³/acre) estimates. Using 10 m wall-to-wall NAIP DHMs, we derive three different strata maps, which are independently used to post-stratify population totals and species group estimates derived with FIA’s standard estimation procedure (described in chapter 4 of [11]). The three strata maps include: 1. the full DHM, which contains height values for both forest and non-forest areas; 2. a CHM, where heights in non-forest areas have been removed using a forest mask derived from a national Landsat-based disturbance map [52]; and 3. a CHM with forest type information (e.g., deciduous, evergreen, and mixed) added from NLCD land cover data. To test the utility of each strata map, we produced two measures of relative efficiency (RE), one comparing the variance from each of the three NAIP stratification maps to the variance from FIA’s current post-stratification approach (referred to as RE_FIA), and another comparing the variance from the best strata map to the variance obtained from simple random sampling (referred to as RE_SRS). Furthermore, we also examine the utility of NAIP DHMs to improve general forest estimate totals for specific hardwood and softwood species groups. We conclude our analysis by making general recommendations for future use of NAIP DAP and derivative DHM/CHM products within FIA and the southeast region.

2. Materials and Methods

2.1. Study Area

Our study area encompasses the states of TN and VA in the southeastern U.S. (Figure 2A). According to FIA’s 2009 state report [53], TN has a land area of 109,150 km², 51.9% of which is covered by forest, with a total growing stock volume of 0.784 billion m³. Most of TN falls in the humid subtropical zone (Köppen climate classification Cfa), except for a few areas in the eastern part of the state that have more of an oceanic (Köppen Cfb) climate. For reporting purposes, FIA divides TN into five survey units: West (T1), West Central (T2), Central (T3), Plateau (T4), and East (T5) (Figure 2B). A large amount of the land area in the West survey unit is comprised of cropland; thus, forest cover is lower than in other survey units. West Central is a mixture of cropland, pasture, woodland, and oak–hickory–pine forest. The Central survey unit includes Nashville, the state’s largest metropolitan area. The East survey unit is located along the Appalachian Mountains. This unit’s elevation difference is approximately 1500 m and Appalachian oak and northern hardwood forests are the dominant forest types. On the other hand, VA has an area of 102,280 km², of which 58.7% was covered by forest in 2016, with a total growing stock volume of 1.116 billion m³ [54]. VA’s climate is humid subtropical (Cfa). FIA divides VA into five survey units for reporting, which include: Southern Mountains (V1), Northern Mountains (V2), Southern Piedmont (V3), Northern Piedmont (V4), and Coastal Plain (V5) (Figure 2B). The Southern and Northern Mountains are located along the ridge of Appalachian Mountains; thus, the geographical and environmental settings of these two units are similar to the East survey unit in TN. The Northern and Southern Piedmont units represent the transitional area between inland mountains and the coastal plain, while the latter survey unit occupies the eastern border of the state. The natural vegetation cover in the coastal plain of VA is dominated by loblolly–shortleaf pine and oak–hickory–pine forest. A basic overview of the forest characteristics, topography, and disturbance rates (based on data from the USFS Landscape Change Monitoring System or LCMS) in each FIA survey unit can be found in Table 1.

2.2. Point Cloud Data

In 2018, FIA and other USFS partners purchased point clouds for the states of TN and VA (among others). In both states the NAIP imagery was acquired using multiple Leica ADS-100 SH100 pushbroom sensors resulting in the collection of 60 cm RGB digital imagery with approximately 30% side lap. Three look angles were collected (25.6° forward, nadir, 19.4° backward) of which only the forward lap was used to derive the stereo imagery. Point clouds were derived using the Leica XPro SGM Module and delivered in .laz 1.4 format with accompanying color infrared (CIR) and RGB values encoded [55]. Contract specifications required the nominal point spacing, or distance between the 3D points, to be double the ground surface distance of the native imagery (or in this case 1.2 m). Both the VA and TN point clouds exceeded these specifications with 0.7 m and 1.0 m point spacings, respectively. Acquisition dates ranged between early August and mid-December (Figure 3D and Table 2A); thus, a portion of the imagery was taken during the leaf-off season for hardwoods. Additionally, because the NAIP point clouds were collected in different seasons and processed by different vendors, there are noticeable differences in point density across the study area (Figure 3A and Table 2B). For example, when TN’s NAIP point cloud was rasterized into a 10 m grid product, about 1% of the pixels had fewer than 30 points per 100 m² while in VA more than 4% of pixels fell in this range. Areas with few if any points were predominately the result of deep shadows caused by low sun angles and steep topography due to the late fall/early winter acquisition period (Figure 3A–C).

As part of our height validation, ALS LiDAR point clouds covering FIA field plot locations were acquired from the USGS 3D Elevation Program (3DEP, [56]). We acquired LiDAR data categorized as Quality Level 2 by the U.S. Geological Survey. Quality Level 2 products have a minimum nominal pulse spacing of 0.7 m and a vertical error of 10 cm, evaluated as the root mean square error (RMSE) in the elevation dimension; in addition, only the data acquired within the last eight years (2011–2019) were included in the Level 2 product [57]. Because the 3DEP 1 m DEMs are derived from the ALS point clouds, their spatial coverage is nearly identical to the coverage shown in red shading in Figure 2A. Note, data comparisons involving ALS are limited to areas with available data.

2.3. Digital Elevation Data

To convert the rasterized NAIP point clouds into normalized heights, we used digital terrain data to represent the ground surface. Here, 1 m and 1/3 arc-second (approximately 10 m) DEMs from USGS 3DEP were used to represent the ground elevation. Both DEM datasets were acquired by several agencies and processed to a common coordinate system and unit of vertical measure. While the 1 m DEM is available for nearly the entire state of TN, large portions of VA are missing 1 m data, including the mountainous highlands and coastal plains located in the west central and eastern parts of the state, respectively (Figure 2A). Unlike the 1 m data, the 10 m DEM is available seamlessly across both states.

2.4. NAIP DHMs and Strata Map Processing

To produce statewide, wall-to-wall DHMs we used the “gridmetrics” tool, part of the FUSION software package developed and distributed by the USFS Pacific Northwest Research Station [58]. To avoid dealing with missing DEMs and the logistical constraints involved with processing and storing numerous high resolution raster datasets, we opted to produce 10 m DHMs for both states. To minimize impacts of missing data and outliers we used the rasterized point density map (Figure 2A) and the North American Forest Dynamics (NAFD) national disturbance map [59,60]. For example, we used the NAFD map to identify forest pixels with five or fewer NAIP DAP points and/or maximum heights greater than 50 m (the height of the tallest measured tree in TN and VA) and assigned them the overall average maximum height value for all forest pixels, which was 18 m. Additionally, areas mapped as water in the NAFD map were assigned a maximum height value of 0.

After processing and filtering the DHMs, we combined them with other remote sensing datasets to form three strata maps for use in post-stratifying the FIA estimates. The three strata maps include: 1. the raw DHM, with surface heights for all objects on the landscape, including vegetation, as well as non-forest objects such as buildings, houses, and other manmade structures (referred to as PS_DHM); 2. a CHM where surface heights of non-forest objects have been removed using a forest/non-forest mask derived from the NAFD disturbance map ([52]; referred to as PS_CHM); and 3. a CHM combined with general forest type information (evergreen, mixed, deciduous) obtained from the 2016 NLCD map ([61,62] referred to as PS_CHM+FT). To develop the three strata maps we first binned the DHM into 8 roughly equal height classes (0–3 m; 3–7.5 m; 7.5–13.5 m; 13.5–18 m; 18–22.5 m; 22.5–27 m; 27–31.5 m; and 31.5–50 m). Then to generate the PS_DHM map we added a ninth class representing areas mapped as water in the NAFD map. Similarly, the PS_CHM map was generated by adding a ninth and tenth class representing non-forest and water areas in the NAFD map. Lastly, for the PS_CHM+FT map we used the NLCD land cover map to assign water and non-forest classes, as well as to further break out the 8 DHM height bins into deciduous (NLCD class 41), evergreen (NLCD class 42), mixed (NLCD class 43) and bottomland hardwood (NLCD class 90) classes. Note that when grouping strata, not all forest types contained pixels from all height bins, thus after iterative recoding, the TN and VA PS_CHM+FT maps contained 14 and 16 forest specific height bins (e.g., deciduous 0–7.5 m; evergreen 0–7.5 m; mixed 0–7.5 m, etc.), plus non-forest and water classes, respectively. Moreover, while the NAFD and NLCD maps both have 30 m pixels, they are only used for masking and filtering the DHMs; thus, the final resolution of all three NAIP strata maps is 10 m.

To test the three NAIP strata maps against the precision gained by FIA’s current post-stratification approach, we also developed a fourth strata map by binning the NLCD TCC map [29] into four separate classes (0–10%; 11–47%; 48–84% and 85–100%). This fourth strata map (referred to as PS_TCC) has a resolution of 30 m and is only used to calculate the relative efficiency gained from using the NAIP height maps instead of tree canopy cover data to post-stratify the FIA estimates (i.e., RE_FIA).

2.5. FIA Data

Field data collected by the USFS FIA program in the states of TN and VA were obtained online from the FIA data mart [63]. Here we use the FIA data at the plot-level to validate the NAIP point cloud heights and as a sample to evaluate the benefit of using NAIP DHMs for post-stratification. The design of the FIA program intends for each permanent plot in TN and VA to be measured every five years [11]; however, due to logistical constraints and other issues, it can take more than 5 years to collect the full cycle of plots. To derive estimates for the time period closest to 2018 (which is the year the NAIP imagery was acquired) we selected the most recent field measurement for each plot, resulting in 2931 forested plots in TN collected from 2012–2019 (with an average n = 366 plots per year) and 3304 plots in VA collected from 2014–2019 (with an average n = 550 plots per year). Each FIA plot is comprised of four, 7.32 m radius circular subplots which are 36.58 m apart and oriented 120°, 240°, and 360° from the plot center. FIA plot coordinates are recorded at the center of subplot 1, which falls in the middle of the 0.042 ha plot area. Here we use the actual FIA plot locations, which are not publicly available; however, since the locations of subplots 2–4 are not precisely known, we only use the FIA tree measurements collected on subplot 1 to validate the point clouds. Since some of the forested plots have no trees on subplot 1 only a subset of measurements from each state are used for the point cloud validation (1936 in TN and 2258 in VA). To determine the tree with the maximum height, we used the FIA tree variable ACTUALHT (m), which is a record of the distance from ground level to the highest remaining part of the tree still present and attached to the bole [64]. Hardwood and softwood species groups are defined using the FIA tree variable SPGRPCD. In the eastern U.S. there are 19 hardwood and 9 softwood species groups (see Table S1 in the Supplementary Materials for individual species groups). Because species group is a tree-level attribute, some degree of simplification is necessary to group plots into similar classes. Here we use the species group of the tallest live tree to assign each subplot to the hardwood or softwood class. While this approach may introduce some errors, we believe they will be minor at the hardwood/softwood level of aggregation because in many cases in southern mixed hardwoods the tallest tree on a plot comes from the most dominant species group, either due to competition (in natural stands) or planting (in managed stands). Volume estimates are based on the FIA tree variable VOLCFGRS (m³), which represents the total cubic volume of wood in the central stem of trees with diameter at breast height (DBH) ≥5.0 inches. Finally, for the post-stratification analysis, information from all 4 FIA subplots is used to derive forest area (in ha), total live volume (in m³) and volume/hectare (in m³/ha) estimates.

2.6. Point Cloud Tree Height Validation

To assess the quality of the NAIP point cloud heights in areas with trees and forest we compared the maximum heights from the ALS and NAIP point clouds to the maximum tree heights taken from the FIA measurements. For this plot-level assessment, we first clip the 3D points falling within a 7.32 m (24 ft) radius circle centered on each FIA plot. Next, using the lidR package in R [65,66] we derive normalized heights by subtracting surface elevation using the corresponding 1 m and 10 m DEMs, then the resulting maximum height values at each resolution are recorded for each plot. This derivation of point cloud maximum height follows the area-based approach [67] and is equivalent to the highest point (in meters) sensed by the point cloud datasets within the 0.02 ha (0.04 acre) area surrounding the center of each FIA plot. Using the FIA tree measurements from this area, which is analogous to the area covered by FIA subplot 1, we present a series of scatter plots to help visualize the impact of increasing the resolution of the DEMs used to derive normalized surface heights, as well as to explore differences in the quality of the NAIP point cloud data by state and species group. The scatter plots include fitted regression lines (R²) and correlation statistics, Pearson r, root mean square error (RMSE), and bias calculated as:

(1)$r=\frac{{\sum }_{i=1}^N\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=1}^N{\left(x_i-\overline{x}\right)}^2}\sqrt{{\sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2}}$

(2)$RMSE=\sqrt{\frac{{\sum }_{i=1}^N{\left(x_i-y_i\right)}^2}{N}}$

(3)$Bias=\frac{{\sum }_{i=1}^N\left(x_i-y_i\right)}{N}$

where y = point cloud height for the ith observation, x = FIA tree height for the ith observation, $\overline{y}$ = average point cloud height and $\overline{x}$ = average FIA tree height. Note for the RMSE and bias estimators x is considered the true value and $\overline{x}$ the true mean.

Since some of the FIA plots were measured up to six years before the NAIP imagery was collected, we used 30 m forest disturbance data from the Landscape Change Monitoring System (LCMS, [68]) to exclude plots that experienced forest cover loss between the year of FIA field measurement and 2018. In addition, a very small number of plots with NAIP maximum heights > 50 m were also excluded, as these were determined to be outliers resulting from noise in the point cloud data.

As 1 m ALS heights are known to be of high quality [69,70], we use them here as part of the validation process to help shed light on the error structure of the NAIP data, as well as to determine if the FIA measurements have more error in certain places or are biased as trees get taller. Given the number of different field crews used to collect FIA data, coupled with challenges associated with measuring tree heights from ground-level, we assume the FIA measurements also contain errors. Therefore, to evaluate the strength of the height relationships we also compare R²’s using reduced major axis (RMA) regression [71], which makes no assumptions about errors in either the X or Y variables. Lastly, to get a sense of how the height errors propagate through time, we present a series of violin plots showing the distribution, median, and direction of height differences between the ALS and NAIP point clouds and the FIA tree heights collected annually from 2012 to 2019.

2.7. Post-Stratified Estimation

To test the utility of the NAIP height maps for post-stratification, we produce a series of FIA estimates using the Forest Inventory ESTimation and Analysis (FIESTA) package available in R [72]. Because the effectiveness of post-stratification is directly related to the correlation between the strata maps and the variable being estimated, we expect the benefit of using the NAIP height strata will fluctuate for different variables and geographic areas of interest. Therefore, estimates for three different variables, forest area (in ha), total forest volume (in m³), and forest volume/area (in m³/ha), are derived independently for each of the survey units found in TN and VA (see Figure 2B for survey unit locations). Using these estimates, we report RE_FIA, which compares the variances obtained from the NAIP strata maps to those from FIA’s current estimation approach (using PS_TCC for post-stratification) and RE_SRS, which compares the variance from the NAIP strata maps to the variance obtained from simple random sampling. For each state and survey unit, RE_FIA is reported for all three NAIP strata maps (PS_DHM, PS_CHM and PS_CHM+FT), whereas RE_SRS is only reported for the strata map producing the highest RE_FIA. For the best strata map we also report the mean value of the variable of interest, its percent standard error, as well as the number of plots gained or lost [derived as (RE_FIA − 1.00) × n], and the apparent sample size after post-stratification (derived as n + the number of plots gained or lost). Note that for certain variables and survey units it is possible the NAIP strata maps will result in higher variance than FIA’s current approach. In these cases, the estimates and standard errors are reported for the tree canopy cover strata map (previously denoted as PS_TCC), whereas the plot gains and apparent sample sizes are reported as losses relative to the initial sample. To help determine which of the FIA variables were most improved by the NAIP height strata, and to get a sense of which of the NAIP maps were most useful for post-stratification, we also provide a summary of these results, showing the overall number of improved estimates and the average number of plots gained/lost by state and survey unit.

A number of helpful functions in FIESTA facilitated importing and clipping the survey unit boundaries and strata maps, as well as assigning the strata map values to the forest inventory plots. After processing the boundaries and strata maps, the modGBpop function is used to define the area units, FIA plot, tree and condition tables, and the map to be used for post-stratification. The characteristics defined by modGBpop are then used in the modGBratio function to generate per-acre estimates by domain and estimation unit using calculations based on chapter 4 of [73]. Here, we use FIA’s base grid plots (i.e., PLOT.INTENSITY = 1) and the ratio of means estimator (described in Section 4.3.4 of [73]) to derive per-hectare estimates of live forest volume using the FIA tree variable VOLCFGRS. To get estimates on a per-hectare basis, the ratio estimator (Equation (4) divides the variable of interest (in this case total live volume) by the estimate of forest area in the population of interest (in this case the FIA survey units in TN and VA), resulting in three separate estimates from which to gauge the reduction in variance that occurs from using the various NAIP height maps for post-stratification.

(4)$\widehat{V{a}_d}=\frac{\widehat{V{t}_d}}{\widehat{F{a}_d}}=\frac{{\sum }_h^H{W}_{h\overline{V{t}_{hd}}}}{{\sum }_h^H{W}_{h\overline{F{a}_{hd}}}}$

Here $\widehat{V{a}_d}$ is the estimated forest volume per hectare in the dth survey unit, $\widehat{V{t}_d}$ is the estimated total forest volume in the dth survey unit, $\widehat{F{a}_d}$ is the estimated forest area in the dth survey unit, ${\overline{Vt}}_{hd}$ is the estimated total forest volume in stratum h of the dth survey unit, ${\overline{Fa}}_{hd}$ is the estimated forest area in stratum h of the dth survey unit, $W_h$ is the weight for stratum h and H is the total number of strata. The variance of the forest volume per hectare estimate $v\left(\widehat{V{a}_d}\right)$ is derived using Equation (5) defined as:

(5)$v\left(\widehat{V{a}_d}\right)=\frac{1}{\widehat{F{a}_d^2}}\left[v\left(\widehat{V{t}_d}\right)+\widehat{V{a}_d^2}v\left(\widehat{F{a}_d}\right)-2\widehat{V{a}_d}cov\left(\widehat{V{t}_d},\widehat{F{a}_d}\right)\right]$

Lastly, to get a sense of how well the NAIP strata perform for hardwood and softwood species, we also develop ratio estimates of live forest volume for the individual species groups found in each FIA survey unit. Note that at this level of estimation, several species had very low sample sizes so only the ones with a percent standard error ≤ 20% are reported. Similar to the broader survey-unit results, we summarize these estimates to determine which of the NAIP strata maps is most effective in reducing variance at the species group-level and to quantify the number of plots gained (or lost) on average for hardwood and softwood species.

3. Results

3.1. Point Cloud Validation

The plot-level validation results comparing the ALS and NAIP point cloud heights and the FIA field-measured tree heights (shown in Table 3) indicate that, regardless of the source or DEM resolution, the point cloud maximum heights underestimated the maximum heights obtained from the field measured trees. Across both states the 1 m ALS point clouds had an average height difference of −0.36 m, which was significantly less than the −1.85 m and −1.27 m average difference observed for the 1 m and 10 m NAIP products, respectively. The results also showed that the average amount of underestimation for both ALS and NAIP is lower for hardwoods (ALS 0.24 m, NAIP −1.34 m) than for softwoods (ALS −0.96 m, NAIP −1.78 m). Despite the overall underestimation, the level of correlation between the point cloud maximum heights and FIA maximum tree heights was moderately to fairly strong, with average across state Pearson r values of 0.82 for the 1 m ALS and 0.76 and 0.72 for the 1 m and 10 m NAIP data, respectively (Table 3). We also found correlations between the point cloud heights and FIA maximum tree heights were much lower for hardwoods (ALS 0.77, NAIP 1 m 0.67, and NAIP 10 m 0.64) than softwoods (ALS 0.88, NAIP 1 m 0.86, and NAIP 10 m 0.79).

The patterns of underestimation of maximum height by the point clouds and the generally better performance for hardwood species revealed in the maximum height comparison are also reflected in the scatter plots and fitted RMA regression lines shown in Figure 4. For example, the scatter plots reveal a much tighter fit and less noise for the point cloud heights derived with the 1 m vs. 10 m DEM data. This is true for softwoods more so than hardwoods. Overall, the 1 m ALS maximum heights have the best model fits, followed closely by the 1 m NAIP (Table 3 and Figure 4). The scatter plots also show a much higher level of noise in the NAIP point cloud data collected in VA. This is apparent by the larger number of significant errors along the bottom part of the x-axis, where NAIP maximum height is at or near zero but FIA maximum heights range from 5–36 m. While the count shading shown in Figure 4 reveals the majority of NAIP heights do fall along the 1:1 line, the fitted relationships show that regardless of DEM resolution, the NAIP point cloud heights are underestimated at the high end (i.e., maximum heights ≥ 30 m) and overestimated at the low end (i.e., maximum heights 0–10 m) compared to the FIA tree measurements. Because ALS LiDAR itself has proven to be an accurate method for measuring the height and structure of trees and vegetation, we also directly compared the 1 m ALS and NAIP maximum heights to help inform the interpretation of their relationships with FIA data. In general, the ALS and NAIP maximum heights are moderately (R² = 0.53) to strongly (R² = 0.87) related, as can be seen by the fitted RMA regression lines shown in Figure 5. These scatter plots show the relationship between NAIP and ALS is stronger in TN than in VA (e.g., R²s in TN are 0.22 to 0.32 higher) and that the relationship improves when the time between acquisition dates is kept within 2 years (Figure 5, right panel). This is especially true in TN where R² improves from 0.75 to 0.87 within this narrower acquisition window. Aside from a few isolated errors at the high end, most of the noise in the NAIP data seems to occur in maximum heights ≤ 25 m. In this range, the NAIP data often significantly underestimates the ALS heights, especially in VA (see lower row, Figure 5). Model fits notwithstanding, the count shading and alignment of the individual errors in relation to the 1:1 line show that above 30 m, the NAIP and ALS maximum heights are generally in good agreement, with little to no systematic bias present (especially in TN, upper right, Figure 5).

Finally, we plotted the differences between the FIA and point cloud maximum heights using a series violin/box plots broken out by FIA measurement year (Figure 6). These plots indicate that, regardless of the DEM resolution or year of FIA measurement, the medians of the point cloud maximum heights are consistently 1 to 5 m less than the FIA measured tree heights. Overall, the box plots in Figure 6 show that across time, the median height differences for the 1 m ALS are lower, and the residual errors are better balanced above and below zero, as opposed to the 1 m NAIP, which has most of its errors consistently below zero. Although the 10 m NAIP contains slightly less canopy detail (Figure 7) its height errors display slightly less systematic underestimation than the 1 m NAIP. Starting in 2016 through the end of 2019 a slight trend of increasing underestimation is observed for both the NAIP and ALS point clouds (Figure 6).

Prior to running wall-to-wall, statewide DHMs and other heights metrics, we used FUSION to process the NAIP data for a number of smaller test areas to help determine the most appropriate grid cell resolution to best meet the needs of our study. Outputs from these initial tests showed that although the 2 m DHMs captured more canopy detail than the 10 m DHMs (e.g., see Figure 7A vs. Figure 7B), the overall differences were not deemed significant enough to justify the additional processing time and storage space required to develop the higher resolution product. Therefore, 10 m statewide DHMs were used as the basis for developing the NAIP strata maps (PS_DHM, PS_CHM and PS_CHM+FT) for TN and VA (Figure 8). In addition, Figure 8 also highlights the difference between the PS_DSM, which includes heights for buildings and other non-forest structures, and PS_CHM, which removes heights in non-forest areas using a forest mask.

3.2. Post-Stratified Estimates

3.2.1. Forest Volume

The NAIP strata maps were used to produce post-stratified volume per area (m³/ha), total volume (m³), and forest area (ha) estimates for the 10 survey units in TN and VA (Table 4). These results show that, for both total forest volume and volume/area estimates, the NAIP strata maps yielded lower variances than the NLCD tree canopy cover map (PS_TCC) currently used by FIA. Although all 20 individual survey unit volume estimates were improved by the use of the NAIP strata maps, the level of improvement varied by survey unit. For example, based on the bold shading in Table 4 (which indicates the best overall performing NAIP strata map per variable in terms of RE_FIA), it is apparent that for forest volume/area and total volume the NAIP strata maps performed best in survey units T4 and T2 in TN and V4 and V2 in VA. The overall highest relative efficiency gains relative to the FIA approach (RE_FIA) were in TN, where both forest volume/area and total volume estimates resulted in maximum RE_FIA’s = 1.46. This level of variance reduction is analogous to a 46% increase in the sample size used to make the estimates, which is analogous to adding 251 plots to the volume/area estimate in T4, and 190 plots to the total volume estimate in T2 (Table 4). In VA, the gains for forest volume were slightly lower than in TN, with maximum RE_FIA’s of 1.38 for total forest volume and 1.33 for forest volume/area. The 38% increase in sample size in V2 is like adding 265 plots to the total forest volume estimate, while the 33% increase in V4 theoretically adds 167 plots to the V4 forest volume/area estimate (Table 4). While the NAIP strata maps resulted in similar levels of RE_FIA for both types of volume estimates, the gains relative to simple random sampling (RE_SRS) were two times higher for total forest volume (which ranged from 2.28 to 2.87) than for forest volume/area (which ranged from 1.24 to 1.58) (Table 4).

To get a sense of the overall performance of the NAIP strata maps across units and to help identify which of the strata maps performed best, we summarized the individual survey unit results in Table 4 to get the average percentage of plots gained/lost resulting from the use of the NAIP strata maps, as well as to sum the number of times each individual map resulted in the highest RE_FIA (denoted by the gray shading in Table 4). This summary (Table 5) indicates that the average percent of plots gained across all the survey units for forest volume ranged between 24% and 36%. The lowest percentage gain for volume, 15% for V3 volume/area, is well above zero indicating there was no penalty for using the NAIP strata maps to post-stratify the forest volume estimates. Table 5 also shows that for forest volume, the PS_CHM strata map produced the lowest variances in 7 of the 10 survey units. This was especially true for total volume, where it had the best results in 4 out of the 5 survey units in both TN and VA (Table 5). Note that because the results in Table 4 allow for ties (with gray shading used for multiple maps) there were two survey units (T4 and V3) where the PS_DHM strata map equaled the performance of the PS_CHM map.

3.2.2. Forest Area

In contrast to the consistently positive efficiency gains observed for the forest volume estimates, the NAIP strata maps were much less effective at improving the estimates of forest area (Table 4). For example, in TN the NAIP DHM strata maps did produce lower variances than the FIA approach in 4 of the 5 survey units; however, in 2 of the improved survey units the RE_FIA’s were barely above 1.00 (see T3 and T5, Table 4). The other two survey units did see slightly higher gains with RE_FIA’s of 1.24 (T2) and 1.11 (T4), which resulted in 100 and 61 additional plots, respectively. In the fifth survey unit (T1), the NAIP strata map generated estimates with higher variance than the current FIA approach, resulting in a theoretical reduction of 21% of its sample plots (Table 4). Overall, on average there was a 4% gain in the number of plots added to the forest area estimates in TN.

On the other hand, in VA, only one of the 5-survey unit forest area estimates was improved by the NAIP strata maps (Table 4). The 8% maximum gain observed in V2 (Table 4) was the maximum achieved in VA, resulting in the theoretical addition of 55 plots to the sample (Table 4). Use of the NAIP height strata in the other 4 survey units yielded forest area estimates with higher variances than the current FIA approach, resulting in a theoretical loss of 22 (V4) to 119 (V1) sample plots (Table 4). Across all the survey units in VA, there was a −8% average reduction in sample size compared to estimates from FIA’s current post-stratification approach (Table 5).

Although the NAIP strata maps were not particularly effective in reducing the variance of the forest area estimates relative to FIA’s current approach, the gains relative to simple random sampling were the highest of all the variables tested (ranging from 2.48 (V2) to 4.05 (T2), Table 4). This level of improvement, which yields a 2.5 to 4× increase in sample size, far exceeds the gains achieved for forest volume (see RE_SRS’s, Table 4).

3.2.3. Species Group Volume

Post-stratified forest volume estimates were also run for the individual species groups found in the 10 survey units in TN and VA. Each of the survey units contained a unique combination of the 9 softwood and 28 hardwood species groups used by FIA in the eastern U.S (Table S1). In total there were 84 (TN) and 80 (VA) species group estimates which produced significant results (i.e., standard errors ≤ 20%). Similar to the general forest volume estimates reported in Table 4, the results for the individual species groups (Table S2) show that in 4 out of 5 survey units in TN and 3 out of 5 survey units in VA, a softwood species group achieved the maximum gain from using the NAIP strata maps. In TN, the volume/area estimate for the other eastern softwood group in unit T3 had the maximum RE_FIA of 1.72, while in VA the maximum RE_FIA of 1.60 was observed for the eastern white pine group in unit V2 (Table S2). To get a sense of performance of the NAIP strata maps across units, and to determine which map was most effective at lowering the variance of the species group volume estimates, we summarized the individual species results by grouping them into hardwood and softwood classes. These results (Table 6) show that across both states, forest volume/area was improved the most by the NAIP height strata, with 85% (TN) and 86% (VA) of its species group estimates improved versus 75% (TN) and 66% (VA) for total forest volume (see the number of improved species, Table 6). Overall, 78% of the significant species group volume estimates were improved by the NAIP strata maps. Table 6 also shows that on average softwood species groups tended to gain the most, with average gains of 10–20% versus 4–7% for hardwoods. Note there were some penalties for using the NAIP height strata, as some of the hardwood and softwood species groups had minimum gains well below zero, indicating a theoretical loss of plots compared to the use of the tree canopy cover strata (PS_TCC).

To determine which NAIP strata map performed best for the species group volume estimates, we tabulated the number of times each map recorded the highest gain (Table 6). Overall, across both states the NAIP PS_CHM+FT map had the highest RE_FIA in 59% of the improved species group volume estimates. The results also show that the PS_CHM+FT map had the highest RE_FIA in 65% of the improved species group models in TN and 52% in VA. The other two NAIP strata maps, improved ~20% (PS_CHM) and ~30% (PS_DHM) of the species group volume estimates in TN and VA, respectively.

4. Discussion

4.1. Validation of NAIP DAP and ALS Point Clouds

To evaluate the quality of the remote sensing point clouds we used a 3-way validation approach to compare ALS and NAIP DAP maximum surface heights with FIA field measured tree heights. This approach not only allowed a robust evaluation of the point clouds ability to accurately resolve the maximum height of trees and forest canopies, but given the high accuracy of ALS LiDAR (e.g., [33,74,75]), it also allowed us to evaluate potential bias in the FIA field measurements. As expected, these comparisons (presented in Figure 4 and Figure 5 and Table 3) show the ALS point cloud maximum heights are in better agreement with the FIA tree heights than the 1 m NAIP data, as evidenced by their higher correlations (r ranging from 0.74–0.88 for ALS vs. 0.63–0.89 for NAIP) and levels of explained variance (R²s ranging from 0.55–0.78 for ALS vs. 0.39–0.79 for NAIP).

Perhaps a more important criterion of data quality is the level of systematic bias in the point cloud datasets. Overall, regardless of the grid cell size of the underlying DEM, the species group, or the state of data collection, the NAIP point clouds underestimated the FIA measured tree heights by an average of −1.56 m (Table 3). Similar levels of underestimation, particularly at the high end of the height range (>15 m), have been reported in other studies (e.g., [76,77]). We also found that the NAIP point clouds had higher levels of bias for softwoods (average bias = −1.78 m) than hardwoods (average bias = −1.34 m), even though in most cases the softwood heights were more highly correlated with the FIA measurements. This occurs because the flaws and artifacts introduced by disturbance and phenology have a lower impact on the NAIP softwood heights (particularly in TN, see Figure 4), so the bias calculations are more impacted by the underestimation that occurs at the high end of the height spectrum (i.e., maximum heights > 30 m). On the other hand, the NAIP hardwood heights displayed a lot more noise that resulted in distinct overestimation at the low end and underestimation at the high end of the height range (see Figure 4), which is much more problematic to account for when making wall-to-wall maps. This phenomenon also explains why the 10 m NAIP heights have less bias than the 1 m NAIP heights (−1.27 m vs. −1.85 m, Table 3), although the 10 m data clearly have more error and residual scatter in their fitted relationships with the FIA measurements (Figure 4). These examples illuminate the challenge of interpreting bias simply by the calculated values for the different point cloud datasets.

4.2. NAIP Point Cloud Anomalies

While ALS performed better overall, a closer look at the scatter plots in Figure 4 reveals that most of the improvement occurred in hardwoods, where ALS had far fewer errors along the lower part of the x-axis (e.g., where FIA maximum height is high and point cloud maximum height is near zero). These errors, which are evident in all the NAIP plots except TN softwoods (top two rows of Figure 4), have two causes. First, some of these errors are caused by gaps of missing data where the NAIP imagery was acquired in leaf-off conditions. In areas with missing leaves the photogrammetric processing algorithms were unable to successfully triangulate enough 3D points, resulting in zero heights above ground. Although most of the NAIP images covering the study area were acquired between October and December (i.e., late fall/early winter, see Figure 3D), the impact of missing leaves was much more pronounced in VA, where a notable increase in the frequency and magnitude of errors along the x-axis can be seen especially for hardwoods (Figure 4). Overall, the effect of missing leaves was mostly confined to eastern TN and western VA (i.e., survey units T5, V1 and V2, Figure 2) along the higher elevations of the Blue Ridge mountains, where cooler temperatures resulted in earlier periods of senescence and leaf-drop compared to lower elevations. The high spatial variability in leaf-off conditions across the entire study area precluded us from using acquisition date alone to identify areas where the quality of the NAIP point cloud was poor. In fact, in many cases (particularly in TN) the late season imagery captured active periods of senescence where leaves were in the process of changing colors but had not yet fallen off the trees. Our analysis found that in these situations the NAIP height estimates still tended to be of high quality.

The late season acquisition window and steep topography also caused isolated patches of deep shadow, which resulted in extremely low point returns. Like the areas of missing leaves, shadowed areas (highlighted in red, Figure 3A–C) tended to occur in isolated patches containing 1–30 points per 10 m grid cell. Unlike areas of missing data, which had 0 points and no height estimate, areas of deep shadow often resulted in erroneously high height values, a finding that has been previously reported in other studies [76]. For the point cloud height validation, we avoided these areas by eliminating plots with estimated height values > 50 m, which is the height of the tallest tree measured by FIA in TN and VA. Overall, fewer than 10 plots were omitted for this reason. Although areas of missing data (due to missing leaves, deep shadows and other anomalies) represented a small percentage of the full NAIP data set, we opted to leave the impacted plots in the validation since these areas are inherent flaws in the data that ultimately limit its effectiveness when used for post-stratification of the FIA estimates (discussed below).

4.3. Other Factors Impacting NAIP Point Cloud Quality

4.3.1. Forest Disturbance

Although the inherent flaws in the NAIP imagery did limit its overall quality, the biggest factor contributing to the large errors in the NAIP point cloud heights is forest management and disturbance. Although we attempted to remove impacts of forest harvesting using LCMS data, there were still many instances where this screening process did not remove all artifacts related to disturbance, resulting in substantial deviations in observed and predicted heights along both the x- and y-axes. This is particularly evident in VA where partial harvesting and thinning operations, which have been shown in other studies to cause errors in NAIP CHMs [77], resulted in additional noise.

Overall, there were far more errors in the VA NAIP data (Figure 4, right two columns) than in TN (Figure 4, left two columns), likely due to better image acquisition dates in TN resulting in higher point densities compared to VA (Figure 3A–D), as well as lower rates of disturbance (10.8% in TN vs. 14.7% in VA based on LCMS data). We also found there was very little impact of using the 1 m vs. the 10 m resolution DEM for height normalization of the NAIP point cloud data, as evidenced by the fact that the amount of residual scatter around the fitted RMA regression models in TN (Figure 4) did not significantly worsen with coarser DEM resolution. This bodes well for the use of 10 m DEMs when running wall-to-wall height maps for use in post-stratification.

4.3.2. Temporal Offset

One aspect that contributes to the noise added by disturbance is the amount of time between the FIA measurements and the date the NAIP images were collected. Figure 6 shows that neither the DEM resolution (1 m vs. 10 m) nor the number of years between field measurement and image acquisition (ranging from 6 years before to 1 year after) significantly impacts the error structure or bias of the NAIP maximum height estimates. For example, the median errors observed across time range between slightly less than 0 and −3.5 m, with the majority (25–75%) of data each year falling within +/− 3.0 m (Figure 6). While Figure 6 does indicate a slight increase in the amount of underestimation as the FIA year approaches the NAIP acquisition year (i.e., 2018), the effect is minimal, with very little of the data exceeding +/− 5.0 m. Overall, the lack of temporal degradation in the NAIP heights is encouraging, as FIA estimates are always made with plots that range backwards in time across a number of years, thus our results suggest NAIP point clouds should continue to provide gains in precision for a number of years after its acquisition date.

4.3.3. GPS Imprecision

Another factor that is difficult to quantify but certainly impacts the level of agreement between the point clouds and field plots is the inherent imprecision of the FIA plot locations. In the southeastern U.S., FIA currently uses recreational grade GPS receivers to collect plot coordinates in the field. Although a broad assessment of FIA plot coordinate precision is beyond the scope of this paper, initial results from a pilot study aimed at collecting High-Precision Global Navigation Satellite System (HPGNSS) coordinates for forested FIA plots in VA reveals that typical FIA plot coordinates collected with recreational grade GPS receivers have an RMSE of ~8.0 m (n = 31) when compared to HPGNSS locations. Although [41] concluded that global positioning error was not a significant factor in determining the level of LiDAR-assisted post-stratification improvement of aboveground biomass estimates in Minnesota, their reasoning was based on the fact that plots measured with HPGNSS and recreational receivers resulted in similar mean estimates of biomass, despite finding that the ratio of standard errors indicated a relative gain in efficiency of around 50% when using the more precise coordinates. This level of gain is equivalent to adding 50% more plots to the FIA sample, which is a significant improvement in both precision and potential cost savings that should not be overlooked. While we did not explicitly examine the effects of plot location accuracy, other studies have shown offsets from imprecise GPS locations can significantly degrade the relationship between point cloud heights and FIA field measurements, leading to a reduced gain in estimated precision compared to plots with HPGNSS locations (e.g., see Table 1 in [27]).

The amount of noise added due to GPS inaccuracy will undoubtedly vary, depending on the resolution of the underlying DEM and the heterogeneous nature of the landscape in which the plot is located. For example, if a 10 m DEM (which covers 100 m²) is used to normalize the point cloud heights, the currently observed offset in VA (+/− 8 m) could cause a plot to be linked to a pixel that mostly lands outside the field measured subplot, which covers 168.16 m² (or roughly a 13 m pixel). On the other hand, if a 1 m DEM is used, and an area-based calculation is performed summarizing all the pixels that are purported to cover the FIA subplot area, the GPS plot locations will result in pixels outside the plot area contributing to the calculation; however, even if a portion of the field area is covered, the inclusion of some of the plot pixels helps smooth the noise and ultimately limits the effect of the spatial location error. This is especially true in dense forests where shifting the spatial location of plot center has little effect due to homogenous canopy conditions. However, in highly managed forests where frequent harvesting and thinning result in a patchier, more structurally diverse landscape, the effect of GPS imprecision may be more significant. For instance, in TN where rates of forest management are lower and the quality of the NAIP imagery is higher, upscaling the NAIP point cloud maximum heights from 1 m to 10 m did not result in an expansion of errors or addition of significant outliers compared to the FIA measurements. On the other hand, in VA, where disturbance rates from harvesting are much higher and the quality of the NAIP imagery is lower, the upscaling from 1 m to 10 m resulted in a proliferation of new errors along the x-axis. While this evidence is far from conclusive, it does warrant further investigation, especially as FIA seeks ways to maximize the efficiency gained from using point clouds in its operational estimation system.

4.3.4. FIA Tree Measurement Uncertainty

Because nearly all the point clouds displayed some level of underestimation for maximum heights > 30 m and given other studies have shown errors in field measured tree heights as trees get taller [78,79], we compared the 1 m NAIP and 1 m ALS maximum heights directly, with the idea of using the proven accuracy of ALS to check for potential bias in the FIA measurements. As with the FIA height validations (presented in Figure 4), we also attempted to remove the effects of disturbance using the LCMS data. Despite our attempts, Figure 5 shows a good amount of residual scatter remained after filtering for disturbance, particularly when the time between measurements exceeded two years (left column, Figure 5). The upper right panel of Figure 5 shows the best overall relationship between ALS and NAIP maximum height occurred in TN, where the lower rate of disturbance and reduced measurement gap of two years yielded a strong linear relationship (R² = 0.87) between the two remote sensing-based height estimates. The fact that this relationship also displays very little residual scatter and more importantly, next to no bias at either end of the height spectrum, suggests the remote sensing point clouds may not underestimate the FIA height measurements as much as Figure 4 might suggest. This is especially true as other studies have also shown better, less biased relationships between LiDAR and NAIP point cloud heights than between NAIP and field measured tree heights (e.g., see [76]). Since FIA tree heights are only recorded to the nearest foot from ground level and considering the challenges associated with locating the highest part of upper canopy trees in dense deciduous forests, there are no doubt errors in the FIA tree heights. Although the ALS and NAIP maximum heights in VA are much noisier (and thus less valuable for evaluating agreement) there is still a noticeable lack of systematic bias evident in these graphs (especially in the lower right panel, Figure 5), which further supports the fact that, despite the effects of disturbance and other noise, the NAIP point cloud heights are in excellent agreement with the best available remotely sensed measures of height from ALS LiDAR.

4.3.5. Post-Stratifying FIA Estimates with Digital Height Maps vs. Tree Canopy Cover Data

Our other main objective, aside from evaluating the quality of the NAIP point cloud heights, was to quantify the level of precision gained from using NAIP statewide height maps to post-stratify FIA forest area and volume estimates. The results in Table 4 clearly show that for forest volume, the height based NAIP strata maps did a better job of reducing variance than the strata maps based on tree canopy cover data. In fact, although the level of relative efficiency gained varied by state (avg RE_FIA in TN = 1.35 vs. VA = 1.28) and survey unit (min/max RE_FIA across TN = 1.28–1.46 vs. 1.15–1.38 in VA), the NAIP height maps outperformed canopy cover for all 20 volume estimates in TN and VA (Table 4). Overall, for forest volume we found that on average, the best NAIP strata map (shaded in gray, Table 4) resulted in a 31% gain in sample size compared to the tree canopy cover maps. This level of improvement resulted in the theoretical addition of 175 sample plots per survey unit on average. Although the least number of plots gained for any one survey unit was 71 (V3 in VA Northern Piedmont, Figure 2), the other survey units gained between 126–265 plots depending on the initial sample size and strength of the NAIP heights relationship with forest volume. The quality of this relationship hinges on the level of noise introduced by GPS imprecision (resulting in spatial mismatch between plots and pixels) and our assignment of plots to strata using a single 10 m pixel (which introduces error because only a portion of the area sampled on the ground is being captured by the remote sensing data).

Because the amount of efficiency gained in the sample is based on the strength of the relationship between the map and the variable of interest (i.e., R², see [27]) substantial errors in the Y direction will serve to flatten the fitted linear relationship, causing the predictions to regress to the sample mean, which ultimately negates the benefit of having wall-to-wall height maps. These errors, which are readily apparent in scatter plots between NAIP maximum height and FIA forest volume (m³/ha) (Figure 9), ultimately reduce the strength of the fitted regression relationships, especially for softwoods, where the outliers along the y-axis capped R² below 0.5. Minus these outliers, the R²’s for softwoods would likely have approached 0.7 or higher given their similarity to the relationships presented in Figure 7 of [27]. Due to the exponential relationship between R² and RE shown in Figure 5 of [27], there is likely a twofold level of improvement that could still be obtained for softwood volume just by eliminating the errors in the NAIP height maps. In light of these potential improvements, Figure 9 also highlights the enormous potential of NAIP height maps to achieve even further gains in efficiency vs. tree canopy cover maps, which in reality have a poor relationship with volume and only serve to reduce estimate variance because of a slightly positive, but mostly spurious relationship with volume (Figure 9, bottom row).

What is also apparent from Table 4 is that, unlike the improvements realized for volume, the NAIP height maps were not a clear improvement over the canopy cover maps when it came to stratifying the FIA forest area estimates. Although 4 of the 5 survey units in TN were improved by the NAIP height maps, the level of relative efficiency gained (avg RE_FIA = 1.10) was 20% lower than was observed for volume. This average level of improvement is somewhat misleading because survey units T3 and T5 only saw 1–2% improvements in sample size, leading to only 6 and 15 plots gained, respectively (Table 4). The fact that only 1 of 5 survey units in VA was improved further draws into question the usefulness of the NAIP height maps for stratifying forest area estimates. Overall, the NAIP strata maps mostly resulted in less precise forest area estimates which at times caused a theoretical loss in sample size in VA compared to the estimates stratified with the tree canopy cover map. This is not terribly problematic, nor is it terribly surprising. Canopy cover is a better indicator of where trees and forests are located spatially (due to the satellite data’s ability to capture information about leaf area), and because the tree canopy cover map is in the public domain and is a product co-produced by the USFS, it can continue to be used by FIA for stratification at no additional cost to the program.

Although not as widely successful as the survey unit estimates, the NAIP strata maps also did a nice job of improving the species group volume estimates. With 66–86% of the significant species group estimates gaining some benefit from the NAIP strata maps it is clear they still add value, even though on average only 5 and 15 plots were gained for hardwoods and softwood models, respectively (Table 6). Although softwood species groups typically had fewer plots, they clearly benefited more from the NAIP height maps than hardwoods, with maximum gains adding between 52–72 plots to their overall sample size (Table 6). Despite this success there were still times when the NAIP strata maps were not an improvement over tree canopy cover, resulting in a small loss of precision (and sample size) for both hardwood and softwood species groups. The same sources of error discussed above (and highlighted in Figure 9) impact the strength of the species group post-stratifications and, because the sample size is much smaller, the presence of any outliers in the NAIP height maps is particularly problematic. For example, our results showed only a few extreme outliers (primarily in the Y direction) were responsible for limiting efficiency gains for the Loblolly/shortleaf pine and Other yellow pine volume estimates (FIA species codes 2 and 3—Table S1).

The last question pertaining to post-stratification is which of the NAIP strata maps is best? While our results are far from conclusive, they do show that, for traditional FIA estimates made at the survey unit level, combining the NAIP DHM with a secondary forest mask (in this case from NAFD) to produce a CHM (or PS_CHM) led to higher efficiency gains for many of the forest area and volume estimates (Table 4 and Table S2). In fact, of the 83% of survey unit estimates that were improved by the NAIP height data across TN and VA, 67% of them involved use of the NAIP CHM (PS_CHM). Although there were times when the other NAIP strata maps (PS_DHM and PS_CHM+FT) yielded better results, we note that in most of these instances the efficiency reported for PS_CHM was only 1–2% lower than these other maps. It was also evident that, at this scale of estimation, adding the extra information about forest type to the NAIP height data (e.g., using NLCD) was not particularly useful.

However, when estimates were made for specific species groups (Table S2), the added forest type information in the PS_CHM+FT map began to take on added importance. For example, of the 78% of the species group estimates that were improved by the NAIP height maps, just over half (53%) were attributed to the NAIP CHM with added forest type information (Table 6 and Table S2). While the percentage improvement tended to be lower at this scale of estimation, we note that in many cases the difference varied between realizing slight positive gains in efficiency vs. small losses in sample size compared to the tree canopy cover strata maps (Table S2). Overall, our results show that strata maps developed from NAIP height information are clearly better at reducing the variance of forest volume estimates than strata derived from tree canopy cover data, and when the NAIP strata are better, information from NAFD (e.g., PS_CHM) and NLCD (PS_CHM+FT) can lead to even bigger gains in efficiency than just using the wall-to-wall NAIP DHM (PS_DHM), which only achieved the highest gains for 20% and 23% of the survey unit and species group estimates, respectively.

4.3.6. Potential for Operational Use of NAIP Point Clouds within FIA

Because the net effect of post-stratification can be quantified in terms of the theoretical number of plots gained (or lost) relative to the initial sample size, perhaps the best way to determine the level of added benefit provided by the statewide NAIP point clouds is to compare the amount it would take to field measure the number of additional plots gained from post-stratification vs. the total cost of the NAIP point cloud data. For example, at a cost of USD 1000 per field visit the 71 plots added to the forest volume/area estimate in unit V3 (Table 4) equates to USD 71,000, which is slightly more than it cost to purchase both the TN and VA NAIP point clouds. In other words, for reasonably sized states, the potential savings achieved in just one survey unit will more than cover the cost of purchasing the NAIP point cloud data. At the high end, stratifying the survey unit total forest volume estimates with the NAIP height maps added 938 and 951 plots to the FIA sample in TN and VA, respectively, which in total nets nearly USD 2 million in potential cost savings from not having to measure 1889 additional field plots. As other studies have reported similar precision gains when using NAIP point clouds with area-level small area estimation approaches [80], there is growing evidence to suggest NAIP DAP is a highly cost-effective option for increasing the precision of NFI forest volume estimates. While these numbers are impressive, they omit the costs associated with processing and storing the NAIP point clouds, which are not insignificant. Based on our experience, the level of effort to produce 10 m statewide DHMs in FUSION can be a bit time consuming but is otherwise a straightforward process that should not impose significant operational or financial burdens; therefore, FIA should consider wider adoption and use of NAIP point clouds in its operational estimation system. We do note that computational limitations and storage requirements still make it difficult to operationally produce wall-to-wall 2 m DHMs, although we envision with continued improvements in computing technology that this barrier will be lessened in the future. Once wall-to-wall DHMs are produced it is not difficult to incorporate the NAIP height data into different stratification routines as we did here with FIESTA. Our results show use of different stratification maps will result in different levels of precision gain, so further testing is needed to fully understand the best way to maximize the use of NAIP point clouds in FIA’s estimation system. As illuminated here, future FIA estimates may well benefit from use of both height and tree canopy cover maps, especially if more advanced model-assisted estimation approaches are used.

5. Conclusions

Our results revealed several interesting findings about the quality of NAIP DAP point clouds and their potential utility for improving forest inventory estimates in southeastern mixed hardwood forests. Overall, regardless of the underlying DEM resolution (1 m or 10 m), the NAIP point clouds underestimated FIA maximum tree heights by an average of −1.56 m but otherwise were in strong agreement with the field measured trees (R²’s = 0.58–0.89). The biggest difference between the NAIP point cloud and field measured maximum tree heights occurred at the state level, where the underestimation in VA was 1 m more than in TN (bias in TN = −0.99 m vs. −2.13 m in VA), likely due to the use of different vendors and processing routines. In terms of post-stratification, using the NAIP DHMs instead of tree canopy cover maps resulted in forest volume estimates which were 15 to 46% more precise than FIA’s current approach. The additional precision gained from using the NAIP height data added 175 plots per survey unit on average to the FIA sample, offering huge financial savings compared to the costs of installing more field plots. Although more testing is needed to determine the full benefit of using NAIP DAP point clouds in other forest types across the U.S., our results (and those of others, e.g., [27,76,77,80]) highlight the benefits, potential, and growing feasibility of using NAIP DAP point clouds to improve the precision of FIA’s operational forest inventory estimates. Moving forward, future work will focus on using NAIP DHMs and tree canopy cover data with more sophisticated model-assisted approaches, such as the generalized regression estimator (GREG), which uses linear regression and calibration techniques such as raking and lasso to improve estimate precision. More work is also needed to improve the quality and consistency of the photogrammetric models used to derive 3D points from the NAIP stereo images so that point clouds from different states and acquisition dates can be used more consistently over time to estimate changes in forest growth. Collecting high precision GPS data at FIA subplot locations would help reduce misalignment errors, leading to cleaner, less noisy models and improved stratification performance.

Author Contributions

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

Data Availability Statement

NAIP imagery is publicly available via U.S. Department of Agriculture. NAIP point clouds were purchased by USFS and are available by contacting the first author. 3DEP LiDAR data are publicly available from the U.S. Geological Survey (USGS). FIA data are publicly available from USFS, but actual plot locations are withheld to protect landowner privacy.

Acknowledgments

The authors appreciate the assistance of FIA staff and partners who have contribute to this research. Special thanks to Tracey Frescino, USDA Forest Service, Rocky Mountain Research Station for assistance with the r FIESTA package and John Pemberton, Virginia Department of Forestry, for collecting the high precision HPGNSS GPS data. This study was subject to agency review and approved for publication. Any use of product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. This paper was written and prepared by U.S. Government employees on official time, and is therefore, in the public domain and not subject to copyright.

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 workflow used to validate the NAIP point clouds and assess the NAIP digital height model (DHM)/canopy height models (CHMs) ability to post-stratify the FIA estimates.

Figure 2. Study area showing (A) the states of Tennessee (TN) and Virginia (VA) with red outline indicating the 1 m DEM coverage and (B) the elevation and survey units in each state. Combinations of letters and numbers represent survey units. T1: West, T2: West Central, T3: Central, T4: Plateau, T5: East, V1: Coastal Plain, V2: Southern Piedmont, V3: Northern Piedmont, V4: Northern Mountains, V5: Southern Mountains.

Figure 3. NAIP point cloud density and acquisition month. (A) Point density raster showing the number of 3D NAIP points in each 10-m grid cell. (B) A zoomed view of the point density raster. (C) A zoomed view of the NAIP imagery natural color composite. (D) NAIP imagery acquisition month.

Figure 4. Scatter plot of FIA maximum tree heights (x-axis) vs. maximum point cloud heights (y-axis). Columns are broken out by state and species group (i.e., hardwood and softwood), while rows show the point cloud heights obtained for NAIP and ALS using different DEM resolutions. Plots are shaded from blue (low) to yellow (high) to indicate the density of observations that occur in places where multiple plots overlap. Blue lines are the reduced major axis regression (RMA) fits and the red dashed lines are the 1:1 lines.

Figure 5. Comparison of the 1 m ALS and 1 m NAIP point cloud heights at the FIA plot locations. Columns split the dataset into the plots where ALS was collected more (left) or less than 2 years from the NAIP acquisitions, which were collected in 2018; rows break out the plots collected in TN (top) and VA (bottom). Shading from blue (low) to yellow (high) indicates the density of observations in places with multiple overlapping plots. Blue lines are the reduced major axis regression (RMA) fits and the red dashed lines are the 1:1 lines.

Figure 6. Differences between FIA maximum tree height and point cloud maximum height by FIA measurement year. Dashed red lines are drawn at +/− 5 m and dashed black lines at zero. For the boxplots, the bottom and top of the rectangles span the 25% to 75% height range, while the thick black bars are the median values. Columns break out the plots by point cloud source (ALS or NAIP) and DEM resolution (1 m and 10 m), while rows highlight the data collected in different states (TN top, VA bottom).

Figure 7. Comparison between 1 m NAIP natural color composite (left), 10 m NAIP DHM (center), and 2 m NAIP DHM (right) for softwood (A) and hardwood (B) areas in the state of TN.

Figure 8. Comparison between 1 m NAIP natural color composite (A), 10 m DHM (B), and the three 10 m NAIP strata maps PSDHM (C), PSCHM (D), and PSCHM+FT (E) in the state of TN (35°6′11″N, 84°22′44″W).

Figure 9. Scatter plots showing the relationship between NAIP maximum height (in meters, top row) and NLCD tree canopy cover (in %, bottom row) vs. FIA volume per area (m3/ha) for hardwood (left column) and softwood (right column) species groups.

Table 1

FIA measurement summary.

Survey Unit	Plots	Total Volume	Forest Area	VPA	Max Height	Mean Height	Forest Loss	Mean Elevation	Major Species Group
		Million m3	%	m3/ha	m	m	%	m
West (T1)	457	0.91	36.91	166.86	24.12	16.22	35	128.56	Loblolly and shortleaf pines
West Central (T2)	412	0.91	66.03	139.4	23.98	16.3	41	198.52	Select white oaks
Central (T3)	541	1.10	43.15	157.51	23.86	16.23	11	246.34	Hickory
Plateau (T4)	546	1.20	66.39	165.49	24.98	16.54	34	480.59	Other red oaks
East (T5)	707	1.51	55.33	184.99	25.76	16.99	17	585.99	Other white oaks
Coastal Plain (V1)	694	1.52	59.36	190.63	23.14	15.7	51	35.49	Loblolly and shortleaf pines
S. Piedmont (V2)	694	1.56	68.36	165.4	22.91	15.78	48	205.78	Loblolly and shortleaf pines
N. Piedmont (V3)	474	1.01	56.19	203.34	25.94	17.52	28	253.92	Yellow-poplar
N. Mountains (V4)	503	1.12	64.51	172.34	22.78	15.66	14	626.32	Other white oaks
S. Mountains(V5)	584	1.24	64	191.05	24.97	17.04	19	757.8	Yellow-poplar

Summary of forest attributes in each FIA survey unit. Total volume, Volume per Area (VPA), and Forest area are estimated from the FIA sample using simple random sampling. Max height represents the tallest measured tree while Mean height is the average height of all the trees measured in each survey unit. Forest loss was estimated using disturbance data from the U.S. Forest Service (USFS) Landscape Change Monitoring System (LCMS). Major species group represents the most frequently observed species group in each survey unit. See Table S1 for the full list of FIA species groups.

Table 2

NAIP point density by acquisition month.

(A) Percentage of NAIP imagery in each acquisition month (%).
	July			August	September	October		November	December
Tennessee	4.13			12.20	6.43	64.25		12.99	0.00
Virginia	0.00			35.71	0.88	38.08		19.14	6.18
(B) Percentage of 10 m pixels in each point density class (%) and the Min, Mean and Maximum point density in each state.
	0–30	31–80	81–130	131–180	181–230	231–280	280+	Min	Mean	Max
Tennessee	0.95	67.92	12.42	15.51	2.63	0.57	0.24	0	86.41	566
Virginia	4.21	64.19	8.38	19.42	3.12	0.68	0.19	0	88.12	447

Table 3

Summary of maximum height comparison between the NAIP point cloud data and FIA field measured trees.

Source	State	Resolution	Species	Plots	Avg. Point Cloud Height	Avg. FIA Height	r	RMSE	Bias
		Meter		n	Meter	Meter		Meter	Meter
NAIP	TN	1	Softwood	224	19.27	20.54	0.82	3.46	−1.28
NAIP	TN	1	Hardwood	1429	24.00	25.05	0.71	3.96	−1.05
NAIP	TN	10	Softwood	224	19.95	20.81	0.81	3.43	−0.85
NAIP	TN	10	Hardwood	1429	24.23	25.02	0.70	4.12	−0.79
ALS	TN	1	Softwood	126	20.08	19.89	0.87	2.56	0.20
ALS	TN	1	Hardwood	794	24.91	24.71	0.79	3.20	0.20
NAIP	VA	1	Softwood	162	17.84	20.76	0.89	3.38	−2.92
NAIP	VA	1	Hardwood	634	23.25	25.40	0.63	4.90	−2.15
NAIP	VA	10	Softwood	462	18.32	20.40	0.77	3.70	−2.08
NAIP	VA	10	Hardwood	1447	23.34	24.69	0.58	5.30	−1.35
ALS	VA	1	Softwood	154	18.62	20.73	0.88	3.42	−2.11
ALS	VA	1	Hardwood	628	25.42	25.15	0.74	3.58	0.28

Table 4

Post-stratified estimation results.

Tennessee
Forest Volume/Area (m3/ha)
Survey unit		REFIA				Volume/Area(m3/ha) *			ApparentSample Size
	Plots	PSDHM	PSCHM	PSCHM+FT	RESRS *		SE% *	Plot Gain/Loss
T1	457	1.19	1.28	1.09	1.49	166.82	2.88	126	583
T2	412	1.30	1.34	1.20	1.58	140.06	2.48	138	550
T3	541	1.27	1.31	1.33	1.33	155.47	2.42	179	720
T4	546	1.34	1.44	1.46	1.39	168.73	2.06	251	797
T5	707	1.26	1.30	1.27	1.32	184.46	1.88	211	918
Total forest volume (million m3)
Survey unit		REFIA				Total volume(million m3) *			ApparentSample Size
	Plots	PSDHM	PSCHM	PSCHM+FT	RESRS *		SE% *	Plot Gain/Loss
T1	457	1.33	1.34	1.16	2.74	156.17	3.02	155	612
T2	412	1.45	1.46	1.27	2.51	127.30	2.60	190	602
T3	541	1.32	1.32	1.33	2.87	165.75	2.63	180	721
T4	546	1.37	1.37	1.33	2.39	203.29	2.14	203	749
T5	707	1.29	1.30	1.28	2.81	275.77	2.02	211	918
Forest area (thousand ha)
Survey unit		REFIA				Forest area(thousand ha) *			ApparentSample Size
	Plots	PSDHM	PSCHM	PSCHM+FT	RESRS *		SE% *	Plot Gain/Loss
T1 **	457	0.66	0.78	0.79	3.93	927.51	1.88	−96	361
T2	412	0.89	1.24	0.69	4.05	908.86	1.38	100	512
T3	541	0.78	0.89	1.01	3.80	1066.18	1.74	6	547
T4	546	0.83	1.11	1.06	3.65	1209.48	1.24	61	607
T5	707	0.80	1.02	0.97	3.80	1495.00	1.29	15	722
Virginia
Forest Volume/area (m3/ha)
Survey unit		REFIA
	Plots	PSDHM	PSCHM	PSCHM+FT	RESRS *	Volume/area(m3/ha) *	SE% *	Plot Gain/Loss	ApparentSample Size
V1	694	1.19	1.18	1.18	1.42	191.04	2.14	130	824
V2	694	1.28	1.22	1.24	1.54	161.97	2.21	196	890
V3	474	1.07	1.15	1.01	1.24	201.46	2.26	71	545
V4	503	1.14	1.33	1.20	1.37	171.40	2.07	167	670
V5	584	1.20	1.26	1.20	1.29	192.15	2.18	150	734
Total forest volume (million m3)
Survey unit		REFIA
	Plots	PSDHM	PSCHM	PSCHM+FT	RESRS *	Total volume(million m3) *	SE% *	Plot Gain/Loss	ApparentSample Size
V1	694	1.34	1.36	1.23	2.39	284.15	2.24	248	942
V2	694	1.36	1.38	1.23	2.33	247.98	2.31	265	959
V3	474	1.29	1.29	1.20	2.37	204.04	2.51	139	613
V4	503	1.31	1.32	1.27	2.46	197.35	2.18	163	666
V5	584	1.23	1.20	1.15	2.28	245.77	2.25	136	720
Forest area (thousand ha)
Survey unit		REFIA
	Plots	PSDHM	PSCHM	PSCHM+FT	RESRS *	Forest area(thousand ha) *	SE% *	Plot Gain/Loss	ApparentSample Size
V1 **	694	0.71	0.83	0.78	3.07	1523.16	1.77	−119	575
V2	694	0.77	1.08	0.77	2.48	1539.16	1.27	55	749
V3 **	474	0.64	0.88	0.85	3.19	1006.63	1.31	−58	416
V4 **	503	0.60	0.96	0.95	3.54	1118.82	1.32	−22	481
V5 **	584	0.71	0.87	0.79	3.70	1240.08	1.88	−73	511

Post-stratified estimation results. Shaded RE_FIA indicates the strata map with the highest variance ratio. Columns with asterisks (*) have values taken from the best post-stratification method while double asterisks (**) indicate survey units where the tree canopy cover map (PS_TCC) resulted in a lower variance than the NAIP strata maps.

Table 5

Summary of post-stratification results by survey unit.

	SurveyUnits	ImprovedSurvey Units	Best PS Map			Plot Gain/Loss (%)
			PSDHM	PSCHM	PSCHM+FT	Avg.	Min	Max
TN Forest Area (acres)	5	4	0	3	1	4	−21	24
TN Total Forest Volume (m3)	5	5	1	4	1	36	30	46
TN Total Forest Vol/Area (m3/acre)	5	5	0	3	2	34	28	46
VA Forest Area (acres)	5	1	0	1	0	−8	−17	8
VA Total Forest Volume (m3)	5	5	2	4	0	32	23	38
VA Total Forest Vol/Area (m3/acre)	5	5	2	3	0	24	15	33
Total	30	25	5	18	4	-	-	-

Improved survey unit represents the number of survey units where RE_FIA is greater than 1. Best PS map indicates the number of times each NAIP strata map was found to be the best method (note in the event of a tie units can have more than 1 best map per survey unit). Plot gain/loss column shows the average percentage of plots gained or lost relative to the initial sample size, as well as the minimum and maximum number of plots gained per survey unit in each state (note losses relative to the NAIP maps are shown as negatives).

Table 6

Post-stratified estimation results by FIA species group.

	TotalSpecies	ImprovedSpecies	Best PS Map			Softwood Plot Gain/Loss (%)				Hardwood Plot Gain/Loss (%)
			PSDHM	PSCHM	PSCHM+FT	Avg.	Min	Max	Avg.	Min	Max
TN Total Forest Volume (m3)	84	63	4	21	44	17	−2	64	4	−7	18
TN Total Forest Volume/Area (m3/acre)	84	71	23	8	44	20	2	72	6	−4	20
VA Total Forest Volume (m3)	80	53	7	19	35	10	−6	52	4	−7	28
VA Total Forest Vol/Area (m3/acre)	80	69	32	21	29	14	−5	60	7	−4	29
Total	328	256	66	69	152	-	-	-	-	-	-

Total species: the number of species groups observed in the FIA sampling plots. Improved species: the number of species groups which have RE_FIA greater than 1.00. Softwood plot gain/loss: the number of plots gained (or lost) for softwood species groups by using the NAIP strata maps instead of tree canopy cover data. Hardwood plot gain/loss: the number of plots gained (or lost) for hardwood species groups by using the NAIP strata maps instead of tree canopy cover data.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/rs14174386/s1, Table S1: The full list of species groups in eastern U.S., Table S2: Total volume and Volume/area estimates by each species group and survey unit.

References

1. FAO. Voluntary Guidelines on National Forest Monitoring; FAO: Rome, Italy, 2017; ISBN 978-92-5-109619-2

2. Nilsson, M.; Nordkvist, K.; Jonzén, J.; Lindgren, N.; Axensten, P.; Wallerman, J.; Egberth, M.; Larsson, S.; Nilsson, L.; Eriksson, J. et al. A Nationwide Forest Attribute Map of Sweden Predicted Using Airborne Laser Scanning Data and Field Data from the National Forest Inventory. Remote Sens. Environ.; 2017; 194, pp. 447-454. [DOI: https://dx.doi.org/10.1016/j.rse.2016.10.022]

3. Tomppo, E.; Gschwantner, T.; Lawrence, M.; McRoberts, R.E.; Gabler, K.; Schadauer, K.; Vidal, C.; Lanz, A.; Ståhl, G.; Cienciala, E. et al. National Forest Inventories; Tomppo, E.; Gschwantner, T.; Lawrence, M.; McRoberts, R.E. Springer: Dordrecht, The Netherlands, 2010; ISBN 978-90-481-3232-4

4. Breidenbach, J.; McRoberts, R.E.; Alberdi, I.; Antón-Fernández, C.; Tomppo, E. A Century of National Forest Inventories—Informing Past, Present and Future Decisions. For. Ecosyst.; 2021; 8, 36. [DOI: https://dx.doi.org/10.1186/s40663-021-00315-x]

5. Kangas, A.; Astrup, R.; Breidenbach, J.; Fridman, J.; Gobakken, T.; Korhonen, K.T.; Maltamo, M.; Nilsson, M.; Nord-Larsen, T.; Næsset, E. et al. Remote Sensing and Forest Inventories in Nordic Countries—Roadmap for the Future. Scand. J. For. Res.; 2018; 33, pp. 397-412. [DOI: https://dx.doi.org/10.1080/02827581.2017.1416666]

6. Mononen, L.; Auvinen, A.-P.; Ahokumpu, A.-L.; Rönkä, M.; Aarras, N.; Tolvanen, H.; Kamppinen, M.; Viirret, E.; Kumpula, T.; Vihervaara, P. National Ecosystem Service Indicators: Measures of Social–Ecological Sustainability. Ecol. Indic.; 2016; 61, pp. 27-37. [DOI: https://dx.doi.org/10.1016/j.ecolind.2015.03.041]

7. Blackard, J.A.; Patterson, P.L. National FIA Plot Intensification Procedure Report; RMRS-GTR-329 U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station: Fort Collins, CO, USA, 2014; 63.

8. Goerndt, M.E.; Monleon, V.J.; Temesgen, H. Small-Area Estimation of County-Level Forest Attributes Using Ground Data and Remote Sensed Auxiliary Information. For. Sci.; 2013; 59, pp. 536-548. [DOI: https://dx.doi.org/10.5849/forsci.12-073]

9. Ståhl, G.; Saarela, S.; Schnell, S.; Holm, S.; Breidenbach, J.; Healey, S.P.; Patterson, P.L.; Magnussen, S.; Næsset, E.; McRoberts, R.E. et al. Use of Models in Large-Area Forest Surveys: Comparing Model-Assisted, Model-Based and Hybrid Estimation. For. Ecosyst.; 2016; 3, 5. [DOI: https://dx.doi.org/10.1186/s40663-016-0064-9]

10. McRoberts, R.E.; Nelson, M.D.; Wendt, D.G. Stratified Estimation of Forest Area Using Satellite Imagery, Inventory Data, and the k-Nearest Neighbors Technique. Remote Sens. Environ.; 2002; 82, pp. 457-468. [DOI: https://dx.doi.org/10.1016/S0034-4257(02)00064-0]

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

12. Tinkham, W.T.; Mahoney, P.R.; Hudak, A.T.; Domke, G.M.; Falkowski, M.J.; Woodall, C.W.; Smith, A.M.S. Applications of the United States Forest Inventory and Analysis Dataset: A Review and Future Directions. Can. J. For. Res.; 2018; 48, pp. 1251-1268. [DOI: https://dx.doi.org/10.1139/cjfr-2018-0196]

13. Hu, H.; Wang, G.G. Changes in Forest Biomass Carbon Storage in the South Carolina Piedmont between 1936 and 2005. For. Ecol. Manag.; 2008; 255, pp. 1400-1408. [DOI: https://dx.doi.org/10.1016/j.foreco.2007.10.064]

14. Brandeis, C.; Hodges, D.G.; Poudyal, N. Forest Resource Economics in Transition: Traditional and Emerging Markets; Gen. Tech. Rep. SRS-247 U.S. Department of Agriculture Forest Service, Southern Research Station: Asheville, NC, USA, 2018; 247, pp. 1-70. [DOI: https://dx.doi.org/10.2737/srs-gtr-247]

15. Bechtold, W.A.; Ruark, G.A.; Lloyd, F.T. Changing Stand Structure and Regional Growth Reductions in Georgia’s Natural Pine Stands. For. Sci.; 1991; 37, pp. 703-717.

16. Reams, G.A. Radial Growth Trends of Loblolly Pine in the Virginia Coastal Plain. Can. J. For. Res.; 1996; 26, 9. [DOI: https://dx.doi.org/10.1139/x26-187]

17. Brady, T.J.; Monleon, V.J.; Gray, A.N. Calibrating Vascular Plant Abundance for Detecting Future Climate Changes in Oregon and Washington, USA. Ecol. Indic.; 2010; 10, pp. 657-667. [DOI: https://dx.doi.org/10.1016/j.ecolind.2009.11.004]

18. Desprez, J.; Iannone, B.V.; Yang, P.; Oswalt, C.M.; Fei, S. Northward Migration under a Changing Climate: A Case Study of Blackgum (Nyssa Sylvatica). Clim. Chang.; 2014; 126, pp. 151-162. [DOI: https://dx.doi.org/10.1007/s10584-014-1207-z]

19. Woodall, C.W.; Oswalt, C.M.; Westfall, J.A.; Perry, C.H.; Nelson, M.D.; Finley, A.O. An Indicator of Tree Migration in Forests of the Eastern United States. For. Ecol. Manag.; 2009; 257, pp. 1434-1444. [DOI: https://dx.doi.org/10.1016/j.foreco.2008.12.013]

20. Patterson, P.L.; Coulston, J.W.; Roesch, F.A.; Westfall, J.A.; Hill, A.D. A Primer for Nonresponse in the US Forest Inventory and Analysis Program. Environ. Monit. Assess.; 2012; 184, pp. 1423-1433. [DOI: https://dx.doi.org/10.1007/s10661-011-2051-5]

21. McRoberts, R.E.; Wendt, D.G.; Nelson, M.D.; Hansen, M.H. Using a Land Cover Classification Based on Satellite Imagery to Improve the Precision of Forest Inventory Area Estimates. Remote Sens. Environ.; 2002; 81, pp. 36-44. [DOI: https://dx.doi.org/10.1016/S0034-4257(01)00330-3]

22. Tomppo, E.; Haakana, M.; Katila, M.; Peräsaari, J. Multi-Source National Forest Inventory: Methods and Applications; Springer: Dordrecht, The Netherlands, 2008; ISBN 978-1-4020-8712-7

23. 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]

24. Hansen, M.H.; Wendt, D.G. Using Classified Landsat Thematic Mapper Data for Stratification in a Statewide Forest Inventory. First Annual Forest Inventory and Analysis Symposium; US Department of Agriculture, Forest Service, North Central Research Station: St. Paul, MN, USA, 2000; pp. 20-27.

25. McRoberts, R.E.; Liknes, G.C.; Domke, G.M. Using a Remote Sensing-Based, Percent Tree Cover Map to Enhance Forest Inventory Estimation. For. Ecol. Manag.; 2014; 331, pp. 12-18. [DOI: https://dx.doi.org/10.1016/j.foreco.2014.07.025]

26. Schroeder, T.A.; Healey, S.P.; Moisen, G.G.; Frescino, T.S.; Cohen, W.B.; Huang, C.; Kennedy, R.E.; Yang, Z. Improving Estimates of Forest Disturbance by Combining Observations from Landsat Time Series with U.S. Forest Service Forest Inventory and Analysis Data. Remote Sens. Environ.; 2014; 154, pp. 61-73. [DOI: https://dx.doi.org/10.1016/j.rse.2014.08.005]

27. 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]

28. 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]

29. Coulston, J.W.; Moisen, G.G.; Wilson, B.T.; Finco, M.V.; Cohen, W.B.; Brewer, C.K. Modeling Percent Tree Canopy Cover: A Pilot Study. Photogramm. Eng. Remote Sens.; 2012; 78, pp. 715-727. [DOI: https://dx.doi.org/10.14358/PERS.78.7.715]

30. Lu, D.; Chen, Q.; Wang, G.; Moran, E.; Batistella, M.; Zhang, M.; Vaglio Laurin, G.; Saah, D. Aboveground Forest Biomass Estimation with Landsat and LiDAR Data and Uncertainty Analysis of the Estimates. Int. J. For. Res.; 2012; 2012, 436537. [DOI: https://dx.doi.org/10.1155/2012/436537]

31. 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]

32. Zhu, X.; Liu, D. Improving Forest Aboveground Biomass Estimation Using Seasonal Landsat NDVI Time-Series. ISPRS J. Photogramm. Remote Sens.; 2015; 102, pp. 222-231. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2014.08.014]

33. Holmgren, J. Prediction of Tree Height, Basal Area and Stem Volume in Forest Stands Using Airborne Laser Scanning. Scand. J. For. Res.; 2004; 19, pp. 543-553. [DOI: https://dx.doi.org/10.1080/02827580410019472]

34. Næsset, E. Determination of Mean Tree Height of Forest Stands Using Airborne Laser Scanner Data. ISPRS J. Photogramm. Remote Sens.; 1997; 52, pp. 49-56. [DOI: https://dx.doi.org/10.1016/S0924-2716(97)83000-6]

35. 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]

36. Baltsavias, E.P. A Comparison between Photogrammetry and Laser Scanning. ISPRS J. Photogramm. Remote Sens.; 1999; 54, pp. 83-94. [DOI: https://dx.doi.org/10.1016/S0924-2716(99)00014-3]

37. Iglhaut, J.; Cabo, C.; Puliti, S.; Piermattei, L.; O’Connor, J.; Rosette, J. Structure from Motion Photogrammetry in Forestry: A Review. Curr. For. Rep.; 2019; 5, pp. 155-168. [DOI: https://dx.doi.org/10.1007/s40725-019-00094-3]

38. Huang, W.; Dolan, K.; Swatantran, A.; Johnson, K.; Tang, H.; O’Neil-Dunne, J.; Dubayah, R.; Hurtt, G. High-Resolution Mapping of Aboveground Biomass for Forest Carbon Monitoring System in the Tri-State Region of Maryland, Pennsylvania and Delaware, USA. Environ. Res. Lett.; 2019; 14, 095002. [DOI: https://dx.doi.org/10.1088/1748-9326/ab2917]

39. Lisein, J.; Pierrot-Deseilligny, M.; Bonnet, S.; Lejeune, P. A Photogrammetric Workflow for the Creation of a Forest Canopy Height Model from Small Unmanned Aerial System Imagery. Forests; 2013; 4, pp. 922-944. [DOI: https://dx.doi.org/10.3390/f4040922]

40. White, J.C.; Wulder, M.A.; Vastaranta, M.; Coops, N.C.; Pitt, D.; Woods, M. The Utility of Image-Based Point Clouds for Forest Inventory: A Comparison with Airborne Laser Scanning. Forests; 2013; 4, pp. 518-536. [DOI: https://dx.doi.org/10.3390/f4030518]

41. McRoberts, R.E.; Chen, Q.; Gormanson, D.D.; Walters, B.F. The Shelf-Life of Airborne Laser Scanning Data for Enhancing Forest Inventory Inferences. Remote Sens. Environ.; 2018; 206, pp. 254-259. [DOI: https://dx.doi.org/10.1016/j.rse.2017.12.017]

42. 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]

43. Vastaranta, M.; Wulder, M.A.; White, J.C.; Pekkarinen, A.; Tuominen, S.; Ginzler, C.; Kankare, V.; Holopainen, M.; Hyyppä, J.; Hyyppä, H. Airborne Laser Scanning and Digital Stereo Imagery Measures of Forest Structure: Comparative Results and Implications to Forest Mapping and Inventory Update. Can. J. Remote Sens.; 2013; 39, pp. 382-395. [DOI: https://dx.doi.org/10.5589/m13-046]

44. White, J.C.; Stepper, C.; Tompalski, P.; Coops, N.C.; Wulder, M.A. Comparing ALS and Image-Based Point Cloud Metrics and Modelled Forest Inventory Attributes in a Complex Coastal Forest Environment. Forests; 2015; 6, pp. 3704-3732. [DOI: https://dx.doi.org/10.3390/f6103704]

45. Gobakken, T.; Bollandsås, O.M.; Næsset, E. Comparing Biophysical Forest Characteristics Estimated from Photogrammetric Matching of Aerial Images and Airborne Laser Scanning Data. Scand. J. For. Res.; 2015; 30, pp. 73-86. [DOI: https://dx.doi.org/10.1080/02827581.2014.961954]

46. Ginzler, C.; Hobi, M. Countrywide Stereo-Image Matching for Updating Digital Surface Models in the Framework of the Swiss National Forest Inventory. Remote Sens.; 2015; 7, pp. 4343-4370. [DOI: https://dx.doi.org/10.3390/rs70404343]

47. Rahlf, J.; Breidenbach, J.; Solberg, S.; Næsset, E.; Astrup, R. Digital Aerial Photogrammetry Can Efficiently Support Large-Area Forest Inventories in Norway. For. Int. J. For. Res.; 2017; 90, pp. 710-718. [DOI: https://dx.doi.org/10.1093/forestry/cpx027]

48. Bohlin, J.; Wallerman, J.; Fransson, J.E.S. Forest Variable Estimation Using Photogrammetric Matching of Digital Aerial Images in Combination with a High-Resolution DEM. Scand. J. For. Res.; 2012; 27, pp. 692-699. [DOI: https://dx.doi.org/10.1080/02827581.2012.686625]

49. Irulappa-Pillai-Vijayakumar, D.B.; Renaud, J.-P.; Morneau, F.; McRoberts, R.E.; Vega, C. Increasing Precision for French Forest Inventory Estimates Using the K-NN Technique with Optical and Photogrammetric Data and Model-Assisted Estimators. Remote Sens.; 2019; 11, 991. [DOI: https://dx.doi.org/10.3390/rs11080991]

50. Jenkins, C.N.; Houtan, K.S.V.; Pimm, S.L.; Sexton, J.O. US Protected Lands Mismatch Biodiversity Priorities. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 5081-5086. [DOI: https://dx.doi.org/10.1073/pnas.1418034112] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25847995]

51. Cohen, W.B.; Yang, Z.; Stehman, S.V.; Schroeder, T.A.; Bell, D.M.; Masek, J.G.; Huang, C.; Meigs, G.W. Forest Disturbance across the Conterminous United States from 1985–2012: The Emerging Dominance of Forest Decline. For. Ecol. Manag.; 2016; 360, pp. 242-252. [DOI: https://dx.doi.org/10.1016/j.foreco.2015.10.042]

52. Schleeweis, K.G.; Moisen, G.G.; Schroeder, T.A.; Toney, C.; Freeman, E.A.; Goward, S.N.; Huang, C.; Dungan, J.L. US National Maps Attributing Forest Change: 1986–2010. Forests; 2020; 11, 653. [DOI: https://dx.doi.org/10.3390/f11060653]

53. Oswalt, C.M.; Oswalt, S.N.; Johnson, T.G.; Brandeis, C.; Randolph, K.C. Tennessee’s Forests, 2009. Resource Bulletin SRS-189; U.S. Department of Agriculture Forest Service, Southern Research Station: Asheville, NC, USA, 2012; pp. 1-165.

54. Brandeis, T.J.; Hartsell, A.J.; Randolph, K.C.; Oswalt, C.M. Virginia’s Forests, 2016. Resource Bulletin SRS–223; U.S. Department of Agriculture Forest Service, Southern Research Station: Asheville, NC, USA, 2018; pp. 1-116.

55. Webb, J.; Schaaf, A.; Davenport, J.; Fisk, H.; Strunk, J.; McGaughey, R.J.; Anderson, H.; Gatziolis, D.; Bell, D.M. NAIP Point Clouds: Evaluating and Processing Statewide High-Resolution Elevation Data; U.S. Department of Agriculture, Forest Service, Geospatial Technology and Applications Center: Salt Lake City, UT, USA, 2020; 11.

56. U.S. Geological Survey 3D Elevation Program. Available online: https://www.usgs.gov/core-science-systems/ngp/3dep (accessed on 10 October 2021).

57. Sugarbaker, L.J.; Eldridge, D.F.; Saghy, D.L.; Stoker, J.M.; Thune, D.R. Status of the 3D Elevation Program, 2015; U.S. Geological Survey: Reston, VA, USA, 2017; 13.

58. McGaughey, R.J. FUSION/LDV: Software for LIDAR Data Analysis and Visualization; United States Department of Agriculture Forest Service: Seattle, WA, USA, 2016; 217.

59. Goward, S.N.; Huang, C.; Zhao, F.; Schleeweis, K.; Rishmawi, K.; Lindsey, M.; Dungan, J.L.; Michaelis, A.R. NACP NAFD Project: Forest Disturbance History from Landsat, 1986–2010; ORNL DAAC: Oak Ridge, TN, USA, 2016; [DOI: https://dx.doi.org/10.3334/ORNLDAAC/1290]

60. USDA Forest Service Rocky Mountain Research Station North American Forest Dynamics (NAFD). Available online: https://www.fs.usda.gov/rmrs/projects/north-american-forest-dynamics-nafd (accessed on 30 January 2022).

61. Jin, S.; Homer, C.; Yang, L.; Danielson, P.; Dewitz, J.; Li, C.; Zhu, Z.; Xian, G.; Howard, D. Overall Methodology Design for the United States National Land Cover Database 2016 Products. Remote Sens.; 2019; 11, 2971. [DOI: https://dx.doi.org/10.3390/rs11242971]

62. Multi-Resolution Land Characteristics Consortium NLCD 2016 Land Cover (CONUS). Available online: https://www.mrlc.gov/data/nlcd-2016-land-cover-conus (accessed on 30 January 2022).

63. USDA Forest Service FIA DataMart. Available online: https://apps.fs.usda.gov/fia/datamart/CSV/datamart_csv.html (accessed on 30 January 2022).

64. Burrill, E.A.; Wilson, A.M.; Turner, J.A.; Pugh, S.A.; Menlove, J.; Christensen, G.; Conkling, B.L.; David, W. The Forest Inventory and Analysis Database: Database Description and User Guide for Phase 2 (Version 7.2); U.S. Forest Service: Washington, DC, USA, 2018.

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

66. Roussel, J.-R.; Auty, D.; Coops, N.C.; Tompalski, P.; Goodbody, T.R.H.; Meador, A.S.; Bourdon, J.-F.; de Boissieu, F.; Achim, A. LidR: An R Package for Analysis of Airborne Laser Scanning (ALS) Data. Remote Sens. Environ.; 2020; 251, 112061. [DOI: https://dx.doi.org/10.1016/j.rse.2020.112061]

67. Gobakken, T.; Næsset, E. Assessing Effects of Laser Point Density, Ground Sampling Intensity, and Field Sample Plot Size on Biophysical Stand Properties Derived from Airborne Laser Scanner Data. Can. J. For. Res.; 2008; 38, 5. [DOI: https://dx.doi.org/10.1139/X07-219]

68. Healey, S.P.; Cohen, W.B.; Yang, Z.; Kenneth Brewer, C.; Brooks, E.B.; Gorelick, N.; Hernandez, A.J.; Huang, C.; Joseph Hughes, M.; Kennedy, R.E. et al. Mapping Forest Change Using Stacked Generalization: An Ensemble Approach. Remote Sens. Environ.; 2018; 204, pp. 717-728. [DOI: https://dx.doi.org/10.1016/j.rse.2017.09.029]

69. Næsset, E. Predicting Forest Stand Characteristics with Airborne Scanning Laser Using a Practical Two-Stage Procedure and Field Data. Remote Sens. Environ.; 2002; 80, pp. 88-99. [DOI: https://dx.doi.org/10.1016/S0034-4257(01)00290-5]

70. Næsset, E. Airborne Laser Scanning as a Method in Operational Forest Inventory: Status of Accuracy Assessments Accomplished in Scandinavia. Scand. J. For. Res.; 2007; 22, pp. 433-442. [DOI: https://dx.doi.org/10.1080/02827580701672147]

71. Cohen, W.B.; Maiersperger, T.K.; Gower, S.T.; Turner, D.P. An Improved Strategy for Regression of Biophysical Variables and Landsat ETM+ Data. Remote Sens. Environ.; 2003; 84, pp. 561-571. [DOI: https://dx.doi.org/10.1016/S0034-4257(02)00173-6]

72. Frescino, T.S.; Patterson, P.L.; Freeman, E.A.; Moisen, G.G. An R-Based Tool for Analysts, to Look at Temporal Trends in Forest Estimates. Moving from Status to Trends; Rocky Mountain Research Station, US Forest Service: Baltimore, MD, USA, 2012.

73. Scott, C.T.; Bechtold, W.A.; Reams, G.A.; Smith, W.D.; Westfall, J.A.; Hansen, M.H.; Moisen, G.G. Sample-Based Estimators Used by the Forest Inventory and Analysis National Information Management System; Gen. Tech. Rep. SRS-80 U.S. Department of Agriculture, Forest Service, Southern Research Station: Asheville, NC, USA, 2005; pp. 53-77.

74. Andersen, H.-E.; Reutebuch, S.E.; McGaughey, R.J. A Rigorous Assessment of Tree Height Measurements Obtained Using Airborne LIDAR and Conventional Field Methods. Can. J. Remote Sens.; 2006; 32, pp. 355-366. [DOI: https://dx.doi.org/10.5589/m06-030]

75. Tesfamichael, S.G.; van Aardt, J.A.N.; Ahmed, F. Estimating Plot-Level Tree Height and Volume of Eucalyptus Grandis Plantations Using Small-Footprint, Discrete Return Lidar Data. Prog. Phys. Geogr. Earth Environ.; 2010; 34, pp. 515-540. [DOI: https://dx.doi.org/10.1177/0309133310365596]

76. Prior, E.M.; Thomas, V.A.; Wynne, R.H. Estimation of mean dominant height using NAIP digital aerial photogrammetry and lidar over mixed deciduous forest in the southeastern USA. Int. J. Appl. Earth Obs. Geoinf.; 2022; 110, 102813. [DOI: https://dx.doi.org/10.1016/j.jag.2022.102813]

77. Ritz, A.L.; Thomas, V.A.; Wynne, R.H.; Green, P.C.; Schroeder, T.A.; Albaugh, T.J.; Burkhart, H.E.; Carter, D.R.; Cook, R.L.; Compoe, O.C. et al. Assessing the utility of NAIP digital aerial photogrammetric point clouds for estimating canopy height of managed loblolly pine plantations in the southeastern United States. Int. J. Appl. Earth Obs. Geoinf.; 2022; in review

78. Sibona, E.; Vitali, A.; Meloni, F.; Caffo, L.; Dotta, A.; Lingua, E.; Motta, R.; Garbarino, M. Direct Measurement of Tree Height Provides Different Results on the Assessment of LiDAR Accuracy. Forests; 2016; 8, 7. [DOI: https://dx.doi.org/10.3390/f8010007]

79. Wang, Y.; Lehtomäki, M.; Liang, X.; Pyörälä, J.; Kukko, A.; Jaakkola, A.; Liu, J.; Feng, Z.; Chen, R.; Hyyppä, J. Is Field-Measured Tree Height as Reliable as Believed—A Comparison Study of Tree Height Estimates from Field Measurement, Airborne Laser Scanning and Terrestrial Laser Scanning in a Boreal Forest. ISPRS J. Photogramm. Remote Sens.; 2019; 147, pp. 132-145. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2018.11.008]

80. Cao, Q.; Dettmann, G.T.; Radtke, P.J.; Coulston, J.W.; Derwin, J.; Thomas, V.A.; Burkhart, H.E.; Wynne, R.H. Increased precision in county-level volume estimates in the United States Forest Inventory with area-level small area estimation. Front. For. Glob. Chang.; 2022; 5, 769917. [DOI: https://dx.doi.org/10.3389/ffgc.2022.769917]