Evaluating the sensitivity of forest structural diversity characterization to LiDAR point density

Abstract

Recent expansion in data sharing has created unprecedented opportunities to explore structure–function linkages in ecosystems across spatial and temporal scales. However, characteristics of the same data product, such as resolution, can change over time or spatial locations, as protocols are adapted to new technology or conditions, which may impact the data's potential utility and accuracy for addressing end user scientific questions. The National Ecological Observatory Network (NEON) provides data products for users from 81 sites and over a planned 30‐year time frame, including discrete‐return light detection and ranging (LiDAR) from an airborne observation platform. LiDAR is a well‐established and increasingly available remote sensing technology for measuring three‐dimensional characteristics of ecosystem and landscape structure, including forest structural diversity. The LiDAR product that NEON provides can vary in point density from 2 to 25+ pt/m² depending on the instrument and acquisition date. We used NEON LiDAR from five forested sites to (1) identify the minimum point density at which structural diversity metrics can be robustly estimated across forested sites from different ecoclimatic zones in the United States and (2) to test the effects of variable point density on the estimation of a suite of structural diversity metrics and multivariate structural complexity types within and across forested sites. Twelve of 16 structural diversity metrics were sensitive to LiDAR point density in at least one of the five NEON forested sites. The minimum point density to reliably estimate the metrics ranged from 2.0 to 7.5 pt/m², but our results indicate that point densities above 7–8 pt/m² should provide robust measurements of structural diversity in forests for temporal or spatial comparisons. The delineation of multivariate structural complexity types from a suite of 16 structural diversity metrics was robust within sites and across forest types for a LiDAR point density of 4 pt/m² and above. This study shows that different metrics of structural diversity can vary in their sensitivity to the resolution of LiDAR data and that users of these open‐source data products should consider the point density of their data and use caution in metric selection when making spatial or temporal comparisons from these datasets.

INTRODUCTION

Recent efforts that make ecological data freely available have created unprecedented opportunities to understand ecosystems across large spatial and temporal scales (Balch et al., 2020; Gewin, 2016; Nagy et al., 2021), but nonstandard data resolution poses a challenge to using these datasets (Michener, 2015). Ecological data sharing networks, such as the National Ecological Observatory Network (NEON), provide a source of large-scale data to anyone with Internet access. Specifically, NEON offers a suite of biological, environmental, and remote sensing data products for a 30-year operational period (2018–2048) within the United States, spanning 81 sites and 20 ecoclimatic domains (Metzger et al., 2019). Despite protocols that were initially standardized for data collection at each site (e.g., plot-level forest woody structure and landscape-level aerial LiDAR), rapid technological advances lead to increased data quality and resolution over time (Ouma, 2016). While these changing data resolutions create benefits, they also pose potential challenges for making comparisons with ecological data across space and time (Reichman et al., 2011; Zipkin et al., 2021). Furthermore, given the rising popularity of open-source data networks, nonspecialists are increasingly accessing these data products and might be unaware of technical considerations for data use and application (Huang et al., 2019; McCord et al., 2021a, 2021b). End users of these data streams require information and analyses focused on the impacts of changing data resolution on research outcomes and the utility of the data for answering scientific questions.

Discrete-return LiDAR is one of the most widely available remote sensing technologies for measuring three-dimensional (3D) aspects of ecosystem and landscape structure (Guo et al., 2021; Lefsky et al., 2002). The NEON airborne observation platform (AOP) is equipped with aerial discrete-return LiDAR sensors, which collect data across the entire spatial footprint of each of the 81 sites that make up NEON (Krause & Goulden, 2015). NEON AOP flies up to three remote sensing payloads that collect LiDAR data each year. Payloads 1 and 2 included Optech ALTM Gemini LiDAR instruments beginning with early test flights in 2012–2015 and supporting operational flights from 2016 to the present. This sensor has a maximum pulse repetition frequency (PRF) of 100 kHz and provides an average LiDAR point density of ~4 pulses/m² (Krause & Goulden, 2015; NEON, 2021b). In 2018, Payload 3 debuted in AOP operations with a RIEGL LMS-Q780, with a PRF of up to 400 kHz and a higher return point density of double or more than that of the older instrument (NEON, 2021b). In 2021, Payload 1 was upgraded with an Optech Galaxy LiDAR, which has a PRF up to 1 MHz. Taking all of the sensor differences into consideration (see Appendix S1 for a detailed technical description of the payload specs that impact point density), older AOP LiDAR data collected with the Optech Gemini will likely have lower point densities influenced more by the wider outgoing pulse width and lower sensitivity. Therefore, the higher resolution LiDAR data products are only available for a subset of sites from 2018 and onward, since the Optech Gemini sensors continue to be used at the remaining sites (see Table 1 for payloads used at our study sites). Previous research has shown that NEON LiDAR data are useful for estimating the stand-level structural diversity of forest ecosystems (i.e., 3D arrangement of vegetation within ecosystems) (LaRue et al., 2019, 2020). However, efforts to characterize ecological patterns across space and time using NEON LiDAR data may be affected by variability in data resolution from these two sensors.

TABLE 1 Airborne observation platform LiDAR payload used at National Ecological Observatory Network study sites each year.

