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

Forest Ecosystems, Jun 2018

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 m2, 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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

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

Vauhkonen Forest Ecosystems Predicting the provisioning potential of forest ecosystem services using airborne laser scanning data and forest resource maps Jari Vauhkonen 0 0 Natural Resources Institute Finland (Luke), Bioeconomy and Environment Unit , P.O. Box 68, Yliopistokatu 6, FI-80101 Joensuu , Finland 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 m2, 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. 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 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 m2 as in Naesset 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 (Naesset 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 ha bitat 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 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 the ALS data may be condensed to a few metrics (Kane et al. 2010; Leiterer et al. 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; Thompson et al. 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 m2. 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) : 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 ¼ nNij ; where vij is the normalized value and nij 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Þ ; where vij 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 ð1Þ ð2Þ 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. 1 or 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) aWhen computing the values for the present study, the following details or exceptions compared to original publications were made: 1The 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). 2Values 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) . 3The estimated carbon was obtained based on conversion factors from species-specific, total stem volumes to carbon contents. bTo 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. both the Eqs. 1 and 2 were 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. The models of Table 1 were 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. 2 resulted 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. 1 and 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. 1 were 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 thresholds of ground, shrub, and suppressed or dominant canopy, respectively. The following categories of the features were considered: obtained as S1 ∪ {xi × xj} ∀ i, j ∈ S1, 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: – Canopy height and density, which are the basic predictors used in ALS analyses (Naesset 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 (ccX_ground, where X is the echo category) or first echoes above the shrub layer threshold (ccshrub), 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 coniferous- or deciduous-dominated forest in the ALS data acquired during the leaf-off period (Liang et al. 2007). These features are denoted by diffx–y, where suffix x–y refers 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 propX/Y_z, where X/Y indicated the ratio of two echo categories X and Y, and z was 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 (runderstory); the standard deviation of the height values of echoes reflected above ground but below the dominant canopy threshold (stdunderstory); and the ratio of the echoes reflected from the shrub layer to all echoes (rshrub). Features xi, i = 1, 2, …, 143, listed above formed the initial set S1 of candidate predictors. To account for useful interactions between the features, the final set S was ð3Þ ð4Þ ^yin ¼ an þ XN n¼1bnxjcnn ; where ^yin is the vector of predicted priority values for the i:th ES, xjn is the j:th feature of S, and an, bn, and cn are model parameters at the n:th round of N = 1, 2, 3, 4 iteration rounds. Parameters an, bn, and cn were estimated using the nls function of R statistical computing environment (R Core Team 2016) . Testing every candidate feature as xjn at every iteration round, the RMSE between the predicted and reference priority values was computed as: RMSE ¼ sPffiffiffiffiffiffiinffi¼ffiffi1ffiffiðffiffi^yffiffiiffiffiffiffiffiffiffiffiffi2ffiffiffiffi −yiÞ ; n where n is the number of observations, and ^yi and yi 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. 1 were 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, 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 yi was the reference priority value of the i:th ES and ^yi its 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 ALS- and MS-NFI-based predictions were compared using the RMSE and coefficient of determination (R2) 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. 1 or 2 applied 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, s2), where E was the expected value and s2 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 X of 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 feature1 × feature2, where feature1 was often an absolute height value (a mean height, percentile, or height difference) and feature2 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 (propfirst/FP_ground). The predictors of BIOD (e.g., propfirst/FP_ground; diffonly–LP; hstdLP) 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 runderstory, diffonly–FP, propfirst/FP_ground, and hstdfirst. 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. 3 and 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 R2 = 0.11 to 0.81. The R2-improvements of the other ESs were close to this magnitude, except for BILB and COWB, which had R2 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. 3 and 4, right column). Decision analyses based on different priority functions, data, and model uncertainties Compared to the use of interval-scaled priority functions, those based on the ratio scale slightly increased the proportion of decisions that had a perfect agreement in both ALS and field data (Fig. 5, left panel). The same observation was made regarding the decisions based on the MS-NFI data, but the differences between the priority function forms were in general minor. When the ratio-scaled priority functions were used in the remaining analyses, the ALS and field data resulted to the same decision in 42% of the plots; the ALS-based ES was at least the second best according to field data in 69%; and among the three best alternatives in 84% of the plots. With MS-NFI data calibrated by the local field sample, the corresponding figures were 46%, 61%, and 72% and slightly lower without the calibration, the distribution of the decision scores of all data sources being shown in Fig. 5 (right panel). Thus, although the calibrated MS-NFI data did better in selecting ESs that matched perfectly with the field data, the proportion of poorest decisions was considerably lower based on the ALS data. The ALS-based models resulted to a better decision in 32%, those based on the MS-NFI data in 27%, and the decision equaled in 41% of the plots. Of the latter group of plots, around 2/3 had either BILB or COWB as the most important ES and the decisions for these plots were generally scored ≥4. Except for data sources, the decision score was highly dependent on the priority difference between the best ESs observed from the field plots: in the plots where the decisions equaled between the data sets, the average priority difference between the ESs prioritized as the first and second was 0.12 (standard deviation 0.08), whereas the corresponding figures for the other plots were 0.07 (0.07). Decision making based on the 5th percentile of the ALS-predicted priority value distributions reduced the proportion of worst decisions (Fig. 6, left panel). Otherwise, the decisions based on either the 5th or 95th percentile did not compare favorably with those based on the expected value in either data (Fig. 6). The confusion matrices for the most important ESs based on the different data sources and either the expected or extreme outcomes are shown as Appendices 2 and 3, respectively. A comparison of the confusion matrices based on either expected or extreme values indicates that many plots with the most important ES as those predicted with the highest or smallest error rates (e.g., CARB or COWB, respectively, based on ALS data) could be prioritized for this ES only by explicitly considering the extreme model outcomes. As a number of plots obtain a different prioritization based on the expected values, it was found interesting to look at the plots for which the prioritization changed depending on the model outcome. Using ALS-based ES proxy models, the prioritization based on the expected, worst, and best outcomes equaled in only 22 plots (in 43 using MS-NFI maps calibrated with the field data). If all decision alternatives based on the three outcomes of the ALS-models were considered, the alternative that matched perfectly with the field data was included in 57% of the plots; an alternative among the two best ones in 85%; and among the three best ones in 93% of the plots. With MS-NFI data, the corresponding figures were 59%, 76%, and 88%. Thus, especially the ALS-models were found useful in confining the decision alternatives to those most feasible according to production possibilities, from which the one preferred most preferred according to the stakeholder preferences could be selected based on further MCDA. Discussion Earlier RS-based decision analyses of multiple ESs have relied on map scales such as 1:25,000 (Roces-Díaz et al. 2017) or pixel sizes such as 500 × 500 m2 (Schröter et al. 2014). Also smaller pixel sizes of 60 × 60 m2 (Lehtomäki et al. 2015) or 48 × 48 m2 (Vauhkonen and Ruotsalainen 2017a) have been used in spatial prioritization analyses. Many of the aforementioned scales are coarse considering the need to formulate management prescriptions at the level of operational units (e.g., forest compartments), which are typically 1.5–2.0 ha in size in Finland (Koivuniemi and Korhonen 2006). Based on this background, the results of the present study are encouraging: the use of ALS in particular allowed deriving accurate predictions for plots of around 250 m2, i.e., in a considerably more detailed resolution than current operational compartments. The RMSEs of predicting the studied ES indicators by ALS were 47%–97% of the RMSEs of corresponding predictions based on benchmark forest resource maps derived using coarse satellite images. 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. In the sub-sections below, the results are discussed from the points of view of using ALS data to predict various ES indicator proxies (Table 1) and using these predictions in decision analyses compared to other data. On the ES proxies and their ALS-based modelling The expert models (Table 1) express the suitability of forests for the ESs as transformations of forest attributes. Applying the expert models yielded the highest biodiversity conservation values for mature, densely stocked forests. The values were weighted by site fertility such that most fertile sites received a considerable weight compared to poorer sites. An increasing basal area and mean diameter increased the soil expectation value of all the species, but otherwise the model included species-specific interactions between these and operational environment related parameters. The value of carbon storage increased according to an increasing stem volume depending on the species-specific biomass expansion factors. An increasing maturity (measured in terms of mean age and diameter) and decreasing stand density (measured in terms of either basal area or stem number) increased the suitability of a site for berry picking and visual amenity. The latter models also included species-specific terms such that especially the presence of pine trees improved the values. Thus, the response variables of all other ESs except CARB were modelled as functions of multiple forest attributes. The models based on the ALS data (Table 2) explained the differences in the provisioning potential of the different ESs with RMSEs of 15%–30%. The selected features and model performances are well in line with the background given in the previous paragraph. The predictions of CARB had the highest accuracies, because those essentially explained the variation in the total stem volume converted to carbon using biomass expansion factors. The ratio of first-of-many to all first echoes (propfirst/ FP_ground) was used as a predictor of CARB in addition to height and density metrics, which are typical to total biomass or carbon models. The aforementioned feature has been found useful for separating species (Ørka et al. 2012; Vauhkonen et al. 2014) and should be studied in a broader forest modeling context. Although the field proxy value for TIMB was computed using species-specific predictors, the model based on the ALS data included only features based on height. The models for the other ESs included ALS-based predictors that were clearly related to forest canopy structure: for example, features indicating the existence or abundance of low vegetation were frequently selected to the models of BIOD or recreation-related ESs, respectively. No clear recommendations for the selection of features can be given, except for using multiple heights and echo types when computing the candidate features. The requirement to estimate a high number of model parameters was likely reduced by including products of all candidate features to account for interactions between them, which is a useful property that has not been reported in earlier ALS studies. The features used here are the most common and easily implementable using available software packages such as R. However, it is acknowledged that not all features presented in the literature were included and it could be possible to improve the results with more experimental ones such as those related to volumetric (Vauhkonen and Ruotsalainen 2017b) or textural (Niemi and Vauhkonen 2016) properties. The proxies listed in Table 1 are measured in different units. In order to use the proxies in decision analyses, those need to be normalized to the same scale, which was done in this paper prior to the modeling. Due to the transformations, the prediction accuracies cannot be easily compared with earlier studies, even in relative scale. Even CARB, which is a species-specific transformation of the total volume, cannot be directly compared to previous studies that predict total biomass or carbon due to the Box-Cox-transformation (Appendix 1) applied to equalize the distribution of the response variables for the decision analyses. If a similar model for CARB had been fit without the Box-Cox-transformation, the RMSE would have risen to around 34%. It is higher than (e.g.) the RMSE of 25.7% obtained by Kankare et al. (2015) for total biomass in the same region, but when comparing the figures the differences in the definition of the response variable and a slightly different plot size must also be considered. For more reasoning on the need for the transformations, please see the next section. On the other hand, the modeling task considered above may have been alleviated by the fact that all response variables resulted from models that had been formulated earlier using common forest mensurational attributes. Actual differences between berry yields of two forests could be much more discrete than those predicted as a function of forest attributes. The performance of the models should thus be tested by re-fitting or validating the models against ES-specific indicators that are not based on models but direct observations made in the field (see also Hegetschweiler et al. 2017; Kohler et al. 2017) . On the other hand, it may be less feasible or even impossible to predict certain ES properties otherwise than using features quantifying the three-dimensional forest structure. Examples in the Scandinavian boreal forest structures include mammal or bird habitats (e.g., Melin et al. 2013, 2016) , which could have realistically occurred in the studied area, but could not be included in the present comparison due to the lack of models based on information obtainable from forest resource maps or field plots. The value of ALS data can be seen to lie especially in applications that require identifying sites with high value for e.g. conservation as ALS data coverage is globally increasing and the method has been proven to provide 3D descriptions of vegetation structure that, in turn, have been long known as primary determinants of habitat quality and biodiversity (MacArthur and MacArthur 1961; Dueser and Shugart Jr, 1978; Brokaw and Lent, 1999) . Nevertheless, Vihervaara et al. (2017) identified only four Essential Biodiversity Variables that could benefit from the use of RS data. According to this study, the obtainable improvements are most likely indicator-specific, with magnitude depending on what RS data are available. Other aspects of RS-based decision analyses: Data, normalization, and uncertainties This study tested normalization resulting to values in either an interval- (Eq. 1) or ratio-scale (Eq. 2). Even though the priority value functions based on the different scales did not differ considerably, those based on the ratio-scale performed slightly better in the priority ranking, which is logical as also the ALS-based features provide information at a ratio scale. Although the purpose is not to exhaustively discuss the implications of different value function forms, one important practical aspect discovered in the exploratory analysis could be brought up: An incorrect selection of a value function form would result in biased decisions by weighting the rank-orderings of the alternatives to an undesired direction, which is exemplified by a comparison of the transformed and non-transformed value function forms (Appendix 1). This discussion highlights the importance of considering data normalization procedures in detail already at the modeling stage, if the final applications aim at decision analyses where all modeled proxies should be measurable at the same scale. The ALS and MS-NFI data sources showed several differences, when predicting the priority value functions. Except having a more deterministic relationship with the forest attributes, one additional improvement of the locally fit ALS models is the ability to predict the 0–1-value ranges directly. When using forest attributes such as those predicted by the MS-NFI, the maximum values used in the normalization should be representative of the entire area. In this study, the minima and maxima predictions based on the MS-NFI were more severely incorrect than with ALS data (Figs. 3–4), which partially explains the poorer performance of this data source. Many multivariate predictions from RS data are based on non-parametric nearest neighbor methods, which could have been considered also in the case of ALS data. However, the aforementioned problem would be seriously present also in those types of predictions. Regardless of the input data or applied scale or mapping technique, it is well recognized that the resulting ES maps will include uncertainties (Eigenbrod et al. 2010; Schulp et al. 2014; Räsänen et al. 2015) . A few earlier spatial prioritization studies (Lehtomäki et al. 2015; Räsänen et al. 2015; Vauhkonen and Ruotsalainen 2017a) used forest attribute maps in a similar resolution as considered here, but without a calibration based on field data. Such a practice cannot be recommended due to the observed uncertainties in the uncalibrated MS-NFI data in particular. However, it should be acknowledged that all analyses carried out in the aforementioned studies might not be sensitive to incorrect pixel-level prioritization decisions. Also, each of the aforementioned studies used either aggregated pixels to somewhat account for these effects. It is not possible to address the uncertainty reduction due to the use of aggregated resolution based on the field data of this study, which was collected from fixed-size plots. Finally, it should be noted that in the actual decision making, the stakeholder preferences may affect or even dictate the allocation of the ESs over suitability of forest structure. Even if the decision maker had no preferences on the ESs, aggregating individual pixels, i.e., deviating from their local optima to compose larger treatment units, could be feasible with respect to the implementations of management prescriptions or achieving ecological or economic objectives that are determined over a larger area (Pukkala et al. 2014) . Importantly, the decisions based on the same data may differ for a risk-avoiding, risk-neutral or risk-seeking decision maker (Pukkala and Kangas 1996) . The risk preferences of the decision maker should clearly be incorporated in the decision making based on the ES maps. Here, a similar technique as in Pukkala and Kangas (1996) was used to account for the uncertainties emerging from different accuracies of the ES proxy models. The results reported in the last paragraph of the Results section can be concluded such that separate forest ES maps should be prepared using the expected, worst, and best outcomes of the model predictions to fully describe the production possibilities of the landscape under the uncertainties in the models. Also those related to the decision makers’ preferences could be accounted for as additional stochastic distributions in the analyses (e.g., Kangas et al. 2007) and incorporated in analyses as pixel-level production constraints. Overall, the prioritization approach should be tested further involving real decision makers. Conclusions The provisioning potential of ESs (biodiversity, timber, carbon, berries, and visual amenity) was modeled as expert model-based proxies that express the suitability of forests for the ESs as transformations of forest attributes. The models based on the ALS data explained the variation in these proxies with RMSEs of 15%–30%. The RMSEs of the ALS-based models were 47%–97% of the RMSEs of corresponding predictions based on the MS-NFI forest resource maps. 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. The RMSE-differences did not fully translate to the accuracies of land use decisions: instead, prioritizing the land use for the ESs with the highest provisioning potential could be done with rather similar accuracies based on both data sources and at a resolution of 250 m2, i.e., in a considerably more detailed scale than current operational forest management units. ALS-based models for the ES proxies can however be recommended based on their better stability regarding the model errors. The results suggest that separate forest ES maps should be prepared using the expected, worst, and best outcomes of the model predictions to fully describe the production possibilities of forest under the uncertainties in the models. Additional file Additional file 1: Empirical cumulative density and probability density function forms and Pearson correlation coefficients of different transformations of the expert model predictions for the ESs considered. (DOCX 1758 kb) Appendix 1 Box-Cox transformation of the reference priority values for the ESs An exploratory analysis revealed a practical comparability issue related to using the expert function values for the ESs directly in Eq. 2, which can be demonstrated by visualizing the empirical cumulative distribution functions and histograms of the data (Additional file 1: Figure S1-Figure S6). The expert functions of Table 1 were fit with many different data sets having different value ranges and resulting to unequal frequencies for the different ESs, when applied in the present data. Using Eq. 2, these differences would have translated directly to the priority values, which was found problematic with respect to their ranking. Look especially at the left-hand columns of Additional file 1: Figure S1-Figure S6 and compare Additional file 1: Figure S6 to the others: the direct use of these values would have resulted to a very high number of plots being prioritized as AMEN only because its distribution more frequently included higher values compared to the distributions of the other ESs, which were more often skewed to the left. For this reason, the expert model values of every ES were transformed to produce frequencies that followed the normal distribution as closely as possible prior to converting the value ranges between 0 and 1. The transformation was obtained as (Box and Cox 1964) : yðλÞ ij ¼ 8 > < yiλj−1 λ >: ln yij ; if λ ¼ 0; ; if λ≠0; ð5Þ where yij is the original value of the i:th observation of the j:th ES proxy and λ is a parameter. The value of λ was selected from an interval of − 3 to 3 as the value that maximized a test statistic on whether the considered sample came from a normally distributed population (Shapiro and Wilk 1965) . This transformation did not affect the order of the observations, but produced approximately equally shaped frequency distributions of every ES. As a result, the ESs could be prioritized using two alternative priority value function forms illustrated in the two rightmost columns of Additional file 1: Figure S1-Figure S6: either according to the interval scale (Eq. 1) based only on the order of the observations or according to the ratio scale (Eq. 2) preserving the ratios between the observations, but thanks to the Box-Cox-transformation, having approximately equal (close to normal) frequency distributions between the ESs. The function forms and frequency distributions of the priority values are shown in Additional file 1: Figure S1-Figure S6 and the correlations between the priority values in Additional file 1: Table S1. During the feature selection process, it was noted that exactly the same features would have been selected to the models, regardless of whether the original or Box-Cox-transformed response variables were used. Therefore, although the Box-Cox-transformation has been rarely used in ALS-based studies, it could potentially aid also other applications by normalizing the response variable to a desired form in the model fitting step. 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. 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: 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 , Naesset 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 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 . 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 , Naesset 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 , Naesset E , Vauhkonen J (eds) Forestry applications of airborne laser Scanning.Managing Forest ecosystems , vol 27 . Springer, Dordrecht, pp 357 - 374 Naesset 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 Naesset E , Gobakken T ( 2008 ) Estimation of above- and 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 . 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 . Ørka HO , Dalponte M , Gobakken T , Naesset 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 , Naesset 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 , Naesset 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 . 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. 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 . 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 . 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 . 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 , Naesset 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 . 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 . 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

This is a preview of a remote PDF:

Jari Vauhkonen. Predicting the provisioning potential of forest ecosystem services using airborne laser scanning data and forest resource maps, Forest Ecosystems, 2018, 24, DOI: 10.1186/s40663-018-0143-1