Assessing Site Productivity via Remote Sensing—Age-Independent Site Index Estimation in Even-Aged Forests

Abstract

Forest productivity is a key driver of forest growth and yield and a critical information need for forest management and planning. Traditionally, this information has come from field plots, but these are expensive to measure and have limited coverage. Remote sensing, on the other hand, can provide forest inventory attributes on landscape scales and with a relatively low cost. A common predictor of forest productivity is site index (SI), traditionally estimated from age and height. In plantations, age can often be treated as a known quantity, but in natural-origin forests (of which Canada has vast swaths), age is often unknown and must be estimated, requiring expensive field work and resulting in a high level of error which, in turn, introduces error into SI estimates. The objective of this study is to generate estimates of SI from two successive LiDAR captures. The 99th percentiles (p99) of LiDAR returns from two successive captures 13 years apart were used along with species-specific SI curves to estimate SI. The results were compared to field-based estimates of SI for two major boreal species, jack pine and black spruce in managed and unmanaged conditions. Overall, the difference between the LiDAR-based SI and the field estimate was 2% with a relative mean squared error of 18%. For the few situations in which the height change was small or negative (less than 0.5%/year), SI was estimated from the average p99 and an assumed age of 100. The advantage of this method is that it does not require field sampling or estimates of age. Using two successive LiDAR captures, wall to wall estimates of SI can be generated at the grid cell level (e.g., 20 × 20 m), a level of detail generally not found in inventories. Overall, our results demonstrate the excellent potential for estimating SI from LiDAR alone, without age, to provide detailed productivity information for forest management and inventory that has been lacking in most large-scale inventories until now.

1. Introduction

Forest site productivity is among the most important stand attributes for forest management planning and is crucial for characterizing growth [1]. Assessments of site quality are needed to identify the productivity potential of forest land and provide context for silvicultural planning. Site quality measures should be quantitative, objective and easily assessed [2]. Effective large-scale forestry planning and management requires these site productivity measures to be available at a fine resolution at the regional level [3].

Site index (SI—tree height at a reference or index age) is the most widely accepted method for estimating forest site quality in North America and Europe [4,5,6,7,8]. SI equations have been calibrated for most of the commercial species in Ontario, Canada [9,10,11,12,13,14,15,16,17,18]. Traditionally, SI has been estimated from tree or stand age and height at a single point in time. Some SI equations can also be reformulated to be age-independent and SI can be estimated from two measurements of height [19].

In Ontario, forest management inventories have been derived from the manual interpretation of aerial photographs. Photo interpreters typically estimate stand height and SI and impute age, or estimate stand height and age and impute SI [20]. More recently, satellite-based imagery and airborne LiDAR have been used to estimate site productivity. For instance, LiDAR-based digital terrain models (DTMs) have been used to generate topographic and wetness variables to predict SI [21]. Other studies have used remote sensing to directly derive age and height to estimate SI [7].

The field estimation of SI is time-consuming and expensive [22,23]. Age is one of the more difficult and expensive inventory attributes to assess or predict and is often one of the less reliable attributes [24]. Measuring age from increment cores is costly and may require specialized equipment to ensure accurate and consistent results. Inventory estimates of age can be obtained from silvicultural records and from performance surveys and time since disturbance (e.g., fire, harvesting), but for natural origin forests (the majority in Canada), this information is often not known.

In Ontario, inventory stand age has typically been estimated from stand height and a reference set of height–age curves [25], but Schumacher et al. [24] have noted that the estimation of age using remote sensing shows promise. Landsat time series have been used to identify the occurrence of stand-replacing forest disturbances [26] and thus estimate the year of stand establishment. Landsat-derived maps of stand age have also been used with LiDAR data to create SI maps [7,8]. This works well if most of the dominant/codominant trees are established immediately following the disturbance. However, when dominant/codominant trees are residuals, established prior to the disturbance, or if regeneration is delayed, the level of error may be unacceptably high.

Height estimates from remote sensing, particularly LiDAR, are generally very reliable. This has led to research into using multi-temporal remote sensing to predict SI [6], to augment SI estimates [27] and to develop SI curves [28,29] at the plot or grid cell level (an area-based approach). All of these methods either require field data for calibration or an estimate of age. Solberg et al. [19] used individual tree classification methods to estimate SI from successive LiDAR captures independent of age and field data. We found only one study [30] that used area-based methods and two successive LiDAR captures to obtain age-independent estimates of SI without the need for field calibration data. The algorithm iterated through a range of SI guesses to find the SI that minimized the residuals between the observed p99 and the SI curve. Our study differs from that study in that we provide an analytic solution for determining SI, rather than a “closest fit” solution.

The objective of the present study was to develop an area-based approach to obtain age- and field-data-independent estimates of SI from two successive LiDAR captures. The results were validated against independent field data.

2. Materials and Methods

2.1. Study Area

Field measured plot data came from the Romeo Malette Forest (RMF) (Figure 1). The RMF is located in the northeast portion of the boreal forest in Ontario, Canada. It is characterised by conifer-dominated stands on poorly drained lowlands and intolerant hardwoods, conifers and mixedwoods on gently rising uplands [31].

2.2. Field Plots