Year	ABBY	GRSM	STEI	UKFS	UNDE
2015	…	Payload1 Optech Gemini	…	…	…
2016	…	Payload1 Optech Gemini	Payload1 Optech Gemini	Payload1 Optech Gemini	Payload1 Optech Gemini
2017	Payload2 Optech Gemini	Payload1 Optech Gemini	Payload1 Optech Gemini	Payload1 Optech Gemini	Payload1 Optech Gemini
2018	Payload3 RIEGL LMS-Q780	Payload3 RIEGL LMS-Q780	…	Payload2 Optech Gemini	…
2019	Payload3 RIEGL LMS-Q780	…	Payload3 RIEGL LMS-Q780	Payload3 RIEGL LMS-Q780	Payload3 RIEGL LMS-Q780
2020	…	…	Payload3 RIEGL LMS-Q780	Payload3 RIEGL LMS-Q780	Payload3 RIEGL LMS-Q780
2021	Payload1 Optech Galaxy	Payload1 Optech Galaxy	…	…	…

LiDAR resolution can influence the ability to resolve fine-scale structural attributes of ecosystems (Pearse et al., 2019; Roussel et al., 2017; Yu et al., 2020), because low-density point clouds have reduced potential to resolve the spatial positioning of ecosystem components (Strunk et al., 2012). For example, the use of lower resolution LiDAR data can lead to discrepancies in the mean and maximum canopy height by at least a meter (Roussel et al., 2017) and other vertical aspects of forested (Roussel et al., 2017; Wilkes et al., 2015; Yao et al., 2014; Yu et al., 2020) and non-forested land cover (Balsa-Barreiro & Lerma, 2014). This presents a challenge when measuring ecosystem structural change through time so that measured values reflect true changes (e.g., geographic variation; Treitz et al., 2012; Wilkes et al., 2015) rather than noise associated with data resolution (Ruiz et al., 2014; Singh et al., 2015; Yao et al., 2014). Indeed, previous work has shown that LiDAR point density can influence values of forest structure metrics that represent cover, density, or volume of forest stands (de Almeida, Stark, Shao, et al., 2019; González-Ferreiro et al., 2012; Jakubowski et al., 2013; Kamoske et al., 2021).

Past work addressing the impact of LiDAR point-cloud resolution for characterizing structural metrics in forests has largely focused on tree height or stem area for forestry applications (González-Ferreiro et al., 2012; Jakubowski et al., 2013; Roussel et al., 2017; Strunk et al., 2012; Treitz et al., 2012). Structural diversity is a promising tool for monitoring biodiversity and ecosystem function (Gough et al., 2019; Hakkenberg et al., 2016; LaRue et al., 2019), but consequently, there is a gap in our understanding of how varying LiDAR point-cloud resolutions impact our ability to describe structural diversity metrics across forest types. Identifying the minimum point densities needed to reliably characterize univariate (Atkins et al., 2018; LaRue et al., 2020) or multivariate (Fahey et al., 2019) suites of structural diversity metrics could help facilitate inter-site and interannual comparisons using data of varying resolution and ensure repeatability in ecological monitoring efforts. Assessing the impact of point density on the characterization of a range of metrics is essential, because metrics are sensitive to different structural features within the forest canopy and thus may be differentially affected by variation in point densities. Furthermore, the degree of sensitivity might be influenced by forest type. Our goal was to use NEON LiDAR data to (1) identify the minimum point density at which structural diversity metrics can be robustly estimated across different forested sites in the United States and (2) to evaluate the effect of point density on the delineation of multivariate forest structural complexity types.

MATERIALS AND METHODS

Structural diversity metrics from LiDAR

For this study, we chose NEON AOP LiDAR data collected with a high enough point density to allow for simulated point-cloud thinning to produce a range of data resolution equivalent to that observed across the NEON sites and across aerial LiDAR acquisitions common in the literature. These data were accessed from the spatial extent of nine (distributed) base plots at each of five different NEON sites (Table 2, Figure 1) and had a point density of at least 25 pt/m² within the 40 × 40 m plot area. Within the NEON sampling design, a series of plots, called base and tower plots, are sampled periodically for woody vegetation structure and other attributes. We focused our analysis on a subset of base and tower plots to capture variability across each site and extracted LiDAR from plots that met the minimum data resolution for this study. We randomly thinned the LiDAR data at these 45 plots to eight different point densities ranging from 2 to 25 pt/m² (Figure 2) and used the thinned data to estimate 16 structural diversity metrics (Table 3). Years of LiDAR data collection for sites included Abby Road (ABBY) 2018, Great Smoky Mountains (GRSM) 2018, Steigerwaldt–Chequamegon (STEI) 2019, University of Kansas Field Station (UKFS) 2019, and University of Notre Dame Environmental Research Center (UNDE) 2020.

TABLE 2 Attributes of the forested National Ecological Observatory Network sites.

Attribute	ABBY	GRSM	STEI	UKFS	UNDE
Domain	Pacific Northwest	Appalachians and Cumberland Plateau	Great Lakes	Prairie Peninsula	Great Lakes
Latitude, longitude	45.762439, −122.330317	35.68896, −83.50195	45.50894, −89.58637	39.040431, −95.19215	46.23391, −89.537254
Forest type	Evergreen	Low elevation—deciduous, high elevation—evergreen	Mixed hardwoods	Mixed hardwoods	Second-growth northern mesic
Mean canopy height (m)	34.0	30.0	20.0	19.0	24.0
Average percent canopy cover	91.7	97.5	97.4	98.6	96.8
Mean annual temperature (°C)	10.0	13.1	4.8	12.7	4.3

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

TABLE 3 Stand-level forest structural diversity metrics estimated from National Ecological Observatory Network aerial LiDAR.

