Optimized cellular automaton for stand delineation

Forest inventories based on remote sensing often interpret stand characteristics for small raster cells instead of traditional stand compartments. This is the case for instance in the Lidar-based and multi-source forest inventories of Finland where the interpretation units are 16 m × 16 m grid cells. Using these cells as simulation units in forest planning would lead to very large planning problems. This difficulty could be alleviated by aggregating the grid cells into larger homogeneous segments before planning calculations. This study developed a cellular automaton (CA) for aggregating grid cells into larger calculation units, which in this study were called stands. The criteria used in stand delineation were the shape and size of the stands, and homogeneity of stand attributes within the stand. The stand attributes were: main site type (upland or peatland forest), site fertility, mean tree diameter, mean tree height and stand basal area. In the CA, each cell was joined to one of its adjacent stands for several iterations, until the cells formed a compact layout of homogeneous stands. The CA had several parameters. Due to high number possible parameter combinations, particle swarm optimization was used to find the optimal set of parameter values. Parameter optimization aimed at minimizing within-stand variation and maximizing between-stand variation in stand attributes. When the CA was optimized without any restrictions for its parameters, the resulting stand delineation consisted of small and irregular stands. A clean layout of larger and compact stands was obtained when the CA parameters were optimized with constrained parameter values and so that the layout was penalized as a function of the number of small stands (< 0.1 ha). However, there was within-stand variation in fertility class due to small-scale variation in the data. The stands delineated by the CA explained 66–87% of variation in stand basal area, mean tree height and mean diameter, and 41–92% of variation in the fertility class of the site. It was concluded that the CA developed in this study is a flexible new tool, which could be immediately used in forest planning.