The data came from a network of permanent sample plots (PSPs) with 3–4 measurements, with the last in 2018. These PSPs are part of the Ontario provincial growth and yield program [12] and are deliberately located in conditions of interest. Plots were selected from the database to cover a range of ages and site productivities in relatively pure species conditions and we targeted natural and planted jack pine (Pinus banksiana Lamb.) and black spruce (Picea mariana [Mill] BSP). The jack pine plots were located on dry to fresh, sandy to coarse loamy soils of glaciofluvial origin. The black spruce plots were located on wet organic soils.

The diameter at breast height (Dbh) was recorded for all live and dead trees on a circular 400 m² plot with Dbh ≥ 7.1 cm, which is the threshold used in Ontario for LiDAR field calibration plots. A subsample was measured for height and crown class (dominant, codominant, emergent, intermediate, overtopped, suppressed, anomaly). Species, status (live or dead) and origin (natural, planted) were recorded for all trees.

2.3. LiDAR Data

The LiDAR data used in this study came from two independent projects to generate operational forest inventories. The 2005 acquisition was one of the first LiDAR-based inventories in Ontario [31], while the 2018 acquisition reflects the current government standard [32]. Leaf-on linear-mode LiDAR data were acquired in 2005 using an upgraded Leica ALS40 sensor flown at an average altitude of 2740 m. Leaf-on Single-Photon LiDAR (SPL) data were acquired in 2018 using a Leica SPL100 sensor flown at an average altitude of 3800 m. Noise filtering was performed by the data provider using the approach detailed in [33]. All LiDAR acquisitions were normalized against the reference 2018 LiDAR-derived DTM for consistency [34]. LiDAR metrics were generated via lidR [35,36], using the full point cloud (all returns) above zero metres. The point densities are summarized in Table 1.

2.4. Site Index

Most SI equations have three variables—age, height and SI. Generally these equations are fit with height as the dependent variable and age and SI as the independent variables.

The most recent SI curves developed in Ontario use the following equation form, where ${bhage}_i$ is breast height age (age 1.3 m above the ground) at time i, ${Ht}_i$ is height at time i and a₀ and a₁ are species- and origin-specific coefficients (Table 2).

(1)${Ht}_2=\frac{a_0}{1-\left(1-\frac{a_0}{{Ht}_1}\right)·{\left(\frac{{bhage}_1}{{bhage}_2}\right)}^{a_1}}$

When ${bhage}_2=50$, then ${{Ht}_2=SI}_{50}$, the SI at reference age 50.

Equation (1) can be rewritten to predict age.

(2)${bhage}_2={\left[\frac{\left(1-\frac{a_0}{{Ht}_1}\right)}{\left(1-\frac{a_0}{{Ht}_2}\right)}\right]}^{\frac{1}{a_1}}·{bhage}_1$

Assuming that age increments with time, Equation (2) can be rewritten as follows, where year_i is the year corresponding to ${Ht}_i$.

(3)${bhage}_1=\frac{\left({year}_2-{year}_1\right)}{{\left[\frac{\left(1-\frac{a_0}{{Ht}_1}\right)}{\left(1-\frac{a_0}{{Ht}_2}\right)}\right]}^{\frac{1}{a_1}}-1}$

${SI}_{50}$ can then be calculated by substituting ${bhage}_i$ from Equation (3) in the following equation.

(4)${SI}_{50}=\frac{a_0}{1-\left(1-\frac{a_0}{{Ht}_1}\right)·{\left(\frac{{bhage}_1}{50}\right)}^{a_1}}$

SI curves can be applied at the tree level, plot level or stand level. When applied at the stand level, the average height of the dominant and codominant trees (CDHt) is generally used along with the stand age.

Two field measures of SI (reference and max) were compared to one SI estimate from LiDAR.

2.4.1. Reference SI

For the PSPs, breast height ages were sampled in the field, usually at the first measurement and sometimes at later measurements. The mean age of the dominant/codominant trees of the leading species was used in the calculation of the field reference SI.

At each measurement, a sample of heights was measured (full tally in 2018). Dominant/codominant height (CDHt) was calculated as the mean height of the dominant and codominant trees of the leading species that were measured for height, based on advice from Dr. Mahadev Sharma, Ontario Forest Research Institute.

SI can be estimated for each measurement and can result in a range of SI estimates. The reference SI was estimated as the SI curve which minimized the squared distances from the age/CDht pairs. The observed age was used for each measurement, if available. The data are summarized in Table 3.

The use of average plot CDHt, age and SI is analogous to using polygon CDHt age and SI in the inventory [25].

2.4.2. Max SI

An alternative to the reference SI was investigated. The max SI was estimated as the SI of the largest Dbh tree of the leading species measured for age for the most recent measurement year. The use of individual tree age and height to calculate SI was analogous to the way that the SI curves used here were developed, from individual tree stem analysis [13].

2.5. LiDAR SI

To estimate SI from LiDAR, the 99th percentile of the LiDAR returns (p99) was used as a surrogate for CDHt. Each field plot had two p99 measurements from different acquisition years. The LiDAR point density varied with acquisition and p99 is much less sensitive to the point density than the reported maximum return. The following algorithm was used to estimate SI for the validation field plots.

• Estimate the breast height age corresponding to time 1 (bhage_i) using Equation (3).

