Mapping Dominant Boreal Tree Species Groups by Combining Area-Based and Individual Tree Crown LiDAR Metrics with Sentinel-2 Data

Abstract

Airborne light detection and ranging (LiDAR) data are increasingly used to inform sustainable forest management practices. Information about species composition is needed for a range of applications; however, commonly used area-based summaries of LiDAR data are limited to accurately differentiate tree species. The objective of this study was to map dominant species groups across a large (>580,000 ha) boreal forest by combining area-based and individual tree metrics derived from single photon LiDAR data with multispectral information derived from Sentinel-2 imagery. Classification of the forest into jack pine (Pinus banksiana), black spruce (Picea mariana), mixed conifer, mixedwood, and hardwood species groups resulted in an overall accuracy of 67.8% compared with field data. A more generalized classification into softwood, hardwood, and mixedwood had an overall accuracy of 83.2%. The reflectance in the red edge region of the electromagnetic spectrum (λ = 740 nm), the area and volume of tree crowns, and the cumulative distribution of LiDAR returns in the canopy were particularly important variables to discriminate between species groups. Wall-to-wall predictions of species groups based on remotely sensed data—as derived herein—could provide a spatially-detailed alternative to stand-level species composition information traditionally derived from photo-interpreted aerial imagery.

RÉSUMÉ

Les données de balayage laser aéroporté (LiDAR) sont de plus en plus utilisées pour soutenir la gestion durable des forêts. Connaître la composition des espèces d’arbres est nécessaire pour de multiples applications; cependant, des limites existent concernant l’utilisation de mesures régionales de données LiDAR pour différencier les espèces de manière précise. L’objectif de cette étude est de cartographier les espèces d’arbres dominantes à travers une grande (>580,000 ha) forêt boréale en combinant à la fois des mesures LiDAR régionales et à l’échelle de l’arbre avec des images satellites multispectrales prises par Sentinel-2. La classification de la forêt en pin gris (Pinus banksiana), épinette noire (Picea mariana), conifères mixtes, forêt mixte, et feuillus mixtes a résulté en une précision de 67.8% en comparaison avec des données de terrain. Une précision de 83.2% a également été obtenue avec une classification plus générale soit conifères, forêt mixte et feuillus. La réflectance dans la région du spectre électromagnétique correspondant au point d’inflexion du rouge (λ = 740 nm), la superficie et le volume des couronnes d’arbres et la distribution cumulative des points LiDAR dans les parties supérieures de la canopée étaient des variables particulièrement importantes pour différencier les espèces. La cartographie des espèces dominantes à partir de données de télédétection—telles qu’elles sont dérivées ici—pourrait être une alternative aux estimations à l’échelle du peuplement issues traditionnellement de la photo-interprétation d’imagerie aérienne.

Introduction

Sustainable operational forest management requires detailed spatial information about forest composition and structure (Wulder et al. 2008). In particular, information related to species composition is needed to assess species-specific growing stocks and yield (e.g., Maltamo and Packalen 2014; Tompalski et al. 2018) and monitor tree species richness (Gamfeldt et al. 2013). Traditionally, assessment of species composition involved manual photo-interpretation of aerial imagery, from which stands of homogeneous forest type can be delineated, and field information obtained from plots sampled across the landscape (Leckie and Gillis 1995). However, the accuracy of these photo-interpreted species compositions is rarely assessed (Tompalski et al. 2021), and when they are, have been found to have a relatively low level of accuracy (54%–56% in Pinto et al. 2007, ∼36% in Thompson et al. 2007). Moreover, the accuracy of photo-interpretation depends on the complexity of the forest environment, as well as the experience of the photo-interpreter, the latter of which are increasingly rare. For these reasons, there has been a growing interest in assessing species composition using remotely sensed data. However, to date, remotely sensed applications have been unable to replicate the detailed species compositions associated with photo-interpreted forest inventories, and scaling up species mapping approaches across large-areas remains a significant challenge (Fassnacht et al. 2016).

Technologies, such as light detection and ranging (LiDAR) data acquired from airborne platforms, also referred to as airborne laser scanning (ALS), provide a 3-dimensional point cloud representation of the landscape from which information about forest structure can be extracted. As a result, the past decade has seen the widespread adoption of ALS information into the forestry inventory cycle to improve the precision and accuracy of forest resources inventory (White et al. 2016; Coops et al. 2021). To date, the most common approach to utilize ALS in a forest inventory is with an area-based approach (ABA), which consists of developing predictive models of forest attributes obtained from ground plot information using plot-level LiDAR summary metrics derived from the vertical distribution of the point cloud as predictor variables (Næsset 2002; White et al. 2013; Wulder et al. 2014). While ALS provides robust area-based estimates of forest attributes, such as Lorey’s height, basal area, stem volume, and above-ground biomass, which are strongly related to the structural information provided by the point cloud, other attributes, such as species composition are more challenging to derive accurately using summary metrics of LiDAR point cloud heights (Fassnacht et al. 2016).

Species-related traits, such as the pigments, water content and internal structure of leaves, branching, foliage density, and crown shape influence the distribution and intensity of returns in ALS point clouds and the spectral response measured by optical sensors, particularly in parts of the visible (λ = 400 nm–700 nm), near-infrared (λ = 700 nm–1300 nm) and shortwave infrared (λ = 1300 nm–2500 nm) regions of the electromagnetic spectrum (Fassnacht et al. 2016). Research examining the prediction of tree species with ALS data has largely focused on the individual tree level by using individual tree detection (ITD) and individual tree crown (ITC) segmentation algorithms from which tree-specific crown shape and density metrics can be derived (Holmgren and Persson 2004; Michałowska and Rapiński 2021). Nevertheless, structural information alone generally has a lower potential than spectral information to classify tree species (Fassnacht et al. 2016). Spectral information is often integrated by combining ALS data and high-resolution multispectral or hyperspectral imagery (Holmgren and Persson 2004; Dalponte et al. 2014; Holmgren et al. 2008; Coops et al. 2004) or by using novel multispectral ALS sensors which emit laser pulses at different wavelengths (Axelsson et al. 2018; Budei et al. 2018; Mikko Kukkonen et al. 2019; Prieur et al. 2021). As tree-level approaches are prone to bias due to issues with segmentation (Kaartinen et al. 2012), approaches that combine ABA and ITC, known as semi-ITC, have also been explored (Breidenbach et al. 2010). In areas where species-specific attribute estimates are required for the inventory (e.g., Finland), non-parametric supervised learning methods, such as k-nearest neighbors algorithms (k-NN) are commonly used for area-based estimates (e.g., Lindberg et al. 2013; Maltamo et al. 2021).

