iForest - Biogeosciences and Forestry

iForest - Biogeosciences and Forestry

Making objective forest stand maps of mixed managed forest with spatial interpolation and multi-criteria decision analysis

iForest - Biogeosciences and Forestry, Volume 6, Issue 5, Pages 268-277 (2013)
doi: https://doi.org/10.3832/ifor0099-006
Published: Jul 01, 2013 - Copyright © 2013 SISEF

Research Articles

The spatial interpolation and multi-criteria decision analysis (MCDA) capabilities of geographic information systems have the potential to create new approaches to forest management. In this study of heterogeneously structured stand maps, the potential use of the regularized spline with tension (RST) interpolation method and of ELECTRE TRI MCDA was investigated. For each species and diameter class, one map of the predicted volume per ha was produced with the RST method. The map used data from a total of 1050 circular sample plots. By repeating the same process for the eight species occurring in the study area, 31 volume maps were produced. The accuracy of these prediction maps was calculated at the pixel (20 x 20 m) level and at the area level (per ha). An accuracy of greater than or equal to 97% was achieved at the pixel level, whereas a minimum accuracy of 86% was achieved for the area-based calculations. In addition, these 31 volume maps were compared with the management report results obtained from the government institute responsible for defining management plans. These comparisons were performed for the total volume of all species with volume ratios greater than 1%. The comparisons showed 21, 14, 4, and 2 % accuracies for Calabrian pine, Oriental beech, black pine and oak species, respectively. Following interpolation, these 31 maps were geo-computed, and a volume-based stand map was produced. The 890 different mixture variations resulting from combinations of the volume of species composition and the stand diameter class in these maps were classified according to expert knowledge. In this classification process, ELECTRE TRI MCDA was used to benefit from the capabilities provided by geographic information systems. Finally, ELECTRE TRI was used to reduce the 890 different mixture combinations to 70 stand type classes.

Regularized Spline with Tension, Multi-Criteria Decision Analysis, Forest Map, Forest Management


Information about extant forest communities is of considerable importance in forest management planning. The stand is the smallest production and silvicultural unit, and the most crucial part of this information. The location of planning units, their size and silvicultural features are used to select management objectives and determine the working cycle, management regulations and silvicultural goals. Thus, the precision of decisions to be implemented and the planning to be realized are both closely related to the accuracy of the stand maps. Moreover, map accuracy depends on the type of forest inventory, the sample size, and the techniques used to evaluate the inventory results.

Stand maps have previously been produced subjectively ([13]) in conjunction with current forest management practices. In conventional management planning, stands are defined through the planning process. The type of of the management objectives is a component of this process ([20]). However, the use of computers has changed this conventional approach, and spatially and temporally dynamic description and treatment units can readily be produced ([25]). Thus, dynamic forest planning can be achieved with the use of geopositioned field plot data and computers ([64]).

The sampling methods used in forest inventories vary according to different countries’ forestry objectives and forest structures. Systematic sampling methods are recommended for large and homogeneous forest areas and have been implemented in all production and conservation areas of Turkey since 1964. However, if this sampling method is implemented with an insufficient number of sample plots and/or in heterogeneous forests, the results will be highly questionable because they may fail to reflect the true forest composition. Although Sherman ([55]), Aurbi & Debouzie ([5]), Flores et al. ([14]) and D’Orazio ([10]) used remote sensing data to improve the results of such inventories, the desired outcome was never achieved in practice.

Spatial interpolation methods, which have been classified by Li & Heap ([37]) as non-geostatistical, geostatistical and combined, came into use at the end of the 1960s and have been investigated for use in forest management. As described in Akhavan et al. ([1]), Guibal ([19]) was the first study to use kriging in the forest inventory process. Jost ([28]) also used geostatistical methods to compare conventional inventory results based on systematic sampling with the results of the kriging method. Geostatistically based methods that utilize textural information are frequently used to analyze remote-sensing (RS) images ([66]). The quality and quantity of these methods have increased through advances in the computer sciences and in the science of geographical information systems. These methods allow the use of a variety of ecological and technical parameters in the assessment of inventory results.

Palmer et al. ([50]) compared the spatial predictions derived from Inverse Distance Weighting (IDW), Partial Least Squares (PLS), Regression Kriging (RK), and Ordinary Kriging (OK) for Pinus radiata mean volume, annual increment and mean top height. Viana et al. ([63]) made estimations of the crown biomass of Pinus pinaster stands and aboveground shrubland biomass using forest inventory data, remotely sensed imagery for auxiliary variables and spatial prediction models. These estimations involved the use of RK and OK, Universal Kriging (UK), IDW and Thiessen Polygon estimations. Magnussen et al. ([40]) examined the advantages of multiresolution spatial models (MRSMs) for large datasets compared with standard algorithms.

Montes et al. ([46]) have proposed a method in which the addition of experienced inventory data areas are combined with environmental information. Nanos et al. ([47]) also used a geostatistical approach for the prediction of diameter distributions in stands of Pinus pinaster Ait. Tang & Bian ([58]) successfully applied a forest site evaluation method which integrates a geographic information system (GIS) and a spatial interpolation method (kriging) in geostatistics.

Köhl et al. ([34]) suggested the use of aerial photographs, satellite images, site data, and thematic maps in the analysis and modeling of survey plans. Mandallaz ([41]) compared different geostatistical methods in forest inventory creation. Höck et al. ([27]) evaluated the kriging method using site index data. Holmgren & Thuresson ([25]) used kriging and remote sensing to predict stand volumes.

Heikkinen ([22]) and many other researchers reported that the determination of sampling errors was a serious problem associated with systematic sampling. To obtain greatly improved results from systematic sampling and inventory measurements, the use of additional explanatory environmental variables has been initiated ([36]).

Recently, the production of stand maps based on GIS has started to develop forest management planning in Turkey. Since 2004, computers have been used to perform stand type delineations and generalizations with classical and subjective operator-oriented methods and to produce digital maps of high cartographic quality. However, as emphasized by Burrough ([7]), the spatial and statistical analysis capabilities of GIS have not been fully explored and exploited yet. In contrast, GIS capabilities, combined with different spatial interpolation methods and multi-criteria decision analysis (MCDA), are being developed and implemented in numerous disciplines.

