Assessing the utility of NAIP digital aerial photogrammetric point clouds for estimating canopy height of managed loblolly pine plantations in the southeastern United States

Remote sensing offers many advantages to supplement traditional, ground-based forest measurements, such as limiting time in the field and fast spatial coverage. Data from airborne laser scanning (lidar) have provided accurate estimates of forest height, where, and when available. However, lidar is expensive to collect


Introduction
Advancements in remote sensing have grown significantly over the past few decades, allowing for new and revolutionary ways to measure forest height.Maintaining an accurate forest height model is essential for forest management (Green et al., 2020;Naesset, 2002;Navarro et al., 2020;Schultz, 1999).Light detection and ranging (lidar) has been widely used for height mapping over the past 2-3 decades; however, several impediments to the use of lidar data exist, including high costs, processing capabilities, and limited availability.Digital aerial photogrammetry (DAP) has been found to be comparable to lidar in forest height estimation when utilizing a precise digital elevation model (DEM) obtained from lidar point clouds (Bohlin et al., 2012;Goodbody et al., 2019;Lebarl et al., 2010;Strunk et al., 2020;White et al., 2015;Zimmermann and Hoffmann, 2017).
In the United States the National Agriculture Imagery Program (NAIP) regularly collects high-resolution aerial imagery for the United States Department of Agriculture's (USDA) Farm Service Agency (FSA) (Davis, 2011).For a given state, NAIP was acquired every-five years from 2003 to 2009 and is now acquired every-three years during the   growing season.NAIP imagery is collected using overlapping transects which are used to produce complete stereo-images for the entire country (Strunk et al., 2019;Strunk et al., 2020).The images are all acquired within one growing season for a state, covering an extensive area in a relatively short timeframe, so there is minimal change between images (Strunk et al., 2019(Strunk et al., & 2020)).In addition, NAIP photos can generate 3D point clouds of the surface of the canopy that are comparable to lidar 3D point clouds of the upper canopy.However, some studies have shown that NAIP imagery and point clouds can be significantly impacted by acquisition parameters including sun position, terrain, and shadow (Schroeder et al., in review;Prior et al., 2022).NAIP also cannot penetrate the forest canopy and only produces a digital surface model (DSM) of the stand.A solution to this shortcoming is the use of lidar-derived DEMs, which are considered the "best available data products" by Goodbody et al., (2019), with the NAIP point clouds to calculate a normalized canopy height model (Michez et al., 2020;Mielcarek et al., 2020;White et al., 2015;Navarro et al., 2020;Bohlin et al., 2012).After using a high-quality DEM (≤1 m ground sampling distance), R 2′ s (coefficient of determination: 1 − sum squared regression (SSR)  total sum squares (SST) ) between lidar and DAP have been reported as high as of 0.9 (Michez et al., 2020;Mielcarek et al., 2020).Lidar-based DEMs are widely available from the USGS site and since there is the assumption that topography is not significantly different over time, the DEM does not need to be coincident with the flight.
Clearly there is potential for DAP to predict pine or other forest heights in Europe and western North America.However, this potential has not been widely explored in the managed pine forests of the southeastern United States.There is a need to explore this potential for managed southern pines, which span a large geographic region (863,778 km 2 ), since they are a vital economic resource for the world and there is growing importance for forest companies to evaluate carbon stocks for certification and for the carbon market (Janowiak et al., 2017).We had two primary objectives, as follows: (1) to develop a model that accurately predicts the top-of-canopy height of managed loblolly pine plantations from NAIP photogrammetric point clouds, and (2) to use the resulting model to create a pine canopy height map for all of Virginia, North Carolina, and Tennessee.

