Contribution of topographic features and categorization uncertainty for a tree species classification in the boreal biome of Northern Ontario

ABSTRACT

Variations within local topography can effectively impact the location of tree species within naturally forested areas. Furthermore, the uncertainty of prediction for classification can vastly differ amongst topography and the overlying tree species groupings. This study investigated the supplementation of a suite of topographic features corresponding to morphometry and hydrological considerations, in addition to multispectral imagery and other LiDAR-derived features, at fine (2 m) spatial resolution for a pixel-based tree species classification of a forested region of the boreal biome in northern Ontario, Canada. The study area conforms to the Abitibi River Forest (ARF) and consists of the tree species of black spruce (Picea mariana), balsam fir (Abies balsamea), trembling aspen (Populus tremuloides), balsam poplar (Populus balsamifera), tamarack (Larix laricina), white spruce (Picea glauca), and eastern white cedar (Thuja occidentalis). Random forest (RF) and support vector machines (SVMs) were implemented for the classification. Topographic features, specifically those conforming to channel base level, valley depth, and multi-resolution valley bottom flatness (MRVBF), were among the most important features for species predictors. The RF and SVM methods were trained on pixels of pure stands (composed of 70%+ of same tree species) for the tree species groupings, which were split by site level. Modelling accuracies for both the pixel and site level were reported, with the best model attaining an overall site level accuracy and corresponding Cohen’s kappa score of 0.79 and 0.69 for classification, respectively. Entropy maps were generated to characterize the uncertainty of prediction, and substantiate that the regions of lowest uncertainty correspond to wetlands, which are dominated by black spruce (Picea mariana). A modified entropy map was calculated from the normalized top two probabilities of tree species groupings predicted per pixel, so as to better highlight regions of prediction uncertainty. A prediction map for the second most-likely tree species groupings was also computed, which supports the presence of balsam fir (Abies balsamea) as a secondary tree species throughout the ARF region.

KEYWORDS:

• Tree species classification
• boreal forest
• topographic features
• entropy maps
• WorldView-2

1. Introduction

The southern transition zone of the boreal biome within Canada is expected to undergo drastic transformations with regard to vegetation, due in large part to climate change resulting from human-induced activities (Boulanger et al. 2017; Damien, Woodcock, and Friedl 2018). There is a vital need to monitor the forest cover of the boreal biome, for purposes of inventory (Beaudoin et al. 2014; Fassnacht et al. 2016) as well as establishing a basis for the change detection of anticipated future land cover conversion (Beaudoin et al. 2018). These activities can involve estimates for biomass (Garcia et al. 2017; Beaudoin et al. 2014), and identifying tree species or ecosystems at risk of deterioration (Boulanger et al. 2017). There has been a growing interest with respect to digital soil mapping for predicting soil properties within the boreal transition zones, in particularly for estimating carbon stocks (Minasny et al. 2019) that could be further released through land cover change. Vegetation has been identified as an important soil formation factor for peatland environments within boreal regions (Minasny et al. 2019; Beguin et al. 2017), with soil carbon concentrations interrelated to tree species (Mueller et al. 2012). Specifically, tree species indicators have been established to be among the most significant covariates for digital soil mapping applications within the boreal biome of northern Ontario (Pittman, Hu, and Webster 2021). To assess the environmental impact of climate change amongst the southern boreal transition zone, there is a pertinent necessity for accurate tree species classification maps.

There exists many challenges in regards to implementing tree species classifications for natural forested areas. It is imperative to take into account study areas large enough for natural forested regions, so as to obtain representation for a more comprehensive set of tree species endemic to a region. This necessitation of larger study areas, and subsequent expense and processing requirements, means that it can be impractical to delineate individual trees. Consequently, the objective is to obtain tree species maps (Modzelewska, Ewald Fassnacht, and Stereńczak 2020; X. Zhu and Liu 2014; Beaudoin et al. 2014, 2018). These studies consider pixel-based approaches, as medium-to-low spatial resolution satellite imagery is ordinarily utilized, as finer resolution renderings can be unavailable. Object-basis approaches have been implemented for studies where individual tree demarcation is of interest (Wang, Wang, and Liu 2018; Fassnacht et al. 2016), and concentrate on modeling characteristics of individual trees. Unmanned aerial vehicles (UAVs), also known as drones, have been employed as a primary means of obtaining very high spatial resolution features necessary for individual tree delineation and classification (Franklin 2018; Natesan, Armenakis, and Vepakomma 2020). However, due in large part to battery capacity constraints and the requirement for UAV pilots to maintain a line of sight, UAVs typically fly over regions spanning less than a few km² in extent. Hence, UAVs are not feasible for the collection of remotely sensed data over larger and remote domains. Correspondingly, pixel-based approaches are befitting for tree species classification over larger areas within natural forested regions.

As an autecological concept, there exists a definitive linkage between tree species and topography within the boreal forest of eastern Canada (Denneler, Bergeron, and Bégin 1999; Kulha et al. 2019). Wetlands are prevalent in this region, with black spruce (Picea mariana) primarily the only tree species to subsist within the wetland environment (Haynes et al. 2021). Other tree species have preference for growth on particular types of soil or exposure considerations, which are strongly influenced by local topography (Denneler, Bergeron, and Bégin 1999). Topographic features can relate to both local and landscape scale morphometry, hydrological characteristics, and landscape exposure considerations (Franklin 2020; Brandon et al. 2016), and are normally derived from a digital terrain model (DTM). For urban forestry, it may not warrant sufficient justification to consider topographic features, as most trees in those environments were planted by humans, and hence topographic considerations render a minimal role in establishing what tree species emerges where. However, for study areas encompassing adequately large spatial scales, such as for natural forested regions, the existent vegetation can be contingent upon the topography, even if there might not exist appreciable variation in terms of elevation within the adjacent topography. This can be the case with wetland environments, low mounds, or shallow banks along rivers or streams. Numerous studies outright ignore topographic features for tree species classification. As described in the next paragraph, although some studies exploited the use of topographic features in tree species classification, or broadly in land cover classification, most of these utilized low-spatial resolution DTM that were not derived from LiDAR. The effects of local topography on the spatial distribution of tree species, and thus for tree species classification, are hardly investigated.