• Estimate the SI using Equation (4).

This is similar to the method of Persson and Fransson [30], except they minimized the residuals between the SI curve and the observed height. As noted, the most common SI equations predict height as a non-decreasing function of age. Consequently, it is not possible to estimate SI when the height change over time is negative. The height growth rate generally decreases at older ages and minor changes in height growth can have a big impact on SI estimates.

One option for older stands where the height change is small or negative is to assume that the height growth is effectively zero and compare the average height to the SI curves at some fixed older age. In this study, this option was investigated by setting the height change threshold to 0.5%/yr. When the height change was less than this threshold, SI was estimated from the average p99 from the separate LiDAR captures and an assumed age of 100.

This algorithm was also used to predict SI for jack pine and black spruce for larger areas based on the p99 raster layers. The raster cell size corresponded to the field plot area (0.04 ha).

2.6. Comparison of LiDAR and Reference SI

The SI predictions were evaluated in terms of accuracy and precision by regressing the reference SI on the LiDAR-predicted SI. High accuracy (low bias) was suggested if the regression had a slope of one and an intercept of zero; precision was characterized by the degree of scatter in the observations around the regression line, as measured via the root mean square error (RMSE) of the regression. A simultaneous test of H0: intercept, b₀ = 0 and slope, b1 = 1, was conducted using the following test statistic [37]:

(5)$F_{2,n-2}=\frac{n{\left(b_0-0\right)}^2+2\sum \widehat{y_i}\left(b_0-0\right)·\left(b_1-1\right)+\sum {\widehat{y}}_i^2{\left(b_1-1\right)}^2}{2\sum {\left(\widehat{y_i}-\overline{y}\right)}^2/\left(n-2\right)}$

One weakness in the above approach is that the power of this test varies with precision. When precision is low, the sensitivity of the test is low.

3. Results

The validation plots cover a range of ages, SI (Table 3) and stand conditions (Table 4).

Results are illustrated for two plots (Figure 2). Plot 254 is a young jack pine plantation. The field measurements are consistent with an SI of 20 m and the LiDAR predictions of SI agree very well. Plot 435 is in an older black spruce stand that shows a decline in CDHt between the last two measurements. Depending on the year of measurement, the reference SI ranges from 12 to 14 m.

The LiDAR-based estimates of SI are highly correlated with the reference SI (Figure 3). However, some plots had a small or negative annual change in p99 and were assigned an SI = 8 m, the lowest observed SI. These plots tended to be older (with an average breast height age of 97 years in 2018). If the annual change in LiDAR heights was less than 0.5%/year, SI was estimated from the average height (the average of the 2005 and 2018 p99) and an assumed age of 100. The results are given in Figure 4 and Table 5.

Overall, there is good agreement between the reference and LiDAR SI (overall bias of 2%), particularly for jack pine (Figure 5). There is more variation for black spruce. The simultaneous F-test in Equation (5) testing that the relationship between the LiDAR SI and reference SI (Figure 4) had a slope = 1 and intercept = 0 was not statistically significant (prob. (F > Fobs) = 0.08).

In general, differences between reference SI and max SI are smaller for plantations and larger for natural origin stands. The age range of trees in a plantation is usually smaller than the range in natural stands, and using a tree level estimate of SI (max reference) rather than a plot level estimate (max) may be preferable for stands with a broader age range. The simultaneous F-test for the LiDAR SI and max SI had prob. (F > Fobs) = 0.024, indicating departure from a 1:1 relationship.

LiDAR estimates of SI depend on the SI curve which depends on the leading species and origin. The SI curves in this study were very similar, except for planted black spruce (Figure 6).

Raster maps of SI can be produced to aid in forest management planning and the identification of areas of higher and lower productivity (Figure 7).

4. Discussion

The results show promise in using successive LIDAR captures to estimate SI. The results are particularly promising for relatively pure species, even-aged conditions. The approach has a number of significant advantages compared to most other methods for deriving landscape-scale maps of SI. It does not require calibration data, does not require estimates of age and has very few assumptions. Field plots were used here for validation rather than calibration. LiDAR height (p99) and time elapsed between two captures replaced age. The approach does require the existence of appropriate SI curves and these were taken from the literature. The approach also requires knowledge of the leading species and stand origin. For some species, particularly natural origin jack pine and black spruce, the differences in the SI curves are small. This may be an advantage in mixed species conditions where the leading species may be less than 50% of the species composition or when the leading species changes over time. If substantial differences in SI among species exist, one option is to develop species-independent SI models. In place of species, Nigh [38] found that models using species shade tolerance (tolerant vs. intolerant) and region (coast vs. interior) rather than species had approximately the same accuracy as species-specific models.

P99 was used as a surrogate for CDHt in estimating SI, as has been used by others (e.g., [30]). Using p99 directly, rather than calibrating a prediction model, eliminates the need for calibration data. P99 was chosen because it generally has a high correlation with dominant/codominant height and, unlike the maximum (p100), should be relatively insensitive to the LiDAR mode (linear vs. SPL), to point density and to plot size. However, Socha et al. [28] found that SI curves generated from p100 were much closer to those generated from the reference trees than SI curves generated from p95 and p90. The two LiDAR acquisitions in this study had substantially different point densities (Table 1). The 2005 acquisition had a point density of 0.56 pts/m² (all returns), which was below the 1 pt/m² recommended in White et al. [39] for area-based approaches (used here), and 4 pts/m² for individual tree identification (not used here). However, Treitz et al. [40], studying the same forest, found that a point density of 0.5 pts/m² was sufficient. Higher point densities increased the processing time but did not present any computational issues in this study.