Category	Metric	Description	Reference	Point-cloud summarization type
Height	MOCH	Mean of maximum height (in meters) in 1-m2 grid of plot	Atkins et al. (2018)	Three-dimensional (3D) pixels (voxels)
Q25	25th quantile of height (in meters)	Roussel et al. (2020)	Entire point cloud
Q50	50th quantile of height (in meters)	Roussel et al. (2020)	Entire point cloud
Q75	75th quantile of height (in meters)	Roussel et al. (2020)	Entire point cloud
Q100	Maximum canopy height (in meters) (100th quantile of height)	Roussel et al. (2020)	Entire point cloud
Density	VAI	Sum of the 1-m horizontal slices (starting 0.5 m above the ground) of leaf area density with the plot (in square meters per cubic meter)	LaRue et al. (2020)	Horizontal slices
Openness	DGF	Fraction of 1-m2 canopy gaps in the plot	Atkins et al. (2018)	3D pixels (voxels)
GFP	Distribution of gaps in the point cloud 0.5 m above the ground	Bouvier et al. (2015)	3D pixels (voxels)
External heterogeneity	Rumple	Ratio of area of canopy surface relative to plot area	Jenness (2004)	3D pixels (voxels)
Top rugosity	SD (in meters) of outer canopy heights in 1-m2 plot	LaRue et al. (2020)	3D pixels (voxels)
Internal heterogeneity	SD(ht)	SD of heights within the plot (in meters)	Roussel et al. (2020)	Entire point cloud
SD(SD(ht))	Plot-level SD of the SD of heights within 9-m2 voxels in the plot area (in meters)	LaRue et al. (2020)	3D pixels (voxels)
VCI	Normalization of diversity and evenness (entropy) of 1-m-height bins within the plot	van Ewijk et al. (2011)	Horizontal slices
CV(ht)	Coefficient of variation of heights within the plot	Roussel et al. (2020)	Entire point cloud
FHD	Foliage height diversity	MacArthur and MacArthur (1961)	Horizontal slices
Gini index	Gini coefficient index	Valbuena et al. (2017)	Entire point cloud

Of the NEON AOP LiDAR data products available, we obtained the Level 1 discrete-return LiDAR (product number DP1.30003.001) (NEON, 2021a), corrected each plot's point cloud for elevation, and randomly thinned the plots to the eight different point densities. The 1-km² tiles of LiDAR that are provided by NEON were downloaded for the selected plots using the neonUtilities R package (NEON, 2020). These data are collected during peak growing-season greenness (Krause & Goulden, 2015) and are collected in a standardized manner across sites. Different LiDAR sensors can result in different vertical transmission probabilities (Shao et al., 2019) when they vary in by their optics, scan angles, and pulse footprints; however, we focused on the high-density datasets collected by the most recent NEON sensor, so our data thus have a uniform LiDAR system setup. The data and specific details about data collection and sensor specifications are found on the NEON Data Portal Website (). Following data acquisition, an 80 × 80 m buffer area was clipped around base plot centroids and each plot was visually checked to ensure that there were no large gaps in the LiDAR data coverage. Clusters of atmospheric and ground outliers were filtered by removing points six SDs above and below the mean height. Isolated point outliers were then identified and filtered using the classify_noise function with the isolated voxel filter algorithm (points that had zero neighboring points in a 3 × 3 × 3 window) in the lidR version 3.1.2 R package (Roussel & Auty, 2018). The 80 × 80 m buffer area was then corrected for elevation using a Delaunay triangulation before being clipped to the 40 × 40 m base plot area. Each plot was then visually inspected to ensure that outliers were successfully removed. Finally, the processed point cloud for each plot was randomly thinned (i.e., points randomly selected until the specified point density was reached) from the original processed point cloud to eight different point densities (Figure 2): 2, 4, 6, 8, 10, 15, 20, and 25 pt/m². The random thinning process was repeated for each point density and plot five times to mimic the random locations that the laser pulse would hit objects during flyover. All returns were used, including in the subsequent calculation of metrics.

Sixteen structural diversity metrics from five categories that describe the height, cover and openness, density, internal heterogeneity, and external heterogeneity of vegetation in the forest canopy were calculated (Table 3) (patterned after LaRue et al., 2020). All metrics, unless specified otherwise, were calculated using functions from the lidR version 1.3.2 R package (Roussel & Auty, 2018). Rumple (rumple_index function) and deep gap fraction (DGF) were calculated from a 1-m² grid canopy height model from all returns (grid_canopy function with the points to raster algorithm that takes the highest point in each cell of the raster) and for which the points below 3 m were converted to a value of zero (i.e., DGF was the fraction of all 1-m² cells below the 3 m height threshold). Mean outer canopy height (MOCH) and top rugosity were calculated from a second 1-m² grid canopy height model for which the points below 0.5 m were removed (ground points). After filtering out points below 0.5 m, we calculated the following metrics from the point cloud using the leafR package and functions listed in parentheses (de Almeida, Stark, Silva, et al., 2019): height quantiles, SD(ht), and CV(ht) (preceding functions all cloud_metrics function); SD(SD(ht)) (grid_metrics function); GFP (gap_fraction_profile function); vegetation area index (VAI) (LAD function); and vegetation complexity index (VCI) (VCI function). Foliage height diversity (FHD) (FHD function) and the Gini index (GC function) were estimated with the leafR package (de Almeida, Stark, Silva, et al., 2019).

Analyses