The implementation of topographic features for forest cover or vegetation mapping can be traced back to previous decades (Strahler, Logan, and Bryant 1978; Stage 1976) and have evolved with the advance of remote-sensing technologies, especially with respect to the availability of high-density LiDAR data. Earlier studies have examined the role of topographic features for tree species classifications, such as shaded hillside (Pitkänen 1998), slope steepness, slope positions or aspect (Thomas and Anderson 1993), or inclination of the slope and fetch considerations along with elevation, exposure, and aspect (Denneler, Bergeron, and Bégin 1999). Nevertheless, within existing studies, topographic features are generally only considered for mountainous regions or areas with appreciable differences with elevation, and then only basic topographic features, such as elevation, slope, convexity, and curvature, are employed. For example, Dalponte, Bruzzone, and Gianelle (2012) utilized a DTM derived from LiDAR data with a point density of 8.6 points per m² for a pixel-based classification of a mountainous environment of approximately 10.9 km² in the southern Alps of Italy, where a distinct interconnection between elevation and tree species was established by the forest types (Dalponte, Bruzzone, and Gianelle 2012). Other studies have also included DTM (Chambers et al. 2013; Suratno, Seielstad, and Queen 2009), or basic topographic covariates such as slope (Qin et al. 2015) or convexity and curvature (X. Zhu and Liu 2014). In a recent study for a montane forest in Mongolia (Chiang and Valdez 2019), a gain in tree species classification accuracy was demonstrated when the topographic variables of slope, curvature, aspects, and topographic wetness index were applied together with surface reflectance. The relationship between topographic features and tree species (vegetation cover types) is more evident in the mountainous regions and thus drives the abovementioned studies. However, minimal research has been conducted on the exploitation of topographic features in tree species classification for areas deemed relatively flat with respect to the landscape scale. High-density LiDAR provides for the generation of detailed topographic features, and research can be pursued to exploit the benefit of LiDAR-derived topographic features for classification. Rudimentary topographic features such as DTM or slope might suffice for feature selection in mountainous environments, but drastic variations within elevation may not be present to yield a sound association between tree species and basic topographic features for other biomes. It is warranted to investigate the utilization of topographic covariates corresponding to other attributes and scales of morphometry and hydrologic characteristics, for the prediction of tree species within natural forested region; particularly for domains with relatively flat terrain at the landscape scale, as was the case in this study.

It is also important to address uncertainty estimates with regard to classification, as these assessments are under-investigated. For many classification studies, only accuracies are reported, in the form of confusion matrices and measures derived from them such as overall accuracy and Cohen’s kappa score (Fassnacht et al. 2016; Dalponte, Bruzzone, and Gianelle 2012; Dalponte et al. 2015; Sheeren et al. 2016). Although these metrics provide the general confidences for generated classification maps, these do not convey information for classification uncertainty on individual pixels. Uncertainty characterization is attracting more attention nowadays; as pointed out by Minasny et al. (2019), supplemental uncertainty evaluations should be expressed for all associated mapping procedures. From most classification methods, such as random forest and support vector machine, confidences can be calculated and reported to indicate the posterior probabilities that a given sample belongs to each class. The class with the highest confidence is assigned to this sample as the prediction. One may argue that the corresponding confidence for the prediction may be regarded as a measure of uncertainty. Nonetheless, it can only provide local information on a single event; that is, how likely the given sample belongs to the predicted class. This measure cannot provide any global information on how different the confidence is amongst other categories. In a previous study (Pittman, Hu, and Webster 2021), entropy was calculated based on the posterior probabilities over all categories to measure the classification uncertainty. This measure was further developed by Hu et al (2021), to consider information measures of correlation with the entropy for categories with the highest posterior probabilities, as well as a difference entropy. In this study, an entropy weighted upon just the top two contending probabilities, here denoted as normalized entropy, is proposed, to which is argued to be a better metric for determining uncertainty for classification prediction.

A diversity of tree species can be present within natural forested regions conforming to a biome. As the number of tree species increase, the level of accuracy attained for classification generally decreases, as it can be difficult to obtain sufficient data comprising each subsisting species for training and testing (Chambers et al. 2013; Fassnacht et al. 2016). Confronting this issue can be challenging, particularly if certain tree species are substantially more widespread in a study area than others. The usage of finer spatial resolution data can alleviate part of this setback, as it can facilitate the collection of pixels corresponding to less prevalent tree species for a study area. MODIS imagery of 250-m spatial resolution (Beaudoin et al. 2014, 2018), or even Landsat (Fassnacht et al. 2016) or Sentinel-2 of 30-m and 10-m spatial resolution, respectively, can be too coarse for obtaining pixels corresponding to pure stands of rarer tree species. Furthermore, there is still the issue of resolving mixed tree stands, as forestry plots can be interspersed of separate tree species within a pixel of rasterized imagery acquired from remote-sensing technology. For the classification of forest scenes at a large scale using remotely sensed data, it is common that low spatial resolution imagery is utilized; thus, broad cover types such as coniferous versus deciduous, or softwood as contrasted with hardwood, are often the focus. With higher spatial resolution imagery, species or genera categories are oftentimes of interest, as individual trees can be better discerned. At the local scale, the coexistence of trees of different species poses challenges, especially within a natural forested area. Within the boreal forest of Ontario, Canada, trees of certain species, such as balsam fir (Abies balsamea), frequently grow together with other tree species, and are rarely present in what can be constituted as pure stands. To address the mixture of stands, dominant species are often identified and considered for classification (Altamash, MacLean, and Hennigar 2019; Dalponte et al. 2013; Lindberg, Holmgren, and Olsson 2021). However, the inclusion of mixed tree species groupings can mitigate misclassification for localities where different tree species are sufficiently interspersed. Designating representative categories for the study area is paramount, which can warrant mixed groupings. The impetus should be to adopt as refined of vegetation cover groupings as practical, that are characteristic of the study area. For this study, species-specific categories were considered where feasible with regard to pure stands, and mixed categories for where noted otherwise, to accommodate the coexistence of tree species.

To enhance the robustness of the classification, numerous separate sites corresponding to ascertained tree species groupings, rather than large adjoining areas of pixels, should be selected. The pixels should be allocated for model training and testing by site, and thus classification accuracies should be assessed at the site level. This would prevent the inflation of accuracy due to pixels from a same contiguous selection (for one plot) being used for both model training and verification.

In summary, exploiting an expanded set of topographic features and finer spatial resolution remotely sensed data, a tree species classification was conducted for a boreal study area within the Abitibi River Forest region of the District of Cochrane in northern Ontario, Canada. Accordingly, the availability of finer (2 m) spatial resolution data from satellite imagery, as well as recent LiDAR data, facilitated the generation of more refined spatial features for a natural forested environment. This amounted to the inclusion of a variety of features corresponding to a more comprehensive array of topographical considerations, to improve classification accuracy for a pixel-based approach.

In addition, the application of entropy with respect to the analysis of classification was explored in this study, to systematically gauge the uncertainty. This was expressed in terms of a conventional entropy metric, as well as a weighted normalized entropy of the top two choices of tree species. The resulting prediction and entropy maps have identified localities for specific tree species groupings, as well as locating areas of greater modeling uncertainty.

2. Materials and methods

2.1. Study area