Study area and field measurements
The overall study area includes thinned and non-thinned evergreen forests of Virginia, North Carolina, and Tennessee, based on the National Land Cover Dataset.For this study, stands that have had rows removed and the bare ground is visible between the rows that remain are considered thinned.Areas where no thinning has been done or the tree crowns have covered back up the ground below are considered nonthinned.Field measured loblolly pine (Pinus taeda L.) heights collected on 534 plots across seven regions in North Carolina and Virginia were used for model calibration.All these regions are in humid temperate climates with USDA plant hardiness zones from 6b to 8a (USDA Agriculture Research Service 2012).The study area ranges from coastal to mountainous, representing both the natural range and non-native areas in which loblolly pine is commonly planted (Table 1).
The seven training studies are distributed over Virginia and North Carolina with three in Virginia state forests, one at the Reynolds Homestead in southwest Virginia, and three in North Carolina (Fig. 1).The three Virginia state forests are Appomattox-Buckingham (ABSF), Cumberland (CUSF), and Prince Edward (PESF), all of which are in the Piedmont physiographic region.These three state forests are managed by the Virginia Department of Forestry and maintain numerous actively managed loblolly pine stands.Plots at ABSF, CUSF, and PRSF were established with a Trimble Geo 7x GPS and were measured in 2018-2019 as part of a study conducted by the Forest Modeling Research Cooperative (FMRC; Green et al., 2020).The data for Patrick, Bladen, Jones, and Brunswick counties were collected as part of other research studies in 2019 (Albaugh et al., 2018;Grover et al., 2020).The stands in Patrick County, VA and Bladen County, NC were planted in 2009.The remaining two training studies were in Jones and Brunswick counties in North Carolina and were planted in 2017.These two studies were each comprised of two plots that were broken into six subplots with 12 plot measurements each.Since the training studies in Jones and Brunswick counties were planted in 2017 and measured in 2018, the trees were short and relatively homogeneous in height.A hypsometer was used to measure tree heights on the older plots and a height pole was used for the young plots.The mean dominant height for each plot was used as the field height for this study.
Forest Inventory and Analysis (FIA) data was used to compare the state-wide model height values to height values collected in the field by the USDA Forest Service.The FIA data repository includes information on forest area and location, species, size, tree health, tree growth, harvests, and many other forest characteristics (USDA Forest Service, 2020).Each year, a fixed-area plot design is used to measure a sample of the over 400,000 permanent plots across the United States (Burkhart et al., 2019).These measurements are recorded across both public and private lands that are classified as forestland use or non-forestland use.

Lidar and DEM
Lidar data were used in this study for (1) DEMs with fine spatial-scale accuracy, needed for normalizing the NAIP point clouds to height above ground, and (2) to estimate 2018 tree heights.Lidar point clouds and associated 1-m DEMs were acquired from publicly available data from the United States Geological Survey (USGS) 3D Elevation Program, which does not cover all the state of Virginia but was in our regions of interest (e.g., USGS 3DEP; USGS, 2021).Since the terrain is unlikely to significantly change over the timespan of interest, older DEMs were used for the plots that did not have recent lidar acquisitions.Patrick County was part of a 2018 lidar campaign, ABSF and CUSF were part of a 2015 lidar campaign, and PESF was from a 2014 campaign (Table 2, (USDA-FPAC-BC-APFO; Green et al., 2020;US Geological, 2019)).However, for the three stands in North Carolina, there were no lidar data available within five years of the field measurements so lidar analysis was only done for the Virginia stands and NAIP analysis was done for both North Carolina and Virginia.Lidar data was accessed form the USGS National Map (USGS, 2017a, 2017b).The DEM for North Carolina was from 3DEP data that was collected in 2013.

DAP canopy height processing
Data analysis began by clipping the NAIP and lidar point clouds and the DEM to the extent of the plots in all three state forests (lasclip).The data were normalized by subtracting the DEM from both the NAIP and lidar point clouds to create a point cloud of heights above the ground (normalize_height) for each plot in the training area.
When comparing various percentiles to find the best fitting model, the 50th, 75th, 85th, 90th, and 95th percentiles were evaluated.These percentiles were chosen because of their common use as height metrics in literature and because they were tested with the state forests data in a preliminary trial of this study (Li et al. 2016;Maimaitijiang et al., 2019Maimaitijiang et al., & 2020)).Due to having the highest R 2 value of the tested percentiles, the 90th percentile was used to train the model to predict the canopy height, which follows the results of previous studies (Naesset & Bjerknes, 2001).After the data were clipped and normalized, the 90th percentile of height for each training plot was calculated to ensure that the top of the canopy was represented (cloud_metrics).
After the lidar data were plotted, it was clear that the data could not serve as a direct comparison in this study due to the time between the acquisition date and the field measurement date (Table 2), but the data did help to highlight the issues associated with NAIP when measuring thinned stands.Lidar data for the four stands for which it was available were still analyzed, but only for general comparison with the NAIP data and its accuracy in estimating height.The lidar data were not used in the model analysis.The 534 mean dominant field heights and respective NAIP 90th percentile heights were then used in a reduced major axis (RMA) (also known as standard major axis (SMA)) regression to develop a model to predict field height from NAIP height using the linear model (lmodel2 and summary) functions in R. RMA was selected for this study because it is assumed that there is error in the NAIP data as well as in the field measured data.The following model was fit (Equation ( 1)) after several models of differing percentile values were assessed and compared using R 2 , adjusted R 2 , predicted residual error sum of squares (PRESS), and RMSE.
Predicted height of the canopy (m) = m*(NAIP 90th percentile of height) + b, Where m (or the slope) is the proportional difference between the NAIP and field measured heights, and b (or the intercept) is the height difference between NAIP and ground measured heights (m).To determine if the canopy heights measured on the 534 field plots were representative of pine canopy heights across the broader three state study area, tree data collected by the U.S. Forest Service, FIA program data were used to develop histograms for comparison with the data collected in this study.This helped to ensure that the training data captured the full range of loblolly pine heights across its natural range for the age span of this study.Thomas et al., (2019) found poor correlation between measured and NAIP data from recently thinned stands.Consequently, we examined our data for outliers where the NAIP estimate was greater than one standard deviation from the field measured height.Six plots met this criterion, and all were dropped form the dataset after examining the NAIP imagery (Fig. 2).