ELECTRE (ELimination Et Choix Traduisant la REalité - Elimination and Choice Expressing the Reality) is family of methods derived from MCDA that associate members to create an alternative set of relationships within the system, and then compare the alternatives obtained to solve the problem under investigation. The main advantage of these methods is that the optimal solution is not subjective, avoiding classical approaches based on individual experience and intuition. Diverse ELECTRE methods differ in minor respects. For instance, ELECTRE TRI was created to solve selection and assignment problems.

A review of more than 250 papers in the forestry literature over the past 30 years showed that 9 different versions of MCDA were applied to 9 different forestry topics (wood production, biodiversity, sustainability, reforestation, regional planning, forest industry, risks, and other subjects - [6]). A variant of the ELECTRE method was recommended by Kangas et al. ([30]) who discussed the application of the ELECTRE III and PROMETHEE II methods in 17 forest activity areas. Mendoza & Martins ([42]) critically reviewed the MCDA methods applied to the management of forests and other natural resources, describing the nature of these models and their inherent strengths and limitations. Jayanath & Gamini ([35]) also provided a comprehensive literature review in forest management evaluating 9 different types of MCDA in terms of their theoretical basis and the risks associated with their use. Currently, many opinions and suggestions exist regarding the use of MCDA in forestry. It has been claimed that MCDA alone may not be sufficient for use in forest management. On the other hand, hybrid approaches has been reported as suitable for application in this field ([29]).

Because each mapping expert has different knowledge and experience, traditional management plans can show spatial and temporal differences. To solve the above discrepancies, management plans produced using objective criteria and suitable methods (like GIS) may help to achieve consistency among maps pertaining to different revision periods and/or maps of adjacent planning units.

The main goal of this study is the creation of semi-automated, objective maps based on forest stand volume calculated from ground survey data. To achieve this goal, the spatial analysis capability of GIS was used. The regularized spline with tension (RST) spatial interpolation method was used to produce volume prediction maps of species, and ELECTRE TRI MCDA was applied to create the final forest map needed for objective decision making.

  Materials and methods 

Study Area

The study was conducted in the Asar and the Yenice forest ranges, two managed forest areas located on Mount Ida in Turkey (Fig. 1). The Asar range includes 11 159 ha of productive and 2473 ha of non-productive area. In the Yenice range the productive and non-productive areas are 10 817 ha and 1816 ha, respectively. The terrain is hilly throughout these forests. Minimum and maximum elevation are 100 m and 1.000 m a.s.l., respectively. The primary tree species found are oak species (Quercus sp.), Anatolian black pine (Pinus nigra Arn.), Calabrian pine (Pinus brutia Ten), Oriental beech (Fagus orientalis Lipsky), common hornbeam (Carpinus betulus L.), Trojan fir (Abies nordmanniana (Stev.) Spach. subsp. equi-trojani (Asch. & Sint. ex Boiss.) Coode & Cullen), Spanish chestnut (Castanea sativa Mill.) and other deciduous species (e.g., silver lime (Tilia sp.), Alnus glutinosa L. - [3], [4]).