The study area is located within the Abitibi River Forest (ARF) in the District of Cochrane in Ontario, Canada, in the vicinity surrounding the community of Smooth Rock Falls. It is composed of representative expanses of boreal forest within northern Ontario, consisting of also a river valley and wetlands. The study area at present is relatively unaltered, but at risk of future land use conversion. The longitudes for this study area range from 81.55° W to 81.80° W, and the corresponding latitudes range from 49.15° N to 49.35° N. This study area consists of four neighboring areas, encompassing a total of 108 km². Elevations above sea level vary from a minimum of 210 m in the north, to 280 m in the south, as the topography of this region is relatively flat at the landscape scale and comprises of wetland environments, with a higher elevation prominence near the center of the study area. The Mattagami river transverses through this region on a northward trajectory, with steep river banks (slopes exceeding 25°) along sections of the river. Clearings corresponding to a highway, railway, utility lines, and other roadways crisscross parts of the study area domain.

Predominantly, the study area conforms to mature tree stands of primary forest. Some secondary forest is present along the confines of the study area, that arose due to logging activity from decades past. This region corresponds to the boreal forest and consists primarily of the tree species of black spruce (Picea mariana) and balsam fir (Abies balsamea), with other subsisting tree species including trembling aspen (Populus tremuloides), balsam poplar (Populus balsamifera), white spruce (Picea glauca), tamarack (Larix laricina), and eastern white cedar (Thuja occidentalis). Black spruce dominates in the wetland environments, but is also prevalent in majority stands in non-wetland expanses. Balsam fir is frequently interspersed with other tree species. Trembling aspen is commonly encountered with balsam fir, or in unmixed groves near balsam poplar in younger forest stands. White spruce grows principally in upslope localities near streams and rivers, and tamarack occurs in younger forest stands or along the edges of wetlands in close proximity to black spruce. The presence of eastern white cedar is relatively scarce for this region, as this coincides to the northern periphery of its range.

2.2. Field data

Species composition data of observed tree species were obtained from the Ontario Forest Resources Inventory (FRI) (Paloniemi 2018). These data corresponded to 200 forested sites for this study area and noted the counts of tree species by prism sweep for each site. The tree species counts were recorded at a location of 90 m into a 200-m-long straight transect, for an approximately 10 m distance in all directions. From these counts, the concentrations of tree species were compiled by FRI from the tally observations. These data were logged during two separate field campaigns by Ontario Ministry of Natural Resources and Forestry (MNRF) and FRI personnel, that occurred during the autumns of 2016 and 2017. Transects for the site plots were selected so that each transect comprised of tree species of a more consistent composition, which was confirmed by counts to ascertain leading tree species at nearby locations of 50 m and 150 m into the 200-m-long transect. Separately, the transects were delineated so as to obtain representation of different tree species compositions, which were recorded from a diversity of topographic settings within the study area. For each site plot, an MNRF forester and an FRI specialist determined the compositions of overstory tree species present by means of counting individual trees that were over 2 m in height. Dead trees, including trees without live buds or green foliage, were omitted from the tallies. Additional attributes such as estimates for stand age and height for leading species sample trees, as well as soil textural properties from soil samples extracted at the 90 m position into the 200 m transect, were also obtained so as to determine ecosite classifications for the sites. Nonetheless, the observed tree species compositions were the details of interest for the tree species classification research.

The threshold for pure stands was set so that the dominant tree species corresponded to at least 70% of the composition of the tree species cover for a site (Bravo-Oviedo et al. 2014). There was a total of 104 sites that corresponded to pure stands; these conformed to 6, 5, 86, and 7 sites for the classes of eastern white cedar (Cw), tamarack (La), black spruce (Sb), and white spruce (Sw), respectively. A category composed of trembling aspen or balsam poplar was denoted populus (P) and consisted of five sites. There were 67 sites that were not pure stands, but where a dominant tree species represented 50% to less than 70% of composition of the tree species; these consisted of 57 sites for balsam fir (Mix-Bf) and 10 sites for black spruce (Mix-Sb). Out of these 57 sites for balsam fir, there were three sites which technically corresponded to pure stands (i.e. 80% or higher composition) but were included with the mixed balsam fir class for practicality reasons. A separate category for pure stands of black spruce was considered versus the analogous mixed category, as for this study area the mixed black spruce category would likely correspond to upland black spruce, in contrast to the stunted lowland black spruce characteristic of the wetland environments. The remaining 24 sites were sufficiently mixed (Mix) so as there were no tree species dominant (i.e. less than 50% composition) at each of these site. The study area, with locations of the FRI sites allotted for training and validation, and a true color composite corresponding to 10 August 2014, are depicted in Figure 1.

Figure 1. Study area for the tree species classification within the ARF, with FRI sites denoted with respect to model training and validation. Background is a true-color composite of the area from WV2 imagery for 10 August 2014. The grid coordinates correspond to the reference system NAD 1983 (2011) UTM zone 17N.

2.3. Data and features used for classification

LiDAR data were obtained via Land Information Ontario from the MNRF. These data were acquired from airborne LiDAR during October of 2016 and September of 2017 (Airborne Imaging 2018), using a Leica ALS70-HP instrument, collected at an average point density of 8 points per m². Up to four returns per retrieval were attained from the LiDAR. The flight height was 1000 m above the ground level, with a single pass swath width of 690 m with 20% overlap. For the LiDAR instrument, the scan angle field of view was an effective 38°, the scan frequency was 53 Hz, and the scan pulse rate frequency was 500 kHz. LiDAR rasterized imagery in terms of a digital surface model (DSM), DTM, and canopy height model (CHM) were each derived from the LiDAR point cloud data, to each correspond to rasters of 2-m spatial resolution. The DSM was generated from the elevations of the first returns of the retrievals from the LiDAR point cloud data, whereas the DTM was computed from the corresponding elevations of the last returns. When calculating the LiDAR-derived attributes per 2-m cell size, both averaged values of the elevations for the first and last returns, as well as the maximum and minimum values, were considered for generating the DSM and DTM, respectively. Upon further investigation, the maximum first return elevations were incorporated for the DSM, and the minimum last return elevations for the DTM, with CHM calculated as the difference resulting from this DSM minus DTM. For quality assurance, the CHM was inspected to ensure that no values exceeded 60 m for the study area, as CHM would very rarely surpass 30 m (the exception would be for old growth mature white spruce along the slopes of the river valley). As well, the elevations from the computed DTM were compared to those of a coarser resolution 30 m DTM obtained from the MNRF.