Appomattox-Buckingham state forests
The model was initially applied to two sites in the ABSF, one without thinning and one with a recently thinned area with some of the identified outliers (Fig. 2).These two sites were used to test the mapping approach and to ensure spatial coherence of the height prediction model across the landscape.The stands were measured on a 5 m × 5 m grid to demonstrate the application of the height prediction model (Table 4: grid_metrics) so that at least one NAIP point (or portion thereof) would be in a typical grid cell.Following the previous steps for normalization, the 90th percentile of height for the NAIP points in each grid cell was calculated and then incorporated into the height prediction model to calculate the predicted 90th percentile of height of the grid cell.

Multistate
After analyzing the canopy height maps for ABSF, the resulting calibration model was applied to all areas likely to contain pines in Virginia, North Carolina, and Tennessee (shown in green, Fig. 1).Because of the large number of point clouds involved, each point cloud's 90th percentile of height was converted into a raster on a 5 m × 5 m grid metric.The calibration model was applied to the area of each raster that fell within the NLCD evergreen forest class shapefile.Once this process was completed a map was produced with predicted height values for all evergreen forests across the three states.

Lidar
The lidar data in this study were inadequate for direct comparison to measured forest height due to the difference in the time of acquisition.In Fig. 3 the lidar measured height falls much lower than the field height (almost a four-year height difference).However, there is still a clear linear relationship between the lidar measured heights and the field measured heights.The minimum (min), mean, and median for the lidar measurements are all lower than the field measurements, as expected, but the maximum (max) lidar is slightly larger than the field max (Table 3).

Naip canopy height processing
The field data was well correlated (R 2 of 0.83) with the NAIP 90th percentile of height values (Fig. 4).The NAIP 90th percentile is strongly related to field-measured tree heights in non-thinned areas, however, the NAIP 90th percentile often underestimates the height of trees in thinned stands.
The model created for the NAIP point cloud demonstrated a strong correlation with field heights, however, there is a large bias for the plots where there has been a recent thinning and the NAIP 90th percentile of height is much less than the mean dominant field height.The intercept and coefficient values for the dataset produced the fitted model (1) used in Fig. 4. The fitted model, shown in red, is as follows: Where X is the NAIP 90th percentile of height.The 95 % confidence interval for the coefficient of X is 0.86 -0.90, and the 95 % confidence interval for the intercept is 0.55 -1.07.Fig. 5 shows the distribution of residuals after applying the PHP model to the training dataset.The PHP model has a correlation coefficient of 0.96, an R 2 of 0.92, and an RMSE of 0.81.
The residual versus fitted plot for the NAIP 90th percentile predicted heights of the training data show a left-to-right increasing trend of the residuals at the 1 to 2 m heights and at the 10 to 15 m heights (Fig. 5).However, the bias of the thinned plots being underestimated is very apparent here.The data are also not well distributed across the range of heights as seen by the separation of the residual points from 3 to 7 m and by the un-equal distribution of the thinned and non-thinned points from 7 to 25 m.
The predicted canopy height models, created by the PHP model, are shown in Fig. 6a and 6b for one of the three non-thinned stands and for the thinned stand, respectively.The point density per grid cell for these stands is shown in Fig. 7a and 7b, respectively, demonstrating what appears to be an equal distribution and count of DAP points across the two stands.Fig. 8a -8d demonstrate the variation in the canopy height model accuracy across the thinned stands where some thinned rows have weak height estimation (Fig. 8a & 8b), and others have a more accurate height estimation (Fig. 8c & 8d).