Introduction
Forest inventory methods based on remote sensing often produce inventory results for small square-shaped pixels.This is the case for instance when the inventory is based on the interpretation of satellite images.When the inventory uses laser scanning and the area-based interpretation approach (Naesset 2002;Maltamo and Packalen 2014;Vauhkonen et al. 2014), the optimal size of the interpretation unit is of the same magnitude as the size of the field plots that are used as ground truth (Pippuri et al. 2013).Commonly used field plot sizes range from 250 to 500 m 2 , corresponding to pixel sizes of about 16 m 9 16 m (256 m 2 ) to 22 m 9 22 m (484 m 2 ).For example, in the Finnish Lidar-based forest inventory, stand variables are predicted for 16 m 9 16 m pixels.Lidar inventories cover most of the privately owned forests of Finland (www.metsakeskus.fi).The results of these inventories are freely available at www.metsaan.fi/paikkatietoaineistot.
Another source of free forest inventory data is the multisource National Forest Inventory of Finland (Tomppo et al. 2008).Also in this inventory, site and stand characteristics are predicted for 16 m 9 16 m pixels, which cover the whole country.The data source is available in the internet for free download (www.paikkatietoikkuna.fi).
These forest inventories are excellent datasets for forest management planning because the data are spatially continuous and include all variables required in planning calculations (Maltamo and Packalen 2014).The essential information for management planning includes a few site characteristics, which in Finland are the main type of the site (mineral soil or peatland), fertility class, and temperature sum (accumulated temperature during the growing season).The required growing stock characteristics are mean tree diameter, mean height and stand density (basal area or number of trees per hectare).
Modern forest planning consists of simulation and optimization steps (e.g., Zubizarreta-Gerendiain et al. 2018).The pixels or grid cells for which the stand variables are predicted, can be used as simulation units in these calculations.However, pixels and cells are too small to be used as treatment units in actual management (e.g.harvest blocks).Fortunately, optimization offers possibilities to aggregate treatments into large enough continuous areas, called as dynamic treatment units (DTUs).Usually, various heuristics are used for this purpose (e.g., Bettinger et al. 1999Bettinger et al. , 2002;;Heinonen and Pukkala 2004).Mathematical programming is often ruled out since the number of cells may be tens or even hundreds of thousands (O ¨hman 2000;O ¨hman and La ¨ma ˚s 2003;Heinonen et al. 2018).
Heuristic methods available for spatial optimization fall into two categories: centralized and decentralized methods (Pukkala et al. 2014).Centralized heuristics may be categorized as global methods, since the effect of changing the solution during the optimization process is evaluated ''globally'', i.e., at the level of the whole planning area.Decentralized or local methods compare treatment alternatives at the pixel level.Therefore, they are faster than global heuristics in solving large spatial forest planning problems (Pukkala et al. 2008(Pukkala et al. , 2014)).
Even when heuristic methods are used, planning calculations (both simulation and optimization) become too time consuming if the number of calculation units is very large.Global heuristics may become impractical already with a few tens of thousands of calculation units (Heinonen et al. 2018).However, real-life planning problems may include much larger datasets.For example, if the planning area is 10,000 ha (10 km 9 10 km) and the size of the pixel is 256 m 2 (16 m 9 16 m), the number of pixels is 390,625.
To alleviate the problems arising from a high number of small spatial units, segmentation methods have been developed to aggregate these units into larger segments (Pekkarinen 2002;Pekkarinen and Tuominen 2003;Mustonen et al. 2008;Wulder et al. 2008;Koch et al. 2009Koch et al. , 2014;;Dechesne et al. 2017).The segments, which are spatially continuous and homogeneous forest patches, are then used as calculation units in forest planning.Segmentation may use interpreted stand variables, ALS (airborne laser scanning) pulse data, intensities of different wavelengths of electromagnetic radiation, external data such and soil and road maps, as well as combinations of different data sources and types.
Many segmentation methods have been suggested in the literature, but only a few are commonly used in forestry (Wulder et al. 2008).The most common category of segmentation methods might be region merging, in which adjacent units are gradually merged into larger areas according to their similarity, until a stopping criterion is reached.The initial units may be small raster cells (Mustonen et al. 2008) or Voronoi polygons, each of which represents an individual tree (Olofsson and Holmgren 2014).The shape of the segments can be taken into account in some of the region-merging algorithms (Mustonen et al. 2008).The size of the final segments can be controlled through the stopping criterion.
Also other types of segmentation methods have been presented in forestry literature.For example, Koch et al. (2009) developed a method in which grid cells of 20 m 9 20 m are first classified based of forest type, roughness of canopy surface and tree height.Then, the adjacent cells of the same class are grouped into forest stands.Wu et al. (2013) describe a method that uses the mean shift algorithm to generate ''raw stands'', which are subsequently refined by using the spectral clustering algorithm.Dechesne et al. (2017) provide a review of segmentation methods for airborne LIDAR data and high resolution multispectral imagery.
The use of segmentation methods needs expertise and special software.An alternative to the commonly used segmentation methods are cellular automata (CA).They are transparent and flexible systems, which are easy to program, use, control and understand.CA are self-organizing systems, which have been used for many purposes (Von Neumann 1966;Wolfram 2002).In forestry, CA have been used for land allocation and planning (Strange et al. 2001;Heinonen and Pukkala 2007), and for simulating the spread of pests and pathogens (Mo ¨ykkynen and Pukkala 2014;Mo ¨ykkynen et al. 2015Mo ¨ykkynen et al. , 2017)).
Although CA are simple, they have parameters, which affect the result of the CA run.Therefore, an analysis on the effects of parameters that guide the CA process is required, in addition to programming the CA tool.When small pixels are aggregated into homogeneous forest segments or stands, common criteria are the shape of the stands, their size, and the homogeneity of forest attributed within stands.The usual aim is to have a stand delineation with small within-stand and large between-stand variation in stand attributes.In addition, irregular stand shapes are avoided, as well as very small stands.
The importance of the above criteria depends on the purpose of aggregation; shape of the segments is not important if the main purpose is to reduce the number of calculation units for scenario analyses, yield predictions, etc.In this case, the segments correspond to strata, which do not need to be continuous.On the other hand, shape is important if the segments are used in forest management to implement treatments.Round, compact, and not too irregular segments are targeted.The segments should not be too small, but not too large either.If the segments are regarded as non-divisible management units, large segments may reduce the possibilities of planning to optimize the use of forest resources.For example, an even-flow constraint and a low number of large segments may be a problematic combination.
This study developed a CA for stand demarcation so that all the main criteria of good stand delineation were taken into account, namely shape, size and homogeneity of the created stands.In addition to describing the CA, the effects of its parameters on stand delineation were analyzed.Since the number of parameters and their combinations is high, numerical optimization was used to find parameters that result in good stand delineation.