We used a segmented regression analysis to determine whether there was a change-point in the structural diversity metric values with increasing LiDAR point density at each plot. The presence of a difference in the slope of the line before and after the break indicated that the value of metrics varied across LiDAR point densities. Each segmented regression analysis included the eight-point densities for a specific plot and the five replicates for each density (N_{replicates/density} = 5 for a total of N_points = 40). When structural diversity metrics do not change as point density increases, the relationship is flat and produces no change-point. When the relationship is nonlinear, we assume that the change-point indicates when the metric is stable with respect to changes in point density. Assuming two segments, we used a pseudo-score statistic test to determine when segmented regression was needed with the segmented R package (Muggeo, 2008). The p value for each pseudo-score test was corrected for multiple comparisons using the false discovery rate (N_tests = 720). If the pseudo-score statistic test had an α < 0.05 after correcting for multiple comparisons, then we employed the segmented regression analysis to determine the density value of the two-segment change-point in the linear slope. We used an initial change-point estimate of 8 pt/m². The mean and SD of the change-point in the linear slope was calculated for each of the nine plots across each site to summarize the average change-point at different forested sites and the variation of change-points in plots within sites (Table 4). If the pseudo-score statistic test determined that segmented regression was not needed (insignificant), then a value of 2.0 pt/m² was assigned to that plot in the calculation of the average and SD of the change-point for that site.

TABLE 4 Average point-cloud density (change-point) at which FSD metrics across different LiDAR plot densities stabilize at individual National Ecological Observatory Network (NEON) forested sites and on average across all sites.

Category	Metric	ABBY	GRSM	STEI	UKFS	UNDE	Average minimum point/m2
Height	MOCH	7.4 (1.0)	7.0 (0.9)	7.1 (0.7)	6.7 (0.1)	6.8 (0.9)	7.0
Q25	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0
Q50	3.2 (2.3)	2.0 (0.0)	2.0 (0.0)	2.7 (2.0)	2.0 (0.0)	2.4
Q75	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0
Q100	2.7 (2.0)	2.6 (1.9)	5.6 (5.3)	2.0 (0.0)	2.0 (0.0)	3.0
Density	VAI	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0
Openness	DGF	5.6 (0.9)	4.8 (0.6)	5.3 (0.8)	5.1 (0.7)	4.9 (0.6)	5.1
GFP	2.0 (0.0)	2.0 (0.0)	2.7 (2.0)	2.0 (0.0)	3.1 (2.2)	2.3
External heterogeneity	Rumple	7.3 (0.8)	7.0 (0.7)	7.0 (0.6)	6.9 (0.6)	6.8 (0.2)	7.0
Top rugosity	8.1 (1.0)	7.8 (1.1)	7.1 (1.1)	7.5 (1.0)	7.0 (1.3)	7.5
Internal heterogeneity	SD(ht)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.0
SDSD(ht)	5.5 (2.2)	6.6 (1.1)	5.6 (0.9)	5.8 (1.1)	5.5 (1.1)	5.8
VCI	2.0 (0.0)	2.0 (0.0)	2.0 (0.0)	2.4 (1.1)	2.0 (0.0)	2.1
CV(ht)	2.2 (0.7)	2.0 (0.0)	2.0 (0.0)	2.5 (1.6)	2.0 (0.0)	2.2
FHD	7.1 (4.2)	4.3 (2.0)	3.5 (3.0)	4.2 (2.3)	4.4 (2.0)	4.7
Gini index	2.2 (0.7)	2.0 (0.0)	2.0 (0.0)	2.5 (1.6)	2.0 (0.0)	2.2

We assessed whether variation in LiDAR point density resulted in variable delineation of plots into multivariate structural complexity types (Fahey et al., 2019) using multivariate analyses with the 16 structural metrics. We used nonmetric multidimensional scaling (NMS) ordination on plot-level structural diversity metrics of one plot replicate from each of the eight LiDAR point densities. NMS ordinations were conducted in PC-ORD version 5.31 (McCune & Mefford, 2006) with Sorensen distance measure, the “slow-and-thorough” setting, and 250 runs of real data and 250 Monte Carlo randomizations for solution robustness testing (McCune & Grace, 2002). Metrics were standardized relative to the maximum value to scale metrics before ordination was conducted. Hierarchical agglomerative clustering in PC-ORD, using Ward method and Euclidean distance measures (McCune & Grace, 2002), was used to delineate plots into canopy structural complexity type groupings following the method of Fahey et al. (2019), but here, clusters (types) were based on structural variation in our own data. The optimal cluster group level was assessed using indicator species analysis and mean p values derived across all metrics for each level (McCune & Grace, 2002); the clustering level with the lowest mean p value was used as the optimal number of groups. We then evaluated the group membership of plots to assess whether the designation of structural type for each plot changed across LiDAR point densities.

RESULTS

Minimum LiDAR point density for estimating structural diversity across forested sites

Four structural diversity metrics, from three categories, remained stable across LiDAR point densities for all plots at each of the five sites. Of these, two were from the height category and the remaining two were from the density and internal heterogeneity categories. These included the height metrics of Q25 and Q75, the density metric of VAI (Table 4, Figure 3; Appendix S1: Table S1), and the internal heterogeneity metrics of SD(ht) (Table 4, Figure 4; Appendix S2: Table S1).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

There were 12 structural diversity metrics that varied as sample point density increased, with seven stabilizing between 2.1 and 4.7 pt/m² and five metrics stabilizing between 5.1 and 7.5 pt/m². The metrics that stabilized at low-to-moderate point density (2.1–4.7 pt/m²) included Q50 and Q100 metrics from the height category; GFP from cover and openness; and VCI, CV(ht), FHD, and Gini index from internal heterogeneity. The metrics that required higher point densities to stabilize (5.1–7.5 pt/m²) were MOCH from the height category; DGF from openness and cover; rumple and top rugosity from external heterogeneity; and SDSD(ht) from internal heterogeneity.