Field, NAIP, and FIA data comparison
The distribution of the data is similar for the Field, NAIP, and FIA data, but the plots used for the model did not capture the full range of heights compared to the range collected by the FIA.After applying the model to the three states, however, the resulting histogram of heights matches the FIA histogram very well (Fig. 9).
Both histograms show a bell like curve with a peak at around 15 m.A smaller peak can also be seen around the 4 m area in the FIA sample and around the 2 m area for the three states.Statistical summaries in Table 4 show the distribution of the three data sources and further demonstrate the normalization of the datasets.FIA data for the three states is shown as a comparison between the state level data and our plots.

Naip canopy height mapping
The PHP equation (1) was applied to evergreen forests across the states of Virginia, North Carolina, and Tennessee, but, due to the issues of calculating a reliable CHM with thinned stands, it is most accurate for the stands that are not thinned (Fig. 10) (Ritz et al., 2022).

Discussion
In this study, NAIP was evaluated for its ability to accurately predict the canopy height of loblolly pine stands in seven training studies throughout Virginia and North Carolina.The point clouds were tested against the mean dominant height measurements to determine the error associated with this remote sensing method.The reduced major axis regression model (PHP) had an R 2 value aligning with the results of Mielcarek et al., (2020) and Maimaitijiang et al., (2020).The strong relationship confirmed that NAIP provides accurate measurements of forest height, as expected from previous studies (Kim et al., 2020;Prior et al., 2022;Strunk et al., 2020).The PHP model was then used to calibrate the NAIP CHM's, so the resulting maps were in better alignment with field measured tree heights.The success of this approach verified the PHP model's ability to capture the predicted height of a pine stand and create an effective canopy height model.The PHP model was then applied to Virginia, North Carolina, and Tennessee where DEMs from lidar were available.Not only is NAIP wall-to-wall coverage, but it can be acquired for 1/3 of the costs of lidar with routine data acquisitions (Goodbody et al., 2019;Michez et al., 2020;White et al., 2015;Navarro et al., 2020).Lidar is known for being highly accurate (Strunk et al., 2020;Bohlin et al., 2012;Masek et al., 2015;Noordermeer et al., 2021;Nelson et al., 2007).However, no lidar data were acquired at or near the time of field measurements.As such, this study could not make a direct comparison in height performance between DAP and lidar.The lack of recent lidar data further supports NAIP as a promising alternative to lidar for height measurements because of its consistent data collection schedules (Strunk et al., 2019;Strunk et al., 2020).
Thinned stands were a challenge in this study as NAIP did not perform well over some of these stands, causing six points to be flagged as outliers.Upon further investigation, it was found that the NAIP point cloud 90th percentile of height in these locations was not capturing the true top of canopy, and in some cases missing it entirely.This issue with the thinned stands was further investigated and found to be present in all the thinned stands in the training dataset.However, none of the plots in these stands were flagged as outliers and that is most likely due to the plot centers being in areas where the NAIP point cloud calculated the height more accurately.This issue was also present in thinned stands outside of the training region.The orientation of the thin rows relative to the sun and time of acquisition is potentially the cause of the significant impact of NAIP image shadows and subsequent identification of the top of the crowns, but further investigation is required to determine the true cause of this issue.This misalignment has led to a failure of identifying the high frequency changes in height from the ground to the top of canopy as a result of the discrepancies associated between autocorrelation and stereo-images (Gruen, 2012).
One limitation to this study is the lack of a well vetted, region-wide loblolly pine map for the southeastern United States within which to apply the model.Although the NLCD evergreen class will contain most of the pines, there will also be other trees species that are not as applicable to the model.Another limitation of this study is that the training data do not capture the full range of heights for loblolly pine in the southeastern United States.Our field plots do not capture the maximum height pines can reach in the three states, with heights in the 30 -40 m range according to the FIA data.Because of this limitation, maps of pine canopy height in the three states should include a cautionary note, particularly for the highly productive areas with taller pine trees.
Future work will include fine-tuning this model and creating canopy height models for the remaining states in the southeastern United States as well as utilizing NAIP imagery to remove non-tree areas and increase point cloud accuracy in thinned stands.In addition, with the repeated acquisition of NAIP imagery every-two to three years, there is the possibility of monitoring growth of pine stands, which may be interesting to managers and important to study longer term changes in productivity due to management, disturbance, stress, and increased atmospheric carbon dioxide.Other applications that build on this study could include creating canopy height models for deciduous tree species and modeling stand productivity over time.