Cellular automaton
Cellular automata are self-organizing algorithms based on the assumption that the interaction between cells decreases rapidly with increasing distance (Strange et al. 2001).Each cell takes one of a limited number of states, which can be management schedules, land uses or, in the case of the current study, stand numbers.In stand delineation, the purpose of CA is to find the optimal stand number for each cell.The cell values evolve in discrete time steps according to a set of rules.These rules only consider the local neighborhood of the cell.
The starting points of the cellular automaton for stand delineation were the studies of Strange et al. (2001Strange et al. ( , 2002) ) and Heinonen and Pukkala (2007).Strange et al. (2001Strange et al. ( , 2002) ) compared six different variants of cellular automata and found that the variant that included the socalled innovation and mutation operations, was the best.Innovation is equal to selecting the best state for the cell, and mutation is equal to drawing the cell state randomly from the set of available alternatives.Both mutation and innovation are applied with certain probability, which means that that this variant of CA is not deterministic.Heinonen and Pukkala (2007) used the CA proposed by Strange et al. (2001Strange et al. ( , 2002) ) for spatial optimization in the context of forest planning.They found that the best results are obtained when the probability of innovation is high, and the probability of mutation is low.Based on this result, the algorithm of the current study was simplified so that the probability of innovation was equal to 1 at every iteration and the probability of mutation was constantly zero.As a result, the CA was deterministic so that the best stand number was assigned for each cell at every iteration.Determinism is an advantage when the parameters of the algorithm are optimized, because the CA needs to be run only once with each tested parameter combination.
In the CA developed in this study, all cells were given an initial stand number.This was accomplished by dividing the area into square-shaped initial stands (two hectares in this study).Then, the best stand number was assigned for each cell for several iterations.The ranking of alternative stand numbers was based on three criteria: (1) common border between the cell and the stand (2) stand area, and (3) similarity of stand characteristics in the cell and the stand (Fig. 1).
The ranking of alternative stand numbers was calculated with the following formula, which has the same form as the general additive utility function (Kangas et al. 2008): where U ij is the score (utility) of assigning the number of stand j to cell i, B ij is the length of common border between cell i and stand j, A j is the area of stand j, and D ij is the difference in stand characteristics between cell i and stand j.
The score was calculated only for the eight immediate neighbors of the cell; the cell took the stand number of one of its neighbors.In this article, immediate neighbor cells to east, west, north and south are called adjacent neighbors.Cells to northeast, northwest, southeast and southwest are called corner cells.A cell has eight neighbors, four of which are adjacent and four are corner cells.
With a square lattice, the corner cells have only one common point with the cell.To make it possible to control the effect of corner cells on the CA, corner cells were given a weight between 0 and 1.The weight of adjacent cells was always 1.The ''border length'' B ij of Eq. 1 was calculated as the sum of the weights of the eight neighbors: where v k is the weight of neighbor k and J is the number of stands.Variable s k indicates whether neighbor cell k belongs to stand j (s k = 1, otherwise s k = 0).The difference of stand characteristics between cell i and stand j was calculated from: where w r is the weight of stand variable r (Rw r = 1), q ri is the value of standardized variable r in cell i and q rj is the mean value of the same variable in stand j.All variables included in the calculations were standardized to mean zero and standard deviation 1.This removed the effect of different units of variables.
Since each cell was joined to one of its neighbors, all one-cell-wide stands disappeared at every iteration.The stand number assigned to the cell depended on the weights (a 1 , a 2 , a 3 ) of the three criteria in Eq. 1.In the example shown in Fig. 1, the evaluated cell (white cell in Fig. 1) will be joined either to Stand 1 or Stand 2, depending of the weights of the criteria.If the weight of common border is high, the cell will be joined to Stand 1, otherwise it will be joined to Stand 2.
The variables used to measure the similarity of stand characteristics were: main site type (mineral soil, spruce mire, pine bog, open bog), fertility class (mesotrophic, herb-rich, mesic, sub-xeric, xeric), mean tree diameter, mean tree height, and stand basal area.These are the variables that are usually assessed in the field in management-oriented forest inventories and imported to forest planning systems.
Each of the three criteria of Eq. 1 had an associated suppriority function (u 1 , u 2 and u 3 ).These functions indicated whether the effect of the criterion on the score was linear or non-linear, and whether the effect was positive or negative.The sub-priority functions were as follows: where B max is the maximum possible value of B ij (4 plus 4 times the weight of corner cells, see Eq. 2).Functions u 1 and u 2 are sigmoid curves (Fig. 2) indicating that small values of B ij (common border) and A j (stand area) do not contribute much to the score, and the effect of increasing B ij and A j levels-off after certain value.Priority function u 3 implies that decreasing difference in stand characteristics between cell i and stand j increases the likelihood of joining cell i to stand j.
A sigmoid type of relationship between B ij and U contributed to the formation of compact, round-shaped stands and acted against the appearance of stands consisting of long and narrow chains of cells.This type of stands may be formed in the transitional zones between two stands, or in the transitional zone between mineral soil site and peatland.In forestry practice, these zones are seldom demarcated as separate stands.
The usual aim in stand delineation is to avoid stands that are too small for implementing treatments.However, the aim is not to delineate as large stands as possible, since large stands may decrease the efficiency of forest utilization when stands are understood as indivisible management units.For example, when the aim is to harvest the same amount of wood every year, a low number of large stands may hamper the achievement of this management goal.Therefore, a sigmoid relationship between stand area and score (Eq.5) was used.
Both sigmoid-type priority functions had two parameters (b 1 , b 2 , c 1 , c 2 ).Since the specific aim of using sigmoid