Several metrics, such as FHD, fluctuated between plots within sites and caused change-point variation between sites (Table 4; Appendix S1: Table S1). The change-point in the segmented regression analysis for a structural diversity metric had an average and SD that varied by site. FHD had the highest SD in change-point across sites. Most of the other 16 metrics had lower SDs at a site, with many sites having no significant change-point at any plot (i.e., metric value was stable across LiDAR point density). In general, for all metrics and plot combinations, when a change-point was found, the change-point was not more than 10 pt/m², except for three plots with Q100 (STEI_046 at 16.7) and FHD (ABBY_025 at 13.7 and ABBY_029 at 13.5) (Appendix S1: Table S1).

The effect of LiDAR point density on delineating structural complexity types

At all LiDAR point densities, five clusters representing structural complexity types were delineated across the NEON study sites (Figure 5), but at the lowest point density, cluster assignment for several plots differed from the consistent assignment pattern observed across the higher point densities. Variation in the first NMS axes (59.7%) was explained predominantly by VCI and Q100 structural diversity metrics, with the structural complexity types split as tall, heterogeneous canopies or short, homogeneous canopies. The variation in the second NMS axis (39.8%) was explained by FHD and Gini index metrics with structural complexity types split as multilayered, less skewness in height distribution, or vice versa. Groupings were identical for densities of 4–25 pt/m² (Figure 5c), but there were eight plots that were placed into different clusters at 2 pt/m², including plots from sites STEI, UKFS, GRSM, and ABBY. Overall, the delineation of structural complexity types was robust within sites and across forest types for a suite of 16 structural diversity metrics and LiDAR point density of 4 pt/m² and above.

[IMAGE OMITTED. SEE PDF]

DISCUSSION

We identified metrics of structural diversity that were not sensitive to LiDAR point density, which LiDAR users could employ to characterize forest structural diversity across time and space with less concern for the impact that data resolution would have on their analyses. Of the 16 metrics evaluated here, four showed little-to-no variation across LiDAR point density. These metrics included descriptors from height, density, and internal heterogeneity categories. Their stability, even at low point densities, suggests that these metrics could be suitable for making accurate multitemporal or cross-site comparisons between LiDAR datasets, regardless of instrumentation within the range of 2–25 pt/m². By focusing on these four metrics, it is probable that any temporal variation detected is largely due to real changes in structural diversity rather than noise. Previous studies relating point density to vertical forest structure found similar results with metrics associated with height and density, with some metrics even improving under lower densities (Jakubowski et al., 2013; Treitz et al., 2012), and recommended that metrics unaffected by low-density point clouds could be calculated with less data.

The remaining three quarters of structural diversity metrics in our study were sensitive to the point density of LiDAR, such that users of low-resolution LiDAR (e.g., some NEON flights with less than 8 pt/m²) should carefully consider the use of these metrics in temporal or spatial comparisons of forest structural diversity. On average, the values of structural diversity metrics across a range of forested NEON sites became stable at or above a LiDAR point density of 8.1 pt/m² when plots were averaged within a site. Our study therefore indicates that for the forest types studied here (i.e., Table 2), LiDAR datasets with a resolution of at least ~8 pt/m² should provide the most robust temporal and spatial comparisons of forest structural diversity. Similarly, Yao et al. (2014) found that once LiDAR achieved 10 pt/m², that any advantage of having a higher data resolution for individual tree detection plateaued. However, other studies that used pulse density to quantify LiDAR resolution found that metrics that require an internal view of the canopy were not well resolved with a low resolution (e.g., LAI and subcanopy cover; de Almeida, Stark, Shao, et al., 2019; Jakubowski et al., 2013). For example, categories of structural diversity metrics that describe the internal heterogeneity of canopies may be most susceptible to low point densities because the subcanopy often has fewer data points due to occlusion (LaRue et al., 2020). Therefore, metrics studied here, such as DGF, should be treated with caution in future work as they may be among the more biased metrics at low LiDAR point densities.

Differences in the stability of structural diversity metrics at low point densities can be explained in part by the way the point cloud is summarized to calculate each metric. The five metrics that had stable change-points at higher densities (average across sites 5.1–7.5 pt/m²) were spread across structural diversity categories but were calculated by splitting the point cloud into voxels (3D pixel; i.e., Table 2). In contrast, the six metrics that had an average change-point on the lower end (2.1–4.7 pt/m²) were calculated using points across the entire point cloud (no spatially explicit separation or splitting the data into voxels; i.e., Table 3). A lower point-cloud resolution likely causes variation in specific locations (Roussel et al., 2017) for height and area in ground surface and sparse cover (Wilkes et al., 2015), such that metric values calculated using a voxelization or grid (canopy height model) have lowered stability. Metrics using the entire point cloud (e.g., maximum height; Roussel et al., 2017) have been shown to be more vertically accurate than those using voxels (e.g., mean height from canopy height model; Roussel et al., 2017). Other studies have categorized their forest structural metrics by similar methodological categories as we have for parsing the point cloud to generate metrics (Ruiz et al., 2014; Yu et al., 2020), but they have not explicitly compared the relative stabilities of these different types of metrics (i.e., Table 3). In summary, vertical stratification of height (e.g., SD height or maximum height) across a plot may be less sensitive to lower point densities than spatial stratification of points across a plot (e.g., metrics such as rumple, mean outer canopy height, and SDSD(ht)).