WorldView-2 (WV2) imagery of less than 2% cloud cover of 2-m spatial resolution for multispectral and panchromatic bands were obtained for the study area for the date of 10 August 2014, which coincided to a period of peak vegetation. This imagery corresponded to surface reflectance acquired from coastal (400–450 nm), blue (450–510 nm), green (510–580 nm), yellow (585–625 nm), red (630–690 nm), red edge (705–745 nm), near infrared I (770–895 nm), near infrared II (860–1040 nm), and panchromatic (450–800 nm) bands (Digital Globe 2016). The optical imagery was preprocessed, and then reprojected and coregistered with the LiDAR rasterized imagery, followed by orthorectification with a 2-m spatial resolution DSM generated from the LiDAR data. The final projection employed for the analysis was NAD 1983 (2011) UTM zone 17N. Normalized difference vegetation index (NDVI) and a modified normalized difference water index (NDWI) were calculated from the near infrared I (NIR1) and red bands, and near infrared II (NIR2) and NIR1 bands, respectively, of the WV2 imagery. A better calculation for NDWI would have involved a shortwave infrared (SWIR) band (Liuzzo et al. 2020), but due to lack of availability for the finer spatial resolution, the NIR2 band was substituted.

System for Automated Geoscientific Analyses Geographic Information System (SAGA GIS) (SAGA Development Team 2020) was utilized to generate topographic features from the DTM. These included covariates for terrain analysis, as well as covariates corresponding to morphometry and topographic openness. The topographic features generated were negative openness, multi-resolution valley bottom flatness (MRVBF), multi-resolution ridge top flatness (MRRTF), closed depressions, channel network distance, aspect, relative slope position, topographic wetness index, total catchment area, and valley depth.

In total, there were 22 features employed, comprising of eight multispectral, one panchromatic, two normalized difference indices, CHM, and 10 topographic features generated from a DTM. Pixel selections of 20 m by 20 m plots were extracted from each of these feature layers, each corresponding to the location along the transect for where these tree species tallies were recorded for each site.

2.4. Classification

A pixel-based classification evaluated at the site level was implemented. Stratified random sampling was applied to sample the pure stands sites by tree species type. This sampling methodology was adopted, so as to better model and account for the tree species groupings that are less prevalent for the study area. The data were split at the site level using stratified random sampling training (70%) and test (30%) sets, which resulted in 143 sites allocated for training and 57 sites for validation, respectively.

Random forest (RF) and support vector machine (SVM) approaches were employed for the tree species classification. Machine-learning approaches such as SVM and RF have been popular methods for tree species classifications (Fassnacht et al. 2016), as these are non-parametric classifiers that are flexible with regard to the number of predictors used with training sample sizes (Wang, Wang, and Liu 2018). From classifications for digital soil mapping research (Brungard et al. 2015; Brandon et al. 2016), RF and SVM approaches have been demonstrated to attain higher accuracies when compared to other machine-learning techniques. The structure of a RF is comprised of an ensemble of decision trees, with each decision tree trained upon a randomized bootstrap sample; node-splitting decisions are based upon randomized subsets of predictors (Sage, Genschel, and Nettleton 2020). The SVM method employs hyperplanes to optimally separate classes in the feature space (Zhang et al. 2006). The default case of a SVM is the linear SVM, which corresponds to the class boundaries of the training set of the feature space being resolved by a linear model; a radial-basis SVM transforms the feature space into a higher dimensional space by a radial basis expansion, by which a linear model is applied onto this transformed feature space for determining the class boundaries (Zhang et al. 2006). The caret package in R (Kuhn 2008) was applied for fitting the RF and SVM methods, and the models were trained by applying 10-fold cross-validation (CV) with three replicates. The RFs were trained with the number of trees (n_tree) set to 1000, and the number of randomly selected predictor variables at each node (m_try) set as 5, which corresponded to the integer-rounded square root of the number of features. Radial-basis SVMs were trained with the tuning cost parameter varying among 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, and 128, with the kernel width parameter (sigma) ranging from 0.005 to 0.05. Parameter values were ranged and tested within the CV process to obtain optimal models with respect to accuracy. Boosting-based methods for classification with the caret package in R were also investigated, but higher accuracies were obtained from RF and SVM approaches.

Variable importance was determined from the RF method, which was evaluated based upon the mean decrease in Gini, which measures the total decrease in node impurity for the nodes of trees for the RF (Sage, Genschel, and Nettleton 2020). The decrease in node impurity was achieved by executing splits on the nodes so as to increase homogeneity within the leaves of the trees of the RF. A higher mean decrease in Gini value signifies a higher variable importance for a predictor variable. It was conjectured that predictors coinciding to topographic features would have high variable importance when compared to predictors relating to vegetation such as surface reflectance, as different species of trees tend to thrive in different localities within a forest that are correlated with topography (Denneler, Bergeron, and Bégin 1999).

As an alternative to variable importance from the RF models, a model-independent approach for ascertaining the most crucial features was resolved by receiver operating characteristic (ROC) curve analysis (Kuhn 2007). For each predictor, ROC analysis was implemented by reducing the classification problem with multiple categories to pairwise problems, where a series of cutoffs were imposed to determine class from the predictor data. Sensitivity and specificity were calculated for each cutoff, and the corresponding ROC curve was generated. The area under the ROC curve was computed as a quantification of variable importance. For each classification category, the maximum area under the ROC curve from the corresponding pairwise problems was regarded as the variable importance measurement. This resultant listing of features, ordered by maximum area under the ROC curve, was compared to the ranked features computed from the variable importance of the RF model.

A feature reduction was implemented so as to reduce effects of multicollinearity, by employing just the top 10 predictors (i.e. features) with the highest mean decrease in Gini for the variable importance from the RF method. This cutoff for the amount of features was implemented so as to keep the amount of features used less than 10% of the number of sites, and to also allow for a sufficient amount of features to capture the bulk of the mean decrease in Gini with the variable importance. The resulting reduced set of features was also utilized for the SVMs.

2.5. Accuracy assessment and quantification of classification uncertainty

Both pixel and site level accuracy were assessed for model evaluation. Accuracy by means of the percent correct classification (PCC), as well as Cohen’s kappa statistic, were the metrics for evaluating the classifications. PCC, here referred to as overall accuracy, reports the correctness portion with prediction. The kappa score considers agreement arising by chance (Kuhn 2008) and takes into account the occurrence of misclassification. Overall accuracy varies between 0 and 1. Kappa score ranges between −1 and 1 (Cohen 1960; McHugh 2012), with 0 denoting agreement equivalent to coincidence for the kappa score, and 1 indicating complete agreement; negative values signify agreement worse than that attributed to chance (Cohen 1960) and point to exceedingly poor agreement (Monserud and Leemans 1992). Kappa scores less than 0.4 suggest poor to at most a fair agreement, and between 0.4 and 0.8 indicate moderate to substantial agreement, whereas kappa scores greater than 0.8 imply strong agreement (Brungard et al. 2015; McHugh 2012). From the confusion matrix of classification results, producer’s accuracy (PA) as well as user’s accuracy (UA) were also computed to evaluate accuracy for the individual classes of the tree species groupings.