Area-based or stand-level predictions of species information have been achieved using spectral data from medium resolution multispectral satellite imagery, such as Sentinel-2 or Landsat exclusively (Astola et al. 2019; Grabska et al. 2019; Persson et al. 2018; Breidenbach et al. 2021; Hemmerling et al. 2021) but the integration of these data with ALS in an area-based context has been less explored (Fassnacht et al. 2016). Therefore, the objectives of this research were: (i) to produce area-based predictions of dominant species groups using a combination of standard LiDAR height metrics, aggregated metrics derived from ABA and ITC methods, and Sentinel-2 imagery across a very large (>580,000 ha) boreal forest; (ii) to explore and assess the ability of these three sets of metrics to differentiate between species groups; and (iii) to discuss implications of using this approach for large-area operational forest management at the stand level.

Data and methods

Study area

The study area was the Romeo Malette Forest, located in Ontario, Canada (48.48° N, 81.33° W). The Romeo Malette Forest consists of 630,000 ha, of which 582,430 ha is forested land. It is part of the Boreal Shield Ecozone of Canada (Figure 1). The area is characterized by long cold winters and short warm summers with monthly average daily temperature ranging between −16.8 °C and 17.5 °C and monthly precipitation ranging between 41.3 mm and 90.9 mm (1981–2010 normals). The terrain is dominated by relatively flat to moderately rolling topography (280 m–450 m above sea level) and dominant tree species include black spruce (Picea mariana), jack pine (Pinus banksiana), white spruce (Picea glauca), trembling aspen (Populus tremuloides) and paper birch (Betula papyrifera).

Figure 1. Map of the study area (A) and location within the province of Ontario and the Canadian boreal shield ecozone (B).

The Romeo Malette Forest is actively managed for timber production and fiber procurement. Management activities are guided by Ontario’s Crown Forest Sustainability Act and related policies, and are based on a 10-year forest management plan, which includes a detailed description and mapping of logging and renewal operations and the development of planned road access. The current plan covers the years 2019–2029 and management activities are reported or projected based on forest units defined from tree species composition groupings (Parton et al. 2006). Annual work schedules identify the location of operations scheduled for implementation during each year and provide information on the areas scheduled for harvest, road construction/maintenance, renewal, tending, and protection operations. The Romeo Malette Forest has been certified under the National Boreal Standard of the Forest Stewardship Council (FSC) since 2006.

Field data

A set of 244 circular sample plots (400 m²; radius = 11.28 m), selected following a structurally guided sampling approach driven by ALS metrics were sampled in 2018 (Queinnec et al. 2022). Plot center coordinates were obtained with the EOS Arrow 100 GPS receiver by acquiring a minimum of 600 points per plot center and the use of differential post-processing to achieve sub-meter position accuracy. The height, species, and diameter at the breast height (DBH) of each tree with DBH ≥ 7.1 cm located within the plot boundary were recorded at each plot. The vertical stratification of tree crowns was also assessed by the field crews to identify suppressed trees or trees located within or above the dominant canopy layer. Additionally, 53 of the plots were stem mapped and the distance and azimuth from the plot center to each stem (DBH ≥ 7.1 cm) located in the plot boundary were recorded using an azimuth and range finder instrument (TruPulse^® 360). For leaning trees, the distance and azimuth between the stem and the vertical projection of the crown top on the ground were also measured. Trees with a crown located inside the plot boundary but with a stem outside were also mapped. The cartesian coordinates of tree tops were determined from the polar coordinates (distance and azimuth of tree tops) measured from the plot center.

For each ground plot, the trees located within—or above—the dominant canopy layer were used to derive dominant species groups as described in Table 1. Two functional groupings were established. The first was more general and grouped plots according to the dominance of forest types: softwood, mixedwood, and hardwood dominated groups. The second functional grouping distinguished the dominance of the two most common commercial tree species: black spruce and jack pine. Plots were assigned to five functional groups: black spruce-dominated, jack pine-dominated, mixed conifer, mixedwood, and hardwood. Similar species groupings were used in previous forest inventories at the Romeo Malette Forest (Murray Woods et al. 2011; Queinnec et al. 2022). The capacity of the developed approach using Sentinel-2 and ALS data to accurately predict each of these groupings was assessed.

Table 1. Stratification of ground plots into 3 or 5 functional groups based on the species of the trees positioned in—or above—the dominant canopy layer.

Species group	Number of plots	Number of stem-mapped plots	Description
3 Functional groups
Softwood	97	26	The proportion of softwood^a tree species is >70%
Hardwood	114	25	The proportion of hardwood^b tree species is >70%
Mixedwood	33	2	Plots not classified as Softwood or Hardwood. The proportion of softwood and hardwood tree species in these plots is between 30 and 70%
5 Functional groups
Black Spruce dominated	43	11	The proportion of softwood tree species is >70 and >50% of the trees are black spruce
Jack Pine dominated	40	10	The proportion of softwood tree species is >70 and >50% of the trees are jack pine
Mixed Conifer	12	4	The proportion of softwood tree species is >70% but the proportion of black spruce or jack pine is not >50%
Mixedwood	35	3	Plots not classified as other species groups. The proportion of softwood and hardwood trees in these plots is between 30 and 70%
Hardwood	114	25	The proportion of hardwood trees is >70%

The number of trees in each plot was used to derive proportions.

^a Softwood species include jack pine, black spruce, balsam fir (Abies balsamea), white spruce, American larch (Larix laricina), and eastern white cedar (Thuja occidentalis).

^b Hardwood species include trembling aspen, paper birch, and red maple (Acer rubrum).

LiDAR data and pre-processing