Forest types that vary in their structural configuration (i.e., maximum height, canopy cover, or dominant species) may influence the point-cloud density dependence of structural diversity metrics. Indeed, the stable point density of structural diversity metrics in different sites varied (2.1–8.1 pt/m²) from the maximum 8.1 pt/m² site average. For example, we observed that ABBY, which is a managed forest dominated by Douglas-fir and western hemlock with the tallest average canopy height (Table 2), usually had the highest average change-point for structural diversity metrics. Most of ABBY's harvests appear to have left behind even-aged stands with a few residual tall trees, and thus high variability in height, with different parcels at the site have been harvested between 1940 and 2016. NEON covers many forested sites across 20 ecoclimatic domains, which vary by structure and dominant species; these include conifer- and deciduous-dominated forests like we have included in our five study sites. However, structural diversity metrics that are calculated using the entire point cloud (no voxels) were not sensitive to forest type in our study and may provide a basis for a less biased comparison of structural diversity in different forest types.

The surface canopy structure (denser canopies or evergreen vs. deciduous trees) is likely to influence how many points reach the subcanopy. There are often more points concentrated at the canopy surface of aerial LiDAR, and high-density products may be more likely to reach the subcanopy than a low-density point cloud. Our approach of artificially reducing point-cloud density may not reflect the true vertical distribution of data points within canopy profiles at low point densities. However, we randomly selected points from point clouds that started at a similar point density and from the same sensor (Table 1), which reduces potential bias for comparing the impacts of starting point-cloud density among plots in our dataset studies (Hansen et al., 2015; Jakubowski et al., 2013). This effect of point-cloud density leading to noise that does not reflect true structural diversity may be more of an issue for metrics that examine the subcanopy and in forests with high cover or dense upper canopies.

There was within-site variation in the change-point identified by the segmented regression analysis for the same metric, which might be explained by heterogeneity in the environment, differences in species composition, and disturbances (both natural and anthropogenic) that affect the structural heterogeneity of the forest across scales. For example, ABBY, located within the western foothills of the Cascade Range, has a long history of logging, resulting in a dynamic forested landscape with a mosaic of age classes and resulting structural signatures of disturbance across the plots within this site (NEON, 2021c). Similarly, UNDE and STEI sites were heavily logged as recently as 1960 and 2005, respectively, with both first and secondary growth patterns found across the sites. The GRSM site located in the Smoky Mountains of southeastern Tennessee had wildfires in 2016 that burned over more than 4,050 ha including areas within the study region. Indeed, structural diversity is known to vary differently in disturbed sites, including at GRSM after the 2016 fire (Atkins et al., 2020). It is likely that these landscape disturbances are driving variation in metric values between plots and sites and that certain structural configurations may be more sensitive to low point density LiDAR data (i.e., some metrics and canopy configurations cannot be measured as accurately with low-resolution data).

To monitor forest structural changes in response to disturbance and succession, a multivariate approach focused on categorizing forest stands into structural complexity types based on variation along axes such as maximum canopy height and internal heterogeneity of vegetation arrangement may be useful (Atkins et al., 2020; Fahey et al., 2019; Franklin & Hemstrom, 1981; Spies & Franklin, 1988). However, results from this study suggest that LiDAR data used for this type of analysis should have a minimum point density to ensure reproducibility of structural-type categorization. Ruiz et al. (2014) found that a multivariate approach of structural indices from aerial LiDAR achieves greater predictive ability of forest inventory variables such as biomass and cover over a point density of 1 pt/m², but that the marginal benefit of increasing point density declined above 5 pt/m² (Ruiz et al., 2014). The identification of forest structural complexity types among plots and sites was robust down to a LiDAR density of 4 pt/m². Therefore, our results suggest that users wanting to describe the structural diversity of canopies in their area of interest along a spectrum of forest structural complexity types using multivariate analyses (Fahey et al., 2019) can do so confidently with NEON LiDAR data of 4 pt/m² and above (for plot sizes of 40 × 40 m). Furthermore, the discrepancies in categorizing plots (8 of 45) into structural complexity types at the lowest data resolution occurred in the middle of the second NMS axis, such that the structural complexity types of these plots were correlated with a more moderate canopy height and internal heterogeneity. LiDAR of 2 pt/m² is likely to have fewer points that reach the subcanopy with considerable noise in point placement. However, our results indicate that relatively low-resolution LiDAR may be suitable for temporal comparisons of changing canopy structural complexity types using a multivariate approach. For example, the use of LiDAR to exhibit changes in structural complexity types that are indicative of successional stages or species compositions in response to disturbance or global change might be important for characterizing forest change over the 30-year time span of NEON monitoring (Dodds et al., 2021; Metzger et al., 2019).

We focused on the standard plot size that would be the focus for most NEON users for forest stand-level processes, but the spatial scale has been shown to influence structural diversity metrics (Fotis et al., 2018; Hardiman et al., 2018). The 40 × 40 m base plot footprint is the most commonly used plot size (although a few sites also use a 20 × 20 base plot size) that NEON uses across its 81 sites for woody vegetation structure and other data products from the terrestrial observation system that are collected within these plots (e.g., plant and tree diversity). NEON data users should also be aware that the plot size (or grid size for rasterization) can interact with LiDAR point density to influence the stability of forest structural metric values, with larger plot sizes, such as those having a grid size of 25 × 25 m being more stable (de Almeida, Stark, Shao, et al., 2019; Kamoske et al., 2019). However, the footprint size of structural diversity metrics is of interest, including when characterizing structural complexity types, for studying the relevance of structural diversity to forest biodiversity and ecosystem and may have a unique interaction with metric stability with footprint and LiDAR pulse density that should be further investigated.

AUTHOR CONTRIBUTIONS

Elizabeth A. LaRue conceived the original idea, conducted the analyses, and wrote the initial draft. Robert Fahey conducted the multivariate analyses. All authors edited the manuscript and approved the final version of the manuscript.

