Airborne-laser-scanning-derived auxiliary information discriminating between broadleaf and conifer trees improves the accuracy of models for predicting timber volume in mixed and heterogeneously structured forests

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY-NC-ND 4.0 https://creativecommons.org/licenses/by-nc-nd/4.0/ Airborne-laser-scanning-derived auxiliary information discriminating between broadleaf and conifer trees improves the accuracy of models for predicting timber volume in mixed and heterogeneously structured forests © 2020 Elsevier B.V. All rights reserved. Accepted version (Final draft) Bont, Leo Gallus; Hill, Andreas; Waser, Lars T.; Bürgi, Anton; Ginzler, Christian; Blattert, Clemens Bont, L. G., Hill, A., Waser, L. T., Bürgi, A., Ginzler, C., & Blattert, C. (2020). Airborne-laser- scanning-derived auxiliary information discriminating between broadleaf and conifer trees improves the accuracy of models for predicting timber volume in mixed and heterogeneously structured forests. Forest Ecology and Management, 459, Article 117856. https://doi.org/10.1016/j.foreco.2019.117856 2020 Page 1 Airborne-laser-scanning-derived auxiliary information discriminating 1 between broadleaf and conifer trees improves the accuracy of models for 2 predicting timber volume in mixed and heterogeneously structured forests 3 4 Leo Gallus Bont a,*, Andreas Hill b, Lars T. Waser c, Anton Bürgi a, Christian Ginzler c, Clemens Blattert d 5 *Corresponding Author 6 a Sustainable Forestry Group, Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), 7 Zuercherstrasse 111, CH 8903 Birmensdorf, Switzerland, 8 b Forest Service Rhineland-Palatinate, Beethovenstrasse 3, 54516 Wittlich, Germany 9 c Department of Land Change Science, Swiss Federal Institute for Forest, Snow and Landscape 10 Research (WSL), CH 8903 Birmensdorf, Switzerland 11 d Department of Biological and Environmental Science, University of Jyväskylä, P.O. Box 35, 40014 12 Jyväskylä, Finland 13 14 Highlights: 15 • Information on forest type proportions improved the accuracy of volume predictions 16 • Forest type proportions weighted by canopy height outperformed area proportions 17 • The best model performance was obtained using forest type maps from leaf-off 18 LiDAR 19 20 Page 2 Abstract 21 Managing forests for ecosystem services and biodiversity requires accurate and spatially 22 explicit forest inventory data. A major objective of forest management inventories is to 23 estimate the standing timber volume for certain forest areas. In order to improve the 24 efficiency of an inventory, field based sample-plots can be statistically combined with 25 remote sensing data. Such models usually incorporate auxiliary variables derived from 26 canopy height models. The inclusion of forest type variables, which quantify broadleaf and 27 conifer volume proportions, has been shown to further improve model performance. 28 Currently, the most common way of quantifying broadleaf and conifer forest types is by 29 calculating the proportions of the corresponding areas of the canopy cover. This practice 30 works well for single-layer forests with only a few species, but we hypothesized that this is 31 not best practice for heterogeneously structured and mixed forests, where the area 32 proportion does not accurately reflect the timber volume proportion. To better represent 33 the broadleaf and conifer volume proportions, we introduced two new auxiliary variables in 34 which the area proportion is weighted by height information from a canopy height model. 35 The main objectives of this study were: (1) to demonstrate the advantage of including forest 36 type (broadleaf/conifer distinction) information in ordinary least squares regression models 37 for timber volume prediction using widely available data sources, and (2) to investigate the 38 hypothesis that including the broadleaf and conifer proportions, weighted by canopy height 39 information, as additional auxiliary variables is favourable over including simple area 40 proportions. The study was conducted in three areas in Switzerland, all of which have 41 heterogeneously structured and mixed forests. Our main findings were that the best model 42 performance can generally be achieved: (1) by deriving conifer and broadleaf proportions 43 from a high-resolution broadleaf/conifer map derived from leaf-off airborne laser scanning 44 Page 3 data, and (2) by using broadleaf/conifer proportions weighted by height information from a 45 canopy height model. Incorporating the so-derived conifer and broadleaf proportions 46 increased the model accuracy by up to 9 percentage points in root mean square error 47 (RMSE) compared with models not using any forest type information, and by up to 2 48 percentage points in RMSE compared with models using conifer and broadleaf proportions 49 based solely on the corresponding area proportions, as done in current practice. Our findings 50 are particularly relevant for mixed and heterogeneously structured forests, such as those 51 managed to achieve multiple functions or to adapt effectively to climate change. 52 53 Keywords: airborne laser scanning, best fit models, canopy height model, forest type map, 54 high-precision forest inventory, image-based point clouds, mixed and heterogeneously 55 structured forest, ordinary least squares regression models, merchantable timber volume 56 57 Page 4 1 Introduction 58 Forest ecosystems provide multiple benefits for humans and are particularly important for 59 the conservation of biodiversity (Bäck et al., 2017; MEA, 2005). These manifold contributions 60 of forests make their management a complex and challenging task requiring accurate and 61 spatially explicit information. The aim of forest inventories in general is to obtain reliable 62 information on the condition and development of the forest (Barrett et al., 2016). As the 63 census of an entire forest area is usually impossible, because of the high costs involved, 64 sampling concepts are used in practice. With these methods, the local hectare density of 65 timber volume, basal area and many other forest attributes are derived from measurements 66 of the trees in randomly or systematically distributed sample-plots in the forest area. This 67 data is then used to estimate mean values and totals for the entire forest area, for example 68 the mean or the total timber stock. 69 More accurate information on an entire forested area or a small area (e.g. parts of a forest 70 enterprise) can be obtained by complementing field based inventories with remote sensing 71 data. Furthermore, this can be a cost-effective alternative to increasing the number of field 72 based sample units. The principle of such two-phase inventories is to use statistical models 73 to predict response variables, such as basal area or timber volume, for the population units 74 where no field data is available. Many studies have already demonstrated the potential of 75 these methods (e.g. Hill et al., 2018; Magnussen et al., 2014; Mandallaz et al., 2013; Næsset, 76 2004, 2002; Steinmann et al., 2013). Usually, such statistical models are based on auxiliary 77 variables derived from a canopy height model (CHM) (Xu et al., 2019). Other sources of 78 information, such as tree species or forest type (we use the term ‘forest type’ for the 79 distinction between either conifer and broadleaf trees or evergreen and deciduous trees, 80 Page 5 see section 2.2) have been added in a few case studies to improve the performance of the 81 models (Gabriel et al., 2018; Hill et al., 2018). 82 For Switzerland, a freely available national forest type map (FTM) based on optical remote 83 sensing data (Waser et al., 2017) exists in the framework of the Swiss National Forest 84 Inventory (NFI). FTMs can also be derived using multi-temporal (i.e. leaf-on, leaf-off) laser 85 scanning data. Such an approach entails collecting two Light Detecting and Ranging (LiDAR) 86 datasets at the same site under leaf-off and leaf-on conditions (Liang et al., 2007). 87 Alternatively, forest types can be derived based on return intensity and ranking distributions 88 of laser scans from a single point in time, preferably under leaf-off conditions for the 89 conifer/broadleaf differentiation (Liang et al., 2007; Ørka et al., 2009; Parkan, 2018). 90 However, deciduous conifers (i.e. conifers which defoliate in autumn, such as larch) cannot 91 be identified as conifers with any of these three approaches (Fassnacht et al., 2016). 92 FTM information has been used in many studies to improve standing timber volume 93 estimations. For example, Breidenbach et al. (2008) used a continuous variable ‘conifer 94 proportion’ and its interaction term with the average canopy height to include FTM 95 information derived from leaf-off LiDAR data. Latifi et al. (2012) formed a FTM based on 96 colour infrared (CIR) orthoimages and included the forest type as a categorical variable in 97 their model. The forest type of the sample-plots was assigned to either broadleaf or conifer 98 if the proportion of the pixels of one particular type exceeded 70%, and it was assigned to 99 the mixed category in all other cases. Straub et al. (2009) used CIR orthoimages to create a 100 FTM. They derived the conifer and broadleaf proportions from the percentages of the 101 corresponding pixels in each sample-plot, and this information was then included in the 102 model to estimate the stem volume of forest stands. Hill et al. (2018) included information 103 on five tree species as a categorical variable to improve timber volume predictions. Finally, in 104 Page 6 Nordic countries such as in Finland, including tree species information as auxiliary data to 105 predict timber volume is quite common, as described by Kukkonen et al.( 2019, 2018), 106 Packalén and Maltamo (2006) and Räty et al. (2016). 107 All of the approaches mentioned above have led to a significant improvement in the 108 accuracy of standing timber volume predictions by including forest type or tree species 109 information in the models. However, in all these studies forest type proportion information 110 was derived based on the area covered by the canopy of the corresponding forest type. This 111 current practice works well for even-aged and single-layer forests. However, we hypothesize 112 that this approach is not best practice for heterogeneously structured and mixed forests, in 113 which different age classes and tree species can occur across a small surface (e.g. within a 114 sample-plot). In such forests the area proportion does not adequately reflect the timber 115 volume proportion. There are two possible reasons for this discrepancy. First, the mean tree 116 size might differ depending on the canopy height. To illustrate this possibility we consider a 117 sample-plot on which conifer and broadleaf trees cover about the same amount of area and 118 where the conifers are all mature whereas the broadleaf trees are much younger. In this 119 case the volume on the half of the sample-plot with the young broadleaf trees is clearly 120 smaller than that on the half with the mature coniferous trees. This point has been 121 confirmed in studies about allometric relationships, such as the work by Reineke, (1933) in 122 establishing the self-thinning rule. In this case, using the area proportion of conifers leads to 123 an underestimation of the proportion of conifer timber volume. Second, stand density (the 124 number of stems in a certain area) also differs for the different species. This point was 125 shown by Pretzsch and Biber (2005) and by Rivoire and Le Moguedec (2012), who 126 generalized the self-thinning relationship of Reineke (1933) for multi-species and mixed-size 127 forests. 128 Page 7 To better represent the volume proportion of different forest types, in this study we 129 introduced two new forest type variables (FTVs), referred to as ‘weighted-canopy-height 130 proportions’, in which the area proportions are weighted by height information from a 131 canopy height model. We hypothesized that the ‘weighted-canopy-height-proportion’ FTVs 132 are favourable over simple area proportion FTVs in mixed and heterogeneously structured 133 forests. 134 The overall objectives of the present study were: (1) to determine the advantage of including 135 forest type information in regression models for timber volume prediction using existing and 136 forthcoming data sources, and (2) to investigate the hypothesis that the ‘weighted-canopy-137 height proportion’ FTVs are favourable over simple area proportion FTVs. Both aims are of 138 high practical interest because the underlying data sources used to derive the FTVs, such as 139 the national broadleaf/conifer map or the national leaf-off airborne laser scanning (ALS) 140 data, are widely available and can be integrated with a small effort into current inventories. 141 However, neither of these sources have been used operationally for inventory purposes. This 142 work was embedded in the implementation of design-based regression estimators 143 (Mandallaz, 2013) for predicting standing timber volume. As these design-based regression 144 estimators rely on ordinary least squares (OLS) regression, we used OLS regression models in 145 our study. 146 We addressed the following specific research questions: 147 (1) Is there a gain in model performance when the new ‘weighted-canopy-height proportion’ 148 FTVs are incorporated into OLS regression models for predicting timber volume, compared 149 with models including the simple area proportion FTVs and models with no forest type 150 explanatory variables? 151 Page 8 (2) Is there a common best practice for integrating FTM information, such as the spatial 152 resolution considered, that is independent of the individual LiDAR and FTM characteristics? 153 2 Materials and Methods 154 2.1 Study areas 155 We studied the effect of incorporating different FTV alternatives in three independent study 156 areas in Switzerland: Bremgarten (Brg), Glâne-Farzin (GF) and Sarine (SA) (Figure 1). All three 157 study areas are heterogeneously structured and mixed temperate forests of conifer and 158 broadleaf species, but they have different properties in terms of age-class structure. 159 2.1.1 Bremgarten (Brg) 160 In Brg, the latest forest inventory data from 2011/2012 was used as input. This data is based 161 on permanent sample-plots on an 80 m x 150 m grid, as described by Schmid-Haas (2003). 162 The location of the sample-plots is shown in Figure 1. A total of 363 sample-plots were 163 measured, with an average timber volume of 274 m3 ha-1, 55% of which is conifer wood 164 (Table 1). The distribution of the timber volumes of the sample-plots indicates that stands 165 with large timber volumes are under-represented compared with a normal age-class 166 structure (Salo and Tahvonen, 2002) (Figure 2). This unbalanced age-class structure of the 167 forest was caused by the storm Lothar in December 1999, which primarily damaged mature 168 stands. Figure 3 shows that sample-plots with a large timber volume are dominated by the 169 conifer volume, whereas sample-plots with a small volume are dominated by the broadleaf 170 volume. LiDAR data was acquired in November 2011, providing excellent temporal 171 synchronization between the field based and the remote sensing data. The density of the 172 LiDAR raw data is at least 8 points m-2 and is comparable with values in the other study areas 173 (Table 1). 174 Page 15 238 Figure 4: (a) – (e) Forest type maps (FTMs) and (f) canopy height model (CHM) for a sub-area around eight 239 field based inventory sample-plots in the study area Bremgarten (coordinate system: EPSG 21781, CH 1903 / 240 LV 03) (Li = LiDAR, Sen = Sentinel NFI, NFI = Swiss NFI Orthoimage, OI = Orthoimage, DLT = Copernicus 241 Dominant Leaf Type). For (a) and (e) the FTM differentiates between evergreen and deciduous foliage; for 242 (b), (c) and (d) the FTM differentiates between conifer and broadleaf trees. 243 Page 16 2.2.2 Swiss NFI Orthoimage (NFI) 244 This map distinguishes between broadleaf and conifer trees (Waser et al., 2017) and was 245 produced for the Swiss National Forest Inventory (NFI). It covers all of Switzerland and has a 246 spatial resolution of 3 m x 3 m. The methodology is based on the classification of more than 247 1700 four-band aerial photo strips. Prior to applying it to all of Switzerland, different 248 classification methods, i.e. Support Vector Machine, Logistic Regression and Random Forest 249 (RF) were tested for selected areas in Switzerland regarding accuracy, computing time and 250 minimum number of required training values. The tests revealed that RF performed best 251 regarding accuracy and processing time (Waser et al., 2017; Figure 4b). 252 2.2.3 Sentinel NFI (Sen) 253 The FTM sentinel NFI is a national FTM that distinguishes between broadleaf and conifer 254 trees at a spatial resolution of 10 m. 255 This map is based on freely available spaceborne Sentinel-1 / -2 data from the European 256 Space Agency’s Copernicus Programme (ESA, 2019). In order to get a cloud-free coverage of 257 all of Switzerland, a total of 50 Sentinel-2 images from June, July and August of 2016–2018 258 were used. The Sentinel-1 SAR data was acquired in summer 2016 and 2017 and processed, 259 i.e. flattening terrain and increasing the spatial resolution, according to Small (2012). This 260 map is a product of the Swiss NFI, and it is a follow-up and improved version of the 261 preliminary dataset (referred to as FTM NFI) and is free from overestimations of conifers 262 caused by topographic and illumination effects. This new approach incorporates a RF 263 classifier. According to Breiman (2001), this is a widely used ensemble classifier that 264 produces multiple decision trees using a randomly selected subset of training samples and 265 variables (Figure 4c). 266 2.2.4 Copernicus Dominant Leaf Type (DLT) 267 Page 17 The FTM Dominant Leaf Type (DLT) is a product of the Copernicus Land Monitoring Service 268 coordinated by the European Environment Agency. The DLT FTM provides information on 269 the dominant leaf type (broadleaf or conifer) at a 20 m x 20 m resolution, and it was derived 270 from multi-temporal satellite image data using Support Vector Machine (Langanke, 2017; 271 Figure 4d). 272 2.2.5 LiDAR (Li) 273 The FTM LiDAR differentiates between deciduous and evergreen foliage and was computed 274 by using the corrected return intensity of leaf-off airborne laser scanning (ALS) data, as 275 proposed in the Digital Forestry Toolbox (Parkan, 2018). The FTM LiDAR has a spatial 276 resolution of 0.5 m and was computed for GF and SA from the same leaf-off data as used for 277 the CHM and for Brg from a leaf-off flight from March 2014 with a density of > 10 points m-2. 278 This map has great potential for application because leaf-off ALS data is widely available and 279 calculation with the Digital Forestry Toolbox is straightforward (Figure 4e). 280 2.3 Canopy height model (CHM) 281 We used the LiDAR raw point data (Table 1) to compute pit-free CHMs. This approach for 282 CHM calculation avoids gaps (‘pits’), following the methods presented by Isenburg (2014) 283 and Khosravipour et al. (2014). This algorithm was implemented in LAStools (Rapidlasso 284 GmbH, Gilching, Germany), but also using the package ‘lidR’ (Roussel and Auty 2017) in R 285 version 3.5 (R Core Team, 2018). An example of a pit-free CHM is given in Figure 4. The 286 CHMs produced in this study have a spatial resolution of 0.33 m (Figure 4f). 287 2.4 Response variable 288 The response variable is the local density of the merchantable timber volume [m3 ha-1] 289 (volume of the stem and branches with a diameter ≥ 7 cm) and is referred to as VOL in the 290 model formulation. The field measurements were carried out according to the sampling 291 Page 18 protocol of Keller (2013) for GF and SA and of Schmid-Haas et al. (1993) for Brg. In all study 292 sites, only the DBH (diameter at breast height, measured at 130 cm above the level of the 293 terrain) was recorded. For the calculation of the timber volume of a single tree (VOL_ST) in 294 Brg, the following one-parameter volume function of Hoffmann (1982) was used: 295 𝑉𝑉𝑉𝑉𝑉𝑉_𝑆𝑆𝑆𝑆(𝐷𝐷𝐷𝐷𝐷𝐷)=𝑒𝑒𝑒𝑒𝑒𝑒�𝑎𝑎1+𝑎𝑎2𝑙𝑙𝑙𝑙(𝐷𝐷𝐷𝐷𝐷𝐷)+𝑎𝑎3�𝑙𝑙𝑙𝑙(𝐷𝐷𝐷𝐷𝐷𝐷)�4�, with the parameters a1, a2 and 296 a3 and with DBH as the only input variable. In GF and SA a volume function based on volume 297 tables was used (Schweizer, 2012), again with DBH as the only input variable. This volume 298 function was chosen to ensure comparability with earlier inventories. While in SA and GF the 299 same volume function was used for trees of both leaf types, in Brg separate volume 300 functions were used for broadleaf and coniferous trees. However, a comparison of the 301 functions showed that the differences in the predicted tree volumes between the volume 302 functions differentiating by forest type (Brg) and the mixed volume functions (GF and SA) 303 were marginal (Figure 5). Clear differences only emerged at ≥ 100 cm DBH. However, this has 304 no relevance for the present study, as only four sample trees were above this DBH threshold. 305 As volume functions are not valid for large areas, they have to be calibrated locally. The 306 parametrization of the volume function in Brg was done during the inventory of 1986. For 307 this purpose, the height (h), DBH and diameter at 7 m height (d7) was measured for all 308 sample trees in one-eighth of each sample-plot area (for trees with a DBH < 20 cm), in one-309 quarter of each sample-plot area (20 cm ≤ DBH < 50 cm) or in the entire sample-plot (DBH ≥ 310 50 cm) (Schmid-Haas et al., 1993). With the help of three-parameter volume functions 311 requiring h, DBH and d7 as input, a reference volume was then calculated for all sample trees 312 that were measured in detail. Finally, a nonlinear curve fitting method was used to 313 determine the locally adapted parameters a1, a2 and a3 (Hoffmann, 1982). 314 Page 19 To ensure comparability with previous inventories, we therefore used the parameters 315 derived in 1986. A validation of the volume functions used in Brg was done by Kaufmann 316 (2001). He reported a coefficient of determination (R2) of 95% and a standard deviation of 317 the residuals (Rs) of 20% (conifer) and 27% (broadleaf) of the mean. The local density of the 318 total timber volume on each terrestrial sample-plot was based on the timber volume and 319 inclusion probabilities of individual trees and was calculated using the Horvitz Thompson 320 Estimator (Mandallaz, 2007). 321 322 Figure 5: Volume functions for Glâne-Farzin and Sarine in canton Freiburg (Frbg) and for Bremgarten (Brg). 323 2.5 Auxiliary variables 324 2.5.1 Auxiliary variables based on canopy height models 325 Auxiliary variables form the basis for model building and were used to derive a model for the 326 response variable (VOL) from the sample data. Auxiliary variables can come from various 327 sources, such as LiDAR raw data, CHMs and FTMs. An overview of auxiliary variables that are 328 potentially interesting for forestry can be found in McCallum et al. (2014) and Saarela et al. 329 (2015). Auxiliaries derived from the CHM, which were used in our models, are listed in Table 330 3. 331 Page 20 Most of the auxiliary variables included in this study are well established for forest inventory 332 purposes. However, to the best of our knowledge, the variable SSq has not previously been 333 described in the literature. It is presented in detail in Figure 6a and represents the sum of 334 the squares of the height values of the individual CHM-pixels within a sample-plot. The idea 335 behind introducing this auxiliary variable was to make it possible to map the allometric 336 relationships more appropriately, as the stand volume is usually not a linear function of the 337 canopy height (Eichhorn, 2013; Pretzsch, 2001). 338 Description Abbreviation Unit Percentile – values Specifies the height percentiles of the pixels in the sample-plot. Example: 95th percentile: P95 PXX XX = {99, 95, 90, 80, 75, 70, 60, … , 30, 25, 20, 10, 05, 01} [m] (same unit as the CHM) Coverage – values Indicates the vegetation cover at a certain height. Example: 2 m coverage: C02 = 0.9 means that 90% of the pixels of the vegetation height model in the plot are higher than or equal to 2 m above the ground. CXX XX = {00, 02, 05, 10, 15, 20, … , 50} [0 – 1] Minimum Minimum pixel value within a plot MIN [m] Maximum Maximum pixel value within a plot MAX [m] Standard deviation Standard deviation of the individual pixels within a plot STD [m] Average Average value (height) of the individual pixels within a plot MEAN [m] Sum of squares Sum of the squares of the individual pixels within a plot SSq Average sum of squares Sum of squares divided by number of pixels in a plot MSSq Coefficient of variation MEAN / STD CV Table 3: Auxiliary variables derived from the canopy height model (CHM). 339 Page 21 340 Figure 6: Principle of the combined evaluation of different grids. The two grids forest type map (FTM) 341 (conifer/broadleaf) and canopy height model (CHM) are superimposed and evaluated for each attribute of 342 the FTM. a) Illustration of the derivation of the auxiliary variable SSq (sum of squares), based on a combined 343 evaluation of the CHM and the FTM. The circle represents the boundary of the sample-plot. b) Computation 344 of the different forest type variable alternatives: 1) ratio of the number of pixels (area_prop), 2) ratio of the 345 sums of canopy height (area_chm_prop), and 3) ratio of the sums of squares of the canopy heights 346 (area_chm_sq_prop). Most pixels in (a) are empty to simplify the comprehensibility of the calculations. If the 347 FTM and the CHM have different spatial resolutions, the analysis is performed on the basis of the grid with a 348 higher spatial resolution. 349 2.5.2 Forest type variables (FTV) 350 In order to include FTM information in the statistical models, auxiliaries (referred to as forest 351 type variables, FTV) were introduced that describe the proportion of conifer trees. There 352 were five FTV alternatives: four continuous and one categorical. The continuous FTV 353 alternatives were: (1) the ratio of the number of pixels of different forest types (area_prop), 354 (2) the ratio of average heights of the forest types (area_chm_prop), (3) the ratio of the sum 355 of squares of the canopy heights of the forest types (area_chm_sq_prop) (see Figure 6b), 356 Page 22 and (4) the variable ground_proportion, which is introduced in section 2.5.3. area_prop is 357 the classical formulation and represents the proportion of the area that is covered by the 358 canopies of conifers. area_chm_prop and area_chm_sq_prop were our new ‘weighted-359 canopy-height proportion’ FTVs, for which the area proportion was weighted by information 360 from a CHM. Therefore, the two grids (CHM and FTM) were superimposed and a separate 361 evaluation was performed for each forest type, as presented in Figure 6a and b. For the 362 categorical variable, we chose a setup that assigns the sample-plot to a main forest type 363 (broadleaf or conifer) if the proportion of the pixels of one particular type exceeds 70% and 364 to the category ‘mixed’ otherwise, as done by Breidenbach et al. (2008) and Latifi et al. 365 (2012). Other thresholds were not explored in this study, even though they may also have an 366 impact on model performance (Hill et al. 2018). However, this study was not focused on the 367 investigation of categorical variables, as categorical variables have the disadvantage of 368 requiring an additional degree of freedom for each level. This effect is amplified when 369 interaction terms between the FTV and other auxiliary variables are included in the model. 370 This is particularly problematic when the number of sample-plots is small, as in this study. A 371 rule of thumb says that 10–20 observations per auxiliary variable are required to be able to 372 detect reasonable-sized effects with reasonable power (Harrell, 2017). 373 2.5.3 Forest type variable ground_proportion and classification accuracy 374 FTVs derived from FTMs (section 2.5.2) are usually not error free, as shown, for example, by 375 Hill et al. (2018) and Straub et al. (2009). We evaluated classification accuracy to understand 376 the precision and the errors of the different FTVs and FTMs. This was done by evaluating the 377 correlation of the continuous FTVs area_prop, area_chm_prop and area_chm_sq_prop with 378 ground_proportion (Figure 7). ground_proportion represented the conifer timber volume 379 proportion for each sample-plot based on field measurements. This variable was derived 380 Page 23 from the volumes of the individual trees, which were predicted by means of the volume 381 functions with DBH as an input variable (see section 2.4). ground_proportion was not an 382 error-free variable, since the underlying volume functions were not error free, but field 383 measurements are usually more precise than variables derived from remote sensing data. 384 The analyses showed that the LiDAR-derived FTVs had the strongest correlation with 385 ground_proportion for all study areas. The FTVs area_chm_sq_prop and area_chm_prop 386 were more reliable than area_prop, except for the FTMs Sen and DLT, which had low 387 resolutions of 10 m and 20 m, respectively. The classification accuracy generally increased 388 with average canopy height in the sample-plot. Therefore, including an interaction term 389 between the FTVs and MEAN had the potential to further improve model accuracy, and it 390 was considered in some model formulations. To explore the potential benefit of having a 391 highly precise FTM (i.e. without classification error) based on field measurements, we 392 considered ground_proportion as a FTV in our statistical models. Although classification 393 errors of forest type are typically small when they are based on accurate field observations, 394 the ground_proportion FTV used in our study was not error free because its derivation was 395 based on volume estimations from volume functions. 396 Page 24 397 Figure 7: Correlation of the proportion of conifer timber volume, derived using the different FTMs and FTV 398 alternatives compared with the variable ground_proportion. The figure is separated into different average 399 classes of the CHM that were considered. min_MEAN = 20 means that only sample-plots with a minimum 400 average canopy height of 20 m were considered; a value of 0 means that all sample-plots were considered 401 for the analysis of the correlation. Brg = Bremgarten, GF = Glâne-Farzin (GF) and SA = Sarine. 402 2.6 Statistical methods 403 For model building, we excluded sample-plots with a timber volume value of zero and 404 sample-plots in which harvests had taken place in the period between the field 405 measurements and the acquisition of the remote sensing data. This step is referred to as 406 ‘cleaning’ in Table 1. The FTMs derived from orthoimage (OI) and LiDAR (Li) differentiated 407 between deciduous and evergreen foliage. To differentiate accurately between conifer and 408 broadleaf trees for these FTMs, we removed sample-plots that include larch (Larix decidua). 409 We labelled these cases with the suffix ‘wo_La’ (without larch). The number of remaining 410 Page 31 FTM Li_wo_La and without FTV. The adjusted r-squared values (Figure A2 in the Appendix) 509 indicate an excellent model performance (high adjusted r-squared value), whereas the high 510 RMSECV value in Figure 10 indicates the opposite. However, this was the only case of over-511 fitting. A further comparison of adjusted r-squared values (Figure A2) demonstrates that the 512 remaining models with the highest adjusted r-squared values (Model 9 and GF_Li_wo_La) 513 had rather low RMSE values. 514 Table 6 shows the model formulations for the best performing models. Although some 515 models had many predictors, they were not over-fitted, as shown through the leave-one-out 516 cross-validated RMSE (Figures 8 and 10). For Brg and SA, Model 10 provided the second best 517 model performance according to RMSE, which was only slightly worse than that of Model 8 518 (Table 7) but included far fewer variables. 519 Page 32 520 Figure 9: Change in RMSECV [%] values under the different model types and forest type variables (FTV) 521 alternatives compared with a model without FTVs in the three case study areas (Brg = Bremgarten, GF = 522 Glâne-Farzin, SA = Sarine). Results are only displayed for variables derived from the LiDAR FTM. The peak in 523 GF and Li_wo_La for Models 7–9 was caused by over-fitted models without FTVs. 524 Page 33 525 Figure 10: Leave-one-out cross-validated RMSECV [%] for each regression model, forest type map and 526 combination of forest type variable alternatives in the case study areas Glâne-Farzin (GF) and Sarine (SA). 527 Study area Model Forest type variable Auxiliary variables used Brg_Li 8 area_chm_sq_prop P40 + P30 + P25 + P01 + C10 + C20 + foresttype + MEAN + foresttype:MEAN Brg_Li 8 area_chm_prop P40 + P30 + P25 + P01 + C10 + C20 + foresttype + MEAN + foresttype:MEAN Brg_Li_wo_La 8 area_chm_sq_prop P99 + P60 + P01 + C10 + C20 + CV + foresttype + MEAN + foresttype:MEAN Brg_Li_wo_La 8 area_chm_prop P99 + P60 + P50 + P30 + P01 + C10 + C20 + CV + foresttype + MEAN + foresttype:MEAN GF_Li 8 categorical P60 + P25 + C02 + C10 + C15 + C25 + C30 + foresttype + MEAN GF_Li_wo_La 8 categorical C02 + C10 + C15 + C25 + C30 + foresttype + MEAN + foresttype:MEAN SA_Li 8 area_chm_sq_prop P20 + C05 + C15 + C20 + C25 + C35 + C50 + foresttype + MEAN + foresttype:MEAN Page 34 SA_Li 8 area_chm_prop P20 + C05 + C15 + C20 + C25 + C35 + C50 + foresttype + MEAN + foresttype:MEAN SA_Li_wo_La 8 area_chm_sq_prop MSSq + P25 + P05 + C02 + C15 + C20 + C25 + C35 + C50 + chm_sq_average + foresttype + MEAN + foresttype:MEAN SA_Li_wo_La 8 area_chm_prop MSSq + P25 + P05 + C02 + C15 + C20 + C25 + C35 + C50 + chm_sq_average + foresttype + MEAN + foresttype:MEAN Table 6: Selected model formulations for the best performing forest type model combinations for Model 8. 528 Interaction terms are indicated by ‘:’. ‘foresttype’ stands for the forest type variable and is either area_prop, 529 area_chm_prop, area_chm_prop, categorical or ground_proportion. 530 531 Study area Model Forest type variable Auxiliary variables used Brg_Li 10 area_chm_sq_prop MEAN + STD + foresttype + MEAN:foresttype Brg_Li 10 area_chm_prop MEAN + STD + foresttype + MEAN:foresttype Brg_Li_wo_La 10 area_chm_sq_prop MEAN + STD + foresttype + MEAN:foresttype Brg_Li_wo_La 10 area_chm_prop MEAN + STD + foresttype + MEAN:foresttype GF_Li 10 categorical MEAN + STD + foresttype + MEAN:STD GF_Li_wo_La 10 categorical MEAN + foresttype SA_Li 10 area_chm_sq_prop MEAN + foresttype + MEAN:foresttype SA_Li 10 area_chm_prop MEAN + foresttype + MEAN:foresttype SA_Li_wo_La 10 area_chm_sq_prop MEAN + foresttype + MEAN:foresttype SA_Li_wo_La 10 area_chm_prop MEAN + foresttype + MEAN:foresttype Table 7: Selected model formulations for the best performing forest type model combinations for Model 10. 532 Interaction terms are indicated by ‘:’. ’foresttype’ stands for the forest type variable and is either area_prop, 533 area_chm_prop, area_chm_prop, categorical or ground_proportion. 534 4 Discussion 535 4.1 Influence of different forest type variables and forest type maps on model 536 performance 537 Our results confirm the findings of previous studies (Breidenbach et al., 2008; Latifi et al., 538 2012) that incorporated FTM information to predict timber volume. Both earlier studies 539 reported an increase in model accuracy (measured as the LOOCV RMSE) of the OLS 540 regression models for their investigated study areas. In doing so, Latifi et al. (2012) found an 541 improvement of 2–4 percentage points in RMSE when stratifying into broadleaf and conifer 542 Page 35 forest types. We found similar improvements when the categorical FTV was included in 543 models (Figure 9). In our study, the best performing FTM in all study areas was the one 544 derived from leaf-off LiDAR data. This is not particularly surprising, as this map had both the 545 highest resolution and the best correlation with the values derived from field measurements 546 (ground_proportion), as visualized in Figure 7. The further order of the best suited FTMs was 547 less clear and varied between the study areas. In general, a map with a higher resolution 548 gave better results than one with a coarser resolution, even if both maps showed the same 549 classification accuracy, as derived in section 2.5.3. This was observed in Brg (Figure 8), where 550 we had the following order of accuracy according to the RMSE: Brg_OI (2 m x 2 m) ≈ Brg_NFI 551 (3 m x 3 m) > Brg_Sen (10 m x 10 m) > Brg_DLT (20 m x 20 m). The NFI Orthoimage map 552 (*_NFI) is biased due to topographic and illumination effects, i.e. shadows, of optical remote 553 sensing data (Waser et al., 2017). 554 We found that the selection of a FTV alternative strongly affected model accuracy, 555 particularly when FTMs with a high resolution (≤ 2 m) were used, such as the LiDAR-derived 556 map (0.5 m). The continuous variables thus outperformed the categorical variable. Among 557 the continuous variables, area_chm_prop and area_chm_sq_prop performed about equal 558 and were associated with higher model accuracies than area_prop. For area_chm_prop and 559 area_chm_sq_prop, the forest type proportion was weighted by canopy height information 560 of a CHM. This representation mapped the volume proportions of conifers and broadleaf 561 trees better, in particular for mixed and heterogeneously structured forests. The difference 562 in RMSE between all the different FTV alternatives was up to 2 percentage points. However, 563 in GF there was no clear difference between the different FTV alternatives, possibly because 564 exactly georeferenced sample-plot positions were not available and the model accuracy was 565 already high without incorporating FTVs (Figure 9). The gain in model accuracy was up to 9 566 Page 36 RMSE percentage points for Brg, and up to 3 percentage points for GF and SA. The 567 improvement in model accuracy, measured as percentage points of the RMSE, was lower in 568 GF and SA compared with Brg, as the results from GF and SA already showed a small RMSE 569 before any forest type information was added to the OLS regression model. 570 4.2 Best performing model formulation 571 Our second research question was how to best include FTM information in a regression 572 model. Model 8 performed best for all study areas and FTV combinations. This model was 573 calibrated using a variable selection procedure based on AIC values, followed by the addition 574 of interaction terms, and finally a second variable selection based on AIC values. Except for 575 GF_Li, all of the best performing models included an interaction term between foresttype 576 and MEAN (Table 6). Therefore, we highly recommend using interaction terms when 577 including FTM information in OLS regression models. Interaction terms are particularly 578 meaningful because, depending on the average vegetation height, the FTMs have a different 579 classification accuracy (Figure 7). 580 Although Model 8 included several variables, over-fitting was not an issue, as shown by the 581 cross-validation results. However, if we were interested in making inferences about the 582 predictors, Model 8 could have failed. In such cases, Model 10 would have been the 583 preferred model formulation because it included remarkably fewer variables (Table 7). Due 584 to its smaller number of variables, Model 10 was also less vulnerable to over-fitting. Further, 585 variance inflation, which occurs as a consequence of including too many correlated variables 586 and is quantified as the variance inflation factor (VIF), could affect the significance of single 587 variables (Fox and Monette, 1992). This was not the case in our study, however, because 588 Model 10 had the best AIC and highly correlated variables were eliminated. The only issue 589 could be that main effects are difficult to interpret in models with interactions. 590 Page 37 4.3 Potential of a high-precision forest type variable 591 The idea behind introducing ground_proportion (conifer volume proportion based on field 592 measurements) as a FTV in the models was to explore the potential benefit of having a high-593 precision FTM without misclassifications. Results from the two study areas Brg and SA show 594 that when the LiDAR-derived FTM is used in combination with the variables area_chm_prop 595 and area_chm_sq_prop, regression model performance can be even better than when this 596 FTM is used in combination with ground_proportion. This result is remarkable and deserves 597 an in-depth discussion. The following issues are relevant for its understanding and 598 interpretation. 599 First of all, it must be emphasized that both the variable ground_proportion and the remote-600 sensing-based FTV contained errors. In the case of ground_proportion, the errors were 601 caused by the volume functions used, which were solely based on the DBH. For example, 602 Kaufmann (2001) explored volume functions for the Swiss NFI. He stated that a tree in a later 603 stage of development has a larger stem volume than a tree with the same DBH in an earlier 604 developmental stage. This behaviour was ignored in the volume functions, but we think that 605 this relationship was mapped in the LiDAR-based ‘weighted-canopy-height proportion’ FTV, 606 as the canopy height is correlated with the developmental stage of a tree. In the case of the 607 remote-sensing-based FTV, the most relevant sources of error were misclassification 608 (confusion of conifer and broadleaf), inadequate temporal synchronization (a time 609 difference between the field measurements and the recording of remote sensing data), and 610 imperfect spatial matching of the sample-plot location with the detail of the remote sensing 611 data. 612 Second, the variables ground_proportion and area_chm_prop / area_chm_sq_prop can only 613 be compared directly if they are derived identically. This would be the case if 614 Page 38 ground_proportion was also derived by weighting the broadleaf/conifer sample tree 615 information by the respective tree heights – analogous to calculating area_chm_prop / 616 area_chm_sq_prop. However, this was not the case because tree heights were not recorded. 617 Third, the response variable, the volume per sample-plot, which was calculated by the same 618 volume functions as used to derive the forest type proportion (ground_proportion), was also 619 not error free and contains the same errors as the field observed FTV. The applied volume 620 functions ignore some real-existing influences on the tree volume and they therefore return 621 the same volume for different values of some influential factors. The evaluation and 622 comparison of different models and auxiliary variables is therefore always limited by the 623 accuracy of the reference values, which were in our case generated by one-parameter 624 volume functions. In order to allow a more precise interpretation of the results, one would 625 have to use more precise, e.g. three-parameter volume functions. Therefore, it is difficult to 626 say whether there was actually a statistically significant difference between the results 627 obtained by ground based forest type information and remote sensing based forest type 628 information. 629 Overall, the following conclusions can be drawn: (1) Although the FTM-based variables are 630 error prone, they can generate a more powerful signal for modelling than the field based 631 information ground_proportion. This supports the statement that the proposed weighting by 632 canopy height information of the FTV generates a very useful signal. 633 (2) Until now, it has always been assumed that field information is the reference in terms of 634 accuracy. However, this concept is limited by the sampling protocol; in this case, no tree 635 heights were measured during the field observation and therefore the corresponding 636 calculation of ground reference including tree height information was not possible. If tree 637 Page 39 height had been recorded, the ground reference would possibly have been as good or better 638 than the FTM-based variables. 639 Nevertheless, we recommend including ground_proportion as a FTV in regression models, 640 particularly during model evaluation, because the variable might serve as a benchmark to 641 assess the effect of including FTM information. Additionally, it might function as an indicator 642 of whether it is worth investing resources into creating an advanced FTM. 643 4.4 Generalization of the Results 644 In order to address our second research question, we studied the effect of incorporating 645 different FTV alternatives in three independent study areas. We considered the results in 646 Brg, with its precisely measured forest inventory sample-plots, as the reference. The main 647 findings were that the best results were generally achieved: (1) by deriving forest type 648 variables from a FTM with a high resolution and (2) by concurrently using forest type 649 variables, such as area_chm_prop and area_chm_sq_prop, derived by a superimposed 650 evaluation of the FTM and the CHM. These findings apply to Brg and SA. For GF, point (2) 651 could not be confirmed, on the one hand because the terrestrial inventory sample-plot 652 centres were not positioned accurately and on the other hand because GF is less 653 heterogeneously structured than the other forests because it is dominated by numerous 654 pure mature conifer sample-plots (Figure 3). To summarize, precisely measured inventory 655 sample-plots and the presence of mixed and heterogeneously structured forests are 656 required for this approach to be beneficial for current inventories. Further, leaf-off LiDAR 657 data is desirable. We worked with a LiDAR data density of 8 and 10 points m-2, which 658 delivered excellent results. Such data is widely available; for example, by 2023 Switzerland 659 will be covered completely with leaf-off LiDAR data with an average density of 15–20 points 660 m-2. 661 Page 40 4.5 Relevance for management 662 Climate change is expected to considerably influence forest ecosystems and their 663 management (Hanewinkel et al., 2013; Schelhaas et al., 2015). In combination with climate 664 change, disturbances are expected to increase and affect forest ecosystems (Seidl et al., 665 2017, 2014). To face these future challenges, Messier et al. (2019) recommended using 666 forest management to increase the number of tree species and the structural diversity of 667 forests at the landscape scale to improve the resilience of forests. Managing forests for 668 higher resilience requires, however, more accurate and more spatially explicit information 669 on the current situation of forests, which is usually assessed during forest inventories. Such 670 information is used to guide forest practitioners in their management aims, which further 671 include providing the demanded ecosystem services and the conservation of forest 672 biodiversity in the best possible way (Bäck et al., 2017; MEA, 2005). Assessment of the 673 standing timber volume therefore provides important information to practitioners because 674 management activities are very often addressed primarily in forest stands with the largest 675 timber volumes, as they contain the highest accumulated timber values and require cautious 676 planning of harvesting activities and stand regeneration. Further, the standing timber 677 volume is often used as an important indicator to quantify the provisioning of timber 678 production (Blattert et al., 2017; Bugmann et al., 2017), which is still one of the most 679 important ecosystem services in forestry. Improved inventory fundamentals have further 680 implications. They lead to improved management solutions, such as more efficient 681 harvesting or road network layouts (Bont et al., 2019, 2012), as timber volume is usually a 682 major component used to form such solutions (Bont and Church, 2018). With our presented 683 method for incorporating FTM information, predictions about the standing timber volume 684 Page 47 Appendix 888 889 Figure A1: Adjusted r-squared values for each regression model, forest type map source and combination of 890 forest type variables in the case study area Bremgarten (Brg). 891 Page 48 892 893 Figure A2: Adjusted r-squared values for each regression model, forest type map source and combination of 894 forest type variables in the case study areas Glâne-Farzin (GF) and Sarine (SA). 895 896 Page 49 897 Study area Variables Brg_OI MSSq + STD + MAX + P40 + P30 + P25 + P05 + P01 + C00 + C02 + C10 + C15 + C20 + MEAN2 Brg_OI_wo_La STD + P99 + P60 + P50 + P30 + P25 + P05 + P01 + C00 + C10 + C20 + CV + STD:MEAN Brg_NFI MSSq + STD + MAX + P40 + P30 + P25 + P05 + P01 + C00 + C02 + C10 + C15 + C20 + MEAN2 GF_NFI P95 + P80 + P60 + P40 + P25 + P10 + P01 + C00 + C02 + C05 + C10 + C15 + C20 + C25 + C30 + C35 + C45 + CV SA_NFI P99 + P95 + P20 + C02 + C05 + C15 + C20 + C25 + C35 + C50 + P60 Brg_DLT MSSq + STD + MAX + P40 + P30 + P25 + P05 + P01 + C00 + C02 + C10 + C15 + C20 + MEAN2 GF_DLT P95 + P80 + P60 + P40 + P25 + P10 + P01 + C00 + C02 + C05 + C10 + C15 + C20 + C25 + C30 + C35 + C45 + CV SA_DLT P99 + P95 + P20 + C02 + C05 + C15 + C20 + C25 + C35 + C50 + P60 Brg_Li MSSq + STD + MAX + P40 + P30 + P25 + P05 + P01 + C00 + C02 + C10 + C15 + C20 + MEAN2 Brg_Li_wo_La STD + P99 + P60 + P50 + P30 + P25 + P05 + P01 + C00 + C10 + C20 + CV + STD:MEAN GF_Li P95 + P80 + P60 + P40 + P25 + P10 + P01 + C00 + C02 + C05 + C10 + C15 + C20 + C25 + C30 + C35 + C45 + CV GF_Li_wo_La STD + P95 + P80 + P60 + P40 + P30 + P10 + P01 + C00 + C02 + C05 + C10 + C15 + C20 + C45 + MEAN2 SA_Li P99 + P95 + P20 + C02 + C05 + C15 + C20 + C25 + C35 + C50 + P60 SA_Li_wo_La MSSq + P25 + P05 + C02 + C05 + C15 + C20 + C25 + C35 + C50 + MEAN2 Brg_Sen MSSq + STD + MAX + P40 + P30 + P25 + P05 + P01 + C00 + C02 + C10 + C15 + C20 + MEAN2 GF_Sen P95 + P80 + P60 + P40 + P25 + P10 + P01 + C00 + C02 + C05 + C10 + C15 + C20 + C25 + C30 + C35 + C45 + CV SA_Sen P99 + P95 + P20 + C02 + C05 + C15 + C20 + C25 + C35 + C50 + P60 898 Table A1: Selected variables for Model 7. In addition to the listed variables, the FTV can be added 899 (area_prop, area_chm_prop, area_chm_prop, categorical or ground_proportion). 900 901