Wall-to-wall ALS data were collected across the Romeo Malette Forest in June and July 2018 with the Leica SPL100 as part of an acquisition across the forest management units of Ontario (Bilyk et al. 2021). The SPL100 is a single-photon LiDAR instrument that emits laser pulses in the green part of the electromagnetic spectrum (λ = 532 nm). Each beam transmitted by the laser is split into a 10 × 10 array of beamlets and the energy returned from each beamlet is recorded with an array of cells that are activated by the arrival of single photons. This differs from conventional discrete-return ALS systems where detector cells are activated when the amount of energy is above a detection threshold (Harding 2018). The SPL100 can be operated at high altitude (typically 2000 m–4000 m above ground level) with a high frequency of low-energy pulses, which gives them an operational advantage over conventional ALS systems when high-density point clouds (typically 30 points/m²–50 points/m²) are required for large ALS acquisition programs (Yu et al. 2020). Previous studies conducted in boreal or mixedwood forests have demonstrated that the SPL100 provided area-based forest attribute estimates with an accuracy comparable to what can be expected from conventional linear-mode LiDAR systems (Queinnec et al. 2022; White et al. 2021; Yu et al. 2020; Wästlund et al. 2018).

For the Romeo Malette Forest acquisition, the SPL100 was flown onboard an aircraft at an average altitude of 3800 m and a flying speed of 90 m/s for an average swath width of 2000 m on the ground. The point cloud was geolocated, denoised, and classified into ground and vegetation returns by the data provider according to the approach outlined in Gluckman (2016). This resulted in an average point cloud density of 40.8 points/m² across the study area.

The point cloud was normalized to heights above ground using a triangular irregular network of the ground returns. A 25-cm resolution canopy height model (CHM) was generated from the normalized point cloud using the highest return within the grid of each cell. To densify the point cloud and reduce the number of pits, each return was considered as a disk of 30 cm diameter, which corresponds to the actual footprint of each beamlet on the ground (beam divergence of 0.08 mrad and altitude of 3800 m above ground). The CHM was then smoothed using a Gaussian filter (σ = 1) on a 5 × 5 pixel moving window.

LiDAR metrics

Two sets of LiDAR summary metrics were calculated from the normalized ALS point cloud and derived CHM for each of the 244 ground plots and across the entire study area at a resolution of 20 × 20 m matching the ground plot area of 400 m² (Table 2). The lidR package (Roussel et al. 2020) implemented in R (R Core Team 2022) was used to generate the ALS metrics. The first set of metrics included standard area-based metrics calculated from the height distribution of the LiDAR returns, as well as the rumple index, which characterizes the surface roughness of the canopy. We also calculated the signal attenuation depth (SAD) of each return (Véga et al. 2016; Irwin et al. 2021) which indicates how far a return travels within the canopy. SAD is calculated for each return as the difference between the CHM and the return height and normalized to the CHM. SAD = 0 indicates a return on the outer canopy layer while SAD = 1 indicates a return on the ground. SAD was calculated only for returns where CHM > 2 m.

Table 2. Description of the predictor variables derived from the ALS data and Sentinel-2 imagery.

	Metric	Description
Area-based predictors	zq5, zq10, …, zq90, zq95	5th, 10th, … 90th, 95th percentiles of returns height >1.3 m classified as vegetation
	zsd, zskew, zkurt	Standard deviation, skewness, and kurtosis of returns height >1.3 m classified as vegetation
	pzabove2, pzabovemean	Percentage of first returns >2 m or > average returns height
	zpcum1, zpcum2, … zpcum9	Cumulative percentage of returns in layers (Woods et al. 2008)
	zentropy	Shannon diversity index quantifying the evenness of the distribution of LiDAR returns height
	zrumple	Ratio of canopy outer surface area to ground surface area
	depth_q25, depth_q50, depth_q75	25th, 50th, and 75th percentile of signal attenuation depth of all returns
ITC predictors	n_itd	Number of detected tree tops
	cr_a_sum; cr_a_med; cr_a_sd	Sum, median, or standard deviation of segmented crown area
	cr_vol_sum; cr_vol_med; cr_vol_sd	Sum, median, or standard deviation of the product of segmented crown area and corresponding tree top height (Duncanson et al. 2015)
	cr_lor	Average of tree top height divided by corresponding segmented crown area (Duncanson et al. 2015)
Sentinel-2 predictors	s2_blue, s2_red, s2_green, s2_redge1, s2_redge2, s2_redge3, s2_nir, s2_siwr1, s2_swir2	Blue (band 2: 490 nm); Green (band 3: 560 nm); Red (band 4: 665 nm); Red-edge 1 (band 5: 705 nm); Red-edge 2 (band 6: 740 nm); Red-edge 3 (band 7: 783 nm); NIR (band 8 A: 865 nm); SWIR1 (band 11: 1610 nm); SWIR2 (band 12: 2190 nm)

The second set of metrics was derived by summarizing the distribution of ITD and ITC metrics at the plot-level or 20 × 20 m pixel-level using the sum, median, and standard deviation of values (Table 2). The ITD was performed using the local maximum filter algorithm implemented in lidR (Popescu and Wynne 2004; Roussel et al. 2020). The algorithm uses a circular moving window to detect local maxima in the CHM and a size parameter controls the diameter of the window. We varied the size parameter between 0.5 m and 5 m at 0.5 m intervals to assess the influence of the window on the ITD performance window shapes. The individual tree tops detected by the ITD algorithm were then used as seed points in the marker-controlled watershed segmentation algorithm implemented in the R package imager to detect ITCs. The marker-controlled watershed algorithm was selected as it is parameter-free, not designed for particular crown shapes, and computationally efficient, which are important considerations for large areas with varying species composition, such as the Romeo Malette Forest. The algorithm inverts the CHM and sets the detected tree tops as local minima. A region is then propagated around each tree top by filling the inverted CHM. To avoid over-propagation, we limited the region-growing algorithm to the pixels of the CHM that were not below 1 meter of the detected tree top height. That way, we expected trees with a conical crown shape to grow in smaller regions than trees with a flatter and broader crown.

Sentinel-2 imagery

We obtained 44 Sentinel-2A and Sentinel-2B level-1C (top-of-atmosphere reflectance) scenes acquired across the RMF between June 1 and September 30, 2018 and processed them to level-2A (bottom-of-atmosphere reflectance) using the R package sen2r (Ranghetti et al. 2020). This R package streamlines the downloading and pre-processing of Sentinel-2 imagery using the atmospheric correction algorithm Sen2Cor (Main-Knorn et al. 2017). All pixels classified as missing values, shadows, clouds, or cirrus by the level-2A scene classification map were masked out and scenes were projected and resampled to align with the 20 × 20 m SPL100 metrics grid using nearest neighbor interpolation. The Sentinel-2 bands used in this study had a spatial resolution of either 10 × 10 m or 20 × 20 m so downscaling to the resolution of the gridded ALS metrics was not required. A cloud-free composite product was then produced by averaging the pixel values across all available cloud-free observations for each pixel. Spectral bands listed in Table 2 were selected as predictor variables for species group classification.