29
Fig. 1 In the cellular automaton developed in this study, a cell is joined to one of its adjacent stands.The white cell will be connected to Stand 1 or Stand 2 if basal area is the only stand variable considered (the boldface italic numbers 10, 15, 25, 28, and 29 are stand basal areas in m 2 ha -1 ).The decision depends on the weights of common border, similarity of stand characteristics in the cell and neighboring stand, and area of the stand (Eq. 1) was to create compact round-shaped stands and avoid creating very small stands, the parameters were constrained to the following ranges (see Fig. 2): The CA described above may lead to splitting stands to two or more parts.If the corner cells have a weight higher than zero, different stands may even cross each other without being discontinuous (Fig. 3).This is not always a problem as stands consisting of different parts may be interpreted as strata, which can be used for yield predictions, scenario analyses, etc.However, if the purpose of stand delineation is to create management plans, spatially continuous non-crossing stands are often pursued.
To overcome the problem arising from crossing and discontinuous stands, renumbering was applied to the layout of stand numbers every now and then so that the cells of each continuous and non-crossing area with the same stand number received a new stand number, different from all other continuous areas.In the example of Fig. 3, the yellow stand would be divided into five new stands, two of which would consist of only one cell.As each cell is always joined to one of its neighbors, these one-cell-wide stands would disappear during the next CA iteration.
CA can be run in parallel and sequential modes (Heinonen and Pukkala 2007).A parallel mode can be implemented so that all the changes in cell state (cell state is equal to the stand number of the cell) are implemented at the end of the iteration.In sequential mode, the change is implemented immediately, which means that a change in one cell may affect the ranking of stand numbers of other cells during the same iteration.According to Heinonen and Pukkala (2007), there is no clear difference between the two modes in the performance of the CA.The CA described in this article used the sequential mode.