ACKNOWLEDGMENTS

Funding was provided by National Science Foundation (NSF) DEB award number 1926538 to Brady S. Hardiman and number 1924942 to Elizabeth A. LaRue and Brady S. Hardiman; number 2106103 to Elizabeth A. LaRue; number 1926442 to Robert Fahey, and number 1926454 to Jaclyn Hatala Matthes. The National Ecological Observatory Network is a program sponsored by the National Science Foundation and operated under cooperative agreement by Battelle. This material is based in part upon work supported by the National Science Foundation through the NEON Program.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

LiDAR data are available from the NEON Data Portal, and instructions for downloading the data through the neonUtilities R package are provided in our archived R code: . Code that was used to process the LiDAR and generate data and analysis results (LaRue, 2022) is available from Zenodo: .

References

de Almeida, D., S. C. Stark, G. Shao, J. Schietti, B. W. Nelson, C. A. Silva, E. B. Gorgens, R. Valbuena, D. de Papa, and P. H. S. Brancalion. 2019. “Optimizing the Remote Detection of Tropical Rainforest Structure with Airborne LiDAR: Leaf Area Profile Sensitivity to Pulse Density and Spatial Sampling.” Remote Sensing 11: 92.

de Almeida, D., S. Stark, C. Silva, C. Hamamura, and R. Valbuena. 2019. “leafR: Calculates the Leaf Area Index (LAD) and Other Related Functions.” CRAN. https://cran.r-project.org/web/packages/leafR/index.html.

Atkins, J. W., G. Bohrer, R. T. Fahey, B. S. Hardiman, T. H. Morin, A. E. L. Stovall, and C. M. Gough. 2018. “Quantifying Vegetation and Canopy Structural Complexity from TLS Data Using the Forestr r Package.” Methods in Ecology and Evolution 9: 2057–66.

Atkins, J. W., B. Bond‐Lamberty, R. T. Fahey, L. T. Haber, E. Stuart‐Haëntjens, B. S. Hardiman, E. LaRue, et al. 2020. “Application of Multidimensional Structural Characterization to Detect and Describe Moderate Forest Disturbance.” Ecosphere 11: [eLocator: e03156].

Balch, J. K., R. C. Nagy, and B. S. Halpern. 2020. “NEON Is Seeding the Next Revolution in Ecology.” Frontiers in Ecology and the Environment 18: 3.

Balsa‐Barreiro, J., and J. L. Lerma. 2014. “Empirical Study of Variation in LiDAR Point Density over Different Land Covers.” International Journal of Remote Sensing 35: 3372–83.

Bouvier, M., S. Durrieu, R. A. Fournier, and J.‐P. Renaud. 2015. “Generalizing Predictive Models of Forest Inventory Attributes Using an Area‐Based Approach with Airborne LiDAR Data.” Remote Sensing of Environment 156: 322–34.

Dodds, W. K., K. C. Rose, S. Fei, and S. Chandra. 2021. “Macrosystems Revisited: Challenges and Successes in a New Subdiscipline of Ecology.” Frontiers in Ecology and the Environment 19: 4–10.

van Ewijk, K. Y., P. M. Treitz, and N. A. Scott. 2011. “Characterizing Forest Succession in Central Ontario Using LiDAR‐Derived Indices.” Photogrammetric Engineering & Remote Sensing 77: 261–9.

Fahey, R. T., J. W. Atkins, C. M. Gough, B. S. Hardiman, L. E. Nave, J. M. Tallant, K. J. Nadehoffer, et al. 2019. “Defining a Spectrum of Integrative Trait‐Based Vegetation Canopy Structural Types.” Ecology Letters 22: 2049–59.

Fotis, A. T., T. H. Morin, R. T. Fahey, B. S. Hardiman, G. Bohrer, and P. S. Curtis. 2018. “Forest Structure in Space and Time: Biotic and Abiotic Determinants of Canopy Complexity and Their Effects on Net Primary Productivity.” Agricultural and Forest Meteorology 250: 181–91.

Franklin, J. F., and M. A. Hemstrom. 1981. “Aspects of Succession in the Coniferous Forests of the Pacific Northwest.” In Forest Succession: Concepts and Application, edited by D. C. West, H. H. Shugart, and D. B. Botkin, 212–29. New York: Springer‐Verlag.

Gewin, V. 2016. “Data Sharing: An Open Mind on Open Data.” Nature 529: 117–9.

González‐Ferreiro, E., U. Diéguez‐Aranda, and D. Miranda. 2012. “Estimation of Stand Variables in Pinus radiata D. Don Plantations Using Different LiDAR Pulse Densities.” Forestry 85: 281–92.

Gough, C. M., J. W. Atkins, R. T. Fahey, and B. S. Hardiman. 2019. “High Rates of Primary Production in Structurally Complex Forests.” Ecology 100: [eLocator: e02864].

Guo, Q., Y. Su, T. Hu, H. Guan, S. Jin, J. Zhang, X. Zhao, et al. 2021. LiDAR Boosts 3D Ecological Observations and Modelings: A Review and Perspective. Piscataway, NJ: Institute of Electrical and Electronics Engineers Inc.

Hakkenberg, C. R., C. Song, R. K. Peet, and P. S. White. 2016. “Forest Structure as a Predictor of Tree Species Diversity in the North Carolina Piedmont.” Journal of Vegetation Science 27: 1151–63.

Hansen, E., T. Gobakken, and E. Næsset. 2015. “Effects of Pulse Density on Digital Terrain Models and Canopy Metrics Using Airborne Laser Scanning in a Tropical Rainforest.” Remote Sensing 7: 8453–68.