Classification of dominant species groups

Random forest models were used to classify the dominant species groups (Table 1) using the area-based, ITC and Sentinel-2 predictors (Table 2). The conditional random forest algorithm implemented in the function cforest of the party package in R (Strobl et al. 2007; 2008; Hothorn et al. 2006) was selected due to its ability to compute an unbiased measure of variable importance (see below). Predictor variables with a pairwise Pearson correlation coefficient >0.9 were identified and from each pair of variables, only the one with the lowest average correlation among the other predictor variables was retained. We also performed different scenario runs, as described in Table 3, to assess the influence of different sets of predictor variables on the accuracy of the species group classification.

Table 3. Description of the predictor variables included in the different scenario runs.

Scenario run	Abbreviation
Sentinel-2 only (s2_only)	s2_blue, s2_green, s2_red, s2_red_edge_2, s2_swir2
LiDAR height only (z_only)	zq5, zq85, zsd, zkurt, zskew, zentropy, pzabove2, pzabovemean, zpcum1, zpcum2, zpcum8, zpcum9, sad_q25, sad_q50, sad_q75
Aggregated ITD and ITC metrics only (itd_crown)	n_itd, cr_lor, cr_a_sum, cr_a_med, cr_a_sd, cr_vol_sum, cr_vol_sd
All LiDAR metrics (all_LiDAR)	Metrics included in scenario runs 2 and 3
All predictor variables (all)	Metrics included in all the other scenario runs

Refer to Table 1 for a description of the metrics.

Accuracy assessment of the ITD and area-based species classification

The performance of the ITD was assessed using the 53 stem-mapped plots. Detected trees were automatically matched to their nearest field-mapped trees based on the 3-dimensional distance between their tree tops. Only pairs of detected and mapped trees with a horizontal distance ≤5 m and height difference ≤3 m were considered as valid matches, or true positives (TP). Since the ITD algorithm only detects trees in the upper canopy layer, only field trees that were in or above the dominant canopy layer and that were not matched to any detected tree top were considered as false negatives (FN). All remaining detected tree tops not matched to any field tree were considered as false positives (FP). The accuracy of the ITD algorithm was assessed for each set of parameters by calculating the recall (Equation 1(1) $$\mathit{recall}=\frac{TP}{TP+FN}$$(1) ), precision (Equation 2(2) $$\mathit{precision}=\mathrm{ }\frac{TP}{TP+FP}$$(2) ), and F-score (Equation 3(3) $$F-\mathit{score}=2\times \frac{\mathit{recall}\mathrm{ }\times \mathit{precision}}{\mathit{recall}+\mathit{precision}}$$(3) ). Recall, or producer’s accuracy, indicates the proportion of correctly detected trees out of all trees mapped on the field, whereas precision (user’s accuracy) indicates the proportion of correctly detected trees out of all trees detected by the ITD algorithm. The F-score combines recall and precision measures by calculating their harmonic mean, which penalizes more low or unbalanced recall and precision values than the arithmetic or geometric means. (1) $$\mathit{recall}=\frac{TP}{TP+FN}$$(1) (2) $$\mathit{precision}=\mathrm{ }\frac{TP}{TP+FP}$$(2) (3) $$F-\mathit{score}=2\times \frac{\mathit{recall}\mathrm{ }\times \mathit{precision}}{\mathit{recall}+\mathit{precision}}$$(3)

The accuracy of the random forest classification models was assessed using repeated k-fold cross-validation (k = 10; repeated 5 times). Due to the imbalance of species groups in the observations, the training folds were randomly sampled so that all species groups had the same number of observations as the least frequent species group. The testing folds remained unchanged. The accuracy of the models was assessed using overall accuracy (OA; Equation 4(4) $$OA=\frac{TP}{TP+FN+FP}$$(4) ), recall (Equation 1(1) $$\mathit{recall}=\frac{TP}{TP+FN}$$(1) ), precision (Equation 2(2) $$\mathit{precision}=\mathrm{ }\frac{TP}{TP+FP}$$(2) ), and F-score (Equation 3(3) $$F-\mathit{score}=2\times \frac{\mathit{recall}\mathrm{ }\times \mathit{precision}}{\mathit{recall}+\mathit{precision}}$$(3) ) from the predictions of the repeated folds. The OA is a global measure of accuracy that indicates the proportion of observations that were correctly classified across all species groups whereas recall, precision, and F-score are calculated for each species group. (4) $$OA=\frac{TP}{TP+FN+FP}$$(4)

The cross-validation was also used to optimize the parameter mtry that controls the number of predictor variables that are randomly selected to split the data in the nodes of each tree. The mtry value of the model with the lowest overall accuracy across the folds was used in the final model for both the 3-group and 5-group classification.

Assessing the importance of predictor variables

The importance of the predictor variables to predict the species groups was assessed using the conditional variable importance measure from the conditional random forest model (Strobl et al. 2008). Variable importance in the random forest is assessed by examining the effect of breaking the relationship between a predictor variable and the response, by permuting the predictor variable values, on the prediction accuracy (measured as OA). The difference between prediction OA before and after permutation of a predictor variable values is used as a measure of variable importance. Strobl et al. (2007) showed that when predictor variables are of different nature or vary in their scale of measurement, which is the case herein, the permutation scheme used in random forest tends to be biased toward certain variables. The conditional variable importance used herein is derived from a novel conditional permutation scheme, which enables to derive unbiased variable importance (Strobl et al. 2008). Conditional importance measures were linearly scaled between 0 (least important variable) and 100 (most important variable).

Random forest variable importance provides a measure of the overall contribution of each predictor variable to the species groups classification model but does not assess if some variables are more useful to discriminate between specific species groups than others. Therefore, we performed a pairwise comparison of the predictor variable values between species groups to gain insight into why some pairs of species groups were more or less separable. For each of the 15 most important predictor variables, we used the Kruskal-Wallis test, a non-parametric equivalent to a one-way analysis of variance, to test the null hypothesis that the mean rank of the variable is equal across the species groups. When a significant difference (p-value < 0.05) was observed across the species group, we used the post-hoc Dunn’s test with Holm’s p-value adjustment method to examine pairwise differences. That way, we were able to examine if some variables had more predictive power for some species groups compared to others.

