Jukuri, open repository of the Natural Resources Institute Finland (Luke) All material supplied via Jukuri is protected by copyright and other intellectual property rights. Duplication or sale, in electronic or print form, of any part of the repository collections is prohibited. Making electronic or print copies of the material is permitted only for your own personal use or for educational purposes. For other purposes, this article may be used in accordance with the publisher’s terms. There may be differences between this version and the publisher’s version. You are advised to cite the publisher’s version. This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Parvez Rana & Jari Vauhkonen Title: Stochastic multicriteria acceptability analysis as a forest management priority mapping approach based on airborne laser scanning and field inventory data Year: 2023 Version: Published version Copyright: The Author(s) 2023 Rights: CC BY 4.0 Rights url: http://creativecommons.org/licenses/by/4.0/ Please cite the original version: Rana P., Vauhkonen J. (2023). Stochastic multicriteria acceptability analysis as a forest management priority mapping approach based on airborne laser scanning and field inventory data. Landscape and Urban Planning 230. 104637. https://doi.org/10.1016/j.landurbplan.2022.104637. Landscape and Urban Planning 230 (2023) 104637 Available online 16 November 20220169-2046/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Research Paper Stochastic multicriteria acceptability analysis as a forest management priority mapping approach based on airborne laser scanning and field inventory data Parvez Rana a,b, Jari Vauhkonen b,c,* a Natural Resources Institute Finland (Luke), Paavo Havaksen tie 3, 90570 Oulu, Finland b University of Eastern Finland, School of Forest Sciences, P.O. Box 111, FI-80101 Joensuu, Finland c University of Helsinki, Department of Forest Sciences, Latokartanonkaari 7 (P.O. Box 27), FI-00014 Helsinki, Finland H I G H L I G H T S  First spatially explicit Stochastic Multicriteria Acceptability Analysis (SMAA).  Potential management alternatives quantified by ALS-based pixel proxies.  A nearest neighbor approach enabled pixel-level SMAA.  SMAA estimated the most acceptable management and the strength of this decision.  Decisions differing from the preferences can identify hotspots for forest structure. A R T I C L E I N F O Keywords: Forest planning Decision-making Remote sensing First-rank acceptability Weight space A B S T R A C T The mapping of ecosystem service (ES) provisioning often lacks decision-makers’ preferences on the ESs pro- vided. Analyzing the related uncertainties can be computationally demanding for a landscape tessellated to a large number of spatial units such as pixels. We propose stochastic multicriteria acceptability analyses to incorporate (unknown or only partially known) decision-makers’ preferences into the spatial forest management prioritization in a Scandinavian boreal forest landscape. The potential of the landscape for the management alternatives was quantified by airborne laser scanning based proxies. A nearest-neighbor imputation method was applied to provide each pixel with stochastic acceptabilities on the alternatives based on decision-makers’ preferences sampled from a probability distribution. We showed that this workflow could be used to derive two types of maps for forest use prioritization: one showing the alternative that a decision-maker with given pref- erences should choose and another showing areas where the suitability of the forest structure suggested different alternative than the preferences. We discuss the potential of the latter approach for mapping management hotspots. The stochastic approach allows estimating the strength of the decision with respect to the uncertainty in both the proxy values and preferences. The nearest neighbor imputation of stochastic acceptabilities is a computationally feasible way to improve decisions based on ES proxy maps by accounting for uncertainties, although the need for such detailed information at the pixel level should be separately assessed. 1. Introduction The Ecosystem Service (ES) concept is widely used to communicate the connection of providing these services for human wellbeing and it has become subject to governance policies and practices to support sustainable provision of ESs (Mann et al., 2022). Many established decision support tools enable an understanding of the impacts of man- agement and land-use changes on ESs in particular in terrestrial land-use policy sectors such as agriculture and forestry (Gre^t-Regamey et al., 2017). Nevertheless, considering ESs in forest management planning is not trivial (see, e.g., Knoke et al., 2021). Generally speaking, it is possible to focus on (1) data on the current provision of ESs, (2) models * Corresponding author. E-mail addresses: parvez.rana@luke.fi (P. Rana), jari.vauhkonen@uef.fi (J. Vauhkonen). Contents lists available at ScienceDirect Landscape and Urban Planning journal homepage: www.elsevier.com/locate/landurbplan https://doi.org/10.1016/j.landurbplan.2022.104637 Received 30 January 2022; Received in revised form 7 November 2022; Accepted 8 November 2022 Landscape and Urban Planning 230 (2023) 104637 2 expressing the change from current levels, or (3) preferences or prior- ities for ES management (Nemec and Raudsepp-Hearne, 2013). We propose Stochastic Multi-Criteria Acceptability Analysis (SMAA; Lah- delma et al., 1998; Lahdelma and Salminen, 2001) as an approach to consider the preferences in spatially-explicit management prioritiza- tions of a Scandinavian forested landscape. In this context, ES-related proxy values are derived by data and models of earlier studies (Vauh- konen and Ruotsalainen, 2017; Vauhkonen, 2018) and used as indirect measures of achieving management objectives, as explained below. We do not analyze the spatial balance between supply and demand of the ESs. Planning a forest landscape involves allocating management prac- tices or silvicultural treatments optimally with respect to the specified management objectives. These choices determine the provision of ESs (Pukkala, 2016), as intensive orientation on timber extraction may result to negative impacts on non-timber ESs and biodiversity (Pohjan- mies, 2018; Pohjanmies et al., 2021). Nevertheless, optimized man- agement can balance between these objectives and produce feasible solutions (Pukkala, 2016), especially if large enough areas are consid- ered (Pohjanmies, 2018). In any scale, the production possibilities of a forest are more efficiently used when making forest management de- cisions at the level of individual map units such as pixels instead of predefined management units such as forest stands (Heinonen et al., 2007). On these grounds, Vauhkonen and Ruotsalainen (2017) proposed prioritizing the landscape tessellated to pixels for ESs based on man- agement. For example, managing areas prioritized to recreation by means of selective cuttings can both improve recreation and other values and produce limited amounts of timber (Silvennoinen et al., 2002; Miina et al., 2016) but, due to the management choice, provide all ESs differently than areas prioritized for timber production or set aside. Overall, it is acknowledged that same areas provide more than just one ES, but the planned management segregates the main ESs that can be provided in the near future in intensively (timber > other ESs) or extensively managed (timber < recreation and other ESs) or set aside areas (no timber production, but higher levels of ESs conflicting with timber production). The management prioritization depends on how forest owners evaluate ESs provided by their forests. If stakeholder preferences and production possibilities of the planning area are known, generic multi- criteria decision analysis (MCDA) or optimization methods can be applied to determine optimal management based on discrete or continuous problem setting, respectively, as reviewed by Pukkala (2008) and Uhde et al. (2015) in forest planning contexts, and Lange- meyer et al. (2016) regarding ESs more generally in land-use planning. In Northern Europe, >70 % of the forest area is privately owned (Weiss and Zivojinovic, 2020), and hundreds of thousands of forest owners make management choices independently. The owners’ management objectives related to ESs have been found challenging to elicitate, and overall require more research (Weiss et al., 2019; Juutinen et al. 2021; Nieminen et al., 2021). Even if exact preferences were unknown, at least relevant criteria or (complete or incomplete) rank order of the criteria can often be assumed (Kangas et al., 2015). Some MCDA tools can be applied with such uncertain, incomplete, or even inaccurate informa- tion. For example, Stochastic Multicriteria Acceptability Analysis (SMAA; Lahdelma et al., 1998; Lahdelma and Salminen, 2001) allowed computing the probability of a certain alternative obtaining a given rank, therefore supporting the selection of the most recommendable forest management plan for unknown or partially known preferences (Kangas et al., 2003; but see Kangas, 2006). In spatial MCDA or optimization integrated to a Geographic Infor- mation System (GIS) framework, 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 (Malc- zewski and Rinner, 2015). Maps indicating the provisioning of forest ESs have been produced from remotely sensed forest inventory, structure, and habitat layers (e.g., Vauhkonen, 2018). Such maps may enable to spatially identify most important areas with respect to ESs (termed ’hotspots’ as in Egoh et al., 2008; Schroter and Remme, 2016; Kangas et al., 2018), thus providing information for forest (Vauhkonen and Ruotsalainen, 2017) or other land-use prioritization (Koschke et al., 2012; Verhagen et al., 2017; Honeck et al., 2020) and subsequent de- cision analyses (Caglayan et al., 2021; Forsius et al., 2021). However, applications of the ES maps often lack aspects applied routinely in multi- attribute forest planning with conventional field inventory data (cf., Pasalodos-Tato et al., 2013; Uhde et al., 2015). First, the applications have mainly related to mapping the present state and not resulting in instructions of forest management that enhance the provisioning of the ESs preferred by the decision-makers. Second, even the inventory ap- plications are highly focused; most recently, carbon-related ESs and biodiversity, while cultural services like recreation received less atten- tion (Knoke et al., 2021; see also Gre^t-Regamey et al., 2017, for similar notes on more general environmental decision making). Third, various uncertainties should be better addressed, particularly in spatially explicit analyses of ES maps aimed at decision making (Boerema et al., 2017; Englund et al., 2017; Barton et al., 2018; Kangas et al., 2018). Taken together, the mapping of ES provisioning often lacks analyzing decision makers’ preferences on the ESs provided, although it can affect or even dictate the allocation of the ESs over the landscape. Many MCDA methods require preferences as numerical weights, reflecting acceptable trade-offs between decision criteria, but for many applications it is infeasible to assume that these values were strictly known. A number of generic methodological possibilities also exist for integrating related uncertainties into MCDA to improve the decision proposals (Kangas et al., 2015; Malczewski and Rinner, 2015). For instance, fuzzy and varying preferences have been taken into account attempting to find alternatives that outperform or dominate (e.g., Yatsalo et al., 2015). ES assessments incorporating uncertainty have used such outranking methods or deterministic MCDA methods augmented by sensitivity an- alyses (Bryant et al., 2018). It is also possible to map locations of robust and more sensitive parameters (Ligmann-Zielinska and Jankowski, 2014). To concretize the problem setting, consider a utility function U ˆ P wi  ui(ci), where w is the weight for criterion i and ui(ci) a (sub-utility) function that transforms the criterion values ci to utility-scale, which can be used to analytically evaluate and compare qualitative, quantitative and, e.g., spatial criteria (Kangas, 1993; Pukkala and Kangas, 1993; Store and Kangas, 2001), also for other than forest-related ESs (Musta- joki et al., 2020; Marttunen et al., 2021). Both wi and ui(ci) can be assumed to be stochastic, and that way used to account for uncertainties (Alho and Kangas, 1997; Kangas et al., 2000). If probability distributions for the uncertain data can be obtained, a preferable approach is to include those in the analyses (Kangas et al., 2007; Convertino et al., 2013). Another way to deal with this uncertainty is to vary weights wi to find out the distribution of decisions and determine the most acceptable alternative under different preferences (Lahdelma and Salminen, 2001). This SMAA methodology is reasoned oppositely to methods ranking alternatives with known wi: It explores the feasible weight space by repeated random sampling of wi to analyze the ranks of alternatives or derive weights that make an alternative the most preferred one. The underlying rank acceptability indices can be used to infer, not only the best (or worst) alternative, but also the strength of this prescription in terms of the volume of the affected weight space. Few earlier studies that developed methods to account for un- certainties in expert maps have considered wi stochastically at the level of pixels for the mapping of robustness and sensitivity of the weights. Ligmann-Zielinska and Jankowski (2014) employed Monte Carlo simu- lation to sample from distributions of criteria weights to produce maps showing both the averages and standard deviations of uncertain criteria. Areas with a high average value and either low or high deviation could then be prioritized with consensus or reservation, respectively; latter also as complementary areas for further analyses as in García Marquez et al. (2017). Vauhkonen (2018) suggested a possibility to consider the P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 3 predictive distributions of ES-related proxies instead of just expected values to produce maps suited for decision-makers with different risk preferences. While possibly enhancing the information content of the maps, depending on implementation, the method may require sampling distributions for each pixel, which can be computationally demanding for large areas. The technique proposed by Ligmann-Zielinska and Jankowski (2014) also involved computations in order of 109 for the whole study site composed of 73,170 pixels. To be operational, the pixel- level uncertainty quantifications seem to need either spatial aggregation or numerical approximation approach to circumvent the need for supercomputing. We propose a new approach for forest management prioritization based on uncertain ES-related proxies at the pixel level over a landscape. The novelty is to make the detailed uncertainty quantification using SMAA for a set of sample plots and impute the results to the map pixels to make the analysis computationally feasible. The decision-makers’ management preferences were totally or partially unknown and sampled from a probability distribution. We focused on 1) estimating the strength of the decisions on forest use prioritization with respect to the uncer- tainty in the preferences and 2) identifying areas where the suitability of the forest structure suggested different decisions than the preferences. 2. Methods 2.1. Study area and experimental data Our study area is located in Evo, southern Finland (61.19N, 25.11E), which belongs to the southern boreal forest zone. The area encompasses managed to semi-natural and natural forests regarding their silvicultural history. Coniferous trees (Scots pine Pinus sylvestris and Norway spruce Picea abies) contribute > 80 % of the total growing stock in the study area, with the remaining contribution from deciduous trees (aspen Populus tremula, alders Alnus spp., birches Betula spp., rowan Sorbus aucuparia and willows Salix spp). The state-owned forest is managed with multiple objectives, including timber production, biodi- versity conservation, recreation, and carbon sequestration. Airborne laser scanning (ALS) and field sample plot data were gathered from previous studies in the same area (see, Vauhkonen, 2018). The ALS data were acquired on May 7, 2012, using a Leica ALS50 scanner. The ALS flight was operated at an altitude of 2,200 m above ground level, which yielded a nominal ALS pulse density of 0.8 m 2. These data were pre-processed and ground-classified to an open data product made available by the National Land Survey of Finland (2015). The data were downloaded as tiles of 3  3 km2, which were further quartered to facilitate the metric extraction. We demonstrate our empirical results with data retained for those quarter-tiles that con- tained the field plots (see below). We removed data from within 20 m of all non-forest areas such as lakes, buildings, and roads using the same mask as Vauhkonen and Imponen (2016) based on an existing forest management planning map. As a result, our study area consisted of altogether 1,750 ha of forest land. The field sampling locations were determined by k-means clustering the area in terms of the Euclidean distance of five ALS metrics and selecting representative areas from each cluster, which was found to efficiently stratify the area according to the forest structural variation (Vauhkonen and Imponen 2016). A total of 102 circular field plots (9 m radius) also studied by Vauhkonen (2018) were considered as the field reference data. In an inventory campaign that took place from June to August 2014, the diameter-at-breast height (DBH) and species were observed for each tree with DBH  5 cm. The median tree of each plot was additionally measured for height. These observations were used in standard inventory equations and methods to derive a set of plot-level attributes as described by Vauhkonen (2018). The plot data are sum- marized in Table 1. We used the derived plot-level attributes as inputs to expert models that predicted numeric scores for the management prioritization as a function of the forest mensurational attributes of the plots. We consid- ered the following six indicators or proxies, which are the same as used by Vauhkonen (2018) and explained in more detail in that publication: 1. Potential of the forest site for picking bilberry and cowberry, 2. obtained as a unitless index-value using the expert models of Iha- lainen et al. (2002). 3. Visual amenity as a unitless index-value based on the expert model of Pukkala et al. (1988). 4. Biodiversity as a unitless index-value, computed as basal area- weighted mean diameter  growing stock volume, scaled using dominant species-specific sigmoidal transformation functions and multiplied by the site fertility specific weights. The expert function of Lehtomaki et al. (2015) was applied as in Vauhkonen (2018). 5. Carbon storage, in which the total stem volume was transformed to carbon (t/ha) using species-specific conversion factors of Karjalainen and Kellomaki (1996). 6. Timber production potential as the soil expectation value (SEV, €/ha) predicted by the models of Pukkala (2005). The SEV is the discounted sum of the projected costs and profits from even-aged forest rotations, starting from bare land and continued in perpetu- ity. The model of Pukkala (2005) uses descriptive statistics of the growing stock (mean diameter, basal area, number of stems, and age) and operational environment (price, interest rate, and temperature) to predict an average SEV of a high number of rotations. We used 1300days as the temperature sum and computed our SEV as an average of applying interest rates of 1–4 % and the same saw-wood/ pulpwood price combinations as in the model fitting (Pukkala, 2005). Many other models for the proxy values above could have been considered (cf., Miina et al., 2020). Nevertheless, the models for berry yields, visual amenity and carbon mentioned above were evaluated against other models by Turtiainen (2015), Silvennoinen (2017), and Lehtonen et al. (2004), respectively, and these studies do not suggest that alternative models would include other predictors or otherwise better account for the modeled phenomena. For biological diversity, proxy values based on forest mensurational data are required in the absence of detailed plant mapping data and we are not aware of alter- native functions to Lehtomaki et al. (2015) compatible to our data. Finally, the SEV model of Pukkala (2005) is used in many studies to approximate the net present value from even-aged rotations repeated to infinity. Using the expert models above produced a different value range for every proxy that was normalized for joint priority ranking. We trans- formed each indicator to vary between 0 and 1 similar to Vauhkonen (2018): we first applied a Box-Cox-transformation to force the indicators to comply with normally distributed values as close as possible, and then scaled the range of the values between 0 and 1 using the well-known min–max scaling. We extracted the ALS points for the 102 field reference plots and 452,727 pixels of 16  16 m2 overlaid to the area as a regular grid with a cell size selected to comply with the size of the field plot. Using the point data from each of these computation units and different threshold values Table 1 Characteristics of 102 field plots used in the analysis. STD Standard deviation. Attributes Mean  STD Range Total volume, m3/ha 200  109 38–673 Pine 107  90 0–359 Spruce 61  93 0–464 Deciduous 31  44 0–228 Total carbon, ton/ha 67  39 12–245 Age, years 66  37 11–186 Diameter, cm 23  7 9–43 Dominant height, m 21  5 12–36 P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 4 for height to distinguish ground, shrub, and canopy layers, we calculated a total of 142 features describing the distributions of first, last, first of many, last of many, and only echoes of the ALS pulses (Table 2). 2.2. Stochastic multicriteria acceptability analysis (SMAA) Our problem was to select the prioritized management for every forest plot. Each plot had n ˆ 6 alternatives, represented for the i:th plot by a row vector of the proxy values xi ˆ [uj(xi1), …, uj(xin)], where uj(.) are the Box-Cox-transformation functions described above to map the proxies j ˆ 1, …, n into the range [0,1]. The priority rank of the l:th alternative was determined as an integer from the best (1) to the worst (6) rank as rank(l, xi, w) ˆ 1 ‡ Pn jˆ1ρ wj  xij > wl  xil  ,where the weight vector w represented the decision maker’s preferences on the alterna- tives and ρ(true) ˆ 1 and ρ(false) ˆ 0. The SMAA method (Lahdelma et al., 1998; Lahdelma and Salminen, 2001) explores the feasible weight space W to discover preferences that would give the specific alternative a certain rank. Realizations of w are sampled among the set of feasible weights such that w 2 W ˆ  w 2 Rn wj  0 and Pn jˆ1wj ˆ 1 o . In addition, the preferences can be assumed to be unknown or partially known. If the preferences are totally unknown, the weight vectors w are sampled from a uniform distribution. A complete or partially known importance order is modeled by omitting the weights that do not meet inequality constraints such as wj 1 > wj 2 ? wj 3 ? …? wjn, where wj 1 is the weight for the most preferred alternative and ? denotes undetermined importance order among the other alternatives (Kangas et al., 2015). In SMAA, rank acceptability is then defined as the expected volume of feasible weights that result in rank r for an alternative proportional to the total volume of the weight space. In practice, this score can be computed by recording the ranks of the alternatives resulting from N ˆ 10,000 draws of the weights from the random distribution (Section 2.3) and dividing the number of times the ranks were obtained by the total number of draws. The obtained values range from 0 to 1 and indicate the proportion of draws the alternative obtained the given rank. 2.3. Nearest neighbor imputation of the SMAA-based rank acceptabilities to the map As motivated in the Introduction, stochastic analyses as described above would become unfeasible if computed at the level of each pixel in big pixel data. Our solution to make the analyses computationally feasible was to run the SMAA at the level of the reference field plots to produce the rank acceptability of the prioritized management (Fig. 1b). We searched for nearest neighbor plots for each pixel and used the nearest neighbor imputation (Fig. 1a) to populate each pixel with the SMAA result. The set of 10,000 random weights to be used in the SMAA of each reference plot was generated as follows. Given the set of n alternatives to be compared, simulate n-1 uniformly distributed random numbers that fell within the range [0,1]. Arrange the random numbers u(1), …, u(n-1) in a descending order, preceded by 1 and succeeded by 0, i.e., u(0) ˆ 1 and u(n) ˆ 0. Compute the weight vector w by concatenating wj ˆ u(j) – u(j ‡ 1), j ˆ 1, …, n. This way, it was possible to obtain weights ranging from low to high for any alternative (Kangas et al., 2015). We stored the plot-specific SMAA result to matrices of n  n, where the matrix elements indicated the number each alternative resulted to the given rank among the random weights. We produced these results assuming either neutral decision-makers’ preferences or that there was one most preferred alternative, but no other information about the importance order. These situations corresponded to sampling weights for the unknown vs complete or partially known importance orders, as described in Section 2.2., and resulted to different types of prioritiza- tions of the map locations as shown in the Results. Using the yaImpute package of the R statistical environment (Crookston and Finley, 2008), we searched for nearest neighbors be- tween the reference plots and grid cells in terms of the Euclidean dis- tance of the 142 ALS features. First, we imputed the SMAA matrix to each of the grid cells from a single nearest neighbor plot. Second, we tested the possibility to employ matrices from more than one neighbor and demonstrated our results with three neighbors. In the latter case, the matrices were obtained as inverse distance weighted, i.e. each element of the SMAA matrix of the m:th nearest neighbor plot for a grid cell was multiplied by the term vm ˆ  1 dm  = X 1 d  where d is the vector of distances to the k nearest neighbors, before summing up the k matrices for the grid cell. As indicated by Fig. 1(c-d), there are many alternative ways to interpret the rank acceptabilities (see also Kangas et al., 2015). We only used the first-rank acceptability values to select the alternative with the highest acceptability and determine the strength of this decision, but we discussed other possibilities. We visualized the results as maps and numerically compared the distributions in the reference and imputed data. 3. Results 3.1. Plot-level prioritization When the decision-makers’ preferences were assumed neutral, 39 % of the reference field plots were prioritized for cowberry, followed by timber (19 %), bilberry (18 %), biodiversity (16), carbon (8 %), and scenic beauty (1 %) (Table 3). When the decision-makers were assumed Table 2 ALS features computed. Amount-column gives the total number of features in each predictor category as (number of features resulting from the described computation)  (number of ALS echo categories to which the computation was applied). The reader is referred to Vauhkonen (2018) for more details of the features. Abbreviation Description Amount Hmax Maximum height of echoes above 0.5 m. 1  5 Hmean Mean height of echoes above 0.5 m. 1  5 Hstd Standard deviation of echoes heights above 0.5 m. 1  5 H5–95 Height of the echoes at the 5 %, 20 %, 30 % … 95 % percentile above 0.5 m. 11  5 D5–95 Proportional density of the echoes at the 5 %, 20 %, 30 % … 95 % percentile above 0.5 m. 11  5 VegRat The proportion of all ALS echoes above 0.5 m used to describe vegetation cover as in Korhonen et al. (2008). 1  5 VegRatshrub The proportion of first echoes above 5 m and the proportion of all echoes above 0.5 m but below 5 m, used to describe the shrub layer (Vauhkonen and Imponen, 2016). 2 VegRatunderstory The proportion of all ALS echoes returned above 0.5 m but below the height of the 60th percentile and its standard deviation used to describe intermediate to understory vegetation (Vauhkonen and Imponen, 2016). 2 DIFFx-y The absolute difference between mean heights of the two echo categories as Diffx–y where x–y denotes the height difference of two echo categories, computed without a height threshold and used to discriminate coniferous or deciduous dominated forests (Liang et al., 2007). 4 PROPx/y The ratio of the number of echoes in different categories used to describe differences in species and size properties as in Ørka et al. (2012) and Vauhkonen et al. (2014). 4 P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 5 to prefer a specific alternative, 88 to 97 % of the plots were prioritized for that alternative. For example, when the decision-maker preferred timber production more than others, 97 % of plots were prioritized accordingly. Scenic beauty had the lowest (88 %) priority rate, which means that 12 % of the plots were prioritized for other alternatives, even if the decision-maker’s preferences were highest for scenic beauty. In the plots where the prioritization did not follow the preferences, the forest structure (evaluated by the ES proxy) was stronger in favor of another alternative than the preference. Example ALS transects of these plots are visualized in Fig. 2. Forests prioritized most closely according to the preference value are shown on the upper-left to the lower-right diagonal of the matrix (Fig. 2). The remaining cells of Fig. 2 represent example forests that were prioritized for other alternatives even if the decision-maker preferred this alternative. For instance, for a decision- maker preferring carbon storage, an “ideal” forest for this preference in terms of the CARB proxy is found from the diagonal (CARB-CARB intersection in Fig. 2). Other plots shown on that row were prioritized for the alternatives indicated in the columns because of their forest structure as implied by the specific proxy. Correspondingly, the other plots shown on the CARB column were prioritized for this alternative due to their forest structure regardless that their decision-makers preferred more the other alternatives shown as the row index. 3.2. Landscape-level prioritization using k-NN (k ˆ 1 or k ˆ 3) The imputation of the SMAA results from the reference plots pro- duced a map of the first rank acceptability index (Fig. 3), which showed the preferred alternative and the strength by which it was preferred at a given location. According to the plot-level results presented above, analyzing those alternatives that were selected for a plot, regardless of the preferences, could produce a way to map locations where the forests should be managed differently from the decision-makers’ general pref- erences due to more favorable forest structure for another alternative. Such hotspot map was produced as the digest of the six alternative maps showing the locations selected for another alternative that was preferred (k ˆ 1) (Fig. 4). The individual maps, from which Fig. 4 was digested, are Fig. 1. A schematic overview of the SMAA approach to analyze landscape composed of pixels. The k-NN search was used to find (here, k ˆ 3) nearest neighbor plots for every pixel (a), to impute the matrix showing how many per cent of the N ˆ 10,000 randomly sampled criteria weights supported the decision alternatives in the reference plots (b). This matrix is typically interpreted by visualizing the rank acceptability indices (c) or sorting them cumulatively (d). Refer to Tables 3–5 for the abbreviations of the ESs considered. Table 3 The proportion of reference plots prioritized for the six alternatives based on the first rank acceptability index when decision-makers were assumed neutral or had the highest preference on the indicated alternative. The numbers give the percentage prioritized for the alternatives based on joint analysis of the preferences and proxies. Bolded values indicate the alternative preferred by the decision-maker. ESs Most preferred alternative and % of plots prioritized for an alternative under these preferences Neutral BILB COWB AMEN BIOD CARB TIMB 1. Suitability for bilberry picking (BILB) 17.6 92.2 0 0 3.9 2.9 1.0 2. Suitability for cowberry picking (COWB) 39.2 3.9 93.1 2.0 1.0 2.0 1.0 3. Scenic beauty (AMEN) 1.0 0 0 88.2 0 0 1.0 4. Biodiversity conservation (BIOD) 15.7 2.9 5.9 2.0 92.2 0 0 5. Carbon storage (CARB) 7.8 1.0 1.0 2.9 0 95.1 0 6. Timber production (TIMB) 18.6 0 0 4.9 2.0 0 97.1 P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 6 shown in the Supplementary material (Appendix A). In addition, Appendices B–D show the results corresponding to Figs. 3, 4, and Ap- pendix A, with k ˆ 3. Similar statistics as for the reference field plots were produced to analyze the maps numerically. When the decision-makers were assumed neutral, the k-NN imputation with k ˆ 1 prioritized 41 % of the land- scape for cowberry, followed by bilberry (17 %), timber (17 %), biodi- versity (15 %), carbon (8 %), and scenic beauty (3 %) (Table 4). When the decision-maker had a specific preference, the prioritized area varied from 84 to 97 % and, as above, the remainder represents the area prioritized for another alternative based on their favorable proxy values. For example, when the decision-maker preferred cowberry more than others, 97 % of the area was prioritized for that alternative. Scenic beauty had the lowest (84 %) prioritized area compared to when decision-makers had other preferences. When increasing the considered numbers of neighboring plots to a k- NN (k ˆ 3) analysis, with neutral decision-makers’ preferences, 37 % of the landscape was prioritized for cowberry, followed by timber (18 %), carbon (16 %), bilberry (15 %), biodiversity (13 %), and scenic beauty (1 %) (Table 5). The area prioritized for the most preferred alternative varied from 94 to 100 %. For example, when the decision-maker prefers cowberry more than others, 100 % of plots were prioritized for that alternative. A comparison of the imputation results based on increasing the number of nearest neighbors from k ˆ 1 to k ˆ 3 suggests changes in the proportion of the area prioritized for the different alternatives (Tables 4 and 5). A map comparison (Fig. 5) indicated a higher degree of aggre- gation of the neighboring pixels with the same alternative occurred in the map produced by k ˆ 3 as compared to k ˆ 1. Especially, the pro- portion of carbon storage increased in the k ˆ 3 imputed map compared to the k ˆ 1 imputation (Fig. 5). There were two big areas where the prioritized management changed completely, which is explored in the Discussion. 4. Discussion Our study is a continuation of Vauhkonen (2018) and Vauhkonen and Ruotsalainen (2017) to use proxy maps for prioritizing forest Fig. 2. Examples of prioritization decisions based on the preferences vs proxy values. The row indices indicate the most preferred alternative assumed in the SMAA analysis. The plots located on the upper-left to the lower-right diagonal represent an “ideal” forest for each alternative as that receiving the highest first rank acceptability value. The remaining cells show forests that were prioritized for the alternatives indicated by the row indices even if the decision-maker preferred the alternative of the column index. The cross () symbol depicts that no plots were prioritized according to the row index when the decision-maker had the preference indicated in the column index. P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 7 Fig. 3. The prioritized alternative for a neutral decision-maker and the first rank acceptability index on which the prioritization was based (k ˆ 1). The scale of each color ramp ranges from 0 (low) to 1 (high). P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 8 management so that their provisioning can be enhanced. Similar to those studies, we considered timber, carbon storage, biological di- versity, and recreational values (composed of visual amenity and suit- ability for berry picking), for which the management was prioritized. As explained in Section 1, this division was reasoned from the perspective of segregating management and not aiming to strictly follow an ES classification system. Similar to Englund et al. (2017), citing Costanza (2008), we find motivated that “there are many useful ways to classify ecosystem goods and services, and that the goal should not be to have a single, consistent system, but rather a pluralism of typologies that can be useful for different purposes”. It is acknowledged that a different set of forest ESs could have been considered in the analyses, especially if it was known that the area in question was under integrated management for the provisioning of ESs typical to rural or urban areas (cf., Tammi et al., 2017). Also, alternative indicators for forest ESs, such as different models for berry yields (Miina et al., 2020) or habitat suitability of Fig. 4. Hotspot map, i.e., the digest of 6 maps showing the locations selected for another alternative regardless of the assumed highest priority (k ˆ 1). P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 9 specific species, could be used (cf., Pukkala, 2016). In all, our example should be sufficiently illustrative on applying a SMAA-based prioriti- zation for multiple objectives, where the exact set of ESs and their in- dicators can be changed depending on the MCDA task at hand. The novelty of the present study was to account for the decision- makers’ preferences as stochastic distributions in the pixel-level ana- lyses. Even the real-world decision-makers are often unaware or may not know exactly how much they value or prefer alternative ESs. Our pro- posed approach utilized the SMAA (Lahdelma et al., 1998; Lahdelma and Salminen, 2001) approach to sample from unknown or partially known weight space. Somewhat resembling approaches based on the fuzzy logic (Biber et al., 2021) and SMAA (Pukkala, 2021) were recently presented for other ES evaluations. While the SMAA approach itself has been used in more than a hundred articles published between 1998 and 2017 (Pelissari et al., 2020), this is, to our best knowledge, the first spatial SMAA approach. The examples by Ligmann-Zielinska and Jankowski (2014), García Marquez et al. (2017) and Vauhkonen (2018) suggested benefits from including uncertainties in the map analyses of the above ESs and pos- sibilities to do this by means of Monte Carlo simulations and similar techniques. These approaches bear a considerable computational burden in processing big data formed by the map units. We circum- vented this problem by imputing the SMAA matrices from the plots to the pixels. We demonstrated two types of maps based on operating the SMAA in the reference plots: (1) prioritization maps (Fig. 3) showing the alternative that a decision-maker with given preferences should choose, and the strength of this decision as the degree of the first rank accept- ability in a pixel (Fig. 3); and (2) a hotspot map that was formed as a digest of maps showing alternatives that should have been selected because their specific proxies surpassed the preferences (Fig. 4). The former type of map can be used similar to those derived by Ligmann- Zielinska and Jankowski (2014) to determine areas with a high value and either low or high deviation to be prioritized with consensus or reservation, respectively. The latter map can be useful in identifying areas where the forest structure favored another alternatives than the preferences and could be therefore checked for potential as hotspots for diversifying the management of the landscape. The stochastic approach allowed estimating the strength of the decision with respect to the un- certainty in both the proxy values and preferences. Producing the maps by means of the nearest neighbor imputation does not suffer from similar computational burden as compared to Ligmann-Zielinska and Jankowski (2014). We emphasize that this is the first study applying SMAA for the above purpose so that the studied setup likely includes alleviations that should be re-examined in later studies. The performance of the k-NN method depends on the parameterization of the method with respect to the available reference data. We note that a very similar distribution of al- ternatives was produced by imputing the rank priorities of the neutral decision-maker from field plots using k ˆ 1 (cf., Tables 3 and 4 of plot- level and landscape-level priorities, respectively). This was expected as Vauhkonen and Imponen (2016) had determined the field sample plot locations by selecting evenly according to the Euclidean distance, which we also used instead of many alternative distance metrics (cf., Crookston and Finley, 2008). Producing a similar landscape from the plots origi- nally selected from that landscape, although using different ALS metrics and approach was, therefore, a sanity check of the method. However, there may be several reasons to weight observations in terms of the distances, for example. We propose that the optimal k-NN parameteri- zation for imputing SMAA matrices is studied in data selected differently and using methods applied for diameter distribution studies (Maltamo et al., 2009; Raty et al., 2018), which require multicriteria consider- ations on the performance of the method and therefore conceptually resemble our analyses. Although the number of our field reference plots (102) can be considered sufficient in the light of findings by Maltamo et al. (2009), for example, the k-NN method cannot extrapolate and evaluating the mapping accuracy of our method is only relevant in broader reference data. Fig. 5 indicates peculiarities that originate from the lack of fully representative reference data for the entire area. Fig. 5 shows two big areas, where the prioritized management changed completely as a result of imputing with k ˆ 3 instead of k ˆ 1. When these areas were evalu- ated against publicly available aerial photographs of the area, both of them were found to be treeless areas: one due to a clear-felling and another was a poorly productive mire forest next to a peatland pond. No close-to-treeless areas were present in the reference data (cf., Table 1), which explains the change and suggests that results for these areas may be artifacts that originate from the lack of representative reference data. Nevertheless, in a mapping exercise, the imputation had to be done to the entire area with the nearest neighbors available in the reference data. In our data, the increased proportion of carbon storage in Table 4 The proportion of land prioritized for the alternatives based on the nearest neighbor imputation (k ˆ 1) of the first rank acceptabilities when decision-makers are neutral or prefer specific alternative. Bolded values indicate the alternative preferred by the decision-maker. Ecosystem services Most preferred alternative and % of pixels prioritized for an alternative under these preferences Neutral BILB COWB AMEN BIOD CARB TIMB 1. Suitability for bilberry picking (BILB) 17.1 88.9 0 0.0 2.0 1.2 0.5 2. Suitability for cowberry picking (COWB) 40.6 6.9 96.5 7.8 0.2 10.1 10.0 3. Scenic beauty (AMEN) 2.5 0.0 0.0 84.4 0.0 0.0 0.3 4. Biodiversity conservation (BIOD) 14.9 2.1 3.3 1.3 95.7 0.0 0.0 5. Carbon storage (CARB) 7.8 2.1 0.2 4.0 0.0 88.6 0.0 6. Timber production (TIMB) 17.0 0.0 0.0 2.6 2.2 0.0 89.3 Table 5 The proportion of land prioritized for the alternatives based on the nearest neighbor imputation (k ˆ 3) of the first rank acceptabilities when decision-makers are neutral or prefer specific alternative. Bolded values indicate the alternative preferred by the decision-maker. Ecosystem services 0Most preferred alternative and % of pixels prioritized for an alternative under these preferences Neutral BILB COWB AMEN BIOD CARB TIMB 1. Suitability for bilberry picking (BILB) 15.2 96.9 0.0 0.0 1.1 0.0 0.0 2. Suitability for cowberry picking (COWB) 37.4 1.3 99.6 1.7 0.4 0.9 1.1 3. Scenic beauty (AMEN) 0.9 0.0 0.0 93.8 0.0 0.0 0.0 4. Biodiversity conservation (BIOD) 13.4 0.5 0.5 0.0 98.1 0.0 0.0 5. Carbon storage (CARB) 15.5 0.6 0.0 0.4 0.0 99.1 0.0 6. Timber production (TIMB) 17.6 0.7 0.0 4.1 0.4 0.0 98.9 P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 10 Fig. 5. A map comparison of the effect of neighborhood in the k-NN imputation for producing the landscape-level results. Two areas, delin- eated by red and blue rectangles from the whole landscape shown in the middle, are highlighted in the top and bottom, respectively. The sub- figures a vs b on the top row and c vs d on the bottom row show results for k ˆ 1 vs k ˆ 3, respectively. The scale of each color ramp ranges from 0 (low) to 1 (high). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 11 imputations with k ˆ 3 (cf., the neutral column of Tables 4 and 5) can be judged as exaggeration, reasoned by the discussion above. It reinforces the indication by Kangas et al. (2018) on how harmful the non-random errors in the data can be for trade-off analyses between different ESs. We mapped the first-rank acceptability but note that SMAA produces many alternative ways to interpret the results (Fig. 1c-d). In particular, acceptability analyses consider the weights that support a certain alternative rather than discovering which alternative is preferred using certain weights (Kangas et al., 2015). We did not analyze the weights resulting to the first-rank acceptabilities but note that spatial optimi- zation approaches (e.g., Heinonen et al., 2007), which require prefer- ence information already for determining the optimization task, could benefit from the feasible weights produced by SMAA. Overall, many additional quantifications could be produced assuming stochastic dis- tributions even for individual pixels (Ligmann-Zielinska and Jankowski, 2014) that enables considering risk preferences (risk aversion) of a decision-maker (Vauhkonen, 2018). These analyses, carried out at a sub- stand-level, may allow more efficiently using the production possibil- ities of the forest (Heinonen et al., 2007). However, the added benefit of these quantifications should obviously be studied as there is a risk of “too big data”, meaning an increased computational cost but a low utility of these quantifications. 5. Conclusions We described incorporating decision-makers’ preferences, quantified as stochastic distributions, into a management prioritization of a Scan- dinavian boreal forest in a realistic, computationally feasible, and spatially explicit manner. The nearest-neighbor imputation was used to derive landscape-level maps from plot-level reference data, where the decision-makers’ preferences were assumed totally or partially un- known and sampled from a probability distribution. We showed the potential to derive two types of maps based on first-rank acceptability indices for the forest use prioritization: one showing the preferred alternative and strength for a decision-maker with neutral preferences and another showing identified areas where the suitability of the forest structure suggested different alternatives than the preferences. When field reference data are available, applying the proposed approach can be done based on publicly available ALS data, which is helpful for the management planning of that landscape. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgments Parvez Rana was supported by internal funding of Natural Resources Institute Finland (Solutions 41007-00183800, RemoResp 41007- 00216200) and Maj and Tor Nessling Foundation (MONESER 201900219). Special thanks to Anne Tolvanen, Natural Resources Institute Finland, for the subfigure (1a). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.landurbplan.2022.104637. References Alho, J. M., & Kangas, J. (1997). Analyzing uncertainties in experts’ opinions of forest plan performance. Forest Science, 43, 521–528. https://doi.org/10.1093/ forestscience/43.4.521 Barton, D. N., Kelemen, E., Dick, J., Martin-Lopez, B., Gomez-Baggethun, E., Jacobs, S., … Lapola, D. M. (2018). (Dis)integrated valuation – Assessing the information gaps in ecosystem service appraisals for governance support. Ecosystem Services, 29, 529–541. https://doi.org/10.1016/j.ecoser.2017.10.021 Biber, P., Schwaiger, F., Poschenrieder, W., & Pretzsch, H. (2021). A fuzzy logic-based approach for evaluating forest ecosystem service provision and biodiversity applied to a case study landscape in Southern Germany. European Journal of Forest Research, 140, 1559–1586. https://doi.org/10.1007/s10342-021-01418-4 Boerema, A., Rebelo, A. J., Bodi, M. B., Esler, K. J., & Meire, P. (2017). Are ecosystem services adequately quantified? Journal of Applied Ecology, 54, 358–370. https://doi. org/10.1111/1365-2664.12696 Bryant, B. P., Borsuk, M. E., Hamel, P., Oleson, K. L., Schulp, C. J. E., & Willcock, S. (2018). Transparent and feasible uncertainty assessment adds value to applied ecosystem services modeling. Ecosystem Services, 33, 103–109. https://doi.org/ 10.1016/j.ecoser.2018.09.001 Caglayan, _I., Yes¸il, A., Kabak, O., & Bettinger, P. (2021). A decision making approach for assignment of ecosystem services to forest management units: A case study in northwest Turkey. Ecological Indicators, 121. https://doi.org/10.1016/j. ecolind.2020.107056 Convertino, M., Baker, K. M., Vogel, J. T., Lu, C., Suedel, B., & Linkov, I. (2013). Multi- criteria decision analysis to select metrics for design and monitoring of sustainable ecosystem restorations. Ecological Indicators, 26, 76–86. https://doi.org/10.1016/j. ecolind.2012.10.005 Costanza, R. (2008). Ecosystem services: Multiple classification systems are needed. Biological Conservation, 141, 350–352. https://doi.org/10.1016/j. biocon.2007.12.020 Crookston, N. L., & Finley, A. O. (2008). yaImpute: An R Package for kNN Imputation. Journal of Statistical Software, 23, 1–16. https://doi.org/10.18637/jss.v023.i10. Egoh, B., Reyers, B., Rouget, M., Richardson, D. M., Le Maitre, D. C., & van Jaarsveld, A. S. (2008). Mapping ecosystem services for planning and management. Agriculture, Ecosystems and Environment, 127, 135–140. https://doi.org/10.1016/j. agee.2008.03.013 Englund, O., Berndes, G., & Cederberg, C. (2017). How to analyse ecosystem services in landscapes—A systematic review. Ecological Indicators, 73, 492–504. https://doi. org/10.1016/j.ecolind.2016.10.009 Forsius, M., Kujala, H., Minunno, F., Holmberg, M., Leikola, N., Mikkonen, N., … Heikkinen, R. K. (2021). Developing a spatially explicit modelling and evaluation framework for integrated carbon sequestration and biodiversity conservation: Application in southern Finland. Science of the Total Environment, 775. https://doi. org/10.1016/j.scitotenv.2021.145847 García Marquez, J. R., Krueger, T., Paez, C. A., Ruiz-Agudelo, C. A., Bejarano, P., Muto, T., & Arjona, F. (2017). Effectiveness of conservation areas for protecting biodiversity and ecosystem services: A multi-criteria approach. International Journal of Biodiversity Science, Ecosystem Services and Management, 13, 1–13. https://doi.org/ 10.1080/21513732.2016.1200672 Gre^t-Regamey, A., Siren, E., Brunner, S. H., & Weibel, B. (2017). Review of decision support tools to operationalize the ecosystem services concept. Ecosystem Services, 26, 306–315. https://doi.org/10.1016/j.ecoser.2016.10.012 Heinonen, T., Kurttila, M., & Pukkala, T. (2007). Possibilities to aggregate raster cells through spatial optimization in forest planning. Silva Fennica, 41, 89–103. htt ps://doi.org/10.14214/sf.474. Honeck, E., Sanguet, A., Schlaepfer, M. A., Wyler, N., & Lehmann, A. (2020). Methods for identifying green infrastructure. SN Applied Sciences, 2. https://doi.org/10.1007/ s42452-020-03575-4 1-25 Ihalainen, M., Alho, J., Kolehmainen, O., & Pukkala, T. (2002). Expert models for bilberry and cowberry yields in Finnish forests. Forest Ecology and Management, 157, 15–22. https://doi.org/10.1016/S0378-1127(00)00653-8 Juutinen, A., Kurttila, M., Pohjanmies, T., Tolvanen, A., Kuhlmey, K., Skudnik, M., … Makipaa, R. (2021). Forest owners’ preferences for contract-based management to enhance environmental values versus timber production. Forest Policy and Economics, 132. https://doi.org/10.1016/j.forpol.2021.102587 Kangas, A. (2006). The risk of decision making with incomplete criteria weight information. Canadian Journal of Forest Research, 36, 195–205. https://doi.org/ 10.1139/x05-243 Kangas, A., Korhonen, K. T., Packalen, T., & Vauhkonen, J. (2018). Sources and types of uncertainties in the information on forest-related ecosystem services. Forest Ecology and Management, 427, 7–16. https://doi.org/10.1016/j.foreco.2018.05.056 Kangas, A., Kurttila, M., Hujala, T., Eyvindson, K., Kangas, J., & (Eds.).. (2015). Decision support for forest management. Berlin: Springer -. Kangas, A., Leskinen, P., & Kangas, J. (2007). Comparison of fuzzy and statistical approaches in multicriteria decisionmaking. Forest Science, 53, 37–44. https://doi. org/10.1093/forestscience/53.1.37 Kangas, J. (1993). A multi-attribute preference model for evaluating the reforestation chain alternatives of a forest stand. Forest Ecology and Management, 59, 271–288. https://doi.org/10.1016/0378-1127(93)90007-A Kangas, J., Hokkanen, J., Kangas, A. S., Lahdelma, R., & Salminen, P. (2003). Applying stochastic multicriteria acceptability analysis to forest ecosystem management with both cardinal and ordinal criteria. Forest Science, 49, 928–937. https://doi.org/ 10.1093/forestscience/49.6.928 Kangas, J., Store, R., Leskinen, P., & Mehtatalo, L. (2000). Improving the quality of landscape ecological forest planning by utilising advanced decision-support tools. P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 12 Forest Ecology and Management, 132, 157–171. https://doi.org/10.1016/S0378-1127 (99)00221-2 Karjalainen, T., & Kellomaki, S. (1996). Greenhouse gas inventory for land use change and forestry in Finland based on international guidelines. Mitigation and Adaptation Strategies for Global Change, 1, 51–71. https://doi.org/10.1007/BF00625615 Knoke, T., Kindu, M., Schneider, T., & Gobakken, T. (2021). Inventory of Forest Attributes to Support the Integration of Non-provisioning Ecosystem Services and Biodiversity into Forest Planning—from Collecting Data to Providing Information. Current Forestry Reports, 7, 38–58. https://doi.org/10.1007/s40725-021-00138-7 Korhonen, L., Peuhkurinen, J., Malinen, J., Suvanto, A., Maltamo, M., Packalen, P., & Kangas, J. (2008). The use of airborne laser scanning to estimate sawlog volumes. Forestry, 81, 499–510. https://doi.org/10.1093/forestry/cpn018 Koschke, L., Fürst, C., Frank, S., & Makeschin, F. (2012). A multi-criteria approach for an integrated land-cover-based assessment of ecosystem services provision to support landscape planning. Ecological Indicators, 21, 54–66. https://doi.org/10.1016/j. ecolind.2011.12.010 Lahdelma, R., & Salminen, P. (2001). SMAA-2: Stochastic Multicriteria Acceptability Analysis for Group Decision Making. Operations Research, 49(3), 444–454. https:// www.jstor.org/stable/3088639. Lahdelma, R., Hokkanen, J., & Salminen, P. (1998). SMAA - Stochastic multiobjective acceptability analysis. European Journal of Operational Research, 106, 137–143. https://doi.org/10.1016/S0377-2217(97)00163-X Langemeyer, J., Gomez-Baggethun, E., Haase, D., Scheuer, S., & Elmqvist, T. (2016). Bridging the gap between ecosystem service assessments and land-use planning through Multi-Criteria Decision Analysis (MCDA). Environmental Science and Policy, 62, 45–56. https://doi.org/10.1016/j.envsci.2016.02.013 Lehtomaki, 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, 10, 1–25. https://doi.org/10.1371/ journal.pone.0135926 Lehtonen, A., Makipaa, R., Heikkinen, J., Sievanen, R., & Liski, J. (2004). Biomass expansion factors (BEFs) for Scots pine, Norway spruce and birch according to stand age for boreal forests. Forest Ecology and Management, 188, 211–224. https://doi. org/10.1016/j.foreco.2003.07.008 Liang, X., Hyyppa, J., & Matikainen, L. (2007). Deciduous-coniferous tree classification using difference between first and last pulse laser signatures. In: P. Ronnholm, H. Hyyppa, J. & Hyyppa (Eds.) Proceedings of ISPRS workshop on Laser Scanning and SilviLaser. ISPRS Archives, XXXVI-3/W52, 253–257. https://www.isprs.org/ proceedings/XXXVI/3-W52/final_papers/Liang_2007.pdf. Ligmann-Zielinska, A., & Jankowski, P. (2014). Spatially-explicit integrated uncertainty and sensitivity analysis of criteria weights in multicriteria land suitability evaluation. Environmental Modelling and Software, 57, 235–247. https://doi.org/ 10.1016/j.envsoft.2014.03.007 Malczewski, J., & Rinner, C. (Eds.). (2015). Multicriteria decision analysis in geographic information science. New York: Springer. Maltamo, M., Næsset, E., Bollandsås, O. M., Gobakken, T., & Packalen, P. (2009). Non- parametric prediction of diameter distributions using airborne laser scanner data. Scandinavian Journal of Forest Research, 24(6), 541–553. https://doi.org/10.1080/ 02827580903362497 Mann, C., Loft, L., Hernandez-Morcillo, M., Primmer, E., Bussola, F., Falco, E., … Winkel, G. (2022). Governance Innovations for forest ecosystem service provision – Insights from an EU-wide survey. Environmental Science & Policy, 132, 282–295. https://doi.org/10.1016/j.envsci.2022.02.032 Marttunen, M., Mustajoki, J., Lehtoranta, V., & Saarikoski, H. (2021). Complementary use of the ecosystem service concept and multi-criteria decision analysis in water management. Environmental Management, 69, 719–734. https://doi.org/10.1007/ s00267-021-01501-x Miina, J., Pukkala, T., & Kurttila, M. (2016). Optimal multi-product management of stands producing timber and wild berries. European Journal of Forest Research, 135, 781–794. https://doi.org/10.1007/s10342-016-0972-9 Miina, J., Kurttila, M., Calama, R., & de-Miguel, S., & Pukkala, T.. (2020). Modelling non- timber forest products for forest management planning in Europe. Current Forestry Reports, 6, 309–322. https://doi.org/10.1007/s40725-020-00130-7 Mustajoki, J., Saarikoski, H., Belton, V., Hjerppe, T., & Marttunen, M. (2020). Utilizing ecosystem service classifications in multi-criteria decision analysis – Experiences of peat extraction case in Finland. Ecosystem Services, 41. https://doi.org/10.1016/j. ecoser.2019.101049 National Land Survey of Finland. (2015). File service of open data. Retrieved from https://tiedostopalvelu.maanmittauslaitos.fi/tp/kartta?langˆen. Accessed October 11, 2015. Nemec, K. T., & Raudsepp-Hearne, C. (2013). The use of geographic information systems to map and assess ecosystem services. Biodiversity and Conservation, 22, 1–15. https://doi.org/10.1007/s10531-012-0406-z Nieminen, E., Kareksela, S., Halme, P., & Kotiaho, J. S. (2021). Quantifying trade-offs between ecological gains, economic costs, and landowners’ preferences in boreal mire protection. Ambio, 50, 1841–1850. https://doi.org/10.1007/s13280-021- 01530-0 Ørka, H. O., Gobakken, T., Næsset, E., Ene, L., & Lien, V. (2012). Simultaneously acquired airborne laser scanning and multispectral imagery for individual tree species identification. Canadian Journal of Remote Sensing, 38, 125–138. https://doi. org/10.5589/m12-021 Pasalodos-Tato, M., Makinen, A., Garcia-Gonzalo, J., Borges, J. G., Lamås, T., & Eriksson, L. O. (2013). Assessing uncertainty and risk in forest planning and decision support systems: Review of classical methods and introduction of innovative approaches. Forest Systems, 22, 282–303. https://doi.org/10.5424/fs/2013222- 03063 Pelissari, R., Oliveira, M. C., Amor, S. B., Kandakoglu, A., & Helleno, A. L. (2020). SMAA methods and their applications: A literature review and future research directions. Annals of Operations Research, 293, 433–493. https://doi.org/10.1007/s10479-019- 03151-z Pohjanmies, T. (2018). Trade-offs among intensive forestry, ecosystem services and biodiversity in boreal forests. Academic dissertation. Jyvaskyla studies in biological and environmental science. Pohjanmies, T., Eyvindson, K., Trivi~no, M., Bengtsson, J., & Monkkonen, M. (2021). Forest multifunctionality is not resilient to intensive forestry. European Journal of Forest Research, 140, 537–549. https://doi.org/10.1007/s10342-020-01348-7 Pukkala, T. (2005). Metsikon tuottoarvon ennustemallit kivennaismaan mannikoille, kuusikoille ja rauduskoivikoille (in Finnish for ”Prediction models for the expectation value of pine, spruce and birch stands on mineral soils”). Metsatieteen Aikakauskirja, 2005, 311–322. https://doi.org/10.14214/ma.5778. Pukkala, T. (2008). Integrating multiple services in the numerical analysis of landscape design. In K. von Gadow, & T. Pukkala (Eds.), Designing Green Landscapes. Managing Forest Ecosystems. Dordrecht: Springer -. https://doi.org/10.1007/978-1-4020-6759- 4_6. Pukkala, T. (2016). Which type of forest management provides most ecosystem services? Forest Ecosystems, 3. https://doi.org/10.1186/s40663-016-0068-5 Pukkala, T. (2021). Measuring the social performance of forest management. Journal of Forestry Research, 32, 1803–1818. https://doi.org/10.1007/s11676-021-01321-z Pukkala, T., & Kangas, J. (1993). A heuristic optimization method for forest planning and decision making. Scandinavian Journal of Forest Research, 8, 560–570. https://doi. org/10.1080/02827589309382802 Pukkala, T., Kellomaki, S., & Mustonen, E. (1988). Prediction of the amenity of a tree stand. Scandinavian Journal of Forest Research, 3, 533–544. https://doi.org/10.1080/ 02827588809382538 Raty, J., Packalen, P., & Maltamo, M. (2018). Comparing nearest neighbor configurations in the prediction of species-specific diameter distributions. Annals of Forest Science, 75. https://doi.org/10.1007/s13595-018-0711-0 Schroter, M., & Remme, R. P. (2016). Spatial prioritisation for conserving ecosystem services: Comparing hotspots with heuristic optimisation. Landscape Ecology, 31, 431–450. https://doi.org/10.1007/s10980-015-0258-5 Silvennoinen, H. (2017). Metsamaiseman kauneus ja metsanhoidon vaikutus koettuun maisemaan metsikkotasolla (in Finnish for ”Scenic beauty of forest stands and impact of management”). Academic dissertation. Dissertationes Forestales, 242. htt ps://doi.org/10.14214/df.242. Silvennoinen, H., Pukkala, T., & Tahvanainen, L. (2002). Effect of cuttings on the scenic beauty of a tree stand. Scandinavian Journal of Forest Research, 17, 263–273. https:// doi.org/10.1080/028275802753742936 Store, R., & Kangas, J. (2001). Integrating spatial multi-criteria evaluation and expert knowledge for GIS-based habitat suitability modelling. Landscape and Urban Planning, 55, 79–93. https://doi.org/10.1016/S0169-2046(01)00120-7 Tammi, I., Mustajarvi, K., & Rasinmaki, J. (2017). Integrating spatial valuation of ecosystem services into regional planning and development. Ecosystem Services, 26, 329–344. https://doi.org/10.1016/j.ecoser.2016.11.008 Turtiainen, M. (2015). Modelling bilberry and cowberry yields in Finland: Different approaches to develop models for forest planning calculations. Academic dissertation. Dissertationes Forestales, 185. https://doi.org/10.14214/df.185. Uhde, B., Hahn, A., Griess, V. C., & Knoke, T. (2015). Hybrid MCDA methods to integrate multiple ecosystem services in forest management planning: A critical review. Environmental Management, 56, 373–388. https://doi.org/10.1007/s00267-015- 0503-3 Vauhkonen, J. (2018). Predicting the provisioning potential of forest ecosystem services using airborne laser scanning data and forest resource maps. Forest Ecosystems, 5. https://doi.org/10.1186/s40663-018-0143-1 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. https://doi.org/10.1093/forestry/cpw011 Vauhkonen, J., & Ruotsalainen, R. (2017). 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. Forest Ecology and Management, 389, 272–284. https://doi.org/10.1016/j.foreco.2016.12.005 Vauhkonen, J., Packalen, P., Malinen, J., Pitkanen, J., & Maltamo, M. (2014). Airborne laser scanning-based decision support for wood procurement planning. Scandinavian Journal of Forest Research, 29, 132–143. https://doi.org/10.1080/ 02827581.2013.813063 Verhagen, W., Kukkala, A. S., Moilanen, A., van Teeffelen, A. J., & Verburg, P. H. (2017). Use of demand for and spatial flow of ecosystem services to identify priority areas. Conservation Biology, 31(4), 860–871. https://doi.org/10.1111/cobi.12872 Weiss, G., Lawrence, A., Lidestav, G., Feliciano, D., Hujala, T., Sarvasova, Z., … Zivojinovic, I. (2019). Research trends: Forest ownership in multiple perspectives. Forest Policy and Economics, 99. https://doi.org/10.1016/j.forpol.2018.10.006 Weiss, G., & Zivojinovic, I. (2020). Indicator 6.1 Forest holdings. In FOREST EUROPE, 2020: State of Europe’s Forests 2020. http://www.foresteurope.org. Yatsalo, B., Didenko, V., Gritsyuk, S., & Sullivan, T. (2015). Decerns: A framework for multi-criteria decision analysis. International Journal of Computational Intelligence Systems, 8, 467–489. https://doi.org/10.1080/18756891.2015.1023586 Further reading Sievanen, T., & Neuvonen, M. (Eds.). (2011). Luonnon virkistyskaytto 2010 (in Finnish for “The recreational use of nature in 2010”). Working Papers of the Finnish Forest P. Rana and J. Vauhkonen Landscape and Urban Planning 230 (2023) 104637 13 Research Institute 212, http://www.metla.fi/julkaisut/workingpapers/2011/ mwp212.htm. P. Rana and J. Vauhkonen