The uncertainty in prediction with the models was resolved by means of entropy maps. To determine the extent that the models assigned predictions to a specific category (i.e. in this case tree species), the probabilities for the categorial assignment (i.e. voting) from the models were utilized to calculate entropy scores (Roulston and Smith 2002; Zhu 1997). The entropy was calculated for each pixel as

(1) $H=\frac{-1}{ln\left(n\right)}{\sum }_{k=1}^n{p}_k\cdot ln\left(p_k\right),$(1)

with $p_k$ the probability of a pixel being assigned to category $k$, for a total of $n$ different categories. These probabilities for categorical assignment were each confined as $0\le p_k\le 1$, with the sum of these probabilities for each pixel adding up to one. For the RF, these probabilities were computed from the fraction of the total number of trees which voted for each respective category (Sage, Genschel, and Nettleton 2020). With the SVM methods, a secondary Platt model resolved these probabilities (Kuhn 2013) by means of a sigmoid function to approximate posterior probabilities (Lin, Jen Lin, and Weng 2007). This entropy metric, as the sum of probabilities multiplied by the natural logarithm of each respective probability, is a form of Shannon’s entropy (Kumar, Kumar, and Kapur 1986). The entropy as defined in this equation is bounded between 0 and 1, with lower values indicating less uncertainty in regards to classification.

For the entropy, similar probabilities amongst categorical assignment for a pixel would signify greater uncertainty and consequently a higher entropy score, whereas a considerably higher probability assigned to one category would result in a low entropy score as there would be minimal uncertainty with prediction. However, it might be problematic to ascertain regions of greater certainty in modeling, if the category with maximal probability was allocated only a slightly higher probability than other respective categories for a pixel. For that reason, the entropies with prediction for each pixel were also calculated from only the top two highest probabilities for prediction; i.e. the two highest probabilities out of the eight tree species groupings. These top two probabilities were normalized so that they summed to one; that is, per pixel,

(2) ${\hat{p}}_1=\frac{p_{max}}{p_{max}+p_{2ndmax}},$(2)

(3) ${\hat{p}}_2=\frac{p_{2ndmax}}{p_{max}+p_{2ndmax}}.$(3)

The normalized entropy for the top two rankings takes the term form of Shannon’s entropy, and was computed for each pixel by

(4) $H_{Top2}=\frac{-1}{ln\left(2\right)}\left[{\hat{p}}_1\cdot ln\left({\hat{p}}_1\right)+{\hat{p}}_2\cdot ln\left({\hat{p}}_2\right)\right].$(4)

From this metric, it is asserted that regions with greater prediction certainty can be identified, when compared to the entropy calculated from the consideration of all probabilities for each respective category. This can be a useful evaluation when only the top rankings for assignment from the corresponding categories are of interest, which would better weight the top rankings, so as to get an idea of the difference between first and second (i.e. contending) rankings. As with the other entropy stated before, this normalized entropy is bounded between 0 and 1, where lower values signify lesser uncertainty. The entropy (1) and normalized entropy (4) were each calculated from the best models for prediction.

3. Results

The results from the variable importance are depicted in Figure 2. It is discerned that the topographic feature of channel network base level had the highest variable importance. This feature was followed in variable importance by the topographic features of valley depth and MRVBF, and then NDVI and relative slope position. Afterwards, those features were followed in contribution by CHM, then by the topographic covariates of MRRTF and negative openness, and then various bands of surface reflectance. The features of channel network correspond to the elevation of the channel network base in terms of interpolated channel base level elevations, or else in the same units as the DTM, respectively, with valley depth denoting the vertical distance above the channel network base level (Etzrodt, Zimmermann, and Conrad 2002; Conrad et al. 2015). MRVBF refines valley bottoms as flat low-lying areas (Gallant and Dowling 2003) and MRRTF analogously resolves ridge tops, respectively; negative openness values are high for below-ground pits or depressions (Ryuzo, Shirasawa, and Pike 2002). The features relating to the topography of the channel or valley (specifically, channel network base level, valley depth and MRVBF), in addition to NDVI and CHM for the surface reflectance and height structure of the corresponding trees, respectively, had the highest variable importance for the prediction of tree species. When ROC curve analysis was performed on the predictors for this tree species classification, the same topographic features as identified by the variable importance from the RF models still attained the highest variable importance.

Figure 2. Variable importance for the tree species classification, ranked by mean decrease in Gini, as ascertained from the RF model.

Accuracies evaluated by percent correct classification and kappa statistic for the models are shown in Table 1. These accuracies are reported for the validation set and are summarized in terms of pixel level and site level, respectively. The highest accuracy SVM radial-basis model was fitted when sigma equaled 0.036, and the corresponding tuning cost parameter was set to 64. Accuracies from the RF and SVM approaches were comparable, with the RF overall attaining the higher accuracies. For the validation results, the accuracies for the RF model exceeded 0.7 and the kappa scores surpassed 0.5, indicating at least moderate agreement for the classification.

Table 1. Accuracies of the tree species classification methods, in terms of percent correct classification (accuracy) and Cohen’s kappa, when models were tested on the validation set. Accuracies are reported at both the pixel level, and the aggregated site level.

Model Accuracies	Tree Species Classification
Validation	Pixels		Sites
	Accuracy	Kappa	Accuracy	Kappa
Random Forest (RF)	0.72	0.58	0.79	0.69
Support Vector Machine Radial-basis (SVM Radial)	0.62	0.46	0.74	0.61

A prediction map for the most probable tree species groupings is presented in Figure 3. This map was generated per pixel for the tree species category which attained the highest probability of prediction from the RF. The most predominant tree species for this region is black spruce, which prevails throughout the wetland regions of the ARF. The second most common tree species is balsam fir, which is typically interspersed with other tree species, and widespread throughout the non-wetland areas. Other tree species each tend to subsist in assorted environments within this region, such as with white spruce along stream or river valley slopes. The populus tree species of trembling aspen and balsam poplar are generally dominant along clearings or in regions with deeper surface soil horizons, and tamarack exists either along the edges of wetland environments or within clearings with compacted soil.

Figure 3. Prediction map for most probable tree species groupings for the ARF study region, as determined from the RF model.