Comparing predicted species group to existing photo-interpreted inventory

We compared the estimates of species groups derived herein to stand-level species information available from the most up-to-date photo-interpreted forest resources inventory information available at the Romeo Malette Forest which was initially established with manual photo-interpretation of multispectral aerial imagery acquired in 2005. About 20% of the stands initially established were updated using new aerial photographs or field-based surveys regularly until 2014 to account for harvesting, natural depletion, silvicultural treatments, regeneration, or changes in stand development stages. The objective of this comparison was to explore the classification results across a broader geographic extent than the cross-validation at the plot-level provided, and at the stand-level which is the management unit used for decision making. However, we did not consider the photo-interpreted data as a reference to assess the accuracy of our classification since their accuracy is poorly defined (Pinto et al. 2007; Thompson et al. 2007).

We selected 43,806 productive forest stands (median area of 5.6 ha) with SPL100 estimates of canopy height over 5 m (zq90 > 5 m) and over 50% canopy cover (pzabove2 > 50%) to remove stands recently harvested or disturbed. In each stand, an estimate of overstory species composition as well as the percentage of crown covers each tree species occupies was available from the photo-interpreted inventory, and this information was used to classify each stand into the five functional species groups listed in Table 1. In each stand, we calculated the most frequently predicted species group and compared it to the species group derived from the photo-interpreted inventory.

Results

Individual tree detection and crown segmentation

The accuracy of the ITD varied depending on the diameter of the local maxima filter window and species groups. The highest median F-score was observed for a circular window diameter of 1.5 m for black spruce dominated (F_median = 0.82) and mixed conifer (F_median = 0.82), 2 m for mixedwood (F_median = 0.72), 2.5 m for jack pine dominated (F_median = 0.85), and 3 m for hardwood (F_median = 0.75) species groups (Figure 2). In all species groups combined, the highest median F-score of 0.72 was obtained for window diameters of 2 m, 2.5 m, and 3 m, and the window diameter of 2 m was retained to apply the ITD algorithm across the Romeo Malette Forest. Median values of recall, precision, and F-score for each species group and window diameter are presented in Supplementary Table 1.

Figure 2. Influence of the local maxima window diameter on the F-score of the individual tree detection algorithm by species group. The window diameter for which the highest median F-score was observed is highlighted in green.

Classification of species groups

The classification of ground plots into the three functional species groups softwood, mixedwood, and hardwood resulted in an overall accuracy ranging between 65.98% for the model including Sentinel-2 bands only, and 83.20% for the model including all predictor variables (Figure 3). When including all predictor variables, the softwood species groups had a precision and recall of 92.13% and 91.75% (F-score = 91.94%), and the hardwood species group had a precision and recall of 90.58% and 82.63% (F-score = 86.42%), respectively. The accuracy of the classification was poorer for the mixedwood species group with a precision of 45.62% and a recall of 60.00% (F-score = 51.83%). The confusion matrix is presented in Table 4 and illustrated in Figure 4A and recall, precision, and F-score for each species group and model scenario are included in Supplementary Figures 1 and 2.

Figure 3. Comparison of the overall accuracy of the random forest classification model of the ground plots into softwood, mixedwood, and hardwood species groups (3 Species Groups) or jack pine dominated, black spruce dominated, mixed conifer, mixedwood, and hardwood species groups (5 Species Groups).

Figure 4. Chord diagram visualizing the confusion matrix of the random forest species classification model into (A) softwood, mixedwood, and hardwood species group; or (B) jack pine dominated, black spruce dominated, mixed conifer, mixedwood, and hardwood species groups for the scenario run that included all predictor variables. Predictions from all cross-validation folds were aggregated and are included in this diagram. The predicted species groups are shown in the bottom half and the reference species groups are shown in the upper half. The predicted species groups are shown in the bottom half and the reference species groups are shown in the upper half. Links illustrate the number of observations that were correctly classified (link between same predicted and reference species groups) and the number of classification errors (link originating from a predicted species group going to a different reference species group).

Table 4. Confusion matrix showing the results of the broader, 3-group classification model using all predictor sets. Diagonal bold values indicate the number of observations for which the predicted species group correspond to the reference species group.

		Reference
		SW	MW	HW	Total	Precision
Predicted	SW	445	22	16	483	92.13%
	MW	35	99	83	217	45.62%
	HW	5	44	471	520	90.58%
	Total	485	165	570	1220
	Recall	91.75%	60.00%	82.63%

SW: softwood; MW: mixedwood; HW: hardwood.

Predictions of the repeated folds are included in the matrix.

When dividing the softwood species group into jack pine dominated, black spruce dominated, and mixed conifer groups, the overall accuracy of the model ranged between 46.22% for the model including Sentinel-2 bands only and 67.89% for the model including all the predictor variables (Figure 3). The hardwood species group remained the most accurate with a precision of 87.34% and a recall of 87.19% (F-score = 87.27%). Jack pine and black spruce species groups had a similar precision of 62.50% and 61.76%, respectively. The recall was, however, higher for black spruce (58.60%) than for jack pine (47.5%) dominated groups, which indicates that jack pine was misclassified as another species group more often than black spruce. Most classification errors in the jack pine and black spruce species groups occurred from the confusion of one with the other or with the mixed conifers species group. Similar to the findings of the classification into softwood, mixedwood, and hardwood, low confusion was observed between the coniferous species and the mixedwood and hardwood species. The confusion matrix is presented in Table 5 and illustrated in Figure 4B and recall, precision, and F-score for each species group and model scenario are included in Supplementary Figures 1 and 2.

Table 5. Confusion matrix for the five species groups classification model using all predictor sets. Diagonal bold values indicate the number of observations for which the predicted species group correspond to the reference species group.

		Reference
		JP	BS	MC	MW	HW	Total	Precision
Predicted	JP	95	47	3	6	1	152	62.50%
	BS	58	126	8	12	0	204	61.76%
	MC	26	34	33	14	22	129	25.58%
	MW	18	8	14	76	50	166	45.78%
	HW	3	0	2	67	497	569	87.34%
	Total	200	215	60	175	570	1220
	Recall	47.5%	58.60%	55.00%	43.43%	87.19%

JP: jack pine dominated; BS: black spruce dominated; MC: mixed conifer; MW: mixedwood; HW: hardwood.