Optimization
The CA described above has 13 parameters: • Weights of the main criteria of Eq. 1: three parameters All these parameters affect the stand delineation result of the CA.Finding the best setting of parameters by testing different values with trial and error would be a very tedious task due to the high number of combinations.For example, if five different values for each parameter are inspected and all their combinations are tested, the number of required CA runs would be 5 13 , which is equal to 1,220,703,125.
Due to the huge number of alternatives, optimization was used to find the best possible values for the parameters of the CA.The method used was particle swarm optimization (Pukkala 2009;Arias-Rodil et al. 2015), which has been found to work well when the number of simultaneously optimized variables is high (Jin et al. 2018).A description of the particle swarm optimization algorithm used in this study is presented in Appendix 1 in ESM.
The criterion used in parameter optimization was the degree of explained variance, i.e., how much the delineated stands explained of the total variation in site and growing stock characteristics.Since five stand characteristics were used, the objective function that was maximized in parameter optimizations was the average degree of determination of these five variables: where h is the set of 13 parameters affecting the functioning of the CA and R2(h) r is the R2 (degree of explained variance) of variable r.The R2 of a stand variable was calculated from where, SSE is variation not explained by the delineation and SST is the total variation of the characteristic in the area.SST and SSE were calculated as follows: where N is the number of stands, n j is the number of cells in stand j, y ij is the value of the variable in cell i belonging to stand j, y is the overall mean of the variable, and y j is the mean value of the variable among the cells that belong to stand j.

Optimization problems
It can be concluded that when R2 is maximized as the only criterion, the optimal weight of the similarity of cell i and stand j will be high (a 3 in Eq. 1) and the weights of the other criteria will be low.The weight of the corner cell would most probably be high.The result of this CA would be a very scattered layout of irregular stands, many of which are small.To prevent this outcome, the weights of the three criteria were constrained to the following ranges: In addition, the maximum weight of the corner cell was set as 0.5.Parameter optimization conducted with these constraints is referred to as the baseline case.Since it was noted that the result was still a rather scattered layout of stands with many very small stands, another parameter optimization was conducted with a higher minimum weight for stand area and a zero weight for corner cells: This optimization is referred to as the modified case.The third case was a penalized optimization, in which every stand smaller than 0.1 ha (smaller than 4 cells, as one cell of 16 m 9 16 m is 0.0256 ha) contributed to penalty, which was subtracted from the objective function value (Eq.7).The penalty was calculated as follows: Penalty = 3 9 P small , where, P small is the proportion of stands smaller than 0.1 ha.
Preliminary runs indicated that the stands stabilized after about 15-20 iterations.In the optimizations of this study, the number of iterations was 17. Renumbering was done after the 5th, 10th and 15th iteration.

Data
The multi-source national forest inventory raster maps of 2015 were used in this study (ÓNatural Resources Institute Finland, 2017).In this inventory, the site and growing stock characteristics are interpreted for 16 m by 16 m cells, which cover whole Finland (Ma ¨kisara et al. 2016).The data are available at www.paikkatietoikkuna.fi.The data are grouped into 36 rectangular sets according to the Finnish TM35 division of maps by the National Land Survey (www.nls.fi).For this study, one of these area (P5), representing eastern part of central Finland (east of Kuopio, north of Joensuu), was selected first (west-east range 500,000-692,000 m; south-north range 6,954,000-7,050,000 m in the ETRS-TM35FIN coordinate system).Then, three rectangular areas of 150 by 200 cells (30,000 cells) were randomly chosen for the CA runs (Grids A, B and C).The CA parameters were optimized separately for these three grids.