Confusion matrices for both the RF training and validation are shown in Table 2. PA as well as UA are also reported for each tree species. Accuracies and confusion matrices were also reported for RF training, so as to compare these metrics to those for the validation data. Note that these confusion matrices also convey the counts of how many sites correspond to each tree species grouping. For the method training, high accuracies were attained for all the tree species, with the exception of eastern white cedar. With the validation data, the RF predicted well for tamarack, white spruce, pure stand black spruce and mixed balsam fir tree species groupings. However, the prediction for eastern white cedar (Cw) and populus species (P) were very poor, with a PA of zero for each. Nevertheless, there were only one validation site each for these tree species groupings. Mixed black spruce also had a PA of zero, but it is worth noting that all three validation sites for mixed black spruce (meaning comprised of 50–70% black spruce) were allocated to pure stand black spruce (i.e. 70%+ black spruce). A majority of sites for the mixed tree species grouping were properly assigned, and 45 out 57 validation sites corresponded to groupings with both PA and UA exceeding 0.80. For this study region in the ARF, balsam fir rarely exists in pure stands and is ordinarily intermixed with other tree species, particularly the deciduous populus species, so poorer prediction for mixed groupings can partially be attributed due to that account. Eastern white cedar is not prevalent in this study region, as for the ARF this tree species subsists at the outer reaches of its natural range, limited in part due to winter hardiness, which could explain the poorer modeling for this tree species. Part of the reasoning for the poor accuracy of the prediction for mixed black spruce could be asserted by that these regions likely correspond to non-wetland regions (where black spruce is not the only tree species present) characterized by the taller upland manifestation of black spruce; so models would have difficulty ascertaining this tree species from a more concentrated grouping (commonly coinciding to lowland black spruce).

Table 2. Confusion matrices for the tree species classification from the RF model, for counts aggregated at the site level. Counts are reported for the training (top) set corresponding to 143 sites, and validation (bottom) set corresponding to 57 sites, respectively. Producer’s accuracy (PA) and user’s accuracy (UA) are also reported for each tree species grouping class.

Confusion Matrices		Tree Species Classification
Training & Validation												Producer’s
		Prediction										Accuracy
Reference		Cw	La	Mix	Mix-Bf		Mix-Sb		P	Sb	Sw	(PA)
Eastern White Cedar	Cw	3			1					1		0.60
Tamarack	La		4									1.00
Mixed	Mix			15	2							0.88
Mixed - Balsam Fir Dominant	Mix-Bf				40							1.00
Mixed - Black Spruce Dominant	Mix-Sb				1		6					0.86
Populus	P								4			1.00
Black Spruce	Sb									61		1.00
White Spruce	Sw										5	1.00
User’s Accuracy (UA)		1.00	1.00	1.00	0.91		1.00		1.00	0.98	1.00	Sites: 143
Eastern White Cedar	Cw									1		0.00
Tamarack	La		1									1.00
Mixed	Mix			4	2					1		0.57
Mixed - Balsam Fir Dominant	Mix-Bf			3	14							0.82
Mixed - Black Spruce Dominant	Mix-Sb									3		0.00
Populus	P				1							0.00
Black Spruce	Sb						1			24		0.96
White Spruce	Sw										2	1.00
User’s Accuracy (UA)		.	1.00	0.57	0.82		0.00		.	0.83	1.00	Sites: 57

Entropies were calculated for the validation sites, to obtain a measure of uncertainty with prediction, as reported in Table 3. The normalized entropy based upon the top two contenders as computed by (4), resulted in a higher value, but likely more accurate assessment for uncertainty, than the entropy as computed by (1). Amongst the different tree species groupings, the categories for balsam fir and black spruce tended to have the lowest predication uncertainties, with these categories more represented due to the higher reported site counts. Note that with the conventional entropy metric, that validation sites comprised of white spruce had a lower entropy score than those for black spruce, but that this situation reversed with the normalized entropy metric. Also observe that the normalized entropy metric resulted in a larger range in values among the various categories, which would help to better distinguish different tree species based upon this entropy consideration. Even for sites that were improperly classified, such as the 1 site each for eastern white cedar (Cw) and populus (P), respectively, the normalized entropy (4) conveyed more uncertainty with prediction than that from the conventional entropy (1) metric.

Table 3. Uncertainty calculations with respect to entropy (1) and normalized entropy from the top 2 choices (4), summarized by tree species grouping, for the validation set sites.

Uncertainty Calculations			Validation Sites
			Entropy	Normalized Entropy
Tree Species Grouping		Site Count		Top 2 Choices
Eastern White Cedar	Cw	1	0.64	0.83
Tamarack	La	1	0.73	0.96
Mixed	Mixed	7	0.64	0.86
Mixed - Balsam Fir Dominant	Mixed - Bf	17	0.55	0.72
Mixed - Black Spruce Dominant	Mixed - Sb	3	0.54	0.66
Populus	P	1	0.38	0.46
Black Spruce	Sb	25	0.52	0.70
White Spruce	Sw	2	0.45	0.86

The entropy map for the study area as generated by (1) is shown in Figure 4. The areas with the lowest prediction uncertainties conform to the wetland environments where black spruce is typically the only tree species which can subsist in those areas. Other areas have higher prediction uncertainties. An entropy map corresponding to the normalized top two probabilities of prediction for the tree species groupings via (4) is depicted in Figure 5. Note that the coloration scale used for both entropy figures are the same, with minimum and maximum intensities matching to 0 and 1, respectively. From the normalized entropy figure, more variation with the entropy is apparent. Areas surrounding the wetlands have higher uncertainties with prediction, than what were presented in Figure 4. From the most probable tree species grouping map in Figure 3, it can be determined that the nearest areas to the wetlands correspond to black spruce. It is probable that many of these localities consist of upland black spruce, which are adjacent to wetland regions that are still dominated by black spruce. Within the upland black spruce stands, other tree species can also subsist, which would account for the higher prediction uncertainties when compared to that for the lowland black spruce in the wetlands. Also note that higher entropies are encountered in areas with mixed tree species groupings (from Figure 3), which is intuitive as mixed categories by definition setup encompass more uncertainty with regard to tree species assignment.

Figure 4. Entropy map for uncertainty of prediction, as calculated by (1), for the ARF study region. Darker color intensity indicates lower prediction uncertainty, whereas lighter color intensity corresponds to higher uncertainty with prediction.

Figure 5. Normalized entropy from the top two choices of tree species groupings, as calculated by (4), for the ARF study region. Darker color intensity indicates lower prediction uncertainty, whereas lighter color intensity corresponds to higher uncertainty with prediction.

Predicting for the second most dominant tree species might convey practical information for a study region. A prediction map for the second most probable tree species for the study area is depicted in Figure 6. Note for this map, that no second most probable tree species was identified for pixels where the most probable tree species was allocated a probability greater than 70%, as these areas would most likely correspond to pure stands (such as the wetlands for black spruce). For this map, the mixed balsam fir grouping is widespread, as balsam fir is a tree species that grows as a secondary species throughout the ARF. The mixed grouping is also more prevalent across the study area, which would be expected for localities that do not conform to pure stands.

Figure 6. Prediction map for second most probable tree species groupings for the ARF study region, as established by the RF model.

4. Discussion

4.1. Features

