Predicting the provisioning potential of forest ecosystem services using airborne laser scanning data and forest resource maps

RESEARCH Open Access Predicting the provisioning potential of forest ecosystem services using airborne laser scanning data and forest resource maps Jari Vauhkonen Abstract Background: Remote sensing-based mapping of forest Ecosystem Service (ES) indicators has become increasingly popular. The resulting maps may enable to spatially assess the provisioning potential of ESs and prioritize the land use in subsequent decision analyses. However, the mapping is often based on readily available data, such as land cover maps and other publicly available databases, and ignoring the related uncertainties. Methods: This study tested the potential to improve the robustness of the decisions by means of local model fitting and uncertainty analysis. The quality of forest land use prioritization was evaluated under two different decision support models: either using the developed models deterministically or in corporation with the uncertainties of the models. Results: Prediction models based on Airborne Laser Scanning (ALS) data explained the variation in proxies of the suitability of forest plots for maintaining biodiversity, producing timber, storing carbon, or providing recreational uses (berry picking and visual amenity) with RMSEs of 15%–30%, depending on the ES. The RMSEs of the ALS-based predictions were 47%–97% of those derived from forest resource maps with a similar resolution. Due to applying a similar field calibration step on both of the data sources, the difference can be attributed to the better ability of ALS to explain the variation in the ES proxies. Conclusions: Despite the different accuracies, proxy values predicted by both the data sources could be used for a pixel-based prioritization of land use at a resolution of 250 m 2 , i.e., in a considerably more detailed scale than required by current operational forest management. The uncertainty analysis indicated that maps of the ES provisioning potential should be prepared separately based on expected and extreme outcomes of the ES proxy models to fully describe the production possibilities of the landscape under the uncertainties in the models. Keywords: Forestry decision making, Spatial prioritization, Light detection and ranging (LiDAR), Remote sensing Background Forestry decision making requires evaluating potential management alternatives with respect to multiple objectives (Kangas et al. 2008). A fundamental decision is related to which goods and services to produce: in addition to conventional timber production, the management objectives may be related to maintaining habitats, providing recreational and aesthetic opportunities, and carbon storage or sequestration (e.g. Pukkala 2016). These goods and services are jointly called “multiple uses”(Kangas 1992) or, following Costanza et al. (1997), Daily et al. (1997) and many others, “ecosystem services” of forest. In the following text, I use ESs to abbreviate “Ecosystem Services”, referring most essentially to indicators of forest-related ESs that can be derived from Remote Sensing (RS) or other digital map data as indirect proxies (Andrew et al. 2014). The mapping of these proxies allows spatial prioritization and other spatially explicit analyses of multiple ESs at various scales (e.g. Schröter et al. 2014; Räsänen et al. 2015; Sani et al. 2016; Roces-Díaz et al. 2017). According to reviews (Martínez-Harms and Balvanera, 2012; Englund et al. 2017) and a collection of case studies (Barredo et al. 2015), however, such analyses can be expected to suffer from the lack of standardized terminology, methodology and data. Increased attention should especially be Correspondence: [email protected] Natural Resources Institute Finland (Luke), Bioeconomy and Environment Unit, P.O. Box 68, Yliopistokatu 6, FI-80101 Joensuu, Finland © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Vauhkonen Forest Ecosystems (2018) 5:24 https://doi.org/10.1186/s40663-018-0143-1 focused on quantifying and communicating the resulting uncertainties to the decision makers in order to make informed decisions (see also Eigenbrod et al. 2010; Schulp et al. 2014; Foody 2015). Accounting for these aspects, the present study examines the robustness of forest land-use prioritization based on maps of the provisioning potential of forest ESs (Vauhkonen and Ruotsalainen 2017a), i.e., the fitness of forest patches to provide goods and services typical to the ESs occurring in the studied area, re-considering the methodological and data workflow proposed in the earlier study. To result in valid conclusions from RS-based decision analyses, the estimates should be accurate already at the level of individual pixels. The use of active RS such as Light Detection and Ranging (LiDAR) is expected to produce more accurate information compared to passive, optical RS (Lefsky et al. 2001; Coops et al. 2004; Maltamo et al. 2006), especially, when using small pixels (e.g., 200 m 2 as in Næsset 2002). Forest structure and habitat related inventories in particular benefit from the ability of LiDAR to provide three-dimensional information, when operated as Airborne Laser Scanning (ALS; Maltamo et al. 2014). Kankare et al. (2015) evaluated the estimation accuracy of biomass attributes based on two different RS setups in an area closely resembling to that presently studied. According to their results, pixel-level predictions based on coarse to medium resolution satellite imagery had a Root Mean Squared Error (RMSE) of 47.7% of the total biomass, which could be reduced to 25.7% using ALS and local field reference data. The ALS data used were acquired by the land survey, and the availability of such data is increasing due to large-area acquisitions for terrain elevation modelling. Such data have also been used to map attributes related to habitat (Melin et al. 2013,2016; Vauhkonen and Imponen 2016), structural (Valbuena et al. 2016b; Vauhkonen and Imponen 2016) and aesthetic (Vauhkonen and Ruotsalainen 2017b) properties of the forest. Overall, when various forest ESs are categorized according to a typology such as the Common International Classification of Ecosystem Services (CICES) as in Englund et al. (2017), the potential of ALS for assessing the suitability of forest areas to provide these ESs can be characterized as: –Regulation and maintenance services: A very high number of studies indicates that the vegetation height and density profiles produced by ALS are useful for a detailed quantification of variations in above-ground biomass (Næsset and Gobakken 2008; Zolkos et al. 2013; Popescu and Hauglin 2014) and, thus, carbon storage (Patenaude et al. 2004). Essentially, ALS produces a three-dimensional description of the forest structure, which can be related to ecological properties such as habitat types (Bässler et al. 2011) or biological diversity in general (Müller and Vierling 2014) and employed to assess suitability of forests to be maintained as habitats for different species (Davies and Asner 2014; Hill et al. 2014; Simonson et al. 2014). –Provisioning services: Several studies carried out especially in boreal forest structures indicate ALS data useful for assessing properties related to wood production. Except that the methods listed in the previous paragraphs can be directly used to assess the production potential of bulk biomass, also more detailed predictions of timber assortments (Korhonen et al. 2008; Kotamaa et al. 2010; Vauhkonen et al. 2014; Hou et al. 2016) or wood fiber-related attributes (Hilker et al. 2013; Luther et al. 2014) are possible. Although the yield studies are mostly related to wood-based biomass, there also are examples of improved assessments of the yield of shrub fruits (Barber et al. 2016) or edible fungi (Peura et al. 2016) based on ALS. –Cultural services: The applicability of ALS highly depends on the cultural service of interest. For example, several archaeological studies indicate the potential to improve the mapping of historical remains in the forest using an ALS-based digital terrain model. Similar techniques to visualize the terrain (Domingo-Santos et al. 2011) or trees (Lämås et al. 2015) could potentially be used to assess the aesthetic properties of the forest. To date, the study of Vauhkonen and Ruotsalainen (2017b), which assessed the preferences on the visual amenity of a forest area based on cuttings simulated to triangulated vegetation point clouds, appears to be the only ALS-based attempt towards this direction. The use of ALS can thus be motivated by the potential to obtain a better correspondence with forest biophysical attributes and these data may be available for some areas in a similar extent as land cover maps and other publicly available data. Despite the high potential, however, also ALS-based information may yield a high degree of uncertainties, if applied in expert models formulated according to conventionally measured field attributes. For example, the suitability index proposed by Pukkala et al. (2012) to map potential habitats of Siberian jay (Perisoreus infaustus L.) would require estimating the availability of Vaccinium myrtillus (L.) berries and epiphytic lichens for food and nests. Although sub-models to estimate these attributes are presented (Pukkala et al. 2012), also those include stand age and site fertility, which are difficult to estimate by ALS. Although some researchers have predicted even understorey-related attributes, the results of Korpela et al. (2012) indicate that Vauhkonen Forest Ecosystems (2018) 5:24 Page 2 of 19 direct measures are difficult to obtain due to transmission losses occurring in the upper canopy (see also Maltamo et al. 2005) and such estimations would be even more unreliable based on passive optical RS methods. Even the recognition of dominant tree species may be challenging in ALS-based inventories: despite promising results based solely on ALS (Ørka et al. 2013; Vauhkonen et al. 2014), the results of Räty et al. (2016) suggest difficulties in detecting species, which dominate a minor proportion of an area otherwise homogeneous in terms of the species. On the other hand, ALS may allow producing other attributes with more relevance from the forest management point of view. For example, forests with multilayered vertical structure can be distinguished based on the data (Zimble et al. 2003; Maltamo et al. 2005), which can be further employed in detecting the prevailing silvicultural system (Bottalico et al. 2014), management intensity (Sverdrup-Thygeson et al. 2016; Valbuena et al. 2016a), or development stage (Valbuena et al. 2016b). Even more detailed indices may be developed based on ecological rationale (Listopad et al. 2015) or a thorough understanding of the properties affecting the ALS response (Valbuena et al. 2013,2014). Earlier studies have suggested that the information in theALSdatamaybecondensedtoafewmetrics(Kane et al. 2010;Leitereretal.2015; Valbuena et al. 2017), the partitioning of which will provide a stratification corresponding closely to the structural complexity observed in the field (Pascual et al. 2008;Thompsonetal. 2016; Vauhkonen and Imponen 2016). Even though properties related to individual ESs have been actively studied, no studies that show how to support management decisions related to the provisioning of multiple forest ESs based on three-dimensional forest structure description obtained by ALS can currently be found from the literature. Barbosa and Asner (2017) and Rechsteiner et al. (2017) derived information from ALS data to prioritize landscapes for ecological restoration and species conservation planning, respectively. Packalén et al. (2011) used ALS data and spatial optimization to derive so called dynamic treatment units to guide the management of pulpwood production in a plantation forest. Although a similar approach could be extended to the decision making of other or multiple ESs (Pukkala et al. 2014), all ALS-based applications are, to date, focused on single ESs. The purpose of this study is to test ALS data for management prioritization of multiple ESs in a boreal forest landscape. Proxies for pixel-wise provisioning potential of biodiversity, carbon, timber, berries, and recreational amenities were formulated using ALS-based features and compared to information obtained from forest resource maps with a resolution of 16 × 16 m 2 .The quality of land use prioritization based on the obtained information was evaluated under two different decision support models: either using the developed models deterministically or in corporation with the uncertainties of the models. Methods A methodological overview Specifically, the ALS data are tested for predicting the provisioning potential of ESs (Vauhkonen and Ruotsalainen 2017a) in a spatial prioritization framework, where land use decisions are based on ranking the set of decision alternatives in the considered location(s) and choosing the best according to the decision makers’ preferences (cf., Malczewski and Rinner 2015). When applied to prioritize forests for single (e.g. Lehtomäki et al. 2015) or multiple uses (e.g. Vauhkonen and Ruotsalainen 2017a) based on ES proxy maps, a simplified workflow for such analyses includes three methodological steps: 1) Data acquisition, feature extraction and/or expert modelling to derive proxy values for the analyzed ESs. 2) Scaling and normalization of the proxy values derived from different sources to the same scale. The resulting values can be called ‘priority’,‘benefit’, or ‘utility’value and used in different ways depending on the literature source (see also Pukkala 2008; Pukkala et al. 2014; Malczewski and Rinner 2015). 3) Decision analyses using the normalized data at selected spatial scale(s). Because the normalized proxy maps resulting from the previous steps ‘measure’the ESs in a same scale and account for the value range of each ES in the entire landscape, they can be used (a) to mutually rank ESs within a spatial unit to subsequently prioritize management to provide most suitable ESs in each unit; and (b) to identify the most important locations of specific ESs in the landscape to be considered as management hot-spots or cold-spots. Because the spatial prioritization is carried out at a sub-stand-level using pixels or other corresponding map units, it is expected to allow a more efficient use of the production possibilities of the forest (Heinonen et al. 2007) and, overall, operationalize the concept of ESs for landscape planning, which is further motivated by de Groot et al. (2010). The present study examines whether changes to each of the three steps listed above could improve pixel-wise analyses of the provisioning potential of forest ESs (cf. the discussion section of Vauhkonen and Ruotsalainen 2017a): Vauhkonen Forest Ecosystems (2018) 5:24 Page 3 of 19 1) What data to use for the expert models of the provisioning potential: A consolidated approach to obtain grid-based, wall-to-wall predictions for the tessellated landscapes would be to use forest resource maps based on generalizing field sample plot measurements to larger areas using coarse to medium resolution RS images and other numeric map data (Tomppo et al. 2008a,2008b, 2014). This approach, referred to as Multi-Source National Forest Inventory (MS-NFI), was used by Vauhkonen and Ruotsalainen (2017a). Even if ALS allows more prediction possibilities, as reviewed above, it is practically reasoned to benchmark the accuracies against the pixel data provided by the MS-NFI approach, because different forest resource maps are readily available in many countries (Tomppo et al. 2008b, Roces-Díaz et al. 2017; Vauhkonen and Ruotsalainen, 2017a). 2) How to scale the ESs originally measured in different units for the joint analyses: Vauhkonen and Ruotsalainen (2017a) used a simple normalization to convert the ES values between 0 and 1: vij ¼nij N;ð1Þ where v ij is the normalized value and n ij is the position of the j:th plot in ascending order of the expert model values for the i:th ecosystem service among altogether N plots. Notably, this normalization produced values in an interval scale, whereas the ratios between the expert model values could also be assumed useful for the priority ranking. An alternative, ratio-scale normalization could be computed as: vij ¼ESij−min ESi ðÞ max ESi ðÞ−min ESi ðÞ ;ð2Þ where v ij is the value (or priority or benefit or utility, depending on literature source; see above) produced by the i:th ES in plot j. 3) How to use the obtained information in decision analyses: Vauhkonen and Ruotsalainen (2017a) deterministically prioritized each pixel to the ES with the highest predicted proxy value, but highlighted the need to consider uncertainties around the predictions. If a quantification of the uncertainties is obtained (e.g., by approximating residual errors of calibration models fitted to the data), the decision analyses can consider distributions of uncertainty in addition to the expected values and produce separate recommendations for different decision makers according to their attitudes towards risk (Pukkala and Kangas 1996). Therefore, in addition to deterministic use of the predicted values, this study considered both the expected and extreme outcomes of the predictions when selecting the most suitable ES for a pixel. The principal idea of this analysis is illustrated in Fig. 1. On this background, the present study tested the data source (ALS or MS-NFI), priority value function form (Eqs. 1or 2), uncertainty management approach, and joint implications of these choices to the predictions of the provisioning potential of forest ESs and subsequent management prioritization decisions. Forest ESs considered were selected based on two criteria: likelihood to occur in the studied landscape and existence of expert models to derive proxies for their provisioning potential based on the field measurements (Table 1). The field and MS-NFI data contained estimates of forest attribute that could be directly inserted to the expert models. Using ALS data, regression analyses were employed to estimate predictive relationships between ALS-features and ES proxy values to fully utilize the different properties of these data (cf., Section “ALS-based models for the priority values of the ESs”below). In the absence of independent, wall-to-wall data for validation, both the predictions and validations were carried out at the level of individual forest plots. The evaluation is therefore limited to the local fitness of the ESs for a specific forest patch in a single point in time and without considering their spatial or temporal continuum. No decision maker was assumed in this study and the values obtained from Fig. 1 A generic example of selecting the best decision alternative based on different outcomes of model predictions (colored curves). The yellow curve yields the highest priority value based on the expected (upper horizontal line) or worst outcome of the model. However, if the decision maker weights best possible outcomes, the alternative depicted by the grey curve should be selected as it produces the highest priority in the right tail accumulation point (the interception of the curve and the lower horizontal line) Vauhkonen Forest Ecosystems (2018) 5:24 Page 4 of 19 both the Eqs. 1and 2were therefore treated with equal weights, even if those could additionally be weighted according to the decision makers’preference structure. Study area and experimental data The study area is located in Evo, Finland (61.19°N, 25.11°E), which belongs to the southern boreal forest zone. The data extended over an area of approximately 3 km × 6 km. The forest stands in the area vary from intensively managed to natural forests in terms of their silvicultural status. Approximately 84% of the growing stock in the studied plots is dominated by coniferous tree species Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies [L.] H. Karst.). Deciduous tree species such as birches (Betula spp. L.), aspen (Populus tremula L.), alders (Alnus spp. P. Mill.), willows (Salix spp. L.), and rowan (Sorbus aucuparia L.) occur in mixed stands and below the dominant canopy. Data sets used were compiled from three earlier studies in the same area (Vauhkonen and Imponen 2016; Niemi and Vauhkonen 2016; Vauhkonen and Ruotsalainen 2017a). Vauhkonen and Imponen (2016) downloaded and processed ALS data acquired by the National Land Survey of Finland to stratify the area according to forest structural properties. The ALS data were acquired from a flying altitude of 2200 m using Leica ALS 50 scanner on 7 May, 2012, to yield a nominal pulse density of 0.8 m −2 . Circular sample plots (9 m radius) were placed by clustering the ALS data with respect to forest structural features, which was found to be an efficient strategy to distribute the sample across the spatial, size, and age distributions of the tree stock (Vauhkonen and Imponen 2016). The field measurements were carried out in June–August, 2014. The species and diameter-at-breast height (DBH) were measured for each tree with a DBH ≥5 cm. For each tree species of the plot, a tree with a DBH corresponding to the median tree was measured for height and used to calibrate height curves for predicting the missing tree heights. Plot-level forest attributes were computed from the tree-level measurements using standard equations and methods, which are described in detail in an open-access article by Niemi and Vauhkonen (2016). Publicly available MS-NFI data (Natural Resources Institute Finland 2017) were included to provide a benchmark for the ALS data. The MS-NFI maps are the same used Vauhkonen and Ruotsalainen (2017a) and details on their pre-processing are given in that paper. These raster maps depicted site fertility, growing stock volume and biomass components by tree species, total basal area and mean diameter and height corresponding to those of the (basal area weighted) median tree, and they were produced using a k-nearest neighbor (k-NN) estimation method based on optimized neighbor and feature selection (Tomppo and Halme 2004; Tomppo et al. 2008a, 2014). The method used various satellite images from 2012 to 2014 and National Forest Inventory (NFI) field plot measurements from 2009 to 2013, which were updated to correspond the situation in mid-2013 using growth models. Altogether 102 field plots were covered by both ALS and MS-NFI data and were included in the analyses. Table 1 The ESs considered in this study and expert models for deriving their reference proxy values Abbr. ES Indicator, unit (citation) a Stand-level forest attributes used as predictors b BIOD Biodiversity Index value based on expert opinion (Lehtomäki et al. 2015) 1 Site fertility, growing stock volume, diameter, dominant species TIMB Timber production Soil expectation value (SEV), €∙ha −1 (Pukkala 2005) 2 Diameter, basal area, age, site fertility, species-specific growing stock volume, number of trees, operational environment (temperature, interest rate, timber prices) CARB Carbon storage Estimated amount of carbon 3 ,t∙ha −1 (Karjalainen and Kellomäki 1996) Growing stock volume BILB Suitability for bilberry picking Index value based on expert opinion (Ihalainen et al. 2002) Age, basal area, height, species-specific growing stock volume, site fertility COWB Suitability for cowberry picking Index value based on expert opinion (Ihalainen et al. 2002) Age, species-specific growing stock volume, diameter, site fertility AMEN Visual amenity Index value based on expert opinion (Pukkala et al. 1988) Diameter, number of trees, species-specific growing stock volume, site fertility a When computing the values for the present study, the following details or exceptions compared to original publications were made: 1 The index values are of form diameter × volume, scaled using dominant-species-specific transformation functions (Lehtomäki et al. 2015) and maximum values of forest attributes in the study area, and multiplied by site fertility specific weights (Lehtomäki et al. 2015). 2 Values of operational environment related parameters were obtained as combinations of effective temperature sum fixed to 1300 degree days, interest rates of 1%–4% and saw-wood/pulpwood prices (units in €∙m −3 ) of 30/15, 30/25, 40/15, 40/25, 40/35, 50/25, and 50/35, and the SEV was obtained as an average of these 28 combinations weighted by the proportions of species. All values were adopted from the study by Pukkala (2005). 3 The estimated carbon was obtained based on conversion factors from species-specific, total stem volumes to carbon contents. b To standardize the computation based on all data sets, the following simplifications or groupings were used: -Species groups: pine, spruce, deciduous trees. -‘Diameter’always referred to the basal-area weighted mean diameter. Vauhkonen Forest Ecosystems (2018) 5:24 Page 5 of 19 The models of Table 1were applied to produce plot-specific reference values for the provisioning potential of the ESs based on field data. According to an exploratory analysis, the expert models of Table 1, fit with many different data sets, had considerably different value ranges over the landscape. As a result, a direct normalization of the expert function values specifically with Eq. 2resulted to emphasizing one ES in the priority rankings only because of the different shape and scale of initial value distributions, as elaborated upon in Appendix 1. For this reason, the expert function values of all ESs were transformed to follow the normal distribution as closely as possible using the Box-Cox-transformation (Appendix 1) prior to applying Eqs. 1and 2. The forest attribute estimates based on the MS-NFI maps were transformed using the same parameter values as with field data. This transformation did not affect the order of the observations, but produced approximately equally shaped frequency distributions of every ES, as detailed in Appendix 1. ALS-based models for the priority values of the ESs Prediction models with independent variables extracted from the ALS data were formulated to predict priority function values of the form of Eq. 2.Priority function values corresponding to Eq. 1were obtained by ordering the aforementioned predictions, i.e., no separate models were constructed for the function form of Eq. 2. As reasoned in the Introduction, the aim was not to model the forest attributes used as the predictors of the expert models, but to identify and quantify such properties of the ALS point clouds that directly explained the variation in the ES proxies. As visualized in Fig. 2, the point clouds of the plots with maximum proxy values did not considerably differ between the ESs in terms of the total distributions. However, when height values or proportions were computed separately according to echo categories, ES-specific differences could be pointed out (Fig. 2). The features were therefore extracted in echo categories, which were “only echoes”(suffix _only), “first of many echoes”(_first), “last of many echoes”(_last), “first echoes”(_FP), and “last echoes”(_LP), where the last two categories included “first of many”and “last of many”echoes, respectively, with “only echoes”duplicated in both. Fixed height values of 0.5 m, 5 m, and below or above an adaptive height value determined as the height of the 60th percentile were used as the Fig. 2 ALS height profiles and descriptive characteristics of the field plots considered to be most important locations of the ESs in the data studied (priority value of 1 based on Eq. 2). For comparison, the lower right panel shows a plot that had low priority values of the considered ESs. The black, green, and blue symbols indicate only, first-of-many, and last-of-many ALS echoes, respectively. Grey horizontal lines indicate the mean heights of these echo categories and all echoes and are drawn to illustrate the differences in terms of these metrics between the ESs Vauhkonen Forest Ecosystems (2018) 5:24 Page 6 of 19 thresholds of ground, shrub, and suppressed or dominant canopy, respectively. The following categories of the features were considered: –Canopy height and density, which are the basic predictors used in ALS analyses (Næsset 2002) and were assumed to discriminate between size-specific attributes of the ESs: the maximum (hmax), the mean (hmean), and the standard deviation (hstd) of the height values above the ground threshold; the 5th, 10th, 20th, ..., 90th, and 95th percentiles (hzz, where zz denoted the percentile value); and the corresponding proportional densities (dzz) were computed according to Korhonen et al. (2008, pp. 502–503). –Proportion of echoes above a given threshold to all echoes, corresponding to a vegetation cover estimate (Korhonen et al. 2011). This proportion was computed in two ways: using echoes of different categories above the ground (cc X_ground ,whereXis the echo category) or first echoes above the shrub layer threshold (cc shrub ), which corresponds to an attempt to quantify the shrub layer thickness (cf., Vauhkonen and Imponen 2016). –Absolute differences between mean heights of different echo categories. These features were computed without height thresholds and assumed to discriminate between properties related to coniferousor deciduous-dominated forest in the ALS data acquired during the leaf-off period (Liang et al. 2007). These features are denoted by diff x–y , where suffix x–yrefers to the height difference of echo categories FP–LP, only–LP, first–only, or first–last. –Proportions of the different echo categories, which were assumed to be affected by the species and size specific ES properties in the canopy similar to the ALS-intensity features (Ørka et al. 2012; Vauhkonen et al. 2014). These features are denoted by prop X/Y_z , where X/Y indicated the ratio of two echo categories Xand Y, and zwas the height threshold employed for computations. –Predictors related to the shrub and understorey layers (Vauhkonen and Imponen 2016): the ratio of the echoes reflected above ground but below the dominant canopy threshold (r understory ); the standard deviation of the height values of echoes reflected above ground but below the dominant canopy threshold (std understory ); and the ratio of the echoes reflected from the shrub layer to all echoes (r shrub ). Features x i ,i=1, 2, …, 143, listed above formed the initial set S 1 of candidate predictors. To account for useful interactions between the features, the final set Swas obtained as S 1 ∪{x i ×x j }∀i,j∈S 1 , which resulted to altogether 10,296 candidate features per plot. Separate models for each ES were constructed by inserting features iteratively into a model template: ^ yin ¼anþXN n¼1bnxcn jn;ð3Þ where ^ yin is the vector of predicted priority values for the i:th ES, x jn is the j:th feature of S, and a n ,b n , and c n are model parameters at the n:th round of N=1,2,3,4 iteration rounds. Parameters a n ,b n , and c n were estimated using the nls function of R statistical computing environment (R Core Team 2016). Testing every candidate feature as x jn at every iteration round, the RMSE between the predicted and reference priority values was computed as: RMSE ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn i¼1^ yi−yi ðÞ 2 n; sð4Þ where nis the number of observations, and ^ yiand y i are the predicted and reference values, respectively. The feature that minimized the RMSE was retained in the model template and the iterations were continued until the model included a maximum of four features. However, more criteria were employed to select the model to be used for the prioritization analyses among the models with one to four features: 1) The final predictor inserted had to improve the RMSE by at least 1%. 2) The residual errors had to satisfy the null hypothesis that the considered sample came from a normally distributed population, which was examined graphically using scatter, residual and QQ-plots, and numerically using the test statistic proposed by Shapiro and Wilk (1965). 3) The model had to pass a “sensitivity of convergence”test, in which the model was fit separately for each plot using Leave-One-OutCross-Validation (LOOCV), i.e., not allowing the plot in question to be available in the training data for model fitting. Implications of including this test are further described in the Results section. Predicting the priority values of the ESs based on the MS-NFI maps Benchmark predictions for those based on ALS were obtained by inserting the forest attribute estimates from the MS-NFI maps to Eq. 2. Priority function values corresponding to Eq. 1were obtained by ordering the aforementioned predictions (cf., previous section). The MS-NFI maps included estimates of all other independent variables except the number of trees per hectare, Vauhkonen Forest Ecosystems (2018) 5:24 Page 7 of 19 which was estimated by dividing the total basal area by the basal area corresponding to the mean diameter, i.e., assuming that the resulting number of average-sized trees existed in a pixel. To compute plot-wise estimates, the pixels of the forest resource maps intersecting with the plot polygons were identified using a spatial query. The estimates of a plot were obtained from the intersecting pixels as weighted averages with the joint areas of the plots and pixels as the weights. Finally, to see if amending the models based on ALS with the MS-NFI layers improved the models, a similar feature selection as with ALS data was run including all MS-NFI-based ES and forest attribute proxies as additional feature candidates. Field calibration and evaluation of the predictions Following the method described in the previous section, potential estimation errors in the MS-NFI maps propagate to the predicted priority values, whereas similar error propagation is avoided in the ALS-based analyses due to local model fitting. An additional calibration step was therefore included to eliminate the contribution of the local field sample to the predictions. Calibration models yi¼fð^ yiÞ, where y i was the reference priority value of the i:th ES and ^ yiits RS-based estimate, of all ESs were fit simultaneously as systems of linear equations. Due to the high inter-correlations (see Additional file 1, Table S1), the models were fit in two steps: first, using Ordinary Least Squares (OLS) to produce model residuals, and second, using Seemingly Unrelated Regression (SUR) to account for the residual error covariance matrices in the final models. The computations were carried out in the LOOCV mode using the systemfit package of R (Henningsen and Hamann 2007). The accuracies of the ALSand MS-NFI-based predictions were compared using the RMSE and coefficient of determination (R 2 ) computed between the reference values and predictions obtained from the LOOCV models. Decision analyses The effects of the aforementioned prediction accuracies to the management decisions were evaluated by comparing the priority ranking of the ESs in each individual plot. The ES with the highest priority value, based on Eqs. 1or 2applied to the field reference data, was assumed to be the most suitable ES for the specific plot. The RS-based decision was considered correct, if the most suitable ES based on the field data and the RS prediction equaled. The degree of incorrect decisions was quantified using two approaches. First, the correctness of every decision was given a numerical score (Gopal and Woodcock 1994): situations where RS and field data resulted in the same decision was given a score of 6; those where the RS-based service was the second best according to the field data a score of 5; and so on, until the situation where the RS-based service was the worst according to field data, which was given a score of 1. The distributions of these “decision scores”were compared between the different data sources. Second, the dispersion in field and RS-data between the services selected as the most suitable for the specific plot was examined using confusion matrices. The priority ranking of the less important ESs was not evaluated. In addition to ‘deterministic’decision making described above, the sensitivity of the decisions was examined by incorporating the uncertainties of the models to the analyses (Fig. 1). Instead of using the expected values of the priority functions, the ranking was carried out assuming the predictions as realized values of a random variable X~N(E,s 2 ), where Ewas the expected value and s 2 was the mean squared error of the model residuals. A similar priority ranking as with the expected values was carried out with predictions that were among the worst and best outcomes of the model, obtained as the values of the 5th and 95th percentiles of the distribution of Xof each ES. Results ALS-based models for the priority values of the ESs The ALS features considered as the predictors of the regression models are listed in Table 2. All feature and echo type categories and a wide range of different height values was employed when building the models, for which reason only a few specific observations on the structure of models can be made. All features selected were products of form feature 1 ×feature 2 ,wherefeature 1 was often an absolute height value (a mean height, percentile, or height difference) and feature 2 a proportion (either a canopy cover proxy or proportional density). This combination was especially frequent among the first features selected to the models. In the models of TIMB, all selected predictors were such combinations employing various height values and echo categories. The models of BIOD and CARB used proportion ×proportion types of interactions and the ratio of first-of-many to only and first returns (prop first/FP_ground ). The predictors of BIOD (e.g., prop first/FP_ground ;diff only–LP ;hstd LP )were most diverse in terms of describing the canopy structure with features from different categories. The models of BILB, COWB, and AMEN differed from those mentioned above in employing low percentile values, last pulse proportions and features such as r understory , diff only–FP ,prop first/FP_ground ,andhstd first . Vauhkonen Forest Ecosystems (2018) 5:24 Page 8 of 19 Overall, the canopy cover proxies were the most frequent feature type, whereas the computing heights and echo categories of all features varied. Although a wide range of different height values was used, a predictor with a percentile value above 70 was selected only once. The graphical assessment of model residuals (detailed results not shown) was mainly in line with the test on residual normality (Table 2): the QQ-plots showed heavy-tailed residuals especially for the models with statistically significant values of the Shapiro-Wilk test statistic. However, the deviations of normality were typically related to one or two plots with the highest or lowest values, and not considered problematic for the further analyses. When examined in the same data used for constructing the models, the performance of every model could be slightly improved by increasing the number of predictors to the maximum number allowed. However, when the models were re-fit using LOOCV, model parameters could not be solved for at least one of the plots in the data, resulting to NA values for this performance factor in Table 2. Although this effect could probably have been avoided by allowing a slightly wider range of initial parameters when fitting the models, it was also considered as a sensitivity issue reflecting an over-parameterization of the initial model to certain types of forest structures. The volume of deciduous trees and the MS-NFI based proxy for BILB would have replaced the last ALS-based features in the models of CARB and BILB, respectively, and in the models of COWB, the corresponding MS-NFI proxy would have been selected as the second feature. However, none of the aforementioned MS-NFI features performed better than the ALS-features of these models in terms of the feature selection criteria. Based on the considerations above, the ALS and MS-NFI data sets were always used separately. Also, a different number of predictors was used in the ALS-based models for the priority ranking: BIOD and COWB were modeled using only one predictor (the one selected first); TIMB, BILB, and AMEN using two predictors (those selected first and second); and CARB using three predictors selected first. Comparison of ALS and MS-NFI for predicting the priority values of the ESs The models based on ALS data always outperformed those based on forest attribute estimates derived from the MS-NFI maps. As shown in Figs. 3and 4,the ALS-based models generally explained more variation in the ES proxies. The regression lines of the MS-NFI data based on the SUR calibration models also differed more severely from the 0–1-lines. Using MS-NFI data, TIMB was predicted most accurately with an RMSE of 30.4%. The RMSEs of other ESs were also close (30.6%–33.2%), except AMEN, which had an RMSE of 40.5% and BIOD, which was predicted worst with an RMSE of 41.6%. Using ALS, the ES predicted worst (COWB) had an RMSE of 29.8%, which is 97% of the RMSE of the corresponding MS-NFI prediction. The RMSEs of all other ESs were in order of 21.7%–27.5% (57%–83% of the RMSEs of MS-NFI predictions), except CARB, which was predicted most accurately with an RMSE of 15.1% (47% of the RMSE of MS-NFI prediction). The degree of determination of CARB also improved most due to using ALS instead of MS-NFI, from R 2 = 0.11 to 0.81. The R 2 -improvements of the other ESs were close to this magnitude, except for BILB and COWB, which had R 2 values close to each other based on both the data sources. The residual errors of models based on ALS and MS-NFI were somewhat correlated for BILB and COWB, but not for the other ESs (Figs. 3and 4, right column). Table 2 The features and performance of ALS-based models for predicting ratio-scaled ES proxy values. W –Shapiro-Wilk test statistic ES Predictor W a RMSE RMSE LOOCV BIOD cc shrub × h40 first 0.965 *** 0.259 0.266 + prop first/FP_ground × d50 first 0.981 0.235 NA + diff only–LP × hstd LP 0.986 0.217 0.226 + h95 first × h10 LP 0.977 * 0.203 NA TIMB cc only_ground × hmean FP 0.977 * 0.220 0.228 + h40 last ×cc LP_ground 0.970 ** 0.202 0.213 + h05 last × d05 first 0.989 0.182 NA +cc only_ground × h10 FP 0.988 0.174 NA CARB cc shrub × h60 first 0.980 0.158 0.163 + h20 last ×cc only_ground 0.988 0.144 0.152 + d60 first × d30 LP 0.987 0.138 0.148 + prop first/FP_ground × d05 first 0.984 0.132 NA BILB h60 first × h70 LP 0.987 0.279 0.286 + d50 only ×cc FP_ground 0.984 0.255 0.274 +r understory × d05 FP 0.982 0.238 NA + d70 first × d50 only 0.984 0.222 0.267 COWB d20 LP × h40 FP 0.966 *** 0.281 0.295 + diff only–FP 2 0.981 0.250 NA + prop first/FP_ground × d05 first 0.983 0.233 NA + d60 first × h05 only 0.971 ** 0.217 NA AMEN h10 first × hmean FP 0.993 0.239 0.246 + d30 last × d70 first 0.991 0.219 0.229 + h20 last ×cc FP_ground 0.994 0.204 NA + hstd first × hmean LP 0.991 0.187 NA a The asterisks refer to the significance of the test statistic at the 90% (*), 95% (**), and 99% (***) confidence level Vauhkonen Forest Ecosystems (2018) 5:24 Page 9 of 19 Appendix 2 Confusion matrices for the most important ESs based on expected outcomes Appendix 3 Confusion matrices for the most important ESs based on extreme outcomes Table 3 Confusion between ESs considered as most important based on the field data (observed) and MS-NFI maps (predicted) using the expected values of the predicted ES proxies Predicted BILB COWB AMEN BIOD CARB TIMB Observed BILB 8 5 2 0 2 4 COWB 6 15 9 1 1 2 AMEN 1 1 0 0 0 0 BIOD 1 0 4 4 5 2 CARB 0 1 4 1 1 1 TIMB 3 4 2 1 3 8 Table 4 Confusion between ESs considered as most important based on the field data (observed) and MS-NFI data calibrated with the local field sample (predicted) using the expected values of the predicted ES proxies Predicted BILB COWB AMEN BIOD CARB TIMB Observed BILB 11 7 0 0 1 2 COWB 9 18 0 2 2 3 AMEN 1 1 0 0 0 0 BIOD 10 0 0 5 1 0 CARB 1 2 0 1 3 1 TIMB 3 4 0 1 3 10 Table 5 Confusion between ESs considered as most important based on the field data (observed) and the expected values of the ALS-based models for ES proxies (predicted) Predicted BILB COWB AMEN BIOD CARB TIMB Observed BILB 12 2 0 3 0 4 COWB 5 17 1 9 0 2 AMEN 1 0 0 1 0 0 BIOD 4 3 0 3 5 1 CARB 0 0 0 2 2 4 TIMB 1 3 0 3 5 9 Table 6 Confusion between ESs considered as most important based on the field data (observed) and MS-NFI data calibrated with the local field sample (predicted) using the worst outcomes of the predicted ES proxies Predicted BILB COWB AMEN BIOD CARB TIMB Observed BILB 9 7 0 0 0 5 COWB 9 18 0 1 1 5 AMEN 0 1 0 0 0 1 BIOD 8 0 0 0 1 7 CARB 1 2 0 1 0 4 TIMB 3 4 0 0 1 13 Table 7 Confusion between ESs considered as most important based on the field data (observed) and the worst outcomes of the ALS-based models for ES proxies (predicted) Predicted BILB COWB AMEN BIOD CARB TIMB Observed BILB 9 0 0 1 10 1 COWB 5 2 3 0 23 1 AMEN 1 0 0 0 1 0 BIOD 1 0 1 0 14 0 CARB 0 0 0 0 7 1 TIMB 1 1 1 0 13 5 Table 8 Confusion between ESs considered as most important based on the field data (observed) and MS-NFI data calibrated with the local field sample (predicted) using the best outcomes of the predicted ES proxies Predicted BILB COWB AMEN BIOD CARB TIMB Observed BILB 8 4 0 7 2 0 COWB 8 9 0 16 0 1 AMEN 0 1 0 1 0 0 BIOD 2 0 0 14 0 0 CARB 0 0 0 8 0 0 TIMB 3 2 0 13 0 3 Vauhkonen Forest Ecosystems (2018) 5:24 Page 16 of 19 Abbreviations ALS: Airborne Laser Scanning; AMEN: Visual amenity (one of the ecosystem services considered, see Table 1); BILB: Suitability for bilberry picking (one of the ecosystem services considered, see Table 1); BIOD: Biodiversity (one of the ecosystem services considered, see Table 1); CARB: Carbon storage (one of the ecosystem services considered, see Table 1); COWB: Suitability for cowberry picking (one of the ecosystem services considered, see Table 1); DBH: Diameter-at-breast-height; ES: Ecosystem service; k-NN: k-Nearest neighbor; LOOCV: Leave one out cross validation; MCDA: Multiple criteria decision analysis; MS-NFI: Multi-Source National Forest Inventory; RMSE: Root mean squared error; RS: Remote sensing; TIMB: Timber production (one of the ecosystem services considered, see Table 1) Funding The acquisition of the studied data was originally supported by the Research Funds of University of Helsinki. Availability of data and materials All data and materials can be obtained by requesting from the author. Author’s contributions JV is the sole author. He carried out all analyses and drafted the manuscript. The author read and approved the final manuscript. Ethics approval and consent to participate Not applicable. Competing interests The author declares that he has no competing interests. Received: 14 March 2018 Accepted: 24 May 2018 References Andrew ME, Wulder MA, Nelson TA (2014) Potential contributions of remote sensing to ecosystem service assessments. Progr Phys Geogr 38:328–353 Barber QE, Bater CW, Braid ACR, Coops NC, Tompalski P, Nielsen SE (2016) Airborne laser scanning for modelling understory shrub abundance and productivity. For Ecol Manag 377:46–54 Barbosa JM, Asner GP (2017) Prioritizing landscapes for restoration based on spatial patterns of ecosystem controls and plant–plant interactions. J Appl Ecol 54:1459–1468 Barredo JI, Bastrup-Birk A, Teller A, Onaindia M, de Manuel BF, Madariaga I, Rodriguez-Loinaz G, Pinho P, Nunes A, Ramos A, Batista M, Mimo S, Cordovil C, Branquinho C, Gret-Regamey A, Bebi P, Brunner SH, Weibel B, Kopperoinen L, Itkonen P, Viinikka A, Chirici G, Bottalico F, Pesola L, Vizzarri M, Garfi V, Antonello L, Barbati A, Corona P, Cullotta S, Giannico V, Lafortezza R, Lombardi F, Marchetti M, Nocentini S, Riccioli F, Travaglini D, Sallustio L, Rosario I, von Essen M, Nicholas KA, Maguas C, Rebelo R, Santos-Reis M, Santos-Martin F, Zorrilla-Miras P, Montes C, Benayas J, Martin-Lopez B, Snall T, Berglund H, Bengtsson J, Moen J, Busetto L, San-Miguel-Ayanz J, Thurner M, Beer C, Santoro M, Carvalhais N, Wutzler T, Schepaschenko D, Shvidenko A, Kompter E, Ahrens B, Levick SR, Schmullius C (2015) Mapping and assessment of forest ecosystems and their services –Applications and guidance for decision making in the framework of MAES. Report EUR 27751 EN, Joint Research Centre, European Union. doi: https://doi.org/10.2788/720519 Bässler C, Stadler J, Müller J, Förster B, Göttlein A, Brandl R (2011) LiDAR as a rapid tool to predict forest habitat types in Natura 2000 networks. Biodivers Conserv 20:465–481 Bottalico F, Travaglini D, Chirici G, Marchetti M, Marchi E, Nocentini S, Corona P (2014) Classifying silvicultural systems (coppices vs. high forests) in Mediterranean oak forests by airborne laser scanning data. Eur J Remote Sens 47:437–460 Box GEP, Cox DR (1964) An analysis of transformations. J Royal Stat Soc Ser B 26: 211–252 Brokaw N, Lent R (1999) Vertical structure. In: Hunter ML Jr (ed) Maintaining biodiversity in Forest ecosystems. Cambridge University Press, Cambridge, pp 373–399 Coops NC, Wulder MA, Culvenor DS, St-Onge B (2004) Comparison of forest attributes extracted from fine spatial resolution multispectral and lidar data. Can J Remote Sens 30:855–866 Costanza R, d’Arge R, de Groot R, Farber S, Grasso M, Hannon B, Limburg K, Naeem S, O’Neill RV, Paruelo J, Raskin RG, Sutton P, van den Belt M (1997) The value of the world’s ecosystem services and natural capital. Nature 387: 253–260 Daily GC, Alexander S, Ehrlich PR, Goulder L, Lubchenco J, Matson PA, Mooney HA, Postel S, Schneider SH, Tilman D, Woodwell GM (1997) Ecosystem services: benefits supplied to human societies by natural ecosystems. Issues Ecol 2:1–16 Davies AB, Asner GP (2014) Advances in animal ecology from 3D-LiDAR ecosystem mapping. Trends Ecol Evol 29:681–691 de Groot RS, Alkemade R, Braat L, Hein L, Willemen L (2010) Challenges in integrating the concept of ecosystem services and values in landscape planning, management and decision making. Ecol Compl 7:260–272 Domingo-Santos JM, de Villarán RF, Rapp-Arrarás Í, de Provens ECP (2011) The visual exposure in forest and rural landscapes: an algorithm and a GIS tool. Landscape Urban Plan 101:52–58 Dueser RD, Shugart HH Jr (1978) Microhabitats in a forest-floor small mammal fauna. Ecology 59:89–98 Eigenbrod F, Armsworth PR, Anderson BJ, Heinemeyer A, Gillings S, Roy DB, Thomas CD, Gaston KJ (2010) The impact of proxy-based methods on mapping the distribution of ecosystem services. J Appl Ecol 47:377–385 Englund O, Berndes G, Cederberg C (2017) How to analyse ecosystem services in landscapes –a systematic review. Ecol Indic 73:492–504 Foody GM (2015) Valuing map validation: the need for rigorous land cover map accuracy assessment in economic valuations of ecosystem services. Ecol Econ 111:23–28 Gopal S, Woodcock C (1994) Theory and methods for accuracy assessment of thematic maps using fuzzy sets. Photogramm Eng Remote Sens 60:181–188 Hegetschweiler KT, Plum C, Fischer C, Brändli UB, Ginzler C, Hunziker M (2017) Towards a comprehensive social and natural scientific forest-recreation monitoring instrument –a prototypical approach. Landscape Urban Plan 167:84–97 Heinonen T, Kurttila M, Pukkala T (2007) Possibilities to aggregate raster cells through spatial optimization in forest planning. Silva Fenn 41:89–103 Henningsen A, Hamann JD (2007) Systemfit: a package for estimating systems of simultaneous equations in R. J Stat Softw 23(4):1–40 Hilker T, Frazer GW, Coops NC, Wulder MA, Newnham GJ, Stewart JD, van Leeuwen M, Culvenor DS (2013) Prediction of wood fiber attributes from LiDAR-derived forest canopy indicators. For Sci 59:231–242 Hill RA, Hinsley SA, Broughton RK (2014) Assessing habitats and organism-habitat relationships by airborne laser scanning. In: Maltamo M, Næsset E, Vauhkonen J (eds) Forestry applications of airborne laser scanning. Managing Forest ecosystems, vol 27. Springer, Dordrecht, pp 335–356 Hou Z, Xu Q, Vauhkonen J, Maltamo M, Tokola T (2016) Species-specific combination and calibration between area-based and tree-based diameter distributions using airborne laser scanning. Can J For Res 46:753–765 Ihalainen M, Alho J, Kolehmainen O, Pukkala T (2002) Expert models for bilberry and cowberry yields in Finnish forests. For Ecol Manag 157:15–22 Kane VR, McGaughey RJ, Bakker JD, Gersonde RF, Lutz JA, Franklin JF (2010) Comparisons between field-and LiDAR-based measures of stand structural complexity. Can J For Res 40:761–773 Table 9 Confusion between ESs considered as most important based on the field data (observed) and the best outcomes of the ALS-based models for ES proxies (predicted) Predicted BILB COWB AMEN BIOD CARB TIMB Observed BILB 11 5 0 4 0 1 COWB 4 20 1 8 0 1 AMEN 0 1 0 1 0 0 BIOD 4 3 0 7 0 2 CARB 0 2 0 4 0 2 TIMB 1 6 0 8 0 6 Vauhkonen Forest Ecosystems (2018) 5:24 Page 17 of 19 Kangas A, Kangas J, Kurttila M (2008) Decision support for forest management. Managing forest ecosystems 16. Springer, Dordrecht Kangas A, Leskinen P, Kangas J (2007) Comparison of fuzzy and statistical approaches in multicriteria decisionmaking. For Sci 53:37–44 Kangas J (1992) Multiple-use planning of forest resources by using the analytic hierarchy process. Scand J For Res 7:259–268 Kankare V, Vauhkonen J, Holopainen M, Vastaranta M, Hyyppä J, Hyyppä H, Alho P (2015) Sparse density, leaf-off airborne laser scanning data in aboveground biomass component prediction. Forests 6:1839–1857 Karjalainen T, Kellomäki S (1996) Greenhouse gas inventory for land use change and forestry in Finland based on international guidelines. Mitig Adapt Strat Glob Change 1:51–71 Kohler M, Devaux C, Grigulis K, Leitinger G, Lavorel S, Tappeiner U (2017) Plant functional assemblages as indicators of the resilience of grassland ecosystem service provision. Ecol Indic 73:118–127 Koivuniemi J, Korhonen KT (2006) Inventory by compartments. In: Kangas A, Maltamo M (eds) Forest inventory: methodology and applications. Managing Forest ecosystems, vol 10. Springer, Dordrecht, pp 271–278 Korhonen L, Korpela I, Heiskanen J, Maltamo M (2011) Airborne discrete-return LIDAR data in the estimation of vertical canopy cover, angular canopy closure and leaf area index. Remote Sens Environ 115:1065–1080 Korhonen L, Peuhkurinen J, Malinen J, Suvanto A, Maltamo M, Packalén P, Kangas J (2008) The use of airborne laser scanning to estimate sawlog volumes. Forestry 81:499–510 Korpela I, Hovi A, Morsdorf F (2012) Understory trees in airborne LiDAR data - selective mapping due to transmission losses and echo-triggering mechanisms. Remote Sens Environ 119:92–104 Kotamaa E, Tokola T, Maltamo M, Packalén P, Kurttila M, Mäkinen A (2010) Integration of remote sensing-based bioenergy inventory data and optimal bucking for stand-level decision making. Eur J For Res 129:875–886 Lämås T, Sandström E, Jonzén J, Olsson H, Gustafsson L (2015) Tree retention practices in boreal forests: what kind of future landscapes are we creating? Scand J For Res 30:526–537 Lefsky MA, Cohen WB, Spies TA (2001) An evaluation of alternate remote sensing products for forest inventory, monitoring, and mapping of Douglas-fir forests in western Oregon. Can J For Res 31:78–87 Lehtomäki J, Tuominen S, Toivonen T, Leinonen A (2015) What data to use for forest conservation planning? A comparison of coarse open and detailed proprietary forest inventory data in Finland. PLoS One. https://doi.org/10. 1371/journal.pone.0135926 Leiterer R, Furrer R, Schaepman ME, Morsdorf F (2015) Forest canopy-structure characterization: a data-driven approach. For Ecol Manag 358:48–61 Liang X, Hyyppä J, Matikainen L (2007) Deciduous-coniferous tree classification using difference between first and last pulse laser signatures. In: Rönnholm P, Hyyppä H, Hyyppä J (eds) Proceedings of ISPRS workshop on laser scanning 2007 and SilviLaser 2007. Int arch Photogramm remote Sens, vol XXXVI, part 3/W52, pp 253–257 Listopad CM, Masters RE, Drake J, Weishampel J, Branquinho C (2015) Structural diversity indices based on airborne LiDAR as ecological indicators for managing highly dynamic landscapes. Ecol Indic 57:268–279 Luther JE, Skinner R, Fournier RA, van Lier OR, Bowers WW, Coté JF, Hopkinson C, Moulton T (2014) Predicting wood quantity and quality attributes of balsam fir and black spruce using airborne laser scanner data. Forestry 87:313–326 MacArthur RH, MacArthur JW (1961) On bird species diversity. Ecology 42:594–598 Malczewski J, Rinner C (2015) Multicriteria decision analysis geographic information science. Advances in Geographic Information Science SpringerVerlag, Berlin Heidelberg Maltamo M, Malinen J, Packalén P, Suvanto A, Kangas J (2006) Nonparametric estimation of stem volume using airborne laser scanning, aerial photography, and stand-register data. Can J For Res 36:426–436 Maltamo M, Næsset E, Vauhkonen J (eds) (2014) Forestry applications of airborne laser scanning - concepts and case studies. Managing Forest ecosystems, vol 27. Springer, Dordrecht Maltamo M, Packalén P, Yu X, Eerikäinen K, Hyyppä J, Pitkänen J (2005) Identifying and quantifying structural characteristics of heterogeneous boreal forests using laser scanner data. For Ecol Manag 216:41–50 Martínez-Harms MJ, Balvanera P (2012) Methods for mapping ecosystem service supply: a review. Int J biodiv Sci Ecosyst Serv Manage 8:17–25 Melin M, Mehtätalo L, Miettinen J, Tossavainen S, Packalen P (2016) Forest structure as a determinant of grouse brood occurrence –an analysis linking LiDAR data with presence/absence field data. For Ecol Manag 380:202–211 Melin M, Packalen P, Matala J, Mehtätalo L, Pusenius J (2013) Assessing and modeling moose (Alces alces) habitats with airborne laser scanning data. Int J Appl Earth Obs Geoinfo 23:389–396 Müller J, Vierling K (2014) Assessing biodiversity by airborne laser scanning. In: Maltamo M, Næsset E, Vauhkonen J (eds) Forestry applications of airborne laser Scanning.Managing Forest ecosystems, vol 27. Springer, Dordrecht, pp 357–374 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 Næsset E, Gobakken T (2008) Estimation of aboveand below-ground biomass across regions of the boreal forest zone using airborne laser. Remote Sens Environ 112:3079–3090 Natural Resources Institute Finland (2017) File service for publicly available data. http://kartta.metla.fi/index-en.html. Accessed 16 Oct 2017 Niemi MT, Vauhkonen J (2016) Extracting canopy surface texture from airborne laser scanning data for the supervised and unsupervised prediction of areabased forest characteristics. Remote Sens. https://doi.org/10.3390/rs8070582 Ørka HO, Dalponte M, Gobakken T, Næsset E, Ene LT (2013) Characterizing forest species composition using multiple remote sensing data sources and inventory approaches. Scand J For Res 28:677–688 Ørka HO, Gobakken T, Næsset E, Ene L, Lien V (2012) Simultaneously acquired airborne laser scanning and multispectral imagery for individual tree species identification. Can J Remote Sens 38:125–138 Packalén P, Heinonen T, Pukkala T, Vauhkonen J, Maltamo M (2011) Dynamic treatment units in Eucalyptus plantation. For Sci 57:416–426 Pascual C, García-Abril A, García-Montero LG, Martín-Fernández S, Cohen WB (2008) Object-based semi-automatic approach for forest structure characterization using lidar data in heterogeneous Pinus sylvestris stands. For Ecol Manag 255:3677–3685 Patenaude G, Hill RA, Milne R, Gaveau DL, Briggs BBJ, Dawson TP (2004) Quantifying forest above ground carbon content using LiDAR remote sensing. Remote Sens Environ 93:368–380 Peura M, Gonzalez RS, Müller J, Heurich M, Vierling LA, Mönkkönen M, Bässler C (2016) Mapping a ‘cryptic kingdom’: performance of lidar derived environmental variables in modelling the occurrence of forest fungi. Remote Sens Environ 186:428–438 Popescu SC, Hauglin M (2014) Estimation of biomass components by airborne laser scanning. In: Maltamo M, Næsset E, Vauhkonen J (eds) Forestry applications of airborne laser scanning. Managing Forest ecosystems, vol 27. Springer, Dordrecht, pp 157–175 Pukkala T (2005) Metsikön tuottoarvon ennustemallit kivennäismaan männiköille, kuusikoille ja rauduskoivikoille (in Finnish for “prediction models for the expectation value of pine, spruce and birch stands on mineral soils”). Metsätieteen Aikakauskirja 3(2005):311–322 Pukkala T (2008) Integrating multiple services in the numerical analysis of landscape design. In: von Gadow K, Pukkala T (eds) Designing Green Landscapes. Managing Forest Ecosystems, vol 15. Springer, Dordrecht, pp 137–167 Pukkala T (2016) Which type of forest management provides most ecosystem services? Forest Ecosyst. https://doi.org/10.1186/s40663-016-0068-5 Pukkala T, Kangas J (1996) A method for integrating risk and attitude toward risk into forest planning. For Sci 42:198–205 Pukkala T, Kellomäki S, Mustonen E (1988) Prediction of the amenity of a tree stand. Scand J For Res 3:533–544 Pukkala T, Packalén P, Heinonen T (2014) Dynamic treatment units in forest management planning. In: Borges JG, Diaz-Balteiro L, McDill ME, Rodriguez LCE (eds) The Management of Industrial Forest Plantations. Managing Forest ecosystems, vol 33. Springer, Dordrecht, pp 373–392 Pukkala T, Sulkava R, Jaakkola L, Lähde E (2012) Relationships between economic profitability and habitat quality of Siberian jay in uneven-aged Norway spruce forest. For Ecol Manag 276:224–230 R Core Team (2016) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org/. Accessed 16 Oct 2017 Räsänen A, Lensu A, Tomppo E, Kuitunen M (2015) Comparing conservation value maps and mapping methods in a rural landscape in southern Finland. Landscape Online 44:1–19 Räty J, Vauhkonen J, Maltamo M, Tokola T (2016) On the potential to predetermine dominant tree species based on sparse-density airborne laser scanning data for improving subsequent predictions of species-specific timber volumes. Forest Ecosyst. https://doi.org/10.1186/s40663-016-0060-0 Vauhkonen Forest Ecosystems (2018) 5:24 Page 18 of 19 Rechsteiner C, Zellweger F, Gerber A, Breiner FT, Bollmann K (2017) Remotely sensed forest habitat structures improve regional species conservation. Remote Sens Ecol Conserv 3:247–258 Roces-Díaz JV, Burkhard B, Kruse M, Müller F, Díaz-Varela ER, Álvarez-Álvarez P (2017) Use of ecosystem information derived from forest thematic maps for spatial analysis of ecosystem services in northwestern Spain. Landscape Ecol Eng 13:45–57 Sani NA, Kafaky SB, Pukkala T, Mataji A (2016) Integrated use of GIS, remote sensing and multi-criteria decision analysis to assess ecological land suitability in multi-functional forestry. J For Res 27:1127–1135 Schröter M, Rusch GM, Barton DN, Blumentrath S, Nordén B (2014) Ecosystem services and opportunity costs shift spatial priorities for conserving forest biodiversity. PLoS One. https://doi.org/10.1371/journal.pone.0112557 Schulp CJE, Burkhard B, Maes J, Van Vliet J, Verburg PH (2014) Uncertainties in ecosystem service maps: a comparison on the European scale. PLoS One. https://doi.org/10.1371/journal.pone.0109643 Shapiro SS, Wilk MB (1965) An analysis of variance test for normality (complete samples). Biometrika 52:591–611 Simonson WD, Allen HD, Coomes DA (2014) Applications of airborne lidar for the assessment of animal species diversity. Methods Ecol Evol 5:719–729 Sverdrup-Thygeson A, Ørka HO, Gobakken T, Næsset E (2016) Can airborne laser scanning assist in mapping and monitoring natural forests? For Ecol Manag 369:116–125 Thompson SD, Nelson TA, Giesbrecht I, Frazer G, Saunders SC (2016) Data-driven regionalization of forested and non-forested ecosystems in coastal British Columbia with LiDAR and RapidEye imagery. Appl Geogr 69:35–50 Tomppo E, Haakana M, Katila M, Peräsaari J (2008a) Multi-source national forest inventory –methods and applications. Managing forest ecosystems, vol 18. Springer, Dordrecht Tomppo E, 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 Tomppo E, Katila M, Mäkisara K, Peräsaari J (2014) The multi-source national forest inventory of Finland - methods and results 2011. Working Papers of the Finnish Forest Research Institute, vol 319. http://www.metla.fi/julkaisut/ workingpapers/2014/mwp319.htm. Accessed 16 Oct 2017 Tomppo E, Olsson H, Ståhl G, Nilsson M, Hagner O, Katila M (2008b) Combining national forest inventory field plots and remote sensing data for forest databases. Remote Sens Environ 112:1982–1999 Valbuena R, Eerikäinen K, Packalen P, Maltamo M (2016a) Gini coefficient predictions from airborne lidar remote sensing display the effect of management intensity on forest structure. Ecol Indic 60:574–585 Valbuena R, Maltamo M, Martín-Fernández S, Packalen P, Pascual C, Nabuurs GJ (2013) Patterns of covariance between airborne laser scanning metrics and Lorenz curve descriptors of tree size inequality. Can J Remote Sens 39(sup1):S18–S31 Valbuena R, Maltamo M, Mehtätalo L, Packalen P (2017) Key structural features of boreal forests may be detected directly using L-moments from airborne lidar data. Remote Sens Environ 194:437–446 Valbuena R, Maltamo M, Packalen P (2016b) Classification of multilayered forest development classes from low-density national airborne lidar datasets. Forestry 89:392–401 Valbuena R, Vauhkonen J, Packalen P, Pitkänen J, Maltamo M (2014) Comparison of airborne laser scanning methods for estimating forest structure indicators based on Lorenz curves. ISPRS J Photogramm Remote Sens 95:23–33 Vauhkonen J, Imponen J (2016) Unsupervised classification of airborne laser scanning data to locate potential wildlife habitats for forest management planning. Forestry 89:350–363 Vauhkonen J, Packalen P, Malinen J, Pitkänen J, Maltamo M (2014) Airborne laser scanning based decision support for wood procurement planning. Scand J For Res 29(Suppl.1):132–143 Vauhkonen J, Ruotsalainen R (2017a) Assessing the provisioning potential of ecosystem services in a Scandinavian boreal forest: suitability and tradeoff analyses on grid-based wall-to-wall forest inventory data. For Ecol Manag 389:272–284 Vauhkonen J, Ruotsalainen R (2017b) Reconstructing forest canopy from the 3D triangulations of airborne laser scanning point data for the visualization and planning of forested landscapes. Ann For Sci 74:9. https://doi.org/10.1007/ s13595-016-0598-6 Vihervaara P, Auvinen AP, Mononen L, Torma M, Ahlroth P, Anttila S, Bottcher K, Forsius M, Heino J, Heliola J, Koskelainen M, Kuussaari M, Meissner K, Ojala O, Tuominen S, Viitasalo M, Virkkala R (2017) How essential biodiversity variables and remote sensing can help national biodiversity monitoring. Glob Ecol Conserv 10:43–59 Zimble DA, Evans DL, Carlson GC, Parker RC, Grado SC, Gerard PD (2003) Characterizing vertical forest structure using small-footprint airborne LiDAR. Remote Sens Environ 87:171–182 Zolkos SG, Goetz SJ, Dubayah R (2013) A meta-analysis of terrestrial aboveground biomass estimation using lidar remote sensing. Remote Sens Environ 128: 289–298 Vauhkonen Forest Ecosystems (2018) 5:24 Page 19 of 19