Results
In the baseline optimization, the weight of the corner cell was the largest allowed (0.5) in all three grids (Table 1).The weight of the common border was the lowest allowed (0.4) in two areas.The values of parameters b 1 and b 2 of the sigmoid function for the effect of common border (Eq.4) were always -30 (b 1 ) and 0.8 (b 2 ), leading to a relationship where the proportion of common border started to increase the score of a stand only when the proportion of common border was higher than 0.6 (see Fig. 2).
The obtained parameters led to stand delineation in which the percentage of small stands (smaller than 0.1 ha) was high, 35-43% of all stands (Table 2).The number of stands was high, and the average stand area was small, 0.23-0.36ha.On the other hand, the degree of explained variance was high for all five stand variables, which indicated low within-stand and high between-stand variation is stand characteristics.Figure 4 shows that the CA produced a stand layout where the stand shape was often irregular, in addition to the presence of many small stands.
When the weight of the corner cell was zero and the weight of stand area had to be at least 0.2 (modified problem), the weight of the common border was still near its smallest allowed value, and the sub-priority function for the effect of common border was similar as in the baseline problem (Table 1).However, the number of stands decreased drastically, and the average stand area increased three-to four-fold (Table 2).The proportion of small stands decreased significantly.This increment in stand size was accompanied with decreased degree of explained variance of stand variables.The R2 of site fertility decreased the most, which means that within-stand variation in site fertility increased.
When the objective function was penalized based on the number of small stands (\ 0.1 ha), the number of small stands and the total number of stands decreased further, and the mean stand area was over 1 ha in all tree lattices for which the CA was optimized (Table 2).The location of the sigmoid curve for the effect of common border changed so that smaller proportions of common border started to contribute to the score of the stand (Eq.1).The degree of explained variance decreased, as compared to the modified problem formulation (Table 2).
Parameters obtained from the penalized problem formulation led to a ''clean'' layout of stands (Fig. 4, right panel), in the sense that the shape of the stands was seldom very irregular and the proportion of small stands was low.This appealing delineation was obtained at the cost of increasing within-stand variation, especially in site fertility (Table 2).The effect of parameters obtained from alternative problem formulations on the functioning of the CA is illustrated in Fig. 5. Parameters found in the baseline optimization led to very small average stand size, already after a few iterations (Fig. 5a).Within-stand variation in stand variables decreased rapidly, which is logical since all heterogeneous areas were divided into many small stands.The effect of renumbering after iterations 5, 10 and 15 can be seen clearly.It can be concluded that, before renumbering, the layout resembled the one shown in Fig. 3, in which the renumbering leads to dividing discontinuous and Base refers to the baseline problem formulation, Mod to modified formulation and Penalty to the formulation in which small stands lead to penalty Corner = weight of corner cells in Eq. 2; a 1 , a 2 , a 3 weights of the criteria in Eq. 1; b 1 -c 2 = parameters of Eqs. 4 and 5; w 1 -w 5 = weights of stand variables in Eq. 3 crossing stands into several continuous non-crossing new stands.
When the CA was run with parameters obtained from the modified and penalized solutions, the development of stand delineation during the CA run was very different from the baseline case in a few respects (Fig. 5).The average stand area was constantly much larger, and the effect of renumbering was smaller, suggesting a smaller presence of crossing and discontinuous stands (see Fig. 4, right panel).Because of larger average stand size, withinstand variation was higher than obtained with the baseline parameters.Most within-stand variation occurred in site fertility, which means that larger and more round-shaped stands were created mainly at the cost of increasing withinstand variation in site fertility.
The parameters obtained from the penalized problem formulation with Grid C were used to delineate stands of a larger area (564 9 334 = 188,376 cells, 4822 ha) shown in Figs. 6 and 7. Grid C is located in the upper left corner of this larger area.The results show that the obtained stand borders very well followed the border between mineral soil sites and peatland (Fig. 6a), which is important for forest planning, because only winter cuttings may be prescribed for peatland sites whereas upland forests can be harvested also in summer.On the other hand, there were many stands, most of which represented upland forest, where the site fertility class alternated between two adjacent classes.In most of these cases, the two fertility classes were mesic and sub-xeric, and in some cases they were herb-rich and mesic.In a few cases, the stand consisted of three different site fertility classes.
The stand delineation obtained with the CA followed the spatial changes in mean tree diameter and stand basal area (Fig. 7).Stand borders were drawn in practically all places where the mean tree size or stand density changed within a short distance, which means that the delineation followed the current stand boundaries of the forest.On the other hand, there were many boundaries in places where no change could be observed in tree size or stand density.Only a part of these boundaries can be explained by a change in stand fertility or main site type (Fig. 6).