Hardiman, B. S., E. A. LaRue, J. W. Atkins, R. T. Fahey, F. W. Wagner, and C. M. Gough. 2018. “Spatial Variation in Canopy Structure across Forest Landscapes.” Forests 9: 474.

Huang, L., W. Xue, and T. Herben. 2019. “Temporal Niche Differentiation among Species Changes with Habitat Productivity and Light Conditions.” Journal of Vegetation Science 30: 438–47.

Jakubowski, M. K., Q. Guo, and M. Kelly. 2013. “Tradeoffs between LiDAR Pulse Density and Forest Measurement Accuracy.” Remote Sensing of Environment 130: 245–53.

Jenness, J. 2004. “Calculating Landscape Surface Area from Digital Elevation Models.” Wildlife Society Bulletin 32: 829–39.

Kamoske, A. G., K. M. Dahlin, S. P. Serbin, and S. C. Stark. 2021. “Leaf Traits and Canopy Structure Together Explain Canopy Functional Diversity: An Airborne Remote Sensing Approach.” Ecological Applications 31: [eLocator: e02230].

Kamoske, A. G., K. M. Dahlin, S. C. Stark, and S. P. Serbin. 2019. “Leaf Area Density from Airborne LiDAR: Comparing Sensors and Resolutions in a Temperate Broadleaf Forest Ecosystem.” Forest Ecology and Management 433: 364–75.

Krause, K., and T. Goulden. 2015. NEON L0‐to‐L1 Discrete Return LiDAR Algorithm Theoretical Basis Document (ATBD). NEON.Doc.001291. Boulder, CO: NEON.

LaRue, E. 2022. “lizlarue/NEONLidarFSDPointDensity: NEON LiDAR Point Density: Forest Structural Diversity (StructuralDiversity).” Zenodo. Code. [DOI: https://dx.doi.org/10.5281/zenodo.6463393]

Roussel, J.‐R., D. Auty, N. C. Coops, P. Tompalski, T. R. H. Goodbody, A. S. Meador, J.‐F. Bourdon, F. de Boissieu, and A. Achim. 2020. “lidR: An R Package for Analysis of Airborne Laser Scanning (ALS) Data.” Remote Sensing of Environment 251: [eLocator: 112061].

Roussel, J. R., J. Caspersen, M. Béland, S. Thomas, and A. Achim. 2017. “Removing Bias from LiDAR‐Based Estimates of Canopy Height: Accounting for the Effects of Pulse Density and Footprint Size.” Remote Sensing of Environment 198: 1–16.

Ruiz, L. A., T. Hermosilla, F. Mauro, and M. Godino. 2014. “Analysis of the Influence of Plot Size and LiDAR Density on Forest Structure Attribute Estimates.” Forests 5: 936–51.

Shao, G., S. C. Stark, D. R. A. de Almeida, and M. N. Smith. 2019. “Towards High Throughput Assessment of Canopy Dynamics: The Estimation of Leaf Area Structure in Amazonian Forests with Multitemporal Multi‐Sensor Airborne LiDAR.” Remote Sensing of Environment 221: 1–13.

Singh, K. K., G. Chen, J. B. McCarter, and R. K. Meentemeyer. 2015. “Effects of LiDAR Point Density and Landscape Context on Estimates of Urban Forest Biomass.” ISPRS Journal of Photogrammetry and Remote Sensing 101: 310–22.

Spies, T. A., and J. F. Franklin. 1988. “Old Growth and Forest Dynamics in the Douglas‐Fir Region of Western Oregon and Washington.” Natural Areas Journal 8: 190–201.

Strunk, J., H. Temesgen, H. E. Andersen, J. P. Flewelling, and L. Madsen. 2012. “Effects of LiDAR Pulse Density and Sample Size on a Model‐Assisted Approach to Estimate Forest Inventory Variables.” Canadian Journal of Remote Sensing 38: 644–54.

Treitz, P., K. Lim, M. Woods, D. Pitt, D. Nesbitt, and D. Etheridge. 2012. “LiDAR Sampling Density for Forest Resource Inventories in Ontario, Canada.” Remote Sensing 4: 830–48.

Valbuena, R., M. Maltamo, L. Mehtätalo, and P. Packalen. 2017. “Key Structural Features of Boreal Forests may Be Detected Directly Using L‐Moments from Airborne LiDAR Data.” Remote Sensing of Environment 194: 437–46.

Wilkes, P., S. D. Jones, L. Suarez, A. Haywood, W. Woodgate, M. Soto‐Berelov, A. Mellor, and A. K. Skidmore. 2015. “Understanding the Effects of ALS Pulse Density for Metric Retrieval across Diverse Forest Types.” Photogrammetric Engineering and Remote Sensing 81: 625–35.

Yao, W., J. Krull, P. Krzystek, and M. Heurich. 2014. “Sensitivity Analysis of 3D Individual Tree Detection from LiDAR Point‐Clouds of Temperate Forests.” Forests 5: 1122–42.

Yu, X., A. Kukko, H. Kaartinen, Y. Wang, X. Liang, L. Matikainen, and J. Hyyppä. 2020. “Comparing Features of Single and Multi‐Photon LiDAR in Boreal Forests.” ISPRS Journal of Photogrammetry and Remote Sensing 168: 268–76.

Zipkin, E. F., E. R. Zylstra, A. D. Wright, S. P. Saunders, A. O. Finley, M. C. Dietze, M. S. Itter, and M. W. Tingley. 2021. “Addressing Data Integration Challenges to Link Ecological Processes across Scales.” Frontiers in Ecology and the Environment 19: 30–8.