Predictions of the repeated folds are included in the matrix.

Variable importance for species groups classification

The variable importance of the 3-groups and 5-group species classification models including all predictor variables are presented in Figure 5. For the most detailed 5 species groups classification model, the second red edge band of Sentinel-2 (s2_red_edge_2; λ = 740 nm), the cumulative percentile of height returns zpcum8 and the sum of segmented crown volumes (cr_vol_sum) showed markedly higher variable importance than other predictor variables for classifying species groups (Figure 5B). Six other variables, including four crown metrics, the 25th percentile of signal attenuation depth and zpcum9, had similarly scaled variable importance values ranging between 23% and 32%. The scaled variable importance values of the remaining variables were markedly lower (0%–14%). The most important metrics for the 5-group species classification model also generally had high variable importance in the 3-group model (Figure 5A). Out of the 9 variables highlighted above that had markedly higher importance in the 5-group model (23%–100%), all except zpcum9 were also ranked in the nine most important variables in the 3-group model.

Figure 5. Conditional random forest variable importance of the 3-group (A) and 5-group (B) species groups classification model including all the predictor variables. Raw importance values were scaled between 0 (least important variable) and 100 (most important variable).

Comparison of predictor variables across species groups

Not all species groups were highly distinct across the predictor variables considered (Figure 6). Hardwood was generally the species group that was most distinct from others, especially when compared to jack pine and black spruce dominated groups. Importance of metrics related to crown area and volume were significantly different across all species groups except when comparing black spruce and jack pine dominated groups with the mixed conifer group as well as when comparing jack pine and black spruce dominated groups. The reflectance values of Sentinel-2 spectral bands were generally significantly different when comparing jack pine and black spruce dominated groups with the other species groups, but not—or less—significantly, different when comparing mixedwood, mixed conifer, and hardwood groups. The SWIR2 band (λ = 2190 nm) reflectance value of Sentinel-2 was significantly different (p < 0.001) when comparing mixed conifer, mixedwood, and hardwood groups with the black spruce dominated group, but not significantly different when these groups were compared to the jack pine dominated group. The only predictor variables displaying significantly different values between jack pine and black spruce dominated groups were zpcum9 (p < 0.05), zpcum8 (p < 0.01), and the Sentinel-2 SWIR2 reflectance (p < 0.05).

Figure 6. Pairwise comparison of significant differences in predictor variable values across species groups. The blue and red color gradients indicate if the pairwise difference between the median values is negative or positive, respectively. JP: jack pine dominated; BS: black spruce dominated; MC: mixed conifers; MW: mixedwood; HW: hardwood.

Comparing predicted species group to existing photo-interpreted inventory

Figure 7 shows the distribution of the most frequently predicted species group in the stands of the photo-interpreted forest resource inventory. The black spruce dominated species group was the most frequently predicted species group in 86.2% of the stands dominated by black spruce according to the photo-interpreted forest resource inventory. In stands classified as jack pine-dominated, jack pine was the most frequently predicted species group in 45.7% of the stands, and black spruce was the most frequently predicted species group in 29.7% of the stands. The most frequently predicted species groups in stands classified as hardwood in the photo-interpreted inventory were hardwood (53.8%) and mixedwood (40.7%). Overall, we observed a negligible number of stands where mixed conifer was the most frequently predicted species group. Black spruce dominated (68.1%) and mixedwoods (28.7%) were the most frequently predicted species groups in stands classified as mixed conifer by the photo-interpreted inventory. Figures 8A–8D show some examples of stands and their species group predictions. A comparison against aerial imagery is also available in Supplementary Figure 3.

Figure 7. Chord diagram showing the distribution of the most frequently predicted species groups in the stands of the photo-interpreted inventory. The most frequently predicted species groups in each stand are shown in the bottom half and the stand-level species groups derived from photo-interpreted species composition are shown in the upper half.

Figure 8. Wall-to-wall predictions of species groups overlaid with forest stands labeled by their inventory (photo-interpreted) species groups. Predictions were only made in pixels where zp90 > 5 m, pzabove2 > 50%, and in productive forest stands. JP: jack pine dominated; BS: black spruce dominated; MC: mixed conifers; MW: mixedwood; HW: hardwood.

Discussion

In this study, we combined area-based and individual tree metrics obtained from ALS data with Sentinel-2 data to predict dominant species groups across a large (>580,000 ha) boreal forest located in Ontario, Canada. Below, we discuss the accuracy of the species group predictions, as well as the main variables that helped discriminate between species groups in the context of findings from previous research. We also discuss the implications of using the such approach in an operational context, where information needs vary depending on the forest management approaches.

Accuracy of species group predictions

Our results indicated that the prediction of species groups was the most accurate (OA = 83.20% for 3-group classification; OA = 67.79% for 5-group classification) when models combined ALS and Sentinel-2 predictors relative to models that included only subsets of these predictors. Kukkonen et al. (2018) also observed that the accuracy of Scots pine (Pinus sylvestris), Norway spruce (Picea abies), and broadleaved species volume improved when using both Sentinel-2 and ALS data (RMSE = 49.8%–69.0%) compared to using ALS data only (RMSE = 53.9%–77.1%). In a 14,000 ha temperate forest, Plakman et al. (2020) also found that the combination of Sentinel-2 data and crown structural metrics derived from ALS data yielded the most accurate classification of eight species (OA = 78.5%) at the individual tree level. However, their results show that models based on Sentinel-2 alone performed better (OA = 73.7%) than models only based on crown structural metrics (OA = 58.4%), which was not the case in our study.

Accuracies comparable to the one obtained in this study with the best performing model combining ALS and Sentinel-2 have been observed in studies using Sentinel-2 data alone. For example, Persson et al. (2018) mapped five species in a 1500 ha boreal forest with an OA of 88.2% and Breidenbach et al. (2021) mapped Scots pine (Pinus sylvestris), Norway spruce (Picea abies) and deciduous forests across Norway with an OA of 90.5%. In these two studies, models using multitemporal Sentinel-2 data, that can capture different phenological dynamics between species, outperformed models based on single-date Sentinel-2 acquisition. This could explain the fact that we observed markedly lower accuracy (OA = 65.98% for 3-group classification; OA = 46.23% for 5-group classification) in our study when using a single composite of Sentinel-2 data acquired between June and September. Moreover, our study considered mixed conifers and mixedwood species groups, which comprise a diverse range of species and are not as clearly distinguishable compared to other forest types dominated by a single species.