Conclusion
In this study, a pine height prediction model was created that adjusts NAIP derived CHM's to better align with mean dominant tree heights of loblolly pine found in the southeastern United States.This study validated that, for areas of non-thinned loblolly pine, NAIP point clouds can be used to produce reliable predictions of canopy height.In addition, we also demonstrate how the pine height prediction model can be used to adjust NAIP CHMs to map pine canopy height across large multi-state areas.This study also shows that NAIP holds promise for mapping pine canopy height in other areas of the United States' southeastern pine region to create a series of canopy height models for loblolly pine plantations that can be updated routinely due to NAIP's repeatable acquisition cycles.The performance of NAIP in this study presents a strong argument for its use in addition to lidar in loblolly pine plantation management across the southeast United States.

Fig. 1 .
Fig. 1.Virginia and North Carolina training studies and the extent of statewide modeling across Virginia, North Carolina, and Tennessee with landcover class code 42 (evergreen forest) from National Landcover Dataset (NLCD) 2016.The 'Loblolly Pine Natural Range' in this figure is the geographic range of loblolly pine for which our calibration model is most applicable (Little 1971).The counties identified are the locations of the training data plots.

Fig. 2 .
Fig. 2. A sample of flagged NAIP points that were visually assessed as having recent harvests or thinning in Appomattox-Buckingham and Cumberland state forests where the black points represent five of the six points that were flagged, and the white points represent the points that were not flagged.The red outline is the extent of the stand.Fig. 2b depicts a zoom-in of one of the black points (point A) and Fig. 2c depicts a zoom-in of one of the white points (point B).It is important to note the shadowed, dark ground making up most of point A compared to mostly light green colored vegetation at point B. The sixth flagged point was in a single plot in CUSF.Coordinates are in projection GRS80 UTM 18 N. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Model of the lidar 90th percentile of height versus the field measured height.The 1:1 line is represented in black.The years noted in the legend are the years the lidar data were acquired.Field measurements were from 2018 to 2019.

Fig. 4 .Fig. 5 .
Fig.4.Model of the NAIP 90th percentile of height versus the mean dominant field height, where the yellow points represent plots that were not recently thinned, and the green points represent plots that had experienced a recent thinning.The 1:1 line is shown in black, and the fitted model (1) is shown by the dashed red line.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Canopy height models (CHM) after running the normalized heights through the PHP model, where (a) non-thinned stand in Appomattox-Buckingham state forest and (b) recently thinned stand in Appomattox-Buckingham state forest.The boxes represent zoom-ins of the boxed areas that are shown in Fig. 8. X and Y axis are coordinates in projection GRS80 UTM 18 N.

Fig. 7 .
Fig. 7. Number of DAP points per 5 m × 5 m cells, where (a) non-thinned stand in Appomattox-Buckingham state forest and (b) recently thinned stand in Appomattox-Buckingham state forest.The boxes represent zoom-ins of the boxed areas that are shown in Fig. 8. X and Y axis are coordinates in projection GRS80 UTM 18 N.

Fig. 8 .
Fig. 8.A zoomed-in comparison of the effects of shadow on height estimation of recently thinned stands in Appomattox-Buckingham state forest, where (a) heavy shadowed and (c) non-heavy shadowed are the canopy height models (CHMs) and (b)heavy shadowed and (d) non-heavy shadowed are the NAIP imagery used to create the photogrammetric point cloud.The heavily shadowed imagery produces an underestimated CHM.

Fig. 9 .
Fig. 9. Histograms of measured tree heights where (a) is the most recent distribution of loblolly tree heights measured on Forest Inventory and Analysis plots in VA, NC, and TN, and (b) is the distribution of the 90th percentile of evergreen tree heights calculated in the 3-state-wide canopy height model using the pine height prediction model on a 5 m × 5 m grid (1).

Fig. 10 .
Fig. 10.The state-wide CHM using PHP equation (1) across the evergreen landcover, where (a) shows the coverage in Virginia, North Carolina, and Tennessee and (b) shows the coverage at a much smaller scale.The red star in Fig. 10(a) depicts the location of Fig. 10(b).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Location and description of the field plots.
+ Total number of plots.++ Total number of thinned plots.*Lidar data available.

Table 2
Lidar collection specifications for the four Virginia training regions.

Table 3
Summary statistics of the data in the lidar and field plots.

Table 4
Summary statistics of height for the data in the training regions and the FIA data.