The inclusion of topographic features improved the classification of tree species for a boreal study area in northern Ontario, as substantiated in Figure 2. For this study, geomorphometry information relating to the depth of valleys, ridge flatness, valley bottom flatness and slope position had high importance. Valleys and channels can offer protection against the effects of wind and excess sunlight that could be detrimental to certain tree species, whereas water can accumulate in low spots and either lead to an increase or decrease in vegetation diversity depending on the context of the topography. Within the boreal biome of northern Ontario, white spruce tends to flourish along river banks that are elevated above the flood plain next to the corresponding waterbodies, whereas black spruce can grow in wetland areas where no other tree species can subsist. Topological features should be selected with regard to the underlying local attributes at a point, textural characteristics of the topographic surface, and contextual characteristics for the larger landscape scale (Franklin 2020). Summarizing (Franklin 2020), local variables include slope, aspect and curvatures, whereas textural features include roughness or homogeneity indices. Contextual covariates would correspond to slope position and hydrological considerations such as catchment areas and topographic wetness indices. Franklin (2020) also writes that employing too few geomorphometric variables leads to inadequate topography characterization, but applying multiple highly intercorrelated covariates can cause inaccurate modeling results. More relevant geomorphometric information can be conveyed from a combination of other topographic covariates, rather than from just employing a DTM or features relating to local characteristics.

It is possible that the WV2 multispectral bands did not contain as much diagnostic information when compared to topographic features, as there are both conifer and deciduous tree species present in the region with similar colorations of foliage, that are related within the same genus, respectively. Specifically, there exists black spruce and white spruce, and then the two populus species of trembling aspen and balsam poplar. There are also categories with mixed species, which would not perform well for differentiating species based on multispectral imagery as pixel intensities would be composited. Nonetheless, NDVI was an exception and thus ranked with higher variable importance, as the normalized difference in surface reflectance between the red and infrared bands was significant. The CHM, corresponding to the structure of the vegetation canopy, had higher variable importance than the individual multispectral bands pertaining to surface reflectance of the foliage.

4.2. Entropy

Entropy metrics can be applied to quantify the measure of uncertainty in prediction with regards to a study area. This process has been common in many other disciplines, more recently in digital soil mapping (Brandon et al. 2016), and can be applied to pixel-based classifications. Entropy metrics are better for assessing prediction uncertainty than the aggregation of output from a variety of modeling approaches, as entropy considers voting probabilities rather than the outright categorical assignments. For this study area, pure stand black spruce in wetlands had higher probabilities of assignment than localities with mixed balsam fir (or other tree species groupings), as seen in Figures 4 and 5. An ensemble of model predictions from equally accurate or valid models would likely predict the same most probable tree species for the overwhelming majority of localities, but this would not explicate the fact that the individual models are allotting more certainty with predicting for pure stand black spruce in wetlands than for areas populated with balsam fir. A substantial discrepancy in prediction maps between various approaches could be attributed to inaccurate models, which in actuality should not be utilized for uncertainty analysis. However, if there are more than a few distinct categories, the entropy metric can become diluted in value if there are multiple categories with negligible probabilities, and correspondingly only a couple of categories with contending probabilities. A diluted entropy score would underestimate uncertainty with prediction. For that reason, a normalized entropy metric should be evaluated, based on a reduced selection of voting probabilities for categorical assignment. Normalized entropies can be calculated from the maximum and minimum voting probabilities for ranking (Yun Chun, Chieh Shih, and Yu Kang 2021) or different formulation metrics (Kumar, Kumar, and Kapur 1986), but a more relevant metric incorporating the decisive candidates would be based on the top two categories by voting probability. With prediction maps, of interest is to determine what is the most probable category for assignment, and the contender would be the second most probable category. Second most probable tree species grouping maps can portray alternative categorical assignment, that can provide additional insight when considered in conjunction with entropy maps, so as to better comprehend factors impacting the vegetation cover within a study area.

Entropy scores tended to be lower for regions populated with black spruce or balsam fir, which was also noted in Table 3 with the validation sites. On the contrary, tamarack is a conifer species that sheds its needles during the autumn, but had high relative entropy scores, in part due to difficulty of modeling this species for the study area; tamarack is present either near clearings with compacted clayey soil, or at the periphery of wetlands in peat with low bulk densities. White spruce also had higher entropy scores, which can relate to the difficulty of distinguishing it from upland black spruce in certain contexts. Due to the prevalence of black spruce in the ARF, the visual implications for this tree species were more apparent from the corresponding prediction and entropy maps.

To elucidate the entropy calculations, we will demonstrate a hypothetical example. Consider the simple case of only at most four classification categories, for four different pixels. The respective posterior probabilities of assignment for the four separate classes are: pixel 1 has 0.90, 0.08, 0.01 and 0.01; pixel 2 has 0.75, 0.23, 0.01, and 0.01; pixel 3 has 0.51, 0.17, 0.16, and 0.16; and pixel 4 has 0.51, 0.47, 0.01, and 0.01. From these posterior probabilities, pixel 1 has the least uncertainty as it has an overwhelmingly majority assigned to one class. Pixel 2 has the second least uncertainty with assignment, as a super majority posterior probability is still allocated to one class. On the contrary, pixels 3 and 4 have higher uncertainty with prediction, with pixel 4 having the highest uncertainty as the top two contenders have almost equal probability. When calculating entropy by (1), the resulting entropy values for pixels 1, 2, 3, and 4 are 0.281, 0.466, 0.888, and 0.570, respectively; thus, yielding the irrational result of pixel 3 being assessed with a higher prediction uncertainty than pixel 4. When considering the normalized entropy of the top two contenders via (4), firstly, the normalized top probabilities are calculated by (2) and the normalized second-most probabilities are computed by (3). These normalized top two probabilities, rounded, for each pixel case here are: 0.92 and 0.08 for pixel 1; 0.77 and 0.23 for pixel 2; 0.75 and 0.25 for pixel 3; and 0.52 and 0.48 for pixel 4. Substituting these values into (4) yields normalized entropies for pixels 1, 2, 3, and 4 of 0.408, 0.786, 0.811, and 0.999, respectively. Hence, pixel 4 is now appropriately evaluated to have the highest uncertainty with prediction. Additionally, one can note that the uncertainty values for the pixels are correspondingly higher with the normalized entropy based upon just the top two choices, than with the conventional entropy metric; thus, preventing a dilution of the measure of uncertainty.

4.3. Accuracy assessment and considerations

Pixel-based classification can yield a more effective evaluation if based on a site level, rather than just pixel level. The objective should be to employ methods for training models corresponding to many sites each of a purer selection of pixels for a tree species grouping, rather than large individual selection blocks of pixels coinciding to a grouping. Ultimately, what is of interest for a pure stand or tree species grouping is to determine whether the tree species class for the site on the whole was correctly predicted, rather than determining the accuracies of the individual pixels corresponding to that site. It is more practical to ensure that a small area (and hence the subsequent selection of pixels) corresponds to a pure stand, rather than a relatively large continuous selection (which in reality would still only constitute as one site). If models are only evaluated at a pixel level, apart from considering sites, there is the risk of overfitting models if pixels utilized for training and testing are intermingled from the same plots. Without aggregating by site level first, there is the potential for greatly inflated pixel level accuracies; thus, accuracy should be evaluated at a site level.