The accuracy of the models using ALS data only (OA = 76.2%–82.6% for 3-group classification; OA = 61.0%–63.8% for 5-group classification) was within the range of accuracy values that Fassnacht et al. (2016) compiled from 17 studies classifying species with ALS data (OA = ∼55%–90%). It is important to note however that these predictions ranged from the individual tree species to species groups, in various forest types and areas as well as with different ALS sensors and acquisition parameters. To our knowledge, only one study by Prieur et al. (2021) has examined tree species classification with the SPL100. Their study area was also located in Ontario, Canada but in a more complex species-environment than the Romeo Malette Forest studied herein. Prieur et al. (2021) obtained accuracies comparable to the ones observed herein, with an OA of 82.9% if grouping predictions into softwood and hardwood species groups and 68.3% if grouping predictions into pine, spruce, other softwood hardwood, and intolerant hardwood species groups.

Variables driving species groups classification

Our study highlighted that the three sets of predictors provided valuable information to discriminate the modeled species groups. As discussed by Fassnacht et al. (2016), species influence the distribution of ALS returns within the canopy as a function of differences in crown shape or foliage density. Two metrics representing the cumulative distribution of returns in the upper layers of the canopy (zpcum8, zpcum9) had a markedly larger conditional variable importance than the other plot-level summary ALS metrics. The signal attenuation depth metrics, which are influenced by foliage density and canopy closure, were also important variables in the random forest classification model, especially to discriminate between coniferous species groups and mixedwood or hardwood species groups. This also confirms results from Irwin et al. (2021) who observed that the signal attenuation depth of the SPL100 was generally lower in closed-canopy deciduous forests, compared to coniferous forest types. It is also important to note here that that point clouds of the SPL100 tend to be dominated by first returns in the upper canopy because of the low energy of the laser beams and the filtering of noise returns that can decrease the amount of returns in the intermediate and lower layers of the canopy (Yu et al. 2020; Brown et al. 2020).

The red edge region of the electromagnetic spectrum is generally found to be of high importance in vegetation mapping because it results from a low reflectance in the red, due to top chlorophyll absorption, and a high reflectance in the near-infrared (NIR), due to the internal structure of leaves (Horler et al. 1983). This high relevance is also reflected in our findings, with the red edge 2 band of Sentinel-2 (λ = 740 nm) being the important spectral band to discriminate species groups, a finding that is supported by other studies carried out in forests with both coniferous and deciduous trees (Astola et al. 2019; Persson et al. 2018). The red edge 2 band was also the most significantly different (p-value < 0.001) when comparing the jack pine and black spruce dominated species groups with the mixedwood or hardwood species groups. Among the blue (λ = 490 nm), green (λ = 560 nm), red (λ = 665 nm), red edge 2 (λ = 740 nm), and SWIR 2 (λ = 2190 nm) Sentinel-2 bands included in the models, the SWIR 2 band was the only one where significantly different reflectance values were observed between the jack pine and black spruce dominated groups in our study. Sentinel-2 red edge 2 band (λ = 740 nm) and narrow-NIR band (λ = 860 nm), which were not included in our models due to its high correlation (r > 0.99) with the red edge 2 band, have also been shown to help discriminate Scots pine from Norway spruce (Mikael Kukkonen et al. 2018). These results are in line with previous research that has observed that the spectral reflectance of black spruce needles was lower than jack pine needles in the NIR (λ = 700 nm–1000 nm) region of the electromagnetic spectrum (Middleton et al. 1997). In a more recent study also examined the shortwave infrared (SWIR) region of the electromagnetic spectrum, Hovi et al. (2017) observed that spectral differences between coniferous boreal tree species needles were larger in the SWIR (λ = 1300 nm–2500 nm) region of the electromagnetic spectrum compared to the visible (λ = 400 nm–700 nm) or NIR (λ = 700 nm–1300 nm) regions. In addition to the reflectance of needles, the site on which tree species occur can also influence the of each species group. For example, black spruce often occurs on wetter sites than jack pine, and the SWIR region of the electromagnetic spectrum is known to be sensitive to moisture (Zhu et al. 2014).

Individual tree crown features derived from ALS can be effective to separate species groups, in particular coniferous from broadleaf species (Fassnacht et al. 2016). On the other hand, Duncanson et al. (2015) found that the accuracy of area-based estimates of above-ground biomass modeled from ITC metrics aggregated at the plot-level, as done in this research, did not outperform the accuracy of models based on plot-level height summary metrics. In our study, including ITC metrics marginally improved the accuracy of our random forest classifications compared to models with area-based metrics alone. ITC metrics were significantly different (p-value < 0.05) when comparing the hardwood species group to the conifer groups, but not among conifer groups. The metrics related to the crown area or proxy of crown volume were more important if they were summed rather than averaged within a plot or a pixel. Summed metrics could possibly capture more information, e.g., in relation to canopy cover, than averaged metrics and explain why they were more important predictors. The fact that we limited the marker-controlled watershed segmentation to the top of each tree was important. Without this limit, delineated crowns tended to occupy the entire plot or pixel area, and summed crown area metrics had limited variability across the species groups, especially when segmenting closed canopies. This approach was suitable to improve the explanatory power of the crown metrics but would not have been as effective if the objective was to isolate individual tree and their entire crowns as best as possible. More complex individual tree crown metrics, such as 3-dimensional alpha shapes that represent the envelope of a tree crown (Korhonen et al. 2013), could have more explanatory power than the ones that were used herein to classify tree species. However, the generation of detailed individual tree metrics across large areas, such as the Romeo Malette Forest is challenging due to the computational requirements and the necessity of tuning of individual tree segmentation algorithm parameters to different forest types.

The measure of the intensity of linear-mode LiDAR instruments, which is related to the amount of energy in the reflected pulse, can be more useful than the structural information of the point cloud to discriminate tree species (Fassnacht et al. 2016; Shi et al. 2018). Equivalent intensity measures cannot be derived from LiDAR instruments based on single photon LiDAR detection since the returned signal is binary instead of being a pulse of energy. A pseudo-intensity measure can be generated by aggregating the number of binary returns detected for each transmitted pulse but its properties and characteristics compared to the intensity of other ALS sensors are not well known (Hartzell et al. 2018). We did not include the SPL100 pseudo-intensity features as predictor variables in this study, and other studies have shown that it provides a limited amount of information for species-specific volume predictions or species classification. However, further research characterizing their properties across different surfaces and suitable standardization methods would be valuable.