The results, an overall bias of 2% and an RMSE of 18%, are comparable to other studies that used an area-based approach to predict SI from LiDAR. Persson and Fransson [30] reported an RMSE% of 25.4% (uncorrected for bias) for Norway spruce based on the change in p99 over three growing seasons. Noordermeer et al. [6] used linear regression to predict SI from metrics from two successive LiDAR captures for spruce and pine in Norway. The models varied with species and district and had RMSE% ranging from 9.8% to 20.0%. Gopalakrishnan et al. [7] predicted SI for young plantation loblolly pine (Pinus taeda L.), using age derived from Landsat and height derived from LiDAR, and reported RMSE% of 15.1%–19.7%. The correlation between the LiDAR-based and reference SI was higher than that reported by Bjelanovic et al. [21] for three boreal species in Alberta (0.85 vs. 0.57), but the RMSE% was higher (19% vs. 8%–12%). The RMSE% for jack pine and natural black spruce were comparable to the results reported by McRoberts et al. [41] for field measurements. Véga and St.-Onge [29], using aerial photography acquired from 1945 to 2003, matched height observations to reference SI curves to predict SI and age for natural jack pine with an average bias of 0.76 m and RMSE of 2.41 m for SI and average bias of 1.86 years and RMSE of 7 years for age. Tompalski et al. [8] used LiDAR-derived height and Landsat time-series-derived age to estimate SI in high-productivity temperate rainforests with an average bias of 6%–16% and RMSE of 47%–58%.

The evaluation of the results is heavily dependent on the reference SI. Most assessments of accuracy and precision assume that the reference data are free from error and that any differences between the reference data and predictions are due to prediction errors. However, in this study, the reference data (field SI) included a number of sources of error. Both tree age and height include measurement error. The plot ages include sampling error. The plot age is the average from a small sample (generally n = 3), usually for one measurement. The plot height includes sampling error, as the trees measured for height are not necessarily the same for each measurement. The SI equations include model error.

In the field estimation of SI, McRoberts et al. [41] identified sources of uncertainty, including the selection of the site trees in the plot, the species selected, the estimation of the tree height and age and the accuracy of the SI curves. In sugar maple stands in Michigan, eight crews measured 16 trees in one plot and nine crews measured 18 trees in another plot. The differences in SI had a standard deviation of 2.90 m. Luoma et al. [42] compared height measurements from four trained mensurationists and found no bias but observed a standard deviation of 0.5 m (2.9%) among the measurements. The relative precision was not affected by tree height but was higher for deciduous trees than conifers. Using the same study site, Wang et al. [43] found that individual tree height estimates from airborne LiDAR were as reliable as or more reliable than field measurements of height. Therefore, it is important to recognize and consider the uncertainty associated with field SI estimates when evaluating alternative estimation methods like remote sensing approaches including LiDAR, which has a much higher level of measurement accuracy.

We selected as a reference SI the SI that minimized the sums of squares. Véga and St.-Onge [29] minimized the absolute deviations. Another option would be to select the measurement closest to the SI reference age. That option offers the advantage of being the least sensitive to SI curves. If the SI reference age is 50, then field measurements at the breast height age of 50 are independent of the SI curves selected.

LiDAR height summaries are generally considered to be very accurate. Sources of error include the reference DTM and the point density. The 99th percentile of LiDAR returns should be relatively unbiased with respect to point density. Matching the field plot to the LiDAR point cloud also includes uncertainty. The GPS location accuracy for the plot centre, although of high quality, still includes error.

Trees generally do not decrease in height except when the main stem is broken or the top is killed, and such trees are generally not used in SI estimation. The SI equations used here were all developed using stem analysis. The equation form selected is a non-decreasing function of age. As a consequence, it is not possible to estimate SI if the height does not change, or if it decreases over time. Individual trees that are dominant in mature stands may not have occupied a similar canopy position at younger ages. While individual trees do not generally decline in height, at the plot or stand level, the dominant height can decline when dominant trees are replaced by shorter trees over time, possibly as a result of mortality or partial harvesting [2]. Also, dominant height may appear to decline when the measurement error is high and height growth is slow. As noted above, the field measurement of individual tree height is subject to measurement error.

In this study, if the height change was small or negative (less than 0.5%/year), SI was estimated assuming an age of 100. The 0.5%/year threshold was selected based on the average height change at an age of 100. For SIs ranging from 10 to 20 m, the annual change in the SI height at an age of 100 was 0.5%/year or less, except for planted spruce, which was 0.6%/year. The choice of an age of 100 was somewhat subjective but corresponded well to the age of onset of old growth for jack pine and black spruce in Ontario, where the age of onset of old growth is defined as when trees reach 75% of their potential Dbh and height growth declines [44].

