TY - JOUR AU1 - Kotivuori, Eetu AU2 - Maltamo, Matti AU3 - Korhonen, Lauri AU4 - Strunk, Jacob L AU5 - Packalen, Petteri AB - Abstract In this study we investigated the behaviour of aggregate prediction errors in a forest inventory augmented with multispectral Airborne Laser Scanning and airborne imagery. We compared an Area-Based Approach (ABA), Edge-tree corrected ABA (EABA) and Individual Tree Detection (ITD). The study used 109 large 30 × 30 m sample plots, which were divided into four 15 × 15 m subplots. Four different levels of aggregation were examined: all four subplots (quartet), two diagonal subplots (diagonal), two edge-adjacent subplots (adjacent) and subplots without aggregation. We noted that the errors at aggregated levels depend on the selected predictor variables, and therefore, this effect was studied by repeating the variable selection 200 times. At the subplot level, EABA provided the lowest mean of root mean square error (⁠|$\overline{\mathrm{RMSE}}$|⁠) values of 200 repetitions for total stem volume (EABA 21.1 percent, ABA 23.5 percent, ITD 26.2 percent). EABA also fared the best for diagonal and adjacent aggregation (⁠|$\overline{\mathrm{RMSE}}$|⁠: 17.6 percent, 17.4 percent), followed by ABA (⁠|$\overline{\mathrm{RMSE}}$|⁠: 19.3 percent, 18.2 percent) and ITD (⁠|$\overline{\mathrm{RMSE}}$|⁠: 21.8, 21.9 percent). Adjacent subplot errors of ABA were less correlated than errors of diagonal subplots, which resulted also in clearly lower RMSEs for adjacent subplots. This appears to result from edge tree effects, where omission and commission errors cancel for trees leaning from one subplot into the other. The best aggregate performance was achieved at the quartet level, as expected from fundamental properties of variance. ABA and EABA had similar RMSEs at the quartet level (⁠|$\overline{\mathrm{RMSE}}$| 15.5 and 15.3 percent), with poorer ITD performance (⁠|$\overline{\mathrm{RMSE}}$| 19.4 percent). Introduction Airborne laser scanning (ALS) based forest inventories are commonly categorized as either Area-Based Approach (ABA) or as Individual Tree Detection (ITD; Vauhkonen et al. 2014). In ABA inventories, stand attribute models are fitted with sample plots using metrics calculated from ALS data. These models are applied in the target area using grid cells as basic inventory units and the predictions are typically aggregated to forest stands. In Nordic countries, ALS data are used together with aerial images to predict tree species (Maltamo et al. 2015; Packalén and Maltamo 2006). In the past decade, different ABA variants have been used in operational stand-level forest management inventories in Norway, Finland and Sweden (e.g. Næsset 2014; Maltamo and Packalen 2014; Nilsson et al. 2017). While ABA is the primary forest inventory approach (McRoberts et al. 2014), like ABA, ITD has also been the subject of considerable research effort. ITD methods can be broken down into four steps: 1) tree detection and delineation; 2) feature extraction; 3) tree species recognition and 4) estimation of stem attributes (Vauhkonen 2010). Most ITD studies have focused on individual steps (Vauhkonen et al. 2014), with a smaller number of ITD studies looking more comprehensively at ITD inventory (e.g. Korpela et al. 2007a; Duncanson et al. 2014; Silva et al. 2016). Combinations of the ABA and ITD approaches have also been considered in literature (e.g. Breidenbach et al. 2010; Xu et al. 2014; Packalen et al. 2015). In ITD, detected crown segments are usually assumed to represent individual trees. However, Breidenbach et al. (2010) examined an ITD approach, where a single crown segment could include one, several or no tree observations. They used nearest neighbour imputation (NN) to predict trees for each crown segment. In turn, Xu et al. (2014) calibrated ABA diameter distributions using ITD. The use of detected tree crowns to improve ABA predictions is an appealing option. Packalen et al. (2015) developed the Edge-tree corrected ABA (EABA) approach. It uses tree positions obtained from the first step of ITD (tree detection and delineation) to adjust plot and grid cell boundaries before ALS metrics are computed. In Packalen et al. (2015), EABA decreased the RMSE of stem volume prediction by two percent points over the traditional ABA approach (17.0 percent vs. 19.1 percent). However, EABA requires high density ALS data, at least 5 pulses m−2, for tree delineation, which is considerably higher and more expensive than the 1 pulses m−2 which is considered sufficient for traditional ABA (e.g. Magnussen et al. 2010; Strunk et.al. 2012). Aggregation is highly relevant to forest inventory because it is the basis for many forest management decisions, which are typically made at some aggregate level, rather than at the tree, plot or grid-cell level. Depending on the approach, variables are predicted for trees or grid cells (e.g. 15 × 15 m) and then aggregated to forest stands (Næsset 2002; White et al. 2013; Maltamo and Packalen 2014). Most studies report prediction errors at the tree or plot level, although in forest management inventories the stand level is the most relevant because treatment decisions are typically made at that level. In stand-level inventories, an established measure of prediction error for continuous variables is RMSE, the square root of the mean squared error, which is estimated for a set of stands. Error in an individual stand can be expressed, for example, with confidence intervals estimated by the model-based variance estimator (Breidenbach et al. 2016). In both cases, we need to account for how the error decreases when the tree- or cell-level predictions are aggregated at units with different size. The motivation for this study arises from our observation that the errors decrease differently when aggregation is performed with different inventory approaches. We henceforth use the word ‘aggregation’ to refer to the process of computing the mean from a group of subplots or grid cells. Technically, in the case of ITD, another aggregation step takes place, which is the summation of ITD objects to the subplot or grid cell level; however, this process is not considered in our use of the term ‘aggregation’. The objective of this study is to demonstrate the behaviour of the aggregate prediction errors for ITD, ABA and EABA. Errors are assessed using different aggregates of 15 × 15 m subplots clipped from larger 30 × 30 m plots. Aggregation scenarios include individual 15 × 15 m subplots (no aggregation), aggregation of two adjacent or diagonal subplots and aggregation of all four subplots. Specific research questions are as follows: 1) How much does aggregation of stem volume from subplots decrease the error rates over predictions without aggregation and does this vary by species? 2) Does the aggregate prediction error behave differently in ABA, EABA and ITD? 3) Does the neighbourhood type (diagonal, edge-adjacent) affect the aggregate prediction error in ABA, EABA and ITD? We repeated the entire analysis 200 times with heuristic variable selection and investigated how the different predictor variables influenced the behaviour of errors in aggregation. Materials and methods Field data The study area is located in southeast Finland (Figure 1), in south-boreal zone (62° 31′ N 29° 23′ E) 90–160 m above the sea level. The area is coniferous dominated (pine and spruce), but a minority of the forests are dominated by deciduous species. The field measurements were made between June and August of 2017. A total of 109 plots were placed in the area following the estimated distribution of the main tree species (pine, spruce and deciduous trees) and development classes. The development classes were estimated using the 2016 ALS data and the main tree species were taken from the publicly available National Forest Inventory rasters (Natural Resources Institute Finland 2013; Tomppo and Halme 2004). Stands with mean height < 7 m were omitted from the study. Figure 1 Open in new tabDownload slide Location of (a) Liperi municipality study area, and (b) the 109 measured 30 × 30 m field plots (red squares). Base maps: Countries 1:1 M (a)/© Euro Geographics 2013; General map 1:1 M (b)/© National Land Survey of Finland 2018. Figure 1 Open in new tabDownload slide Location of (a) Liperi municipality study area, and (b) the 109 measured 30 × 30 m field plots (red squares). Base maps: Countries 1:1 M (a)/© Euro Geographics 2013; General map 1:1 M (b)/© National Land Survey of Finland 2018. Trees were measured on 109 square 30 × 30 m sample plots, which were placed entirely within individual forest stands. Forest stands were delineated by the maturity of the trees and by tree species composition using the pre-information from areal images and ALS-based Canopy Height Models (CHM). Tree measurements included tree species, diameters at breast height, heights and positions for all trees with a diameter ≥5 cm. Stem positions were mapped using the procedure described by Korpela et al. (2007b). The stem-mapping procedure uses (CHM) derived ITD tree objects as reference locations, and other trees are mapped relative to the distances and azimuths of ITD reference trees. The stem map of individual tree locations was then used to assign individual tree records to 15 × 15 m subplots (i.e. total of 436 subplots). Individual stem volumes were predicted using the Laasasenaho (1982) stem volume equations. While conifer tree volumes were aggregated by species, deciduous tree species were grouped together, which is the current practice in Finnish stand-level management inventories. Summary statistics of stand attributes for the 30 × 30 m plots are shown in Table 1. Table 1 Summary statistics of stand attributes in the 30 × 30 m plots. . Hg . Dg . N . G . V . V pine . V spruce . V decid . Min 10.2 9.8 188.9 7.8 67.0 0.0 0.0 0.0 25% quartile 15.0 16.5 566.7 18.7 138.2 0.0 10.5 2.4 Median 18.7 21.8 916.7 22.5 205.2 44.9 45.3 11.8 Mean 19.1 22.0 1068.2 23.6 220.8 82.6 97.8 40.4 75% quartile 22.4 25.6 1405.6 27.5 264.3 146.0 149.5 44.6 Max 30.9 42.2 3866.7 45.5 593.7 384.8 439.9 216.4 SD 5.0 6.7 619.7 7.7 108.5 96.1 113.4 57.9 . Hg . Dg . N . G . V . V pine . V spruce . V decid . Min 10.2 9.8 188.9 7.8 67.0 0.0 0.0 0.0 25% quartile 15.0 16.5 566.7 18.7 138.2 0.0 10.5 2.4 Median 18.7 21.8 916.7 22.5 205.2 44.9 45.3 11.8 Mean 19.1 22.0 1068.2 23.6 220.8 82.6 97.8 40.4 75% quartile 22.4 25.6 1405.6 27.5 264.3 146.0 149.5 44.6 Max 30.9 42.2 3866.7 45.5 593.7 384.8 439.9 216.4 SD 5.0 6.7 619.7 7.7 108.5 96.1 113.4 57.9 Note: Dg = basal area weighted mean diameter (cm), Hg = basal area weighted mean height (m), N = stem number (trees ha−1), G = basal area (m2 ha−1), V = total stem volume (m3 ha−1), V pine = stem volume of pine trees (m3 ha−1), V spruce = stem volume of spruce trees (m3 ha−1) and V decid = stem volume of deciduous trees (m3 ha−1). Open in new tab Table 1 Summary statistics of stand attributes in the 30 × 30 m plots. . Hg . Dg . N . G . V . V pine . V spruce . V decid . Min 10.2 9.8 188.9 7.8 67.0 0.0 0.0 0.0 25% quartile 15.0 16.5 566.7 18.7 138.2 0.0 10.5 2.4 Median 18.7 21.8 916.7 22.5 205.2 44.9 45.3 11.8 Mean 19.1 22.0 1068.2 23.6 220.8 82.6 97.8 40.4 75% quartile 22.4 25.6 1405.6 27.5 264.3 146.0 149.5 44.6 Max 30.9 42.2 3866.7 45.5 593.7 384.8 439.9 216.4 SD 5.0 6.7 619.7 7.7 108.5 96.1 113.4 57.9 . Hg . Dg . N . G . V . V pine . V spruce . V decid . Min 10.2 9.8 188.9 7.8 67.0 0.0 0.0 0.0 25% quartile 15.0 16.5 566.7 18.7 138.2 0.0 10.5 2.4 Median 18.7 21.8 916.7 22.5 205.2 44.9 45.3 11.8 Mean 19.1 22.0 1068.2 23.6 220.8 82.6 97.8 40.4 75% quartile 22.4 25.6 1405.6 27.5 264.3 146.0 149.5 44.6 Max 30.9 42.2 3866.7 45.5 593.7 384.8 439.9 216.4 SD 5.0 6.7 619.7 7.7 108.5 96.1 113.4 57.9 Note: Dg = basal area weighted mean diameter (cm), Hg = basal area weighted mean height (m), N = stem number (trees ha−1), G = basal area (m2 ha−1), V = total stem volume (m3 ha−1), V pine = stem volume of pine trees (m3 ha−1), V spruce = stem volume of spruce trees (m3 ha−1) and V decid = stem volume of deciduous trees (m3 ha−1). Open in new tab Remotely sensed data The ALS data for the Liperi study area were collected between 2 July and 10 July 2016 under leaf-on conditions using an Optech Titan multispectral ALS instrument that provides measurements at three different wavelengths: channel 1 (ch1) = 1550 nm, channel 2 (ch2) = 1064 nm and channel 3 (ch3) 532 nm. The flying height was 850 m above ground level, the half scan angle was 20° and the pulse repetition frequency was 250 000 Hz. The nominal swath width was 650 m, the lateral overlap between flight lines was 55 percent, and footprints for channels 1, 2 and 3 were 0.3, 0.3 and 0.6 m, respectively. Average pulse densities were 13.9, 14.0 and 8.3 points m−2 for channels 1, 2 and 3, respectively. Optech Titan channels are aligned such that the 1064 nm channel points at nadir, the 1550 channel points 3.5° forward of nadir and the 532 nm channel points 7° forward of nadir. We did not calibrate LiDAR intensities because an earlier study with the same data showed that it did not improve tree species discrimination (Kukkonen et al. 2019). The ALS echoes from channel 2 were classified as ground and other echoes using the approach presented by Axelsson (2000). The echo heights from channels 1, 2 and 3 were then scaled to above ground level using the digital terrain model interpolated from ground echoes by Delaunay triangulation. Negative echo heights were set to zero. Hereafter, ALS height refers to height above ground level. Kukkonen et al. (2019) showed that the fusion of multispectral ALS data and spectral values from aerial images provided lower error rates in species specific predictions compared with the use of multispectral ALS data only. Therefore, the aerial image data were also included in our study. Aerial images were captured for the Liperi study area on 23 May and 24 May 2016 using a DMC Z/I Intergraph (01-0128) digital aerial camera. The flying altitude was 4100 m above ground level and the nominal side and end laps were 30 and 80 percent, respectively. The camera had a 300 mm focal length lens and a CCD sensor of 1920 × 3456 pixels (12 × 12 μm pixel size) for each colour band (red, green, blue and near infrared). This setting resulted in a ground sample distance of ~150 cm. Multispectral images were not pan sharpened or orthorectified, i.e. the panchromatic band was not used. External orientation parameters were determined using the bundle block adjustment techniques presented by Mikhail et al. (2001). ALS metrics We calculated two sets of ALS metrics based on echo type including first echo metrics (F = ‘first of many’ + ‘only’) and last echo metrics (L = ‘last of many’ + ‘only’). Specific metrics were means (hmeanF and hmeanL), medians (hmedF and hmedL) and standard deviations (SD) (hsdF and hsdL) of height, height quantiles (h40/50/60/70/80/90F and h40/50/60/70/80/90L), covers (veg2/5/10/15F and veg2/5/10/15L) and echo proportions (echo%F and echo%L). Quantiles were calculated as proposed by Hyndman and Fan (1996). Cover metrics were computed as the ratios of the numbers of echoes under height thresholds, divided by the total number of echoes in the same echo category. Echo proportions were calculated as the ratios of first (F) and last (L) echoes by the number of all echoes (‘first of many’ + ‘intermediate’ + ‘last of many’ + ‘only’). Corresponding metrics were also computed for intensity values including means (intmeanF and intmeanL), medians (intmedF and intmedL), SDs (intstdF and intstdL) and quantiles (int40/50/60/70/80/90F and int40/50/60/70/80/90L). Metrics were computed separately for each LiDAR channel (ch1, ch2 and ch3). In ITD (see chapter ‘Individual tree detection’), a 2 m cut-off was used with all metrics other than covers. In ABA and EABA (see chapters ‘Area-based approach’ and ‘Edge-tree corrected area-based approach’), no cut-off for the ALS metrics was used. Aerial image metrics Spectral values from unrectified aerial images were extracted to ALS echoes using photogrammetry, as proposed in Packalén et al. (2009) and examined in Valbuena et al. 2011. Because of overlap in the images, each ALS echo was seen in multiple images. For each ALS echo, we calculated a mean pixel value from all of the images in which the ALS echo was observed. Aerial image metrics were then calculated from these averaged pixel values and assigned to first return ALS echoes (F) located ≥2 m above ground level (other echoes were not considered). Calculated metrics were maximum, minimum, mean and SD by band (R/G/B/NIR), separately for each ALS channel (ch1, ch2 ch3). Area-based approach Species-specific plot volumes (pine, spruce and deciduous trees) were predicted simultaneously with k-NN imputation using most similar neighbor (MSN) distance (Moeur and Stage 1995). The method has been used in similar forests (Packalén and Maltamo 2007; Breidenbach et al. 2010) and is also used operationally in Finland (Maltamo and Packalen 2014). In this study, the predicted values for the volumes of pine, spruce and deciduous trees were calculated as weighted means of k-NNs (Haara et al. 1997). After initial tests (not shown), the number of NNs (k) was set to five. In ABA, the ALS and aerial image metrics used as predictor variables were extracted for the subplots (Figure 2). Figure 2 Open in new tabDownload slide Plot footprints for the three inventory approaches for the same plot where the red lines represent the original square plot (subplot) footprint. For ABA approach, the area from which ALS metrics are computed is the same as plot footprint. For EABA approach, the green and striped areas were included but in striped areas the heights of echoes were set to zero before computing ALS metrics. For ITD approach, only blue areas were considered. Figure 2 Open in new tabDownload slide Plot footprints for the three inventory approaches for the same plot where the red lines represent the original square plot (subplot) footprint. For ABA approach, the area from which ALS metrics are computed is the same as plot footprint. For EABA approach, the green and striped areas were included but in striped areas the heights of echoes were set to zero before computing ALS metrics. For ITD approach, only blue areas were considered. Predictor variables for k-NN models were selected using a metaheuristic algorithm VSSA proposed by Packalén et al. (2012). In VSSA, the premise is to minimize the squared error loss by solving the NN model repeatedly with a slightly different set of predictor variables. VSSA is a multivariate optimization method (simulated annealing; Kirkpatrick et al. 1983) that minimizes the mean relative RMSE values over all response variables. In our analysis, we used 4000 iterations in the VSSA to select seven predictor variables from 174 candidate variables. Square root, natural logarithm, inverse and third-order polynomial transformations for candidates were also tested. The decision to use seven predictor variables was based on preliminary tests (not shown) and corresponds to the number of variables used in a related study by Packalén et al. (2012). Variable selection was repeated 200 times, and each repetition was used as described in chapter ‘Model validation’. Individual tree detection Tree identification was based on the protocols described by Pitkänen et al. (2004) and Packalen et al. (2015). An initial CHM raster was created by setting a pixel value either to the maximum ALS echo height within each 0.5 m cell, or as NoData if there were no ALS echoes within a pixel. A median filtering with a 3 × 3 neighbourhood was then run repeatedly on the NoData pixels until none remained. We next identified ‘pit’ pixels; a pixel was considered as a pit pixel if at least six of its neighbours in the 3 × 3 neighbourhood were >5 m higher than the pixel itself. Pit pixels were then replaced with the medians of the 3 × 3 neighbourhood, excluding itself. The CHM was then low-pass filtered using a Gaussian kernel. The scale of Gaussian filtering depended on the 80th height percentile of first echoes (h80F). We used a kernel with a sigma (σ) of 0.5 when h80F < 20 m and 0.8 when h80F ≥ 20 m. Trees were identified from CHM using watershed segmentation. Watershed segments were delineated using a drainage direction following an algorithm (Narendra and Goldberg 1980; Gauch 1999) on the inverted CHM (filtered as explained earlier). To maintain separation between the crown and ground, pixels with heights of <2 m were masked from the original segments. Finally, small segments having <4 pixels were merged to the most similar in height neighbour segment. An ALS-detected tree segment was considered to be in a plot if its local maximum was located inside the plot borders (Figure 2). Before the modelling, we linked each ALS tree to one field measured tree using the following conditions: 1) the local maximum was located <1 m from the field tree in the horizontal plane and 2) the height difference between the field measured height and the CHM local maximum was ≥0 and <2 m. If the conditions were not met, ALS tree was not linked. The objective of linking trees was to enable modelling field measured attributes from auxiliary remote sensing measurement for ITD objects. The models were then used to predict forest attributes for all ITD trees. For all ITD objects we calculated the same suite of ALS and aerial image metrics as were computed for ABA plots (see chapters ‘Airborne laser scanning metrics’ and ‘Aerial image metrics’). The only difference was that we used a 2-m cut-off when computing ALS metrics other than cover ratios. Prediction of individual ITD tree object attributes was performed with k-NN imputation using MSN distance as in ABA. The tree-level k-NN model used species-specific volumes (dm3) of linked trees as response variables. Because a tree cannot have more than one species, we set the number of NNs to one (k = 1) in the k-NN imputation. Predictor variables for k-MSN were chosen by VSSA in the same manner as explained in the context of ABA. At the tree level, we selected five variables and repeated variable selection 200 times. Square root, natural logarithm, inverse and third-order polynomial transformations for candidate predictor variables were also tested. Edge-tree corrected area-based approach The motivation for EABA is that some of the foliage for field measured trees extends beyond the plot border, and some of the foliage from un-measured trees extends into field plots. This results in edge tree effects due to errors of omission and commission for ABA. In EABA, predicted tree positions and crown shapes are used to adjust plot and grid cell boundaries to better encompass only measured trees (Packalen et al. 2015). Remote sensing metrics are then computed for ALS echoes within the ITD-adjusted plot and grid cell boundaries. EABA is an inventory variant that contains components of both ABA and ITD methods. Detection of individual trees for EABA was conducted as described in the previous chapter. The edge-tree correction in one of the subplots is shown in Figure 2. The plot edges were adjusted using the detected tree segments. If the local maximum of a segment was located inside the original plot, the border of the plot was extended to include the whole segment. Respectively, if a local maximum was located outside of the plot, the heights of the ALS echoes within the segment were set to zero (i.e. to ground level). Otherwise the EABA was implemented the same as ABA. Figure 3 Open in new tabDownload slide Levels of aggregation. Shaded squares indicate 15 × 15 m plots (subplots) and their combinations. Figure 3 Open in new tabDownload slide Levels of aggregation. Shaded squares indicate 15 × 15 m plots (subplots) and their combinations. Levels of aggregation We used individual plot level as the baseline and aggregated the predictions to two different levels using two different neighbourhood types (Figure 3). The levels are termed as follows: individual = one 15 × 15 m plot (no aggregation). diagonal = two 15 × 15 m plots connected diagonally. adjacent = two 15 × 15 m plots connected on a side horizontally or vertically. quartet = four 15 × 15 m plots i.e. the whole 30 × 30 m plot. While our investigation focuses on the empirical properties of aggregation, our understanding of aggregation is improved by looking at theoretical variance properties. This proves critical in understanding our results when we compare the results for aggregation of diagonal and adjacent pixels. The behaviour of variance when two random variables are added can be defined with respect to their correlation (⁠|$\rho$|⁠): $$\begin{equation} \mathrm{V}\left(\mathrm{X}+\mathrm{Y}\right)=\mathrm{V}\left(\mathrm{X}\right)+\mathrm{V}\left(\mathrm{Y}\right)+2\rho \sqrt{\mathrm{V}\left(\mathrm{X}\right)\times \mathrm{V}\left(\mathrm{Y}\right)\ } \end{equation}$$(1) $$ \mathrm{V}\left(\mathrm{X}\right)=\mathrm{the}\ \mathrm{variance}\ \mathrm{of}\ \mathrm{X} $$ $$ \mathrm{V}\left(\mathrm{Y}\right)=\mathrm{the}\ \mathrm{variance}\ \mathrm{of}\ \mathrm{Y} $$ $$ \rho =\mathrm{the}\ \mathrm{correlation}\ \mathrm{between}\ \mathrm{X}\ \mathrm{and}\ \mathrm{Y}. $$ For our purposes, X and Y can represent errors from two adjacent or diagonal subplots. There should be no difference between |$V(X)$| and |$V(Y)$| for adjacent versus diagonal subplots, since they belong to the same pool of observation. This means that whether diagonal or adjacent subplots have lower variance when aggregated is determined by their respective correlations |$\rho$|⁠. If the correlation between diagonal subplots is less than the correlation between adjacent subplots, then the aggregation of diagonal subplots will have lower variance than the aggregation of adjacent subplots. In this paper, the correlations between diagonal subplots and adjacent subplots were calculated for ABA, ITD and EABA as follows: $$\begin{equation} {\rho}_{\mathrm{diagonal}}=\mathrm{cor}\left(\left[{e}_{ul},{e}_{ur}\right],\left[{e}_{lr},{e}_{ll}\right]\right), \end{equation}$$(2) $$ {\rho}_{\mathrm{adjacent}}=\mathrm{cor}\left(\left[{e}_{ul},{e}_{ll},{e}_{ul},{e}_{ur}\right],\left[{e}_{ur},{e}_{lr},{e}_{ll},{e}_{lr}\right]\right) $$ where |${e}_{ul}$| = vector of errors for upper left predictions, |${e}_{ur}$| = vector of errors for upper right predictions, |${e}_{ll}$| = vector of errors for lower left predictions, |${e}_{lr}$| = vector of errors for lower right predictions and [] denotes that vectors are concatenated. The correlations are presented as mean values of 200 repetitions (see chapter ‘Model validation’ for details). Another thing that we can learn from the theoretical variance properties is that the SD of a sum is generally smaller (in percent), than the sum of component SDs. This holds as long as the correlation between X and Y is <1.0. We demonstrate this property for the simplified case where: $$\begin{equation} \mathrm{V}\left(\mathrm{X}\right)=\mathrm{V}\left(\mathrm{Y}\right), \end{equation}$$(3) which yields $$ \mathrm{V}\left(\mathrm{X}+\mathrm{Y}\right)=\mathrm{V}\left(\mathrm{X}\right)+\mathrm{V}\left(\mathrm{X}\right)+2\rho \sqrt{\mathrm{V}\left(\mathrm{X}\right)\times \mathrm{V}\left(\mathrm{X}\right)\ } $$ |$ =2\times \mathrm{V}\Big(\mathrm{X}\Big)+2\rho \times \mathrm{V}\Big(\mathrm{X}\Big).$| This means that $$ \mathrm{SD}\left(\mathrm{X}+\mathrm{Y}\right)=\sqrt{2+2\rho}\times \mathrm{SD}\left(\mathrm{X}\right) $$ and $$ \mathrm{SD}\left(\mathrm{X}+\mathrm{Y}\right)<\mathrm{SD}\left(\mathrm{X}\right)+\mathrm{SD}\left(\mathrm{Y}\right) $$ $$ \mathrm{when}\ \rho <1. $$ For our case, this means that $$ \mathrm{SD}\left(\mathrm{Quartet}\right)<\mathrm{SD}\left(\mathrm{Adjacent}\right),$$ $$ \mathrm{SD}\left(\mathrm{Diagonal}\right) <\mathrm{SD}\left(\mathrm{Individual}\right)\!. $$ Model validation The performances of ABA, EABA and ITD were evaluated using the same leave one out cross-validation (LOOCV) approach. For each repetition of LOOCV, we omitted an entire 30 × 30 m plot, trained k-NN models with the remainder of the plots (ABA, EABA) or trees (ITD) and predicted the species-specific stem volumes for the plot (ABA, EABA) or the trees (ITD) in the plot that were omitted. The ABA and EABA imputations were conducted for the subplots, and the ITD imputations (dm3) were aggregated to the subplots (m3 ha−1) as well. The total stem volumes were calculated as a sum of the species-specific volumes. Due to the heuristic nature of the VSSA algorithm, a different combination of predictors is found almost every time. We noted that the errors at aggregated levels depend on the selected predictor variables, and therefore, this effect was studied by repeating the variable selection 200 times. Hereafter, these 200 repetitions are called LOOCV distributions. The repetitions were analysed for deeper understanding of the error behaviour in aggregation, although only one realization would be selected in a practical application case. Errors between observed and predicted values from each LOOCV repetition were used to compute performance statistics. Model performances were assessed using RMSE and bias: $$\begin{equation} \mathrm{RMSE}=100\times \sqrt{\sum_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\left({\hat{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{\mathrm{n}}}/\overline{\mathrm{y}} \end{equation}$$(4) $$\begin{equation} \mathrm{bias}=100\times \sum_{\mathrm{i}=1}^{\mathrm{n}}\frac{\left({\hat{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}\right)}{\mathrm{n}}/\overline{\mathrm{y}} \end{equation}$$(5) where ŷi is the predicted value of attribute y on (sub)plot i, yi is the observed value of attribute y on (sub)plot i, y--- is the observed mean of attribute y and n is the number of (sub)plots. In the first chapter of the results section, we report the mean RMSE (⁠|$\overline{\mathrm{RMSE}}$|⁠) and mean absolute bias (MAB) values of the LOOCV distributions: $$\begin{equation} \overline{\mathrm{RMSE}}=\frac{\sum_j^{200}{\mathrm{RMSE}}_j}{200} \end{equation}$$(6) $$\begin{equation} \mathrm{MAB}=\frac{\sum_j^{200}\mid{\mathrm{bias}}_j\mid }{200} \end{equation}$$(7) Absolute values were used instead of signed bias because otherwise negative and positive biases cancelled out. Results Prediction errors at the different aggregation levels Performances (⁠|$\overline{\mathrm{RMSE}}$| and MAB) for the different aggregation levels are provided in Table 2. EABA had the lowest |$\overline{\mathrm{RMSE}}$| values for total, pine and spruce volumes (not deciduous volume) across all levels of aggregation. ABA fared second best (volume |$\overline{\mathrm{RMSE}}$|⁠) after EABA for total, pine and spruce volumes. ABA and EABA differences in |$\overline{\mathrm{RMSE}}$| for total, pine and spruce volumes were all <3 percentage points, with quartet level aggregations generally have the smallest difference in performances and individual subplots having the greatest differences. For deciduous volumes, ABA out-performed EABA |$\overline{\mathrm{RMSE}}$| values; however, the maximum difference was only 1.2 percentage point across all aggregation levels. ITD consistently had the highest |$\overline{\mathrm{RMSE}}$| values across all species groupings and aggregation levels. Table 2 |$\overline{\mathrm{RMSE}}$| and MAB of 200 variable selections with the ABA, ITD and EABA at different levels of aggregation. . . V total . V pine . V spruce . V deciduous . |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB ABA Individual 23.5 0.7 65.7 1.8 57.5 1.1 64.9 2.0 Diagonal 19.3 0.7 55.3 1.8 47.1 1.1 54.0 2.0 Adjacent 18.2 0.7 54.0 1.8 46.7 1.1 53.4 2.0 Quartet 15.5 0.7 47.7 1.8 40.4 1.1 46.9 2.0 ITD Individual 26.2 2.4 69.1 1.5 69.9 6.3 83.7 2.2 Diagonal 21.8 2.4 60.3 1.5 62.0 6.3 69.4 2.2 Adjacent 21.9 2.4 61.1 1.5 62.2 6.3 70.4 2.2 Quartet 19.4 2.4 56.2 1.5 57.9 6.3 62.1 2.2 EABA Individual 21.1 0.9 63.8 2.1 56.2 1.1 65.0 2.1 Diagonal 17.6 0.9 53.0 2.1 46.1 1.1 55.2 2.1 Adjacent 17.4 0.9 53.2 2.1 45.8 1.1 53.7 2.1 Quartet 15.3 0.9 46.9 2.1 39.8 1.1 47.8 2.1 . . V total . V pine . V spruce . V deciduous . |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB ABA Individual 23.5 0.7 65.7 1.8 57.5 1.1 64.9 2.0 Diagonal 19.3 0.7 55.3 1.8 47.1 1.1 54.0 2.0 Adjacent 18.2 0.7 54.0 1.8 46.7 1.1 53.4 2.0 Quartet 15.5 0.7 47.7 1.8 40.4 1.1 46.9 2.0 ITD Individual 26.2 2.4 69.1 1.5 69.9 6.3 83.7 2.2 Diagonal 21.8 2.4 60.3 1.5 62.0 6.3 69.4 2.2 Adjacent 21.9 2.4 61.1 1.5 62.2 6.3 70.4 2.2 Quartet 19.4 2.4 56.2 1.5 57.9 6.3 62.1 2.2 EABA Individual 21.1 0.9 63.8 2.1 56.2 1.1 65.0 2.1 Diagonal 17.6 0.9 53.0 2.1 46.1 1.1 55.2 2.1 Adjacent 17.4 0.9 53.2 2.1 45.8 1.1 53.7 2.1 Quartet 15.3 0.9 46.9 2.1 39.8 1.1 47.8 2.1 Note: Predictions were obtained by LOOCV. Bold values represent the inventory approach with the best performance for a given aggregation level and response variable. Open in new tab Table 2 |$\overline{\mathrm{RMSE}}$| and MAB of 200 variable selections with the ABA, ITD and EABA at different levels of aggregation. . . V total . V pine . V spruce . V deciduous . |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB ABA Individual 23.5 0.7 65.7 1.8 57.5 1.1 64.9 2.0 Diagonal 19.3 0.7 55.3 1.8 47.1 1.1 54.0 2.0 Adjacent 18.2 0.7 54.0 1.8 46.7 1.1 53.4 2.0 Quartet 15.5 0.7 47.7 1.8 40.4 1.1 46.9 2.0 ITD Individual 26.2 2.4 69.1 1.5 69.9 6.3 83.7 2.2 Diagonal 21.8 2.4 60.3 1.5 62.0 6.3 69.4 2.2 Adjacent 21.9 2.4 61.1 1.5 62.2 6.3 70.4 2.2 Quartet 19.4 2.4 56.2 1.5 57.9 6.3 62.1 2.2 EABA Individual 21.1 0.9 63.8 2.1 56.2 1.1 65.0 2.1 Diagonal 17.6 0.9 53.0 2.1 46.1 1.1 55.2 2.1 Adjacent 17.4 0.9 53.2 2.1 45.8 1.1 53.7 2.1 Quartet 15.3 0.9 46.9 2.1 39.8 1.1 47.8 2.1 . . V total . V pine . V spruce . V deciduous . |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB |$\overline{\mathrm{RMSE}}$| MAB ABA Individual 23.5 0.7 65.7 1.8 57.5 1.1 64.9 2.0 Diagonal 19.3 0.7 55.3 1.8 47.1 1.1 54.0 2.0 Adjacent 18.2 0.7 54.0 1.8 46.7 1.1 53.4 2.0 Quartet 15.5 0.7 47.7 1.8 40.4 1.1 46.9 2.0 ITD Individual 26.2 2.4 69.1 1.5 69.9 6.3 83.7 2.2 Diagonal 21.8 2.4 60.3 1.5 62.0 6.3 69.4 2.2 Adjacent 21.9 2.4 61.1 1.5 62.2 6.3 70.4 2.2 Quartet 19.4 2.4 56.2 1.5 57.9 6.3 62.1 2.2 EABA Individual 21.1 0.9 63.8 2.1 56.2 1.1 65.0 2.1 Diagonal 17.6 0.9 53.0 2.1 46.1 1.1 55.2 2.1 Adjacent 17.4 0.9 53.2 2.1 45.8 1.1 53.7 2.1 Quartet 15.3 0.9 46.9 2.1 39.8 1.1 47.8 2.1 Note: Predictions were obtained by LOOCV. Bold values represent the inventory approach with the best performance for a given aggregation level and response variable. Open in new tab As with |$\overline{\mathrm{RMSE}}$|⁠, the performances of ABA and EABA with respect to MAB (Table 2) exceeded the performance of ITD. In ABA and EABA, MABs were similar, always within 0.3 percentage points of each other. The MAB across all methods was generally <2.5 percent. The only exception was for ITD spruce volumes which had a much larger MAB (6.3 percent). In a closer examination of ITD spruce volumes, we found that spruce volumes were nearly always underestimated (199 times out of 200). This was because shade tolerant spruce trees growing under larger trees were generally not detected by ITD, and there was a tendency for adjacent spruce crowns to be merged and remain un-differentiated by ITD. The effect of adjacent vs. diagonal plots In ABA, the |$\overline{\mathrm{RMSE}}$| values computed for adjacent plots were lower than for diagonal plots (Table 2). This effect is illustrated using results of 200 LOOCV distributions in Figure 4 (ABA, top-left), where the RMSE values of total volume predictions were consistently lower in adjacent plots than in diagonal plots. This phenomenon was caused by edge-tree effects along the shared borders between adjacent subplots. Branches and foliage of measured trees near the edge of a subplot commonly extend into the adjacent subplot. For example, 100 percent of the volume of a border tree is associated with a single subplot based on bole position, but the ALS could detect a significant proportion of the foliage in the adjacent subplot, even if the adjacent subplot is bare of trees. This would cause an over-estimation of volume in the bare subplot, and an under-estimation of volume in the subplot where the tree was measured. The edge-tree effect is weaker between diagonal subplots, which are only connected by their corners. In ABA, the error correlation was |$\rho =$| 0.20 between adjacent subplots and |$\rho =$| 0.35 between diagonal subplots (Table 3). From Formula 1 we know that the reduced correlation causes lower aggregate RMSE values for adjacent subplots than for diagonal subplots, as shown empirically in Table 2. Figure 4 Open in new tabDownload slide Relationship of RMSE associated with total stem volume predictions in adjacent and diagonal plots. Each dot represents an individual variable selection, totally 200 repetitions. Predictions were obtained by LOOCV. Relationships are presented separately for the ABA, EABA and ITD. Solid black line is the 1:1 line and the solid red line is the regression line. Figure 4 Open in new tabDownload slide Relationship of RMSE associated with total stem volume predictions in adjacent and diagonal plots. Each dot represents an individual variable selection, totally 200 repetitions. Predictions were obtained by LOOCV. Relationships are presented separately for the ABA, EABA and ITD. Solid black line is the 1:1 line and the solid red line is the regression line. Table 3 Mean correlations between adjacent and diagonal subplots on 200 variable selections in different inventory approaches. . |${\rho}_{\mathrm{adjacent}}$| . |${\rho}_{\mathrm{diagonal}}$| . ABA 0.20 0.35 ITD 0.39 0.39 EABA 0.36 0.40 . |${\rho}_{\mathrm{adjacent}}$| . |${\rho}_{\mathrm{diagonal}}$| . ABA 0.20 0.35 ITD 0.39 0.39 EABA 0.36 0.40 Open in new tab Table 3 Mean correlations between adjacent and diagonal subplots on 200 variable selections in different inventory approaches. . |${\rho}_{\mathrm{adjacent}}$| . |${\rho}_{\mathrm{diagonal}}$| . ABA 0.20 0.35 ITD 0.39 0.39 EABA 0.36 0.40 . |${\rho}_{\mathrm{adjacent}}$| . |${\rho}_{\mathrm{diagonal}}$| . ABA 0.20 0.35 ITD 0.39 0.39 EABA 0.36 0.40 Open in new tab In ITD, the shared plot border had no discernible effect on total volume RMSE values (Table 2), which is consistent with the fact that adjacent and diagonal subplots had the same error correlation (Table 3). In EABA, the differences in RMSE and correlation values between diagonal and adjacent aggregations were minimal (Table 3; Figure 4). Variation of RMSE due to different predictor variables The LOOCV distributions of RMSE were examined to see whether they differ between ABA, EABA and ITD. The quartet-level distributions associated with total volume are shown in Figure 5. The most common RMSE value was clearly lower in EABA than in ABA, although the |$\overline{\mathrm{RMSE}}$| values were similar in both approaches. This observation indicates that with the appropriate variable selection, EABA results in lower RMSE values also at the aggregated-level. However, the EABA distribution had a clear lower limit of 13.0 percent, but there were ABA models having lower RMSE values than the most accurate EABA model (ABA minimum 11.9 percent). The RMSE distributions of ITD were slightly positively skewed, and all of the RMSE values were >16 percent. Figure 5 Open in new tabDownload slide Distributions of RMSE for total stem volume with the ABA, ITD and EABA at the quartet level. Predictions were obtained by LOOCV and variable selection was repeated 200 times. Red vertical line represents the mean and blue line represents the median of the RMSE values. Note: skew = skewness, med = median. Figure 5 Open in new tabDownload slide Distributions of RMSE for total stem volume with the ABA, ITD and EABA at the quartet level. Predictions were obtained by LOOCV and variable selection was repeated 200 times. Red vertical line represents the mean and blue line represents the median of the RMSE values. Note: skew = skewness, med = median. Discussion Our results showed that the improvement in |$\overline{\mathrm{RMSE}}$| associated with total volume predictions from using EABA instead of ABA for individual 225 m2 plots (21.1 vs 23.5 percent) was similar to the RMSE values reported by Packalen et al. (2015) for 254 m2 plots (17.0 vs 19.1 percent). However, the EABA and ABA had similar performances at the quartet level, suggesting that EABA and ABA would provide similar stand level performances. Differences in species-specific |$\overline{\mathrm{RMSE}}$| values were marginal between ABA and EABA at all levels, although EABA typically provided lower |$\overline{\mathrm{RMSE}}$| values. We also found that for ABA, edge-tree related errors reduced the spatial correlation of adjacent plots, resulting in improved performance for adjacent subplot aggregation relative to diagonal aggregation. Since the intention of EABA is to specifically mitigate edge trees as a source of error, there was not an obvious difference between adjacent and diagonal plots. The effect of edge-tree errors on aggregation with ABA was surprisingly large: at the quartet level the difference in |$\overline{\mathrm{RMSE}}$| values between EABA and ABA was negligible (<1 percentage point). In ITD, there was no difference in aggregation performances from adjacent and diagonal plots because predictions are assembled from individual tree objects that are not impacted by omission and commission of trees along the plot border. The LOOCV distributions of RMSE differed considerably between ABA and EABA at the quartet level. The RMSE distribution of EABA was greatly right skewed, whereas ABA distribution was fairly symmetric. The median and mode of the EABA distribution were smaller than for ABA, although with EABA the right skewness means that a particular variable combination has a greater chance of performing poorly (e.g. exceeding 20 percent RMSE) than with ABA. In practice there does not appear to be an advantage to using EABA over ABA at the quartet level, which likely also applicable to the stand level. In this study, ITD had the largest (worst) |$\overline{\mathrm{RMSE}}$| value in every case. This was partly caused by the incorrect prediction of tree species in the k-NN imputation. The use of stem volume from an incorrect tree species leads to erroneous stem volumes when the predictions are aggregated from the tree to the plot level, similarly with total volume predictions. However, except for spruce, the MAB in ITD were small and similar to EABA and ABA MABs (<2.5 percent). The smaller bias observed here relative to what has been seen in other ITD studies (Vauhkonen et al. 2010) is likely due to the use of tree level k-NN imputation by tree species. Most other ITD studies used separate models for each step, such as species classification, tree diameter prediction and stem volume prediction that can cause great cumulative errors. We used multispectral ALS data, aerial images, field measurements and methods that we consider applicable in Nordic boreal forests. Our aim here was not to assess the contributions of different data sources, but to provide sufficient data for a fair comparison between each of the three methods that vary in their data requirements. Both ITD and EABA require higher-density ALS data than ABA, because individual trees need to be delineated from the data. In the past this was a deterrent to using ITD-related methods, although most recent ALS datasets are often sufficiently dense (>5 echoes per m2) for ITD-related approaches. Field data requirements also varied between methods. ITD as implemented here requires exact tree locations, because the field trees and the trees recognized from ALS data need to be linked to model tree attributes from remote sensing attributes. Tree positioning increases the inventory costs of an ITD inventory considerably, especially if the required number of trees is large. While EABA shares some requirements of ITD (high-density ALS data), it does not require field tree locations. In EABA the delineated tree objects are used to adjust where metrics are calculated, but there is no attempt to link field measured trees with ITD tree objects. Knowing the field tree locations can be, however, useful in EABA (or even in ABA) because it enables further adjustment of plot locations after positioning with Global Navigation Satellite Systems (GNSS; commonly GPS and GLONASS). Plots positioned with GNSS are not always accurate, especially when a plot is under a closed canopy (Mauro et al. 2011). An error of 2 m, for instance, likely cancels out the benefit of the edge-tree correction used in EABA methods if tree crowns on a plot are only a few metres wide. Our study has demonstrated the theoretical property that RMSE decreases with aggregation (Formula 3). However, the amount that aggregation can reduce variance is limited by bias. The magnitude of bias does not change with aggregation, assuming that bias is constant across an area of aggregation. This was apparent when we observed constant MAB values within an inventory method in Table 2, regardless of aggregation type. This property can also be seen theoretically from a bias-variance RMSE decomposition (not provided here; see e.g. Chapter 7.3 in Hastie et al. 2009). Assuming a constant bias, this means that systematic bias contributes a greater proportion to RMSE values as larger areas are aggregated. It is important to understand how the errors in ALS inventory behave at the aggregated level. Our results show that prediction errors associated with different ALS inventory methods decrease differently in aggregation. The forest stand is currently the primary level of aggregation in operational forestry. Therefore, the next step is to develop models to estimate stand-level accuracy when the particular ALS inventory method is used. Conclusions In this study, EABA was the best performing method at the individual plot level. However, after aggregation to the quartet level, the prediction errors associated with EABA and the traditional ABA were similar. In contrast, ITD provided the largest prediction errors for total stem volume and by species at each level of aggregation. Aggregation to adjacent and diagonal plots showed that edge-tree errors in ABA compensate for each other in adjacent plots. Minor compensation was also observed in edge-corrected EABA, whereas, the edge-tree effect between plots in ITD was nonexistent. EABA is a potential method for applications where plot-level accuracy is important. However, at stand-level forest management planning, the benefit of the current version of EABA over the traditional ABA approach seems to be limited. We conclude that the behaviour of prediction errors in aggregation is not the same with every method, which needs to be taken into account when estimating uncertainty at the aggregated level. Data Availability Statement The data (X and Y variables) of this article will be shared on reasonable request to the corresponding author. Original data are partly owned by a third party (not shared). Acknowledgements The authors thank the Finnish Forest Centre and the Strategic Research Council of the Academy of Finland (FORBIO project, Grant No. 314224) for supporting data collections in the Liperi test area. The authors also thank the anonymous reviewers for their comments that helped to improve the manuscript. Conflict of interest statement None declared. Funding Strategic Research Council of the Academy of Finland (FORBIO project, Grant No. 314224). References Axelsson , P. 2000 DEM generation from laser scanner data using adaptive TIN models . Int. Arch. Photogramm Remote Sens 33 , 110 – 117 . Google Scholar OpenURL Placeholder Text WorldCat Breidenbach , J. , Næsset , E., Lien , V., Gobakken , T. and Solberg , S. 2010 Prediction of species specific forest inventory attributes using a nonparametric semi-individual tree crown approach based on fused airborne laser scanning and multispectral data . Remote Sens. Environ. 114 , 911 – 924 . doi: 10.1016/j.rse.2009.12.004 . Google Scholar Crossref Search ADS WorldCat Breidenbach , J. , McRoberts , R.E. and Astrup , R. 2016 Empirical coverage of model-based variance estimators for remote sensing assisted estimation of stand-level timber volume . Remote Sens. Environ. 173 , 274 – 281 . doi: 10.1016/j.rse.2015.07.026 . Google Scholar Crossref Search ADS PubMed WorldCat Duncanson , L.I. , Cook , B.D., Hurtt , G.C. and Dubayah , R.O. 2014 An efficient, multi-layered crown delineation algorithm for mapping individual tree structure across multiple ecosystems . Remote Sens. Environ. 154 , 378 – 386 . doi: 10.1016/j.rse.2013.07.044 . Google Scholar Crossref Search ADS WorldCat Euro Geographics . 2013 Countries 1:1 M. Euro Geographics for the administrative boundaries, GISCO, Eurostat. http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/countries#countries [Cited 4 May 2018] Gauch , J.M. 1999 Image segmentation and analysis via multiscale gradient watershed hierarchies . IEEE Trans. Image Process. 8 , 69 – 79 . Google Scholar Crossref Search ADS PubMed WorldCat Haara , A. , Maltamo , M. and Tokola , T. 1997 The k-nearest-neighbour method for estimating basal-area diameter distribution . J. Scand. J. For. Res. 12 , 200 – 208 . doi: 10.1080/02827589709355401 . Google Scholar Crossref Search ADS WorldCat Hastie , T. , Tibshirani , R. and Friedman , J. 2009 The Elements of Statistical Learning . 2nd edn. Springer . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Hyndman , R.J. and Fan , Y. 1996 Sample quantiles in statistical packages . Am. Stat. 50 , 361 – 365 . Published by: Taylor & Francis, Ltd. on behalf of the American Statistical Association . doi: 10.2307/2684934 . Google Scholar Crossref Search ADS WorldCat Kirkpatrick , S. , Gelatt , C.D. and Vecchi , M.P. 1983 Optimization by simulated annealing . Science 220 , 671 – 680 . doi: 10.1126/science.220.4598.671 . Google Scholar Crossref Search ADS PubMed WorldCat Korpela , I. , Dahlin , B., Schäfer , H., Bruun , E., Haapaniemi , F., Honkasalo , J., et al. 2007a Single-tree forest inventory using LiDAR and aerial images for 3D treetop positioning, species recognition, height and crown width estimation. ISPRS Workshop on Laser Scanning 2007 and Silvi Laser 200 35 (Part 3/W5) : 227 – 233 . Korpela , I. , Tuomola , T. and Välimäki , E. 2007b Mapping forest plots: An efficient method combining photogrammetry and field triangulation . Silva Fennica 41 , 457 – 469 . doi: 10.14214/sf.283 . Google Scholar Crossref Search ADS WorldCat Kukkonen , M. , Maltamo , M., Korhonen , L. and Packalen , P. 2019 Comparison of multispectral airborne laser scanning and stereo matching of aerial images as a single sensor solution to forest inventories by tree species . Remote Sens. Environ. 231 . doi: 10.1016/j.rse.2019.05.027 . Google Scholar OpenURL Placeholder Text WorldCat Crossref Laasasenaho , J. 1982 Taper curve and volume functions for pine, spruce and birch . Communicationes Instituti Forestalis Fenniae 108 . 74 p. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Magnussen , S. , Næsset , E. and Gobakken , T. 2010 Reliability of LiDAR derived predictors of forest inventory attributes: a case study with Norway spruce . Remote Sens. Environ. 114 , 700 – 712 . doi: 10.1016/j.rse.2009.11.007 . Google Scholar Crossref Search ADS WorldCat Maltamo , M. , and Packalen , P. 2014 Species-specific management inventory in Finland. In: Maltamo , M., Næsset , E., and Vauhkonen , J. Forestry Applications of Airborne Laser Scanning - Concepts and Case Studies. Managing Forest Ecosystems . Springer, Springer Science + Business Media , pp. 241 – 252 . doi: 10.1007/978-94-017-8663-8 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Maltamo , M. , Ørka , H.O., Bollandsås , O.M., Gobakken , T. and Næsset , E. 2015 Using pre-classification to improve the accuracy of species-specific forest attribute estimates from airborne laser scanner data and aerial images . Scand. J. For. Res. 30 , 336 – 345 . doi: 10.1080/02827581.2014.986520 . Google Scholar Crossref Search ADS WorldCat Mauro , F. , Valbuena , R., Manzanera , J.A. and García-Abril , A. 2011 Influence of global navigation satellite system errors in positioning inventory plots for tree-height distribution studies . Can. J. For. Res. 41 , 11 – 23 . doi: 10.1139/X10-164 . Google Scholar Crossref Search ADS WorldCat McRoberts , R.E. , Andersen , H-E., and Næsset , E. 2014 Using airborne laser scanning data to support forest sample surveys. In: Maltamo , M., Næsset , E., and Vauhkonen , J. 2014. Forestry applications of airborne laser scanning - concepts and case studies. Managing Forest Ecosystems . Springer, Springer Science + Business Media , pp. 269 – 292 . doi: 10.1007/978-94-017-8663-8 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Mikhail , E.M. , Bethel , J.S. and McGlone , J.C. 2001 Introduction to Modern Photogrammetry . Wiley , p. 479 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Moeur , M. and Stage , A.R. 1995 Most similar neighbor: an improved sampling inference procedure for natural resource planning . Forest Science 41 , 337 – 359 . doi: 10.1093/forestscience/41.2.337 . Google Scholar OpenURL Placeholder Text WorldCat Crossref Narendra , P.M. and Goldberg , M. 1980 Image segmentation with directed trees . IEEE Trans. Pattern Anal. Mach. Intell. 2 , 185 – 191 . doi: 10.1109/tpami.1980.4766999 . Google Scholar Crossref Search ADS PubMed WorldCat National Land Survey of Finland . 2018 General map 1:1 M. https://tiedostopalvelu.maanmittauslaitos.fi/tp/kartta?lang=en. [Cited 4 May 2018] Natural Resources Institute Finland . 2013 Volume of pine, spruce and birch. http://kartta.luke.fi/opendata/valinta-en.html. [Cited 16 April 2018] Næsset , E. 2002 Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data . Remote Sens. Environ. 80 , 88 – 99 . doi: 10.1016/S0034-4257(01)00290-5 . Google Scholar Crossref Search ADS WorldCat Næsset E . 2014 Area-based inventory in Norway–From innovation to an operational reality. In: Maltamo , M., Næsset , E., and Vauhkonen , J. Forestry Applications of Airborne Laser Scanning - Concepts and Case Studies. Managing Forest Ecosystems . Springer, Springer Science + Business Media Dordrecht , pp. 215 – 240 . doi: 10.1007/978-94-017-8663-8 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Nilsson , M. , Nordkvist , K., Jonzéna , J., Lindgren , N., Axenstena , P., Wallerman , J., et al. 2017 A nationwide forest attribute map of Sweden predicted using airborne laser scanning data and field data from the National Forest Inventory . Remote Sens. Environ. 194 , 447 – 454 . doi: 10.1016/j.rse.2016.10.022 . Google Scholar Crossref Search ADS WorldCat Packalen , P. , Strunk , J.L., Pitkänen , J.A., Temesgen , H. and Maltamo , M. 2015 Edge-tree correction for predicting forest inventory attributes using area-based approach with airborne laser scanning . IEEE J. Sel. Top. Appl Earth Obs. Remote Sens. 8 , 1274 – 1280 . doi: 10.1109/JSTARS.2015.2402693 . Google Scholar Crossref Search ADS WorldCat Packalén , P. and Maltamo , M. 2006 Predicting the volume by tree species using airborne laser scanning and aerial photographs . Forest Science 52 , 611 – 622 . doi: 10.1093/forestscience/52.6.611 . Google Scholar OpenURL Placeholder Text WorldCat Crossref Packalén , P. and Maltamo , M. 2007 The k-MSN method for the prediction of species-specific stand attributes using airborne laser scanning and aerial photographs . Remote Sens. Environ. 109 , 328 – 341 . doi: 10.1016/j.rse.2007.01.005 . Google Scholar Crossref Search ADS WorldCat Packalén , P. , Suvanto , A. and Maltamo , M. 2009 A two stage method to estimate species-specific growing stock . Photogram. Eng. Remote Sens. 75 , 1451 – 1460 . doi: 10.14358/pers.75.12.1451 . Google Scholar Crossref Search ADS WorldCat Packalén , P. , Temesgen , H. and Maltamo , M. 2012 Variable selection strategies for nearest neighbor imputation methods used in remote sensing based forest inventory . Can. J. Remote Sens. 38 , 557 – 569 . doi: 10.5589/m12-046 . Google Scholar Crossref Search ADS WorldCat Pitkänen , J. , Maltamo , M., Hyyppä , J. and Yu , X. 2004 Adaptive methods for individual tree detection on airborne laser based canopy height model. International archives of photogrammetry . Remote Sens. Spatial Inform. Sci. 36 , 187 – 191 . Google Scholar OpenURL Placeholder Text WorldCat Silva , C. , Hudak , A.T., Vierling , L.A., Loudermilk , E.L., O’Brien , J.J., Hiers , J.K., et al. 2016 Imputation of individual longleaf pine (Pinus palustris mill.) tree attributes from field and LiDAR data . Can. J. Remote Sens. 42 , 554 – 573 . doi: 10.1080/07038992.2016.1196582 . Google Scholar Crossref Search ADS WorldCat Strunk , J. , Temesgen , H., Andersen , H.-E., Flewelling , J.P. and Madsen , L. 2012 Effects of LIDAR pulse density and sample size on a model-assisted approach to estimate forest inventory variables . Can. J. Remote Sens. 38 , 644 – 654 . doi: 10.5589/m12-052 . Google Scholar Crossref Search ADS WorldCat Tomppo , E. and Halme , M. 2004 Using coarse scale forest variables as ancillary information and weighting of variables in k-NN estimation: a genetic algorithm approach . Remote Sens. Environ. 92 , 1 – 20 . doi: 10.1016/j.rse.2004.04.003 . Google Scholar Crossref Search ADS WorldCat Valbuena, R., Mauro, F., Arjonilla, F.J., Manzanera, J.A. 2011 Comparing airborne laser scanning-imagery fusion methods based on geometric accuracy in forested areas . Remote Sens. Environ. 115 , 1942 – 1954 . doi: 10.1016/j.rse.2011.03.017 . Crossref Search ADS WorldCat Vauhkonen , J. 2010 Estimating single-tree attributes by airborne laser scanning: methods based on computational geometry of the 3-D point data . Dissertationes Forestales 104 , 44 . doi: 10.14214/df.104 . Google Scholar Crossref Search ADS WorldCat Vauhkonen , J. , Korpela , I., Maltamo , M. and Tokola , T. 2010 Imputation of single-tree attributes using airborne laser scanning-based height, intensity, and alpha shape metrics . Remote Sens. Environ. 114 , 1263 – 1276 . doi: 10.1016/j.rse.2010.01.016 . Google Scholar Crossref Search ADS WorldCat Vauhkonen , J. , Maltamo , M., McRoberts , R.E., and Næsset , E. 2014 Introduction to forest applications of airborne laser scanning. In: Maltamo , M., Næsset , E., and Vauhkonen , J. Forestry Applications of Airborne Laser Scanning - Concepts and Case Studies. Managing Forest Ecosystems . Springer, Springer Science + Business Media , pp. 1 – 16 . doi: 10.1007/978-94-017-8663-8 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC White , J.C. , Wulder , M.A., Varhola , A., Vastaranta , M., Coops , N.C., Cook , B.D., et al. 2013 A best practices guide for generating forest inventory attributes from airborne laser scanning data using an area-based approach. In Information Report FI-X-10. Natural Resources Canada, Canadian Forest Service . Canadian Wood Fibre Centre , p. 50 . doi: 10.5558/tfc2013-132 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Xu , Q. , Hou , X., Maltamo , M. and Tokola , T. 2014 Calibration of area based diameter distribution with individual tree based diameter estimates using airborne laser scanning . ISPRS J. Photogramm. Remote Sens. 93 , 65 – 75 . doi: 10.1016/j.isprsjprs.2014.03.005 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2021. Published by Oxford University Press on behalf of Institute of Chartered Foresters. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © The Author(s) 2021. Published by Oxford University Press on behalf of Institute of Chartered Foresters. TI - Prediction error aggregation behaviour for remote sensing augmented forest inventory approaches JF - Forestry DO - 10.1093/forestry/cpab007 DA - 2021-08-09 UR - https://www.deepdyve.com/lp/oxford-university-press/prediction-error-aggregation-behaviour-for-remote-sensing-augmented-h0Gdx0VZ58 SP - 576 EP - 587 VL - 94 IS - 4 DP - DeepDyve ER -