Implication for operational forest management

Whether to use tree-level, area-based or stand-level forest inventory ultimately depends on the information needs End users have come to expect detailed species composition estimates, such as the ones provided by manually photo-interpreted inventories. However, in some contexts, knowledge of dominant species may be sufficient to support information needs, such as estimating timber volume (Tompalski et al. 2014). Area-based species group predictions, such as the ones derived herein can inform extensive forest management practices, which prevail in Canada where forests cover large extents, are mainly publicly owned and managed for multiple objectives (Wulder et al. 2007). In nations with smaller forest landbase and larger share of private ownership, intensive forest management practices are more common and require more detailed and precise species composition information than the ones derived herein. Changes in silvicultural approaches and practices, associated with changing climate broader sustainable use of forest resources (Achim et al. 2022), could also require a more precise inventory of forest resources across large scales.

The operational uptake of forest inventories based on individual tree attributes has partly been limited by the bias induced by over-segmentation or under-segmentation errors associated to the ITC delineation process (Vastaranta et al. 2012; Breidenbach et al. 2010; White et al. 2016). Semi-ITC approaches, in which detected segments can be matched to one or multiple trees have been used to reduce this bias (Breidenbach et al. 2010; Breidenbach and Astrup 2014). Although the range of forest types on which the development and accuracy of ITC algorithms are being assessed is increasing (Kaartinen et al. 2012), complex forest types in which a range of species and conditions coexist are limiting the application of ITC methods across large areas (Zhen et al. 2016; Fassnacht et al. 2016). Regardless, in the context of species mapping, there is no demonstrated gain in the accuracy of forest attribute prediction when using ITC or semi-ITC approaches compared to an ABA. For example, Ørka et al. (2013) used ALS and high-resolution multispectral and hyperspectral imagery to estimate the proportion of spruce, pine and deciduous tree species in a boreal forest of Norway with an R² = 0.66–0.96 when using an ITC or semi-ITC approach and an R² = 0.83–0.89 when using an ABA. The three approaches provided estimates with relatively similar accuracy, but the ITC particularly outperformed the ABA to estimate the proportion of hardwood species that were not as frequent as pine and spruce in their study area. There are few studies however that have compared the relative performance of species predictions generated with area-based and ITC approaches and additional research quantifying differences in a range of forest environments would be informative.

Our area-based predictions of species groups had more spatial detail than the stand-level estimates from the photo-interpreted inventory. The overall stand-level patterns of species distributions were also relatively well captured, except in stands of the mixed conifers species groups where misclassifications of black spruce or mixedwood dominated. It is important to consider that the photo-interpreted data should not be necessarily considered as reference data. Differences between the predicted species groups and the ones available in the forest inventory polygons could be explained by errors in the area-based predictions, photo-interpretation errors, or time span between the photo-interpretation and the SPL100 acquisition. By using consistent data sources and methods, species predictions with remotely sensed data also avoid issues associated with subjective interpretation, particularly when many photo-interpreters are involved as is common when inventories span large forest management units.

The derivation of the aggregated crown metrics developed herein is more computationally intensive than the standard area-based summary metrics, which is an important consideration for operational uptake. However, the ITD and ITC algorithms that we used were applied on the raster-based CHM, which is faster than applying these algorithms directly on the point cloud. Moreover, no predictions were made at the individual tree level since ITC metrics were aggregated, which is an important implication for processing time in large areas, such as the Romeo Malette Forest. Area-based species information, together with area-based forest inventory attributes could be used to automatically delineate homogeneous forest stands similar to photo-interpreted stands as part of the forest inventory. In addition, the wall-to-wall species group predictions could be used to select which optimal set of ITC methods and parameters should be used for a stand of interest.

In this study, we observed a relatively marginal gain in species group classification OA when combining ALS and Sentinel-2 predictor variables compared to using ALS metrics only, particularly for the 3-groups species classification (Figure 3). The Romeo Malette forest is actively managed and some of the structural differences captured by the ALS summary metrics could be related to past silvicultural treatments targeting specific stands or species groups. In unmanaged forests, ALS metrics, in particular those related to height, could be less important and the gain in prediction accuracy obtained by combining ALS and Sentinel-2 data could be larger than the one we observed herein. Sentinel-2 data are global, free, and publicly available and their integration into the generation of enhanced forest inventory products is relatively straightforward. As discussed above, few studies have examined the integration of ALS and Sentinel-2 imagery for forest inventory products, and further research, particularly in areas with more diverse species composition and different management regimes, would be interesting to better understand the benefits of integrating such data into forest inventory frameworks.

Conclusion

ALS data are being increasingly integrated into forest inventory workflows to provide detailed and spatially extensive area-based estimates of key forest attributes for forest managers. However, information related to species composition is often obtained from other data sources, such as photo-interpreted imagery, which does not match the spatial resolution of ALS-based inventories and can have relatively poor accuracy. Previous research has examined individual tree species classification with ALS data, often in combination with spectral data from high-resolution imagery, but scaling this information up to large extents remains a challenge. In this study, we demonstrated across a large forest area (>580,000 ha) that some of the species composition information provided by photo-interpreted inventories can be derived with relatively high accuracy using a combination of area-based and individual tree crown summary metrics of forest structure from ALS data as well as moderate resolution multispectral data from Sentinel-2. We found that the spectral information of Sentinel-2, particularly in the near-infrared region of the electromagnetic spectrum (λ = 740 nm), the area and volume of upper tree crowns, and the cumulative distribution of returns in the upper sections of the canopy were the most important metrics in the classification models. Area-based species grouping products, such as the ones derived herein could be integrated operationally to inform regional forest management and monitoring as well as used to guide the derivation of more precise information (i.e., tree-level) where needed. Future research improving the amount of detail in species composition from remotely sensed products in an operationally viable manner would be particularly useful for forest management.

Acknowledgments

We thank the Ontario Ministry of Natural Resources and Forestry for providing us with the LiDAR and field data used in this research.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Forestry Futures Trust Ontario eFRI Accuracy Assessment and Change Update Approaches project under grant 18-2018; and the Natural Sciences and Engineering Research Council of Canada under grant CRDPJ 462973-14.