Total forest acreage for Turkey in 2004 was 21 188 747 ha (1 831 536 ha of deciduous forests, 11 403 791 ha of coniferous forests, 2 204 268 ha of mixed forests and 5 749 152 ha of coppice ([16]).

The study area is mainly covered by mixed forests, whose structural heterogeneity is known to hinder the application of interpolation methods, thus representing a challenging test for the assessment of the applied methodology in such environments.

Vegetation Data

The 2010 management year inventory data for the Asar and Yenice forest range were used in this study. The inventory data consisted of 1 050 circular sample plots (size: 200, 400, 600, and 800 m2) systematically placed at 300 m intervals. Within each plot, measurements were taken for diameter, mean height, top height, age, health status and stem quality of each tree. The most frequent sampling plot size was chosen to achieve a pixel resolution for the interpolation maps equal to the size of the sample plots of development age c (400 m2 - see below). Due to the inventory method adopted, the forest structure within approximately 50 m of the center of the sample plot was checked visually. Any distinctive, observable characteristics within the plot and along the path from plot to plot were recorded.

The classification of even-aged stands was performed according to the 2008 Turkish forest management instructions. Diameter ranges were grouped in the following classes: (a) 1.0-7.9 cm; (b) 8.0-19.9 cm; (c) 20.0-35.9 cm; (d) 36.0-51.9 cm; and (e) > 52.0 cm. In diameter class a, no measurements were taken for the calculation of stand parameters.

Forest inventory and the production of stand maps in Turkey include the following steps:

  1. The management-planning unit is divided into homogeneous segments based on remote sensing images.
  2. In the planning units, measurements and determinations are performed in sampling plots distributed according to a systematic sampling scheme.
  3. The sampling plots are grouped by the similarity of tree species and mixture, diameter class, site class, and canopy density.
  4. The averages calculated for these groups are taken as representative of the parameters for the whole stand.
  5. Each different group is placed within a common stand type.
  6. Stand boundaries are defined based on the information recovered from the sampling points located on the stand map.
  7. A stand map based on the defined boundaries and features is produced.

The vector and attribute data for the sample plots were entered into the GRASS GIS software database ([18], [59], [49]) in accordance with the methodological steps to produce an objective forest map (Fig. 2). The calculated volumes of tree species per ha were combined with the vector data. Eight different tree species occur in the study area. Of these eight species, the most common are Anatolian black pine, oak species, Oriental beech, and Calabrian pine (Tab. 1). Oak was present in 670 of 1050 plots in the b diameter class, and black pine was present in 564 of 1050 plots in the c diameter class. Black pine had the highest mean volume per ha (124.26 m3 ha-1) in the e diameter class, and the maximum volume per ha (549.65 m3 ha-1 - Tab. 1). Based on the data for the presence of species by diameter class in the survey plots, the seven most common species according to diameter class (Qb, Pnc, Pnb, Qc, Pnd, Pne, Foc) were selected to determine the spatial interpolation tension and the smoothing parameters of the RST interpolation method. These seven species and diameter classes were chosen because potential error in the class with the largest volume and number of sample points could potentially influence the entire evaluation.

Fig. 2 - Synopsis of the methodological steps carried out in this investigation.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Tab. 1 - Summary statistics for inventory data and error results. (MAE): Mean absolute error; (MAD): Mean absolute error/Mean.

Min Max Mean Standard
400 m2
400 m2
1 ha
1 ha
b 555 0.28 142.43 14.5 21.09 0.37 0.03 1.65 0.11
c 564 3.16 348.4 64.14 61.03 1.61 0.03 7.55 0.12
d 463 13.03 500.85 91.88 80.55 2.56 0.03 10.31 0.11
e 299 30.56 549.65 124.26 84.84 2.54 0.02 11.07 0.09
Calabrian pine b 145 0.37 77.08 10.13 14.99 0.07 0.01 0.35 0.03
c 167 2.31 208.8 43.11 39.01 0.29 0.01 1.5 0.03
d 131 1.3 373.28 79.04 67.08 0.51 0.01 2.47 0.03
e 75 25.31 413.68 101.83 79.41 0.35 0 1.84 0.02
b 166 0.44 23.63 4.94 4.74 0.06 0.01 0.3 0.06
c 49 3.06 51.1 11.29 9.96 0.05 0 0.24 0.02
d 3 22 29.42 25.71 5.25 0 0 0.02 0
e 2 34 34 34 0 0 0 0.02 0
b 60 0.4 31.1 8.06 7.27 0.03 0 0.16 0.02
c 54 4.47 114.9 28.36 25.9 0.11 0 0.51 0.02
d 29 22.73 191.53 63.73 41.75 0.14 0 0.67 0.01
e 10 30.34 105.33 56.95 24.62 0.05 0 0.25 0
b 12 0.58 27.63 9.38 9.8 0.01 0 0.05 0.01
c 6 4.5 34.83 17.41 11.5 0.01 0 0.04 0
d 3 18 27 22.5 6.36 0 0 0.02 0
Oriental beech b 202 0.47 44.93 11.69 9.83 0.13 0.01 0.53 0.05
c 179 3.88 202.78 54.29 41.53 0.68 0.01 2.2 0.04
d 98 14.4 360.27 60.56 51.64 0.48 0.01 1.74 0.03
e 36 29.7 338.8 92.37 70.37 0.37 0 1.33 0.01
Spanish chestnut b 31 0.6 19.3 5.53 4.17 0.02 0 0.08 0.01
c 25 6.3 103.6 23.48 22.27 0.06 0 0.27 0.01
d 7 17.63 60.1 37.33 14.3 0.02 0 0.11 0
e 5 17.77 136.6 74.61 51.2 0.02 0 0.13 0
Oaks sp. b 670 0.44 70.85 11.89 11.35 0.35 0.03 1.61 0.14
c 514 2.71 148.75 30.55 25.19 0.72 0.02 3.2 0.10
d 289 8.73 173.85 43.62 31.53 0.75 0.02 3.39 0.08
e 158 17.85 240.25 65.46 45.29 0.71 0.01 3.43 0.05

  Enlarge/Reduce  Open in Viewer

Determining spatial interpolation parameters

To minimize the effects of prediction error and to determine the prediction parameters, different tension and smoothing parameters were tested using the “v.surf.rst” module (spatial approximation and topographic analysis from a given point or isoline data in vector format to floating point raster format using regularized spline with tension) of the GRASS software ([18]). The cross-validation and deviation values were calculated for the seven most common species by diameter class.

Cross-validation ([11]) was performed by dividing the data into two groups. The first group was used for interpolation using the method under evaluation. The second group was then used for comparison ([53], [54]). “Leave-one-out” cross-validation is a very common procedure and is frequently implemented in GIS packages. By minimizing the root-mean-square error (RMSE), cross-validation ensures that optimal interpolation parameters are found ([24]).

Prediction of forest map

Spatial interpolation methods transform data representing a continuous phenomenon in a vector point-and-line. This transformation involves mathematical functions allowed to pass through or near given points ([48]). Data collected for forest inventories are generally pinpoint-based because of the cost associated with the inventories. Pinpoint data are converted into continuous surfaces with the spatial interpolation functions implemented into the GIS package. Although many interpolation methods exist, GRASS software provides the analytic capability needed for methods such as Voronoi polygons, IDW and RSWT ([48]).

The changeable approach to intermediary estimation depends on the solution to the following problem: does the intermediary estimation function pass through or close to the data points, or can it be smoothed as much as possible? These two requirements are combined into a single condition to decrease the total deviation at the measured points and to smooth the spline function semi-norm ([44], [43]). Splines including the selection of the covariance function, as specified by numerous authors, are equivalent to universal kriging ([23], [48]).

Predictions based on the volume values per ha were obtained with the estimation parameters determined from the cross-validation and deviation method. This procedure (RST) is implemented within GRASS’s “v.surf.rst” module ([18]). Subsequently, the maps obtained were processed with particular GIS commands to produce the species contribution ratio maps and the mixture map.

Error assessment

The accuracy of the prediction maps computed with RST was calculated on both point and area basis. The MAE (Mean Absolute Error) and the MAD (Mean Absolute Deviation) were determined during the evaluation of the errors in the volume-estimation maps. The “v.sample” command was used to determine the point-wise error, and the “v.rast.stats” command was used to determine the area-wise error. The volumes (in ha) derived from the measured values from the inventory sampling plots in the field were compared with the predicted values at the applicable point on the maps obtained by interpolation. The area-based evaluation considered that each sample point represented precisely 1 ha. In addition, the total area-based volume was calculated and compared with the value obtained from conventional inventory methods.

Production of stand map from interpolation maps

The computed raster maps of the volume per ha for each species (Pb, Pn, Ab, Fo, Q, Cb, Cs, Od) were converted to vector maps and then overlaid to produce a vector map displaying the mixture ratios of the 8 species. The combinations acquired from these maps were evaluated with expert knowledge of forest management. This evaluation was then used to determine the criteria and classes for making the final decision map for the managed forest. For example, in Tab. 2, 8 species and diameter class combinations were categorized into a single category by defined expert knowledge criteria. All other similar combinations were categorized with this approach. This procedure used the ELECTRE TRI “Bouyssou-Marchant” model in Quantum GIS ([52], [57]).

Tab. 2 - ELECTRE TRI categorization sample.

1. 6Pn2Fo1Q1Cb de 6Pn2Fo1Q1Cb de
2. 6Pn2Fo1Q1Cb de
3. 6Pn3Fo1Q de
4. 6Pn2Fo2Q de
5. 6Pn2Fo2Q ed
6. 6Pn2Fo1Q1Cb ed
7. 6Pn1Fo2Q1Cb d
8. 6Pn2Fo2Cb de

  Enlarge/Reduce  Open in Viewer


Determining spatial interpolation parameters

Based on the cross-validation results obtained in this study, the interpolation parameters showing the lowest RMSE value varied among the seven diameter classes considered (Tab. 3). Therefore, a map showing the deviation of the seven most common species by diameter class was calculated. All diameter classes gave the minimum mean absolute-deviation error, with a value of 80 for the tension parameter and a value of 0.2 for the smoothing parameter (Tab. 3). The use of these values for interpolation would likely decrease the error rate.

Tab. 3 - Deviation values for the seven most common species by diameter class.

Cross validation Deviation Smoothing
20 40 60 80 20 40 60 80
Pnb 6.64 6.6 6.6 6.6 0.39 0.26 0.22 0.2 0.2
6.63 6.6 6.6 6.6 0.72 0.49 0.42 0.38 0.4
6.62 6.6 6.6 6.6 1 0.71 0.61 0.55 0.6
6.61 6.6 6.61 6.61 1.26 0.91 0.78 0.72 0.8
Pnc 30.36 30.29 30.35 30.4 1.83 1.23 1.04 0.93 0.2
30.32 30.3 30.36 30.42 3.36 2.32 1.98 1.79 0.4
30.29 30.31 30.38 30.44 4.69 3.32 2.85 2.59 0.6
30.28 30.33 30.4 30.46 5.87 4.23 3.66 3.35 0.8
Pnd 39.66 39.4 39.37 39.37 2.38 1.59 1.34 1.2 0.2
39.58 39.38 39.36 39.37 4.38 3 2.55 2.31 0.4
39.52 39.37 39.36 39.38 6.12 4.3 3.68 3.35 0.6
39.47 39.36 39.36 39.39 7.65 5.49 4.73 4.32 0.8
Pne 43.38 42.92 42.84 42.82 2.73 1.8 1.51 1.35 0.2
43.23 42.88 42.83 42.81 4.98 3.38 2.86 2.58 0.4
43.13 42.86 42.81 42.8 6.92 4.82 4.11 3.73 0.6
43.05 42.84 42.8 42.79 8.63 6.14 5.27 4.8 0.8
Foc 8.5 8.49 8.53 8.57 0.53 0.35 0.29 0.27 0.2
8.48 8.5 8.55 8.58 0.98 0.67 0.57 0.51 0.4
8.47 8.51 8.56 8.59 1.37 0.96 0.82 0.75 0.6
8.46 8.52 8.57 8.61 1.71 1.22 1.06 0.97 0.8
Qb 6.5 6.42 6.41 6.4 0.4 0.26 0.22 0.2 0.2
6.47 6.41 6.4 6.4 0.73 0.5 0.42 0.38 0.4
6.46 6.41 6.4 6.4 1.01 0.71 0.61 0.55 0.6
6.44 6.41 6.4 6.41 1.26 0.9 0.78 0.71 0.8
Qc - - - - 0.78 0.52 0.43 0.39 0.2
- - - - 1.44 0.98 0.83 0.75 0.4
- - - - 2 1.4 1.2 1.09 0.6
- - - - 2.5 1.78 1.54 1.4 0.8

  Enlarge/Reduce  Open in Viewer

Prediction of forest maps

A total of 31 independent interpolation maps were obtained for the diameter classes of each species (e.g., black pine in Fig. 3a, b, c, and d). The total volumes of particular species per ha (e.g., black pine: Pnb + Pnc + Pnd + Pne = Pnt) were computed by summing these species-diameter-class interpolation maps (Fig. 4a). The percentage contributions of each diameter class to the mixture were obtained by calculating the ratios of the diameter class volume of the particular species and the total volume of the species (e.g., Pnb/Pnt) with the help of this map (Fig. 4b). Similarly, the summation of the total volumes of each species per ha (Pnt + Pbt + Odt + Abt + Cbt + Fot + Cst + Qt = Vt) yielded the total volume per ha of the trees in the planning units (Fig. 4c). The percentage contributions of each species to the mixture were obtained by calculating the ratios of the total volume of the species and the overall total volume (e.g., Pnt/Vt) with the help of this map (Fig. 4d).

Fig. 3 - Interpolation maps of black pine. (A) diameter class; (B) diameter class ; (C) diameter class ; (D) diameter class

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 4 - (A) Total volume map for black pine; (B) volume ratio map for black pine (diameter class b); (C) total volume map for all the species; (D) volume ratio map for black pine.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Error assessment of predictions

If the predictions were evaluated on the basis of the size of one pixel (400 m2), the results were very satisfactory (MAD < 4%). However, on the basis of a sample plot that was approximately one ha in area, the maximum MAD increased to 12% (Pnb) and 14% (Qb). The maximum standard deviations for the black-pine diameter classes d and e, were 81% and 85%, respectively (Tab. 1). High standard deviations were also found for the Calabrian-pine diameter classes d and e (67% and 79%, respectively), Oriental beech (51% and 70%, respectively), and oak spp. (31% and 45%, respectively).

The total volume by species as determined by conventional methods ([15]) and was compared to that determined by the methods applied in our study for the entire study area (Tab. 4). The data over the entire area showed that four species had volume ratios greater than 1%. In descending order of volume ratios, these species were Pn > Q > Pb > Fo, with volume ratio values of 55%, 19%, 16%, and 8%, respectively. These four species represented a total volume of 98%. The maximum relative errors were calculated for the Calabrian pine diameter classes b and c. These relative error values were +33% and +23%, respectively. The minimum relative errors for these four species were calculated for diameter classes b, c, and d of the oak species. These relative error values were 2%, 4%, and 2%, respectively.

Tab. 4 - Accuracy of total volume calculated with spatial interpolation (SI), compared with conventional forest inventory results (CFIR).

Type of
Diameter Total
Total partecipant
Class b Class c Class d Class e
Pb SI 36226 165430 248434 220339 670428 0.16
CFIR 27332 134388 211765 180807 554292
% 33 23 17 22 21
Pn SI 146309 611316 695638 559590 2012853 0.55
CFIR 160647 529222 614516 625805 1930190
% -9 16 13 -11 4
Ab SI 5507 15973 20363 5326 47169 0.01
CFIR 4273 12838 16481 6034 39626
% 29 24 24 -12 19
Fo SI 34856 141304 95689 52129 323978 0.08
CFIR 31873 128509 80135 44008 284525
% 9 10 19 18 14
Q SI 132299 233543 176357 150154 692353 0.19
CFIR 132519 228881 170362 147514 679276
% 0 2 4 2 2
Cb SI 1111 1023 710 0 2845 0
CFIR 1101 871 576 0 2548
% 1 17 23 - 12
Cs SI 2367 7859 3324 4979 18529 0
CFIR 2092 7351 2759 2035 14237
% 13 7 20 145 30
Od SI 13044 9291 1333 785 24452 0.01
CFIR 12228 7363 2414 257 22262
% 7 26 -45 205 10

  Enlarge/Reduce  Open in Viewer

Black pine had a maximum volume ratio of 55% for the entire study area and a total relative error of 4%. The relative error rates for the individual diameter classes were 16% for class c, 13% for class d and -11% for class e.

Oriental beech has the lowest volume ratio (8%) of the four highest-volume species. The relative error for the total volume was 14%. The relative errors for the diameter classes were 10%, 19%, and 18% for classes c, d, and e, respectively.

Generally, the highest relative errors in volume estimation occurred for diameter classes c and d.

The pixel-based (sampling area-based) MAE and the area-based MAE followed similar trends. The pixel-based MAD values ranged between 0% and 3%.

When producing the forest stand map from the interpolation maps, the diameter class for each value of species volume derived from the raster maps was uploaded with the “v.rast.stats” command onto the vector map as attribute data. The resulting species mixtures, including both volume-based dominants and accompanying secondary species, were then considered. In the next stage of the process, the 890 different mixture variations were reduced to a smaller number of categories based on expert knowledge. Criteria were established based on the mixture rates of the tree species and the average volumes for the diameter classes. These ELECTRE TRI criteria and the values for the weights were then entered, and 70 different profiles were generated. The ELECTRE TRI module of QGIS was used to create a final decision map. This map represents 70 different profile categories based on species mixture rates and diameter classes (Fig. 5).

Fig. 5 - Stand type maps obtained using the ELECTRE TRI approach by combining the 31 prediction maps of the volume per ha.

  Enlarge/Shrink   Download   Full Width  Open in Viewer


Diameter classes d and e showed the larger standard devistions for species representing 98% of the total volume (Tab. 1). In general, the observed error values were as low as expected from the associated standard deviations, with occasional exceptions. Minor differences were found between MAD errors at the pixel level and at the 1-ha level.

The observed differences in MAE among diameter classes for the same species could be interpreted as due to variation in site index between neighboring plots, that could cause the volume to increase by approximately 15%. Moreover, differences in canopy closure between neighboring stands could have also contributed to the observed differences in MAE described above. Therefore, MAE could be reduced by including the above parameters in the volume estimation, although the number of maps to be produced would also increase.

In the species diameter class with approximately 300 observations distributed around a value of about 300, the MAD values calculated on the basis of 1 ha were 8% and 9% for Qd and Pne, respectively. Increasing the number of species found at the sampling points (Qb, Pnc, Pnb, Pnd, and Qc), the same trend was observed. This increase in the error values is not affected by standard deviation or mean volume (Tab. 1). These maximum MAD values could result from the distribution of black pine and oak species in the study area. Indeed, data heterogeneity may fairly affect the performance of the methods ([38]).

The design of the RST method adopted in this study ensures that the predicted values at a given points remain very close to the observed values at that point. Thus, the MAD of the sampling plots (400 m2 = 1 pixel) did not increase at the ha level (10 000 m2/400 m2=25 pixels). This result reflects the characteristics of the inventoried sample plots, whose structure was fairly similar to the forest structure in the surrounding area.

The volumes obtained by the RST method were slightly higher than the corresponding values in the forest management plan data (Tab. 4). Black pine showed the highest volume ratio (55%), with a relative error of 4% only. The largest error value was detected for diameter class b of Calabrian pine. The observed differences can be accounted for by differences between the two methods. The classical method adopted in forest management plans evaluates the data groups separately to create separate stand types. Instead, the spatial interpolation approach used here evaluates the data relative to the overall configuration of the observations and the spatial scales. Different results are also expected because of the subjective personal evaluations involved in conventional approaches.

According to the 2008 Turkish Forest Management Regulation Handbook ([15]), an error ≤ 10% in volume assessment would be expected for production forests. However, unlike production forests, an error rate > 10% is accepted for protected forests because each systematic sampling point is 600 x 600 m = 36 000 m2. Thus, based on the above criteria, the results of this study could be considered as acceptable.

In this study, a final decision map has been produced representing 70 different profile categories based on species mixture rates and diameter classes in the area analyzed (Fig. 5), obtained combining the prediction maps produced. The method adopted has facilitated the reductions performed for adjacent diameter classes. In addition, the small forest fragments were integrated into a general framework by the analysis. Therefore, this approach facilitated the reduction of small fragments until the structural and statistical features of the desired stand types were attained, in accordance with management goals. To perform this reduction process in a objective way, the ELECTRE TRI method was used. Aktas & Yilmaz ([2]) also used the same methodology to produce a stand map from interpolation maps. The above process would require the use of expert knowledge on forest management to make an improved decision map.

The area considered in this study enlarged upon both forested and non-forested zones (e.g., agricultural areas, settlements, and grasslands). To avoid the effect of the latter areas on total volume assessment, adequate masks were used to exclude non-forested area from the analysis. An additional problem encountered is the accuracy of predictions near the boundaries of the forest management units. The incorporation of neighboring plan units in these areas is an appropriate way of improving the accuracy of predictions. Nevertheless, the resulting map should be refined through aerial photography or satellite images.

Spatial interpolation methods (especially kringing) have been previously argued as unsuitable for assessment purposes in managed forests ([1]). On the other hand, several studies have reported successful application of such methods in homogenous natural forest areas ([20], [62]).

The RST method adopted here allows the user to create stand maps very efficiently. In order to test the accuracy of volume predictions, cross-validation and resampling methods have been applied ([47], [66], [50], [58], [33], [60], [63]), aimed at comparing the observed value at a given point with its prediction obtained from the interpolated continuous surface. However, in areas characterized by highly heterogeneous forest structures, the above approach may not be adequate, and additional control points in the field may be required. On the contrary, a subset of sampling points may be used in areas densely sampled for the above comparisons, keeping the points excluded from analysis for control purposes.

Further improvement of the map obtained may be achieved by the inclusion of auxiliary data (e.g., soil, climate, topography etc.) and the use of multivariate interpolation methods or mixed models, thus combining the most important parameters of the forest stands. For example, several environmental parameters have been used to model the spatial variability in Pinus radiata productivity in New Zealand ([32]). Moreover, remote sensing data may be incorporated into the analysis to obtain an improved accuracy of predictions. Dubayah et al. ([9]) found that the LiDAR remote sensing has a vast potential for the direct estimation of several key forest characteristics, such as canopy heights, stand volume, basal area and aboveground biomass ([8]). Lovell et al. ([39]) have performed studies on optimal LiDAR acquisition parameters for forest height retrieval. Estornell et al. ([12]) have made estimations of shrub biomass by airborne LiDAR in small forest stands, obtaining accurate measures of the most important plant parameters (volume, average height, diameter and canopy) with minimal costs. Haywood ([21]) adopted a semi-automated stand delineation process in the mapping of natural eucalypt forests using remote sensing, multi-spectral imagery, lasers and stand characteristics. Furthermore, the application of LiDAR allows the creation of accurate three-dimensional relief models of forest stands. For example, Mitasova et al. ([45]) used a simultaneous spline approximation and topographic analysis for LiDAR elevation data in open-source GIS to create accurate relief models. LiDAR, in combination with inventory data, has also been used to evaluate the forest habitats of birds and other wildlife ([17], [61], [65], [56]).

Forest management activities depends on international, national, and regional forestry strategies. Moreover, regional planning should include community participation and the analysis of multiple functions of forest resources ([51]). For this reason, multi-function criteria are often used in forest planning in addition to conventional management criteria, and MCDA has been carried out in different decision-making contexts in forestry. Khadka & Vacik ([31]) used multi-criteria analysis to support community forest management including 6 criteria and 40 profiles. Huth et al. ([26]) used multi-criteria analysis to support decisions on harvest volume based on the criteria of rotation, harvest intensity, and target diameter.

In many forestry studies, MCDA was applied using multiple criteria for the assessment of a single parameter (e.g., total volume, site class). Similarily, the MCDA map obtained in our study has also multiple criteria and allows the evaluation of multiple functions. The conventional methods used in Turkey require professional experience and a background of appropriate knowledge. The use of an MCDA method, such as the ELECTRE TRI adopted here, may bypass the above problem.


Akhavan R, Amiri GhZ, Zobeiri M (2010). Spatial variability of forest growing stock using geostatistics in the Caspian region of Iran. Caspian Journal Environmental Sciences 8: 43-53.
Online | Gscholar
Aktas S, Yilmaz OY (2011). Producing of stand draft maps with spatial prediction methods. Journal of the Faculty of Forestry, Istanbul University, Turkey 62(2):131-146.
Anonymous (2011a). Forest Management Plan Forestry Asar 2010-2019. Regional Directorate of Forests in Çanakkale, General Directorate of Forests, Ministry of Environment and Forests, Ankara, Turkey.
Anonymous (2011b). Forest Management Plan Forestry Yenice 2010-2019. Regional Directorate of Forests in Çanakkale, General Directorate of Forests, Ministry of Environment and Forests, Ankara, Turkey.
Aurbi P, Debouzie D (2000). Geostatistical estimation variance for the spatial mean in two dimensional systematic sampling. Ecology 81: 543-553.
CrossRef | Gscholar
Balteiro LD, Romero C (2008). Making forestry decisions with multiple criteria: a review and an assessment. Forest Ecology and Management 255: 3222-3241.
CrossRef | Gscholar
Burrough PA (2001). GIS and geostatistics: essential partners for spatial analysis. Environmental and Ecological Statistics 8: 361-377.
CrossRef | Gscholar
Dubayah RO, Drake JB (2000). Lidar remote sensing for forestry applications. Journal of Forestry 98 (6): 44-52.
Dubayah R, Knox R, Hofton M, Blair JB, Drake J (2000). Land surface characterization using Lidar remote sensing. In: “Spatial information for land use management” (Hill MJ, Aspinall RJ eds). International Publishers Direct, Singapore, pp. 25-38.
D’Orazio M (2003). Estimating the variance of the sample mean in two-dimensional systematic sampling. Journal of Agricultural, Biological and Environmental Statistics 8: 280-295.
CrossRef | Gscholar
Efron B, Gong G (1983). A leisurely look at the bootstrap, the jackknife, and cross-validation. The American Statistician 37: 36-48.
CrossRef | Gscholar
Estornell J, Ruiz LA, Velázquez-Martí B, Fernández-Sarría A (2011). Estimation of shrub biomass by airborne LiDAR data in small forest stands. Forest Ecology and Management 262: 1697-1703.
CrossRef | Gscholar
Feng Y, Tang S, Li Z (2006). Application of improved sequential indicator simulation to spatial distribution of forest types. Forest Ecology and Management 222: 391-398.
CrossRef | Gscholar
Flores LA, Martinez L I, Ferrer CM (2003). Systematic sample design for the estimation of spatial means. Environmetrics 14: 45-61.
CrossRef | Gscholar
Forest Management Regulation Handbook (2008). General Directorate of Forests, Ministry of Environment and Forests in Turkey. Official Gazette of the Republic of Turkey No: 26778, pp. 18.
General Directorate of Forests (2006). Forests in Turkey. Ministry of Environment and Forests in Turkey, Ankara, pp. 160.
Graf RF, Mathys L, Bollmann K (2009). Habitat assessment for forest dwelling species using LiDAR remote sensing: Capercaillie in the Alps. Forest Ecology and Management 257: 160-167.
CrossRef | Gscholar
GRASS Development Team (2011). Geographic resources analysis support system (GRASS) Software. Open Source Geospatial Foundation Project.
Online | Gscholar
Guibal D (1973). Estimating oukoumés Gabon. Internal Memo #333, Center Morphology, Mathematics, Fontainbleau, France.
Gunnarsson F, Holm S, Holmgren P, Thuresson T (1998). On the potential of kriging for forest management planning. Scandinavian Journal of Forest Research 13: 237-245.
CrossRef | Gscholar
Haywood A (2011). Semi-automating the stand delineation process in mapping natural eucalypt forests. Australian Forestry 74 (1): 13-22.
CrossRef | Gscholar
Heikkinen J (2006). Assement of uncertainty in spatially systematic sampling. Forest inventory, methodology and applications. Managing Forest Ecosystems 10 (1): 155-176.
CrossRef | Gscholar
Hengl T (2007). A practical guide to geostatistical mapping of environmental variables. EUR 22904 EN Scientific and Technical Research Series no. 143, Office for Official Publications of the European Communities, Luxemburg, Belgium.
Hofierka J, Parajka J, Mitasova H, Mitas L (2002). Multivariate interpolation of precipitation using regularized spline with tension. Transactions in GIS 6: 135-150.
CrossRef | Gscholar
Holmgren P, Thuresson T (1997). Applying objectively estimated and spatially continuous forest parameters in tactical planning to obtain dynamic treatment units. Forest Science 43: 317-326.
Huth A, Drechsler M, Köhler P (2005). Using multi-criteria decision analysis and a forest growth model to assess impacts of tree harvesting in Dipterocarp lowland rain forests. Forest Ecology and Management 207: 215-232.
CrossRef | Gscholar
Höck BK, Payn TW, Shirley JW (1993). Using a geographic information system and geostatistics to estimate site index of Pinus radiata for Kangaroo Forest, New Zealand. N. Z. J. Forrest. Science 23: 264-277.
Jost A (1993). Geostatistical analysis of sampling error for systematic sampling. Ph.D. thesis, University of Freiburg,Breisgau, Germany, pp. 113. [In German]
Kangas J, Kangas A (2005). Multiple criteria decision support in forest management-the approach, methods applied, and experiences gained. Forest Ecology and Management 207: 133-143.
CrossRef | Gscholar
Kangas A, Kangas J, Pykäläinen J (2001). Outranking methods as tools in strategic natural resources planning. Silva Fennica 35 (2): 215-227.
Online | Gscholar
Khadka C, Vacik H (2012). Use of multi-criteria analysis (MCA) for supporting community forest management. iForest 5: 60-71.
CrossRef | Gscholar
Kirschbaum MU, Watt MS (2011). Use of a process-based model to describe spatial variation in Pinus radiata productivity in New Zealand. Forest Ecology and Management 262: 1008-1019.
CrossRef | Gscholar
Kovacs JM, Liu Y, Zhang C, Flores-Verdugo F, de-Santiago FF (2011). A field based statistical approach for validating a remotely sensed mangrove forest classification scheme. Wetlands Ecology and Management 19: 409-421.
CrossRef | Gscholar
Köhl M, Magnussen S, Marchetti M (2006). Sampling methods, remote sensing and GIS multiresource forest inventory. In: “Remote Sensing” (chapt. 4). Springer Verlag, Heidelberg, Germany.
Jayanath A, Gamini H (2009). A critical review of multi-criteria decision making methods with special reference to forest management and planning. Ecological Economics 68(10): 2535-2548.
Lappi J, Kangas A (2006). Use of additional information. In: “Forest inventory - Methodology and Application” (chapt. 7). Springer, The Hague,the Netherlands. pp. 107-117.
Li J, Heap AD (2008). A review of spatial interpolation methods for environmental ecientists. No. Record 2008/23, Geoscience Australia, Canberra, Australia.
Li J, Heap AD (2011). A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors. Ecological Informatics 6: 228-241.
CrossRef | Gscholar
Lovell JL, Jupp DJB, Newnham GJ, Coops NC, Culvenor DS (2005). Simulation study for finding optimal Lidar acquisition parameters for forest height retrieval. Forest Ecology and Management 214: 398-412.
CrossRef | Gscholar
Magnussen S, Næsset E, Wulder MA (2009). Efficient multiresolution spatial predictions for large data arrays. Remote Sensing of Environment 109: 451-463.
CrossRef | Gscholar
Mandallaz D (1991). A unified approach to sampling theory for forest inventory based on infinite population and super population models. Chair of forest management and planning, Swiss Federal Institute of Technology, Zurich, Switzerlands.
CrossRef | Gscholar
Mendoza GA, Martins H (2006). Multi-criteria decision analysis in natural resource management: a critical review of methods and new modeling paradigms. Ecology and Forest Management 230: 1-206.
CrossRef | Gscholar
Mitas L, Mitasova H (1999). Spatial interpolation. In: “Geographical information systems: principles, techniques, management and applications” (Longley P, Goodchild MF, Maguire DJ, Rhind DW eds). GeoInformation International, Wiley, pp. 481-492.
Mitasova H, Mitas L, Brown WM, Gerdes DP, Kosinovski I, Baker T (1995). Modeling spatially and temporally distributed phenomena: new methods and tools for GRASS GIS. International Journal of Geographical Information Systems 9: 433.
CrossRef | Gscholar
Mitasova H, Mitas L, Harmon RS (2005). Simultaneous spline approximation and topographic analysis for lidar elevation data in open-source GIS. Geoscience and Remote Sensing Letters IEEE 2: 375-379.
CrossRef | Gscholar
Montes F, Alicia F, Ledo A (2010). Incorporating environmental and geographical information in forest data analysis: a new fitting approach for universal kriging. Canadian Journal of Forest Research 40 (9): 1852-1861.
CrossRef | Gscholar
Nanos N, Calama R, Montero G, Gil L (2004). Geostatistical prediction of height/diameter models. Forest Ecology and Management 195: 221-235.
CrossRef | Gscholar
Neteler M, Mitasova H (2008). Open Source GIS: a GRASS GIS approach. Kluwer. Academic Publishers/Springer, Boston, USA. pp.406.
Neteler M, Bowman MH, Landa M, Metz M (2012). GRASS GIS: a multi-purpose open source GIS. Environmental Modelling and Software 31: 124-130.
CrossRef | Gscholar
Palmer DJ, Höck BK, Kimberley MO, Watt MS, Lowe DJ, Payn TW (2009). Comparison of spatial prediction techniques for developing Pinus radiata productivity surfaces across New Zealand. Forest Ecology and Management 258: 2046-2055.
CrossRef | Gscholar
Perez-Soba M, Edwards D, Tabbush P (2006). Assessing the impacts of EU policies on sustainable forest management: the sensor approach. In: EFI Scientific Seminar “The role of Forestry in Integrated Environmental Assessment”. , The Netherlands.
Quantum GIS Development Team (2012). Quantum GIS geographic information system. Open Source Geospatial Foundation Project.
Online | Gscholar
Schaffer C (1993). Selecting a classification method by cross-validation. Machine Learning 13 (1): 135-143.
CrossRef | Gscholar
Shao J (1993). Linear model selection by cross-validation. Journal of the American Statistical Association 88: 486-494.
CrossRef | Gscholar
Sherman M (1996). Variance estimation for statistics computed from spatial lattice data. Journal of the Royal Statistical Society 58: 509-523.
Online | Gscholar
Smart LS, Swenson JJ, Christensen NL, Sexton JO (2012). Three-dimensional characterization of pine forest type and red-cockaded woodpecker habitat by small-footprint, discrete-return Lidar. Forest Ecology and Management 281:100-110.
CrossRef | Gscholar
Sobrie O (2010). Integration of ELECTRE TRI in a GIS coupling with a XMCDA web service for inference. In: Proceedings of the 73 Meeting of the European Working Group “Multiple Criteria Decision Aiding”. 8 Decision Deck Workshop. University of Mons (Belgium) 13 April 2011.
Online | Gscholar
Tang D, Bian F (2009). Forest site evaluation based on GIS and kriging. In: ICISE ’09 Proceedings of the “2009 First IEE International Conference on Information Science and Engineering”. Nanjing, Jiangsu China, December 26-December 28, IEEE Computer Society Washington, DC, USA, pp. 2063-2067.
Tattoni C, Ciolli M, Ferretti F, Cantiani MG (2010). Monitoring spatial and temporal pattern of Paneveggio forest (northern Italy) from 1859 to 2006. iForest 3 (1): 72-80.
CrossRef | Gscholar
Tattoni C, Ciolli M, Ferretti F (2011). The fate of priority areas for conservation in protected areas: a fine-scale Markov chain approach. Environmental Management 47: 263-278.
CrossRef | Gscholar
Tattoni C, Rizzolli F, Pedrini P (2012). Can LiDAR data improve bird habitat suitability models? Ecological Modelling 245: 103-110.
CrossRef | Gscholar
Tuominen S, Fish S, Poso S (2003). Combining remote sensing, data from earlier inventories and geostatistical interpolation in multi-source forest inventory. Canadian Journal of Forest Research 33: 624- 634.
CrossRef | Gscholar
Viana H, Aranha J, Lopes D, Cohen WB (2012). Estimation of crown biomass of Pinus pinaster stands and shrubland above-ground biomass using forest inventory data, remotely sensed imagery and spatial prediction models. Ecological Modelling 226: 22-35.
CrossRef | Gscholar
Wallerman J, Joyce S, Vencatasawmy CP, Olsson H (2002). Prediction of forest stem volume using kriging adapted to detected edges. Canadian Journal of Forest Research 32: 509-518.
CrossRef | Gscholar
Wilsey CB, Joshua J, Lawler JJ, David A, Cimprich DA (2012). Performance of habitat suitability models for the endangered black-capped vireo built with remotely-sensed data. Remote Sensing of Environment 119:35-42.
CrossRef | Gscholar
Zawadzki J, Cieszewski CJ, Zasada M, Lowe RC (2005). Applying geostatistics for investigations of forest ecosystems using remote sensing imagery. Silva Fennica 39 (4): 599-617.
Online | Gscholar

Authors’ Affiliation

S Destan
Department of Forest Management, Faculty of Forestry, Istanbul University, 34473 Bahçeköy, Istanbul (Turkey)
O Yilmaz
Department of Forest Engineering, Faculty of Forestry, Istanbul University, 34473 Bahçeköy, Istanbul (Turkey)
A Sahin
Department of Forest Management, Regional Forest Directorate of Istanbul, Istanbul (Turkey)

Corresponding author


Destan S, Yilmaz O, Sahin A (2013). Making objective forest stand maps of mixed managed forest with spatial interpolation and multi-criteria decision analysis. iForest 6: 268-277. - doi: 10.3832/ifor0099-006

Academic Editor

Agostino Ferrara

Paper history

Received: May 18, 2012
Accepted: Mar 18, 2013

First online: Jul 01, 2013
Publication Date: Oct 01, 2013
Publication Time: 3.50 months

© SISEF - The Italian Society of Silviculture and Forest Ecology 2013

  Open Access

This article is distributed under the terms of the Creative Commons Attribution-Non Commercial 4.0 International (https://creativecommons.org/licenses/by-nc/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Creative Commons Licence

Breakdown by View Type

(Waiting for server response...)

Article Usage

Total Article Views: 50982
(from publication date up to now)

Breakdown by View Type
HTML Page Views: 42353
Abstract Page Views: 2632
PDF Downloads: 4672
Citation/Reference Downloads: 23
XML Downloads: 1302

Web Metrics
Days since publication: 4033
Overall contacts: 50982
Avg. contacts per week: 88.49

Article citations are based on data periodically collected from the Clarivate Web of Science web site
(last update: Feb 2023)

Total number of cites (since 2013): 7
Average cites per year: 0.64


Publication Metrics

by Dimensions ©

List of the papers citing this article based on CrossRef Cited-by.


iForest Similar Articles


This website uses cookies to ensure you get the best experience on our website. More info