If a plot height decrease occurs when the stand is mature/overmature, it is likely that older trees are dying from age-related mortality and being replaced by shorter trees that are possibly of a younger age cohort and possibly of a different species. The height growth rate decreases with age. When the height growth is small, our method is very sensitive to small changes in height measurements.

When the height decrease occurs at younger ages, it is unlikely to be related to age-related mortality but more likely to be a result of disturbance (e.g., wind, insect, diseases) killing dominant/codominant trees. Referring to LIDAR-based height estimates, Wulder et al. [45] cautioned that the magnitude of measurement error should be considered when estimating growth.

The methods here can be used to develop a raster layer for SI to replace the current polygon-based SI layer. A higher-resolution SI layer should lead to improved forest management through the application of site-specific silviculture and growth projections.

We focused here on estimating SI. Most studies that have looked at predicting productivity using remote sensing have included a relatively large number of attributes. Rahimzadeh-Bajgiran et al. [3] used 21 variables derived from Sentinel-2 imagery along with more than 7700 field plots to derive a productivity map. Rather than remote sensing, some have used geospatial attributes along with forest attributes to generate large-scale forest productivity maps. Parresol et al. [46] predicted SI with an average bias of 11% and an RMSE of 22%. Bjelanovic et al. [21] used almost 200 field plots and site variables derived from a LiDAR-derived DTM to predict SI and obtained RMSE% of 7%–9% for aspen, 10%–15% for lodgepole pine and 10%–13% for white spruce.

As noted above, one of the limitations of predicting SI is obtaining reliable reference data. Although the methods described here do not require field calibration data, field data were used to validate the models. The LiDAR-derived DTM was used to normalize the LiDAR data. The DTM can also be used to derive topographic variables including slope, aspect, topographic wetness index and depth to water table [21] to predict SI.

The approach used here has relatively few assumptions. It does require two successive LiDAR acquisitions normalized to a common DTM as well as spatial coverage of the dominant species and associated species-specific SI curves. Our approach is relatively insensitive to the LiDAR mode (linear vs. SPL) and point density.

A future step is to combine the methods here based on the canopy height with SI prediction methods based on the DTM (e.g., Bjelanovic et al. [21]) to potentially further improve predictions. In addition, the application of the methods should be tested for additional tree species and forest environments.

5. Conclusions

The results here demonstrate that successive LiDAR captures can be used to estimate SI accurately and precisely for relatively pure forest conditions with positive height growth. The method does not require calibration data or an estimate of age. Calibration data are costly to acquire and age is one of the more expensive and less reliable inventory attributes.

The results reported here are for relatively operational conditions. LiDAR technology has improved rapidly in the last two decades. Despite changes in LiDAR sensors and increasing point density over the time interval, the accuracy and precision of the SI estimates were comparable to other studies. An option for estimating SI when the change in p99 is small or negative showed good results. The methods can be used to efficiently map SI over large areas given repeated LiDAR captures.

Author Contributions

Conceptualization, M.P., M.W. and A.B.; methodology, M.P. and M.W.; software, M.P. and M.W.; validation, M.P. and M.W.; formal analysis, M.P.; writing—original draft preparation, M.P.; writing—review and editing, M.P., M.W. and A.B.; project administration, A.B.; funding acquisition, M.P., M.W. and A.B. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Restrictions apply to the availability of these data. Data were obtained from the Ontario Ministry of Natural Resources and Forestry and are available from the authors with the permission of the Ontario Ministry of Natural Resources and Forestry.

Acknowledgments

Thank you to the Ontario Ministry of Natural Resources and Forestry for access to the LiDAR and field plot data. We are grateful to all of the crews who acquired the field data used in our analysis. Thank you to the editor and reviewers for their thorough and constructive comments.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The Romeo Malette Forest location is given.

Figure 2. The reference (field) and LiDAR heights are compared for plot 254, a young jack pine plantation, and plot 435, an older black spruce stand. SI curves are also given. For plot 254, the reference SI = 20.1 m and the LIDAR SI = 20.1 m. For plot 435, the reference SI = 13.6 m and the LiDAR SI = 14.0 m.

Figure 3. SI estimates based on the p99 values in 2005 and 2018 are compared to reference SI (least squares) based on the field measurements. Plots with negative height change were assigned a minimum SI = 8 m. The correlation, r, was 0.8188, prob. (H0: r = 0) < 0.0001.

Figure 4. The same as Figure 3 except that if the annual p99 change is less than 0.5%, SI is estimated from the average p99 and assumed age of 100. Eight out of twenty-seven plots had an annual height change less than 0.5%; three had a negative change. The correlation, r, was 0.8534, prob. (H0: r = 0) < 0.0001.

Figure 5. SI estimates are given by species and origin.

Figure 6. The height predictions associated with SI = 20 m are given by species and origin. There is little difference except for planted black spruce. The maximum breast height age associated with planted black spruce was 46 years.

Figure 7. The jack pine site index (top) and digital terrain model (bottom) for a forest polygon are given.

Table 1

The point density for the field plots is given by acquisition year.

Year	LiDAR Type	Point Density (Pts/m2)
		All Returns	Last Returns
2005	Linear mode	0.56	0.23
2018	Single-Photon LiDAR (SPL)	41.69	35.45

Table 2

The coefficients associated with Equation (1) are given by species and origin (P = planted, N = natural) along with the maximum age used in calibrating the equations.