A primary advantage of utilizing finer scale features is that 2-m spatial resolution would result in pixels that are more likely to correspond to individual trees, which is evident with WV2 imagery. Finer spatial resolution scale predictors would enable for more representation of tree species endemic for a region, rather than just the most dominant tree species for the region on the whole. For the ARF, black spruce is the most prevalent tree species, followed by balsam fir, which is also observed from site counts in Tables 2 and 3. From the National Forestry Institute (NFI) model prediction output (Beaudoin et al. 2014, 2018), black spruce is over predicted for the ARF, whereas other tree species have minimal predicted concentration. The NFI analysis was based on MODIS imagery with a minimum spatial resolution of 250 m, which is too coarse to generalize for certain localities within the ARF. However, issues other than spatial resolution can impact a tree species classification, which would include data limitations and class imbalance, especially with regard to tree species that are rarer or of less presence.

Random stratified sampling can be effective when there is an unequal distribution of tree species groupings in a region; random stratified sampling can be employed to ensure a more representative portion of tree species for both training and testing purposes (Gimaret-Carpentier et al. 1998). For this study area, black spruce was the dominant tree species for a considerable proportion of sites, followed by balsam fir. Employing random stratified sampling permitted a more equitable split for tree species groupings of lower incidences, such as for tamarack and white spruce. The adoption of random stratified sampling would extent this research to other areas within the boreal forest where tree species populations are highly unequal in distribution.

The accuracy of the GPS receivers operated for field campaigns might be problematic if finer scale predictors are utilized. For commercial grade GPS receivers, accuracy can be limited to the spatial scale of 10 m (Curtis and Wing 2012). With Landsat imagery of 30-m spatial resolution, employing lower accuracy GPS receivers may not be a concern, but this definitely would be of consequence for WV2 imagery of 2-m spatial resolution. To minimize the occurrence of noisy labels for a tree species classification, it is important to select pixels coinciding to pure stands (or characteristic groupings) of tree species for the region. Choosing pixels that correspond to a station near the center along a transect, particularly if the same leading tree species can be ascertained for the neighboring stations of that transect, can better guarantee a selection conforming to the tree species grouping of interest. However, for larger study areas the incidence of noisy pixels may not necessarily be averted. The issue of noisy labels can be addressed by hierarchical constrained energy minimum (HCEM) methods (Bing et al. 2020) or by weakly supervised multiple instance learning (MIL) (Carbonneau et al. 2018). Pixels with noisy labels can also arise due to tree species observation data that are no longer recent or the most relevant; so time span considerations in conjunction with geospatial accuracy need to be taken into account.

A 70% or greater composition belonging to the same tree species for a site, was specified as the threshold to correspond to a pure stand. This level was considered because entirely pure stands can be rare (particularly for relatively larger areas) in a natural forested environment; for the ARF study area, the overwhelming majority of sites that comprised of pure stands existed in the wetlands where black spruce is the only tree species which can endure for the long term. For eastern white cedar, this limit was lowered to 50%, as its occurrence is scarce within the ARF. Bravo-Oviedo et al. (2014) write that in the recent past, a site was considered a pure stand if at least 80% of the tree species in the top story were composed of the same species. Even if a region is a pure stand, there can exist gaps in the top canopy layer, resulting in surface reflectance for optical imagery coinciding to the underlayer. There can also be misclassification of tree species, especially if differentiating between tree species belonging to the same genus; the error of misclassification warrants the relaxation of stringent thresholds for pure stands. To reduce this occurrence of error, it is essential to employ field data collected by forestry experts and specialists.

Pixels for black spruce pure stands were selected from both wetlands and non-wetland environments. Although for this region, it is presumed that mixed black spruce localities would likely correspond to upland black spruce, for a future study it might be warranted to explicitly differentiate between lowland black spruce and upslope (i.e. upland) black spruce (El Abidine et al. 1994) with field observation data. This distinction could also alleviate the effect of the overrepresentation of one characterization of black spruce; lowland black spruce exists in the wetlands and is usually stunted, whereas upland black spruce can grow to taller heights comparable with other endemic tree species. Separating black spruce into categories for lowland and upslope variants would extent the applicability of this classification research for a larger domain of the boreal forest within North America where the tree species in this study are prevalent.

5. Conclusions

A pixel-based classification was carried out for a natural forest region within the boreal biome of northern Ontario in Canada. RF and SVM approaches were trained from plots corresponding to pure stands and mixed groupings from numerous sites, and classification accuracy was ultimately assessed by site rather than by the pixel level. These methods were able to obtain reasonable accuracies for the prediction of tree species groupings within the ARF, with accuracies exceeding 0.7 for correctness and 0.5 for kappa score for the best models.

Features relating to topography were determined to have high variable importance. Specifically, the features of channel network base level, valley depth and MRVBF had the highest relevance, with mean decrease in Gini of 1929.13, 1426.82, and 1140.33, respectively. NDVI and CHM also attained high variable importance, with a mean decrease in Gini of 1124.26 for NDVI and 763.63 for CHM.

Entropy maps discerned regions with regard to prediction uncertainties. In this study, the areas of lowest uncertainties correspond to the wetland environments where black spruce is generally the only tree species present. An entropy map generated from the top two normalized prediction probabilities can assist to better resolve areas that have higher uncertainties with prediction, from amongst regions with lower uncertainties. The normalized entropy map aided in ascertaining areas of mixed tree species groupings, but also in potentially differentiating upslope black spruce from the lowland black spruce characterization, in tandem with tree species prediction maps, within the ARF. The regions of higher uncertainties can coincide with greater tree species diversity, but consequently with where mixed pixels can be of concern; thus, indicating where further approaches might have to be investigated for classification improvement.

Acknowledgments

LiDAR data were obtained from the Ontario Ministry of Natural Resources and Forestry (MNRF) via Land Information Ontario, and contain information licensed under the Open Government license – Ontario. Site data for tree species concentrations were obtained from the Ontario Forest Resources Inventory (FRI). WorldView-2 imagery was purchased from Maxar Technologies.

Disclosure statement

Both authors declare that they have no known competing financial or non-financial interests to report that could have appeared to influence the research presented in this paper.

Data availability statement

The LiDAR-derived feature data that support the findings of this study is openly available in figshare at https://dx.doi.org/10.6084/m9.figshare.21282675. The satellite imagery data are not publicly available.

Additional information

Funding

This research was funded by the Ontario Ministry of Agriculture, Food and Rural Affairs (OMAFRA) [Funding Number ND2017-3179] and the Natural Sciences and Engineering Research Council of Canada (NSERC) [Discovery Grant].