Discussion
This study described a cellular automaton for delineating homogeneous stands from inventory data where stand variables are interpreted for grid cells.The parameters that guide the CA run were optimized for three different problem formulations, using particle swarm optimization.Since the CA can be understood as an algorithm that creates optimal stand delineations, the method used in this study is in fact an example of two-stage optimization, which may be called as meta optimization, super optimization, or hyper optimization (Jin et al. 2018).The upper stage optimized the parameters of the CA while the lower stage optimized the stand delineation.The objective function of the upper level optimization was Eq. 7 while the lower level optimization maximized Eq. 1.
However, when the CA is used in forestry practice, the CA parameters do not need to be optimized every time.If the purpose of CAs to create homogeneous and continuous  stands for forest management planning, the results of the penalized optimization can be used.Running the CA with these parameters leads to stands fully applicable as calculation units in forest planning and treatment units in the implementation of the plan.If further aggregation of treatments or forest features is required, spatial optimization can be used, as shown in previous studies (O ¨hman and La ¨ma ˚s 2003; Heinonen and Pukkala 2004).Several methods are available for spatial optimization: mathematical programming, centralized heuristics such as simulated annealing, tabu search or genetic algorithms (Borges et al. 2002;Bettinger et al. 2002), or decentralized heuristics such as CA and the spatial version of the reduced cost method of Hoganson and Rose (1984), suggested by Pukkala et al. (2008).Spatial optimization can organize a pixel forest into large enough stands also during the forest planning calculations, i.e., without any prior-planning segmentation.However, a large forest consisting of small pixels leads to so large spatial optimization problems that they cannot be solved in reasonable time (Heinonen et al. 2018).For these problems, the most capable category of spatial optimization methods are the same decentralized heuristics which can be used in pre-planning segmentation (Strange et al. 2001;Mathey et al. 2007;Heinonen and Pukkala 2007).Therefore, CA and other decentralized heuristics can be used in two different stages of the planning process: first, to In addition to creating harvest blocks, spatial optimization may also be used to form large continuous habitat patches (Heinonen et al. 2007), minimize the risk of wind damage (Heinonen et al. 2009;Zeng et al. 2010) or create fuel breaks (Gonza ´lez-Olabarria and Pukkala 2011).
When it was required that the CA should delineate compact, continuous and large enough stands, the result was a delineation where the stands were homogeneous with respect to stand variables other than fertility class.Figure 6 shows that it was impossible to create a delineation in which there were no small stands and the stands were homogeneous with respect to site fertility.Therefore, a compromise was necessary.Areas where site fertility changes within a short distance may represent sites where fertility is near the limit of two adjacent fertility classes.Therefore, the main reason for the seemingly heterogeneous stands might be that a continuous feature (site fertility) was classified into discontinuous classes.
The advantage of CA, as compared to several other segmentation methods, is that CA is transparent and flexible one-step method that is easy to program and control.The same CA can be used with several types of lattice data (e.g.Lidar, multisource, or satellite inventory).In addition, the shape of the interpretation unit does not need to be square.Hexagonal shape, for example, has the advantage that it removes the effect of corner cells.Therefore, if the shape of the interpretation unit can be decided in forest inventory, hexagon is preferable to square (Heinonen et al. 2007).
The shape of the interpretation unit can also be irregular (Pukkala et al. 2014).This requires that the length of the common border between adjacent units is calculated, which can be easily done with GIS tools.The stand characteristics used to measure the similarity of cells can be chosen freely.In this study, they were the same variables which are imported to Finnish forest planning systems.However, basal area, mean tree height and mean diameter are usually assessed separately for different tree species.Using species-specific values was not regarded necessary in stand delineation.However, after the delineation, the growing stock variables can be calculated separately for different tree species, and these variables can be imported to the planning systems.
The much used region-merging method (Mustonen et al. 2008;Wulder et al. 2008;Mora et al. 2010;Olofsson and Holmgren 2014) is ''one-directional'' in the sense that the area of the segments always increases when adjacent units are merged.New stand boundaries are never drawn as larger segments are formed by removing boundaries between existing adjacent units.In CA, the segments may enlarge or shrink, and the boundaries may change location.The starting segments may be large or small, and the final segments may be smaller or larger than the initial segments.For these reasons, the CA can be regarded as a more dynamic and flexible method than the region-merging algorithm.An interesting topic for future research would be to compare the performance of the region-merging method and the CA algorithm presented in this study.
One feature of the stand layouts created by the CA is that large homogeneous forest areas (in terms of site and growing stock characteristics) were divided into several stands, which is partly due to the shape of the sub-priority function for stand area (Fig. 2).This is not a disadvantage in forest planning since a higher number of stands increases the possibilities to optimize management.Planning may combine similar adjacent stands into larger units if further aggregation is regarded useful to reduce harvesting costs (Lu and Eriksson 2000), prevent forest fragmentation, reduce the risk of wind damage (Heinonen et al. 2009;Zubizarreta-Gerendiain et al. 2018) or create large enough continuous habitat patches for forest-dwelling species.
As a summary, the CA described and tested in this article is a useful new tool, which could be immediately used in forest planning to reduce the number of calculation units and simplify forest planning calculations.The parameters of the CA can be adjusted based on the purpose of the aggregation.If the purpose is to create stands for forest planning, parameters obtained from the penalized problem formulation of this study can be used.The need for aggregation tools similar to the CA described here is increasing because forest variables are increasingly interpreted for units that are too small to serve as treatment units in the implementation of forest management plans.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 2 Fig. 3
Fig.2Effect of the parameters of Eqs. 4 and 5 on the contribution of common border between cell and stand (a) and stand area (b) to the score of the stand (Eq.1).The score describes the likelihood of joining the cell to the stand

Fig. 4
Fig. 4 Evolution of stand delineation in a grid of 200 9 200 cells during a CA run after the first (a, d), 5th (b, e) and 15th (c, f) iteration when the parameters of the CA are obtained from the baseline (a-

Fig. 5
Fig. 5 Development of stand area (left), within-stand standard deviation (middle) and degree of explained variance (right) during a CA run in Grid C with parameters obtained from baseline (top), modified (middle) or penalized (bottom) parameter optimization

Fig. 6
Fig. 6 Stand delineation when the CA was used with parameters obtained from the penalized problem formulation for Grid C. Grid C (200 9 150 pixels) is located in the northwest corner of the area (consisting of 564 9 334 pixels).a Main site type (red = mineral soil, green = spruce mire, yellow = pine bog, light blue = open bog).b Fertility class (red = herb-rich, green = mesic, blue = sub-xeric, yellow = xeric)

Table 1
Optimal parameter values of the cellular automaton in three areas (Grids A, B and C) for three problems(base, modified, penalty)

Table 2
Results calculated for stand delineations obtained with optimized CA parameters in different problem formulations for three grids