Species	Origin	a 0	a 1	Reference	Maximum Age (Breast Height)
Jack pine	N	32.2567	1.2156	Table 4 of [14]	157
Jack pine	P	30.7690	1.1103	Table 2 of [13]	59
Black spruce	N	35.3570	1.1233	Table 4 of [14]	180
Black spruce	P	63.8046	1.1638	Table 2 of [13]	46

Table 3

The breast height age at the last measurement and the reference SI are given by species and origin. The mean is followed in brackets by the range.

Species	Origin	Number of Plots	Breast Height Age(Years)	Reference SI (m)
Jack pine	N	4	74.6 (54.3–98.0)	17.5 (14.4–21.7)
Jack pine	P	7	42.2 (35.0–55.0)	19.4 (17.9–21.9)
Black spruce	N	13	88.9 (55.7–116.3)	10.9 (8.0–15.1)
Black spruce	P	3	40.4 (27.7–56.8)	13.2 (13.0–13.6)

Table 4

The plot basal area and gross total stem volume are given by species and origin as well as the average Dbh and height of the dominant/codominant trees in the plot. The mean is followed in brackets by the range.

Species	Origin	Number of Plots	Plot		Dominant/Codominant Trees
			Basal Area(m2/ha)	Gross Total Stem Volume (m3/ha)	Dbh(cm)	Height (m)
Jack pine	N	4	32.0 (23.0–43.4)	297 (205–453)	22.3 (18.7–28.5)	20.7 (17.9–22.5)
Jack pine	P	7	33.5 (28.8–38.3)	274 (214–357)	16.3 (13.4–20.8)	17.7 (14.4–20.3)
Black spruce	N	13	30.2 (3.1–42.0)	217 (13–320)	17.8 (10.0–23.1)	16.0 (8.5–20.1)
Black spruce	P	3	20.2 (9.7–30.6)	115 (37–204)	11.3 (8.4–14.8)	10.7 (7.6–14.2)

Table 5

The results are summarized by species and origin. If the change in p99 is less than 0.5%/yr, the LiDAR SI is based on the average p99 and assumed age of 100. The bias and root mean squared error (RMSE) are expressed as a percent of the reference SI.

Species	Origin	N		SI		Ref. vs. LiDAR		Max vs. LiDAR
			Reference	Max	LiDAR	% Bias	% RMSE	% Bias	% RMSE
Jack pine	N	4	17.5	19.9	19.2	10%	17%	−3%	22%
Jack pine	P	7	19.2	19.4	19.1	−2%	14%	−2%	17%
Black spruce	N	13	10.5	12.8	11.8	9%	17%	−8%	23%
Black spruce	P	3	13.3	15.1	9.6	−28%	30%	−37%	45%
All	All	27	14.1	15.8	14.4	2%	18%	−8%	24%

References

1. Tompalski, P.; Coops, N.; White, J.; Goodbody, T.; Hennigar, C.; Wulder, M.; Socha, J.; Woods, M. Estimating changes in forest attributes and enhancing growth projections: A review of existing approaches and future directions using airborne 3D point cloud data. Curr. For. Rep.; 2021; 7, pp. 1-24. [DOI: https://dx.doi.org/10.1007/s40725-021-00135-w]

2. Burkhart, H.; Tomé, M. Modeling Forest Trees and Stands; Springer: Dordrecht, The Netherlands, 2012; 457.

3. Rahimzadeh-Bajgiran, P.; Hennigar, C.; Weiskittel, A.; Lamb, S. Forest potential site productivity mapping by linking remote-sensing-derived metrics to site variables. Remote Sens.; 2020; 12, 2056. [DOI: https://dx.doi.org/10.3390/rs12122056]

4. Carmean, W. Forest site-quality estimation using forest ecosystem classification in northwestern Ontario. Environ. Monit. Assess.; 1996; 39, pp. 493-508. [DOI: https://dx.doi.org/10.1007/BF00396164]

5. Hemingway, H.; Kimsey, M. Estimating forest productivity using site characteristics, multipoint measures, and a nonparametric approach. For. Sci.; 2020; 66, pp. 645-652. [DOI: https://dx.doi.org/10.1093/forsci/fxaa023]

6. Noordermeer, L.; Gobakken, R.; Næsset, E.; Bollandsås, O. Predicting and mapping site index in operational forest inventories using bitemporal airborne laser scanner data. For. Ecol. Manag.; 2020; 457, 117768. [DOI: https://dx.doi.org/10.1016/j.foreco.2019.117768]

7. Gopalakrishnan, R.; Kauffman, J.; Fagan, M.; Coulston, J.; Thomas, V.; Wynne, R.; Fox, T.; Quirino, V. Creating Landscape-scale site index maps for the southeastern US is possible with airborne LiDAR and Landsat imagery. Forests; 2019; 10, 234. [DOI: https://dx.doi.org/10.3390/f10030234]

8. Tompalski, P.; Coops, N.; White, J.; Wulder, M.; Pickell, P. Estimating forest site productivity using airborne laser scanning data and Landsat time series. Can. J. Rem Sen.; 2015; 41, pp. 232-245. [DOI: https://dx.doi.org/10.1080/07038992.2015.1068686]

9. Carmean, W.; Hazenberg, H.; Deschamps, K. Polymorphic site index curves for black spruce and trembling aspen in northwest Ontario. For. Chron.; 2006; 82, pp. 231-242. [DOI: https://dx.doi.org/10.5558/tfc82231-2]

10. Sharma, M.; Parton, J. Climatic effects on site productivity of red pine plantations. For. Sci.; 2018; 64, pp. 544-554. [DOI: https://dx.doi.org/10.1093/forsci/fxy013]

11. Sharma, M.; Parton, J. Analyzing and modelling effects of climate on site productivity of white spruce plantations. For. Chron.; 2018; 93, pp. 173-182. [DOI: https://dx.doi.org/10.5558/tfc2018-026]

12. Sharma, M.; Parton, J.; Woods, M.; Newton, P.; Penner, M.; Wang, J.; Stinson, A.; Bell, W. Ontario’s forest growth and yield modelling program: Advances resulting from the Forestry Research Partnership. For. Chron.; 2008; 84, pp. 694-703. [DOI: https://dx.doi.org/10.5558/tfc84694-5]

13. Sharma, M.; Subedi, N.; Ter-Mikaelian, M.; Parton, J. Modeling climatic effects on stand height/site index of plantation-grown jack pine and black spruce trees. For. Sci.; 2015; 61, pp. 25-34. [DOI: https://dx.doi.org/10.5849/forsci.13-190]

14. Sharma, M.; Reid, D. Stand height/site index equations for jack pine and black spruce trees grown in natural stands. For. Sci.; 2017; 64, pp. 22-40. [DOI: https://dx.doi.org/10.5849/FS-2016-133]

15. Sharma, M. Climate effects on black spruce and trembling aspen productivity in natural origin mixed stands. Forests; 2022; 13, 430. [DOI: https://dx.doi.org/10.3390/f13030430]

16. Subedi, N.; Sharma, M. Evaluating height-age determination methods for jack pine and black spruce plantations using stem analysis data. North. J. Appl. For.; 2010; 27, pp. 50-55. [DOI: https://dx.doi.org/10.1093/njaf/27.2.50]

17. Subedi, N.; Sharma, M. Individual-tree diameter growth model for black spruce and jack pine plantations in northern Ontario. For. Ecol. Manag.; 2011; 261, pp. 2140-2148. [DOI: https://dx.doi.org/10.1016/j.foreco.2011.03.010]

18. Subedi, N.; Sharma, M. Climate-diameter growth relationships of black spruce and jack pine trees in boreal Ontario, Canada. Glob. Change Biol.; 2013; 19, pp. 505-516. [DOI: https://dx.doi.org/10.1111/gcb.12033] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23504788]

19. Solberg, S.; Kvaalen, H.; Puliti, S. Age-independent site index mapping with repeated single-tree airborne laser scanning. Scand. J. For. Res.; 2019; 34, pp. 763-770. [DOI: https://dx.doi.org/10.1080/02827581.2019.1616814]

20. OMNR. Ontario Forest Resources Inventory Photo Interpretation Specifications. Ont. Min. Nat. Res.; 2017.

21. Bjelanovic, E.; Comeau, P.; White, B. High resolution site index prediction in boreal forests using topographic and wet areas mapping attributes. Forests; 2018; 9, 113. [DOI: https://dx.doi.org/10.3390/f9030113]

22. Coops, N. Characterizing forest growth and productivity using remotely sensed data. Curr. For. Rep.; 2015; 1, pp. 195-205. [DOI: https://dx.doi.org/10.1007/s40725-015-0020-x]

23. Racine, E.; Coops, N.; St-Onge, B.; Bégin, J. Estimating forest stand age from LiDAR-derived predictors and nearest neighbor imputation. For. Sci.; 2014; 60, pp. 128-136. [DOI: https://dx.doi.org/10.5849/forsci.12-088]

24. Schumacher, J.; Hauglin, M.; Astrup, R.; Breidenbach, J. Mapping forest age using National Forest Inventory, airborne laser scanning, and Sentinel-2 data. For. Ecosyst.; 2020; 7, 60. [DOI: https://dx.doi.org/10.1186/s40663-020-00274-9]

25. Wylie, R.; Woods, M.; Dech, J. Estimating stand age from airborne laser scanning data to improve models of black spruce wood density in the boreal forest of Ontario. Remote Sens.; 2019; 11, 2022. [DOI: https://dx.doi.org/10.3390/rs11172022]

26. Zhu, Z. Change detection using Landsat time series: A review of frequencies, preprocessing, algorithms, and applications. ISPRS J. Photogramm. Remote Sens.; 2017; 130, pp. 370-384. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2017.06.013]

27. Tompalski, P.; Coops, N.; White, J.; Wulder, W. Augmenting site index estimation with airborne laser scanning data. For. Sci.; 2015; 61, pp. 861-873. [DOI: https://dx.doi.org/10.5849/forsci.14-175]

28. Socha, J.; Hawrylo, P.; Sterenczak, K.; Miscicki, S.; Tyminska-Cazbanska, L.; Mlocek, W.; Gruba, P. Assessing the sensitivity of site index models developed using bi-temporal airborne laser scanning data to different top height estimates and grid cell sizes. In. J. Appl. Earth Obs. Geoinf.; 2020; 91, 102129. [DOI: https://dx.doi.org/10.1016/j.jag.2020.102129]

29. Véga., C. St-Onge. Mapping site index and age by linking a time series of canopy height model with growth curves. For. Ecol. Mgmt.; 2009; 257, pp. 951-959. [DOI: https://dx.doi.org/10.1016/j.foreco.2008.10.029]

30. Persson, J.; Fransson, J. Estimating site index from short-term TanDEM-X canopy height models. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2016; 9, pp. 3598-3606. [DOI: https://dx.doi.org/10.1109/JSTARS.2016.2563158]

31. Woods, M.; Pitt, D.; Penner, M.; Lim, K.; Nesbitt, D.; Etheridge, D.; Treitz, P. Operational implementation of a LiDAR inventory in Boreal Ontario. For. Chron.; 2011; 87, pp. 512-528. [DOI: https://dx.doi.org/10.5558/tfc2011-050]

32. Bilyk, A.; Pulkki, R.; Shahi, C.; Larocque, G. Development of the Ontario Forest Resources Inventory: A historical review. Can. J. For. Res.; 2020; 51, pp. 198-209. [DOI: https://dx.doi.org/10.1139/cjfr-2020-0234]

33. Gluckman, J. Design of the processing chain for a high-altitude, airborne, single-photon lidar mapping instrument. Proceedings of the Laser Radar Technology and Applications XXI; Baltimore, MD, USA, 17–21 April 2016; Volume 9832, 983203. [DOI: https://dx.doi.org/10.1117/12.2219760]

34. Riofrio, J.; White, J.; Tompalski, P.; Coops, N.; Wulder, M. Harmonizing multi-temporal airborne laser scanning point clouds to derive periodic annual height increments in temperate mixedwood forests. Can. J. For. Res.; 2022; 52, pp. 1334-1352. [DOI: https://dx.doi.org/10.1139/cjfr-2022-0055]

35. Roussel, J.; Auty, D.; Coops, N.C.; Tompalski, P.; Goodbody, T.R.; Meador, A.S.; Bourdon, J.; 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]

36. Roussel, J.; Auty, D. Airborne LiDAR Data Manipulation and Visualization for Forestry Applications; R Package Version 4.0.3; CRAN: 2023. Available online: https://cran.r-project.org/package=lidR (accessed on 1 June 2023).

37. Yang, Y.; Monserud, R.A.; Huang, S. An evaluation of diagnostic tests and their role in validating forest biometric models. Can. J. For. Res.; 2004; 34, pp. 619-629. [DOI: https://dx.doi.org/10.1139/x03-230]

38. Nigh, G. Species-independent height-age models for British Columbia. For. Sci.; 2001; 47, pp. 150-157.

39. White, J.; Wulder, M.; Varhola, A.; Vastaranta, M.; Coops, N.; Cook, B.; Pitt, D.; Woods, M. A Best Practices Guide for Generating Forest Inventory Attributes from Airborne Laser Scanning Data Using an Area-Based Approach (Version 2.0). Natural Resources Canada Information Report FI-X_010. 2013; Available online: https://cfs.nrcan.gc.ca/publications/download-pdf/34887 (accessed on 1 June 2023).

40. Treitz, P.; Lim, K.; Woods, M.; Pitt, D.; Nesbitt, D.; Etheridge, D. LiDAR sampling density for forest resources inventories in Ontario, Canada. Remote Sens.; 2012; 4, pp. 830-848. [DOI: https://dx.doi.org/10.3390/rs4040830]

41. McRoberts, R.; Hahn, J.; Hefty, G.; Van Cleve, J. Variation in forest inventory field measurements. Can. J. For. Res.; 1994; 24, pp. 1766-1770. [DOI: https://dx.doi.org/10.1139/x94-228]

42. Luoma, V.; Saarinen, N.; Wulder, M.; White JVastaranta, M.; Holopainenen, M.; Hyppä, J. Assessing precision in conventional field measurements of individual tree attributes. Forests; 2017; 8, 38. [DOI: https://dx.doi.org/10.3390/f8020038]

43. Wang, Y.; Lehtomäki, N.; Liang, X.; Pyörälä, J.; Kukko, A.; Jaakkola, A.; Liu, J.; Feng, Z.; Chen, R.; Hyppä, 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]

44. Uhlig, P.; Harris, A.; Craig, G.; Bowling, C.; Chambers, B.; Naylor, B.; Beemer, G. Old Growth Forest Definitions for Ontario; Queen’s Printer for Ontario: Toronto, ON, Canada, 2001; 53p.

45. Wulder, W.; Bater, C.; Coops, N.; Hilker, R.; White, J. The role of LiDAR in sustainable forest management. For. Chron.; 2008; 84, pp. 807-826. [DOI: https://dx.doi.org/10.5558/tfc84807-6]

46. Parresol, B.; Scott, D.; Zarnoch, S.; Edwards, L.; Blake, J. Modeling forest site productivity using mapped geospatial attributes with a South Carolina Landscape, USA. For. Ecol. Mgmt.; 2017; 406, pp. 196-207. [DOI: https://dx.doi.org/10.1016/j.foreco.2017.10.006]