Generating a Raster Map Presentation of a Forest Resource by Solving a Transportation Problem

Necessary tools for raster map generation, for the approach based on the calibration estimator, were developed and implemented. The allocation of the area weight of each pixel to sample plots was formulated as a transportation problem, using a spectral distance measure as a transportation cost, and solved using the transportation simplex algorithm. Pixel level accuracy was calculated for the methods based on the calibration estimator so that the results could be compared with the results of the nearest neighbour estimation, the reference sample plot method (RSP) at pixel level. Local averaging in a 3 × 3 window was performed for each generated raster map as a postprocessing phase to smooth the map. Test plot results were calculated both for the unfiltered raster map and the filtered raster map. RSP produced the smallest RMSE in the pooled test data. Local averaging with a 3 × 3 filter decreased the pixel level error – and the bias – and the differences between the methods are smaller. Without local averaging, the pixel level errors of the methods based on solving the transportation problem were high. Raster map generation using the methods of this study forms an optional part – followed possibly by the classification of the pixel level results – of the whole computation task, when the area weight computation is based on the calibration estimation. For larger areas than in the present study, such as municipalities, the efficiency of the method based on the transportation model must be improved before it is a usable tool, in practice, for raster map generation. For nearest neighbour methods, the area size is not such a problem, because the inventory area is processed pixel by pixel.


Introduction
Forest inventory applications using satellite remote sensing data are needed to achieve necessary input data for forest planning calculations for large areas.Tokola (2000) reviewed and categorised the main methods for forest inventory applications as stratification (or classification) and regression analysis.Further, the applications of nearest neighbour regression have been of great interest -the reference sample plot method (RSP) (Kilkki and Päivinen 1987, Kilkki 1988, Muinonen and Tokola 1990, Tokola et al. 1996, Tokola 2000) and the Finnish multisource forest inventory method (k-NN) (Tomppo 1996, Katila and Tomppo 2001, Katila and Tomppo 2002, Tomppo et al. 2002).McRoberts et al. (2002), who used nearest neighbour approach as a part of stratified estimation, gave a comprehensive analysis, with precautions, for applications based on the nearest neighbour technique.
The concept of area weight, or sample plot weight, refers to interpreting weight as the amount of similar forest the plot represents in the inventory area (i.e., the region for which the inventory results are calculated).A method, denoted with SAE in this study, for computing the sample plot weights for small area estimation of forest parameters has been presented by Lappi (2001).Small areas, in this context, refer to areas of a size comparable to a municipality, which is usually 200-1000 km 2 in land area, or smaller ones (see also Lappi 2001).This method is based on the calibration estimation, where satellite data form a basis for the weight calibration acting as calibration variables, because the satellite data are known for every pixel.The initial sample plot weights (later also referred to as prior weights) are calibrated to satisfy the known area level totals or mean values, when computed from the satellite data and weights of the sample plot pixels (i.e., pixels each of which cover a field sample plot).This is performed by formulating calibration equations including the weights of the sample plot pixels and the known totals, or mean values of the calibration variables (see Lappi 2001).The known totals of the calibration variables are computed by analysing all the pixels in the inventory area.
When the sum of the calibrated weights equals the number of pixels, the estimated mean value of a forest parameter can be calculated as a weighted average of the forest variables obtained from the sample plots.Besides the original satellite data, channel transformations, that is, volume estimates for each pixel based on the sample plot pixels, can be incorporated into the calibration as auxiliary calibration variables.A variance estimator, based on a spatial model, is also presented (Lappi 2001).
The SAE method, based on calibration estimation, includes the computation of the area level totals of the calibration variables (satellite data) for the inventory area, designing a set of initial weights for the sample plots and calculating the calibrated weights.In the methods, based on nearest neighbour approach, the weights for the sample plots are the cumulative weights from all the pixel level estimates of the inventory area.
According to Lappi (2001), there is no guarantee in RSP that the chosen nearest neighbours add up to unbiased or statistically optimal estimates for the whole region.In the SAE approach, the total area that the sample plot represents (the weight of the plot) is estimated, first, starting from a set of initial weights.The initial weights cannot be set according to the inclusion probabilities if the plots outside the inventory area are to be used, hence the set of initial weights becomes a user-defined parameter (Lappi 2001).As a conclusion, the RSP-based weighting of the sample plots is one alternative for having the initial weights for SAE.
Both of the methods, RSP and calibration estimation, offer means to utilise the plot weights in forest planning systems.For forest planning, the input data should be accurate and localised.The methods, based on nearest neighbour analysis, have the capability to produce results for a forest variable of interest in a raster map form.This cannot be carried out directly in the case of SAE which offers only a set of sample plot weights as a computational result.Thus, for generating a raster map of the estimates, the origin of the weight, allocated to each sample plot pixel, should be known at a single pixel level (RSP), or, if only the estimated weights are known (SAE), the allocated weights should be distributed to the pixels of the inventory area in a separate process.Lappi (2001) suggests a formulation of a linear programming transportation problem to find such a weight allocation, where an overall spectral distance between plot pixels and destination pixels is as small as possible.In the case of using the calibration estimation-based approach, this would combine the small-area calibration estimator and the basic principle of the nearest neighbour method in map production.
In the RSP method, the estimation of forest variables is carried out separately for each pixel in the inventory area, and the raster map can be generated as a by-product, if needed.In the SAE-based methods, the map must be generated in a separate analysis, because the calibration is based on area level means or totals.In the raster map generation, each pixel, in the inventory area, receives weight in such a way that an overall distance is minimised in the whole inventory area.Area level accuracy is not changed in a map generation as the sample plot weighting does not change in distributing the weight over the area.However, the accuracy at pixel level can then be estimated if there is located information on the forest variables available.However, due to the nature of the SAE method, the pixel level results are, then, dependent on the inventory area in which it belongs, and this has to be taken into account when comparing the pixel level results of the method, after map generation, with the results of the RSP method.
The objective of this study is to present and test the feasibility of the approach of solving a transportation problem for the raster map generation, as a part of a forest inventory system based on calibration estimation.For comparison, the results are calculated using a method (RSP) based purely on the nearest neighbour approach.The accuracy of results, at pixel level, is calculated to evaluate the raster map.The calibration estimator-based method, together with initial weights from RSP (i.e., a combined method) is also applied.The framework of the study is to present a map production tool for the SAE-based technique (Lappi 2001), and to present a methodology of using the RSP in the prior weight production.

Material
The satellite data consist of a Landsat 7 ETM image (WRS 187/17; a floated scene) from eastern Finland, acquired on August 2, 1999.The image was rectified to the Finnish National Coordinate system using first order regression models and 43 control points which were located from the base maps of the area.A 25-m pixel size was produced using nearest neighbour resampling in which the root mean square error in the rectification was 0.4086 pixels.
Digital map data were used both in the rectification of the image and also in the selection of the field plots to the analysis.A numerical map data mask, processed in the NFI remote sensing laboratory, included the following categories: forestry land, agricultural areas, water areas, and other non-forestry land (e.g., roads, buildings).For details of the map data see, for instance, Katila et al. (2000).The map data were used to mask out categories other than forestry land.The selection of the sample plots was based on this map data to mimic the true inventory situation.
Field data consist of the sample plots measured in 1999 for the 9th National Forest Inventory in the forest district of Etelä-Savo, which fall in the forestry land category based on the map data and are located in the satellite image window.The field sample is a systematic cluster sample, where a cluster consists of 14 temporary relascope plots, or 10 permanent plots (every 4th cluster).The distance between the clusters is 6 km in both North-South and East-West directions (Valtakunnan metsien… 1999).The resulting data set comprised 2863 sample plots.The main forest variable of interest is the mean volume of the growing stock (m 3 /ha) (Table 1).
For the pixel level checking and evaluation of the techniques, two 15 km × 10 km rectangular inventory areas (test areas TPA and TPB) were generated.The sample plot data were used inside a 12-km inclusion zone around TPA and TPB, respectively.The pixel level accuracy of the estimation was computed for each raster map generation method using randomly selected test data sets of the sample plot pixels.The proportion of the sample plot pixels in TPA and TPB selected to the test data was 75%.The number of outside plots (plots in the inclusion zone) of the areas TPA and TPB was 271 and 287, respectively.Thus, the number of plots available in the estimation for TPA is 287, and there are 298 plots available for TPB.The pixels and sample plots, for the analysis, were selected from the category 'forestry land', based on the digital map data.In TPA, there were 191 895 pixels, and TPB included 193 498 pixels.
The relascope sample trees (tally trees), located on the land use classes, forest land or other wooded land, are measured in the sample plot.In this study, the sample plots on the forestry land map stratum with no trees measured, were given a 0.0 m 3 /ha value of mean volume in the analysis.These plots most usually fall into other land use classes in a field inventory, and, thus, no trees are sampled (Table 2).

Nearest Neighbour Estimator
The method applied in the estimation of the forest variable, at pixel level, is the reference sample plot method (RSP) (Kilkki and Päivinen 1987, Kilkki 1988, Muinonen and Tokola 1990, Tokola et al. 1996, Tokola 2000).The spectral distance (d ij ) between pixels i and j is calculated using a distance function where n c = number of bands, i,j = pixels, b ih = the spectral value of a pixel i at band h; p h = empirical constant for band h.
After computing the distance function values from the target pixel i (i.e., the pixel for which the forest variables are estimated) to each sample plot pixel j, m spectrally nearest sample plot pixels are sought and the weight (w ij ) for each neighbour pixel j = 1,…,m of the pixel i is calculated as After computing the weights w ij , the estimator of a forest variable y, for a pixel i, can be formulated as in Eq. 3.For a pixel i, it is a weighted average of the field measured forest variable values for the spectrally m nearest sample plot pixels: where y ij denotes the forest variable value of the jth nearest neighbour of the target pixel i.
The empirical coefficients p h , in the distance function, describe the relative importance of the variables.They were selected by testing all the combinations of the values p h ∈{0.1, 5.0, 10.0} for the variables in the distance function.The pixel under estimation was left out from the analysis when selecting the m spectrally nearest neighbours for the pixel.The analysis was carried out using the so-called leave one out method, that is, a forest variable estimate for each sample plot pixel i in turn was computed using the other sample plot pixels j as reference.

Calibration Estimation
Calibration estimation (Eqs.4 and 5), (see Lappi 2001), can be applied in the estimation of the sample plot weights (w k , in Eq. 5), In Eqs. 4 and 5, y k denotes a forest variable for a sample unit k, drawn from a finite population, where each unit is associated with a vector x k of spectral variables or auxiliary variables; n is the number of sample plots in the analysis.In the case of this study, the sample units are pixels associated with y variables -forest variables.Eq. 4 describes an estimator (denoted using subscript d) for forest variable y, based on the prior weights d k .Eq. 5 is the calibration estimator for the forest variable y (see Lappi 2001).The calibrated weights w k are obtained by calibrating the prior weights (see Eq. 4), respecting the calibration equation: where t x is the known population total vector of the known x variables (i.e., satellite image pixel spectral values, or other known auxiliary variables for each pixel).In the SAE application, sample plot pixels, outside the inventory area, are included in the estimation (Lappi 2001).This means that the prior weights -originally based on the inclusion probabilities -had to be defined in another way in this study.They were defined using a geographical inclusion zone, where the prior weight is constant, inside the inventory area, and decreases linearly, in the inclusion zone, as the distance increases.The other method applied, in the present study, is to use the RSP sample plot weights as prior weights in the calibration.The computation of the area weight estimation was carried out as in the study of Lappi (2001), (see also Deville and Särndal 1992), incorporating the following distance function: Parameters L and U define the range of the weights, Ld k < w k < Ud k .To have nonnegative weights, L and U were given values 0 and 6.5, respectively (see Lappi 2001).The nonlinear optimisation problem, needed in the search for the weights, was solved by using the routines for the Newton method presented by Press et al. (1992).
In the estimations for this study, the calibration was based on the mean value vector (instead of total values) and t x was computed using satellite pixel data in the test area under examination.Pixel values were scaled, by dividing them by their mean value in the data, to obtain a more stable convergence (see Lappi 2001).

Transportation Problem
Pixel level accuracy is calculated for the methods based on the calibration estimator to allow these results to be compared with the results of the nearest neighbour estimator (RSP) at pixel level, also taking into consideration that in SAE, the pixel level results depend on the properties of the whole inventory area it belongs to.The prerequisite for the pixel level accuracy estimation is that the weight allocated to each sample plot has to be assigned to the pixels in the inventory area.This means that a raster map generation technique has to be implemented.
The weight allocation was formulated as a transportation problem which, in general, concerns transporting goods or services from supply centres to multiple demand centres in an optimal manner (see Dykstra 1984).The problem can be viewed as a network problem, where supply and demand centres are represented by nodes, and the transportation system joining them is represented by a set of directed links.Each link has a value of any relevant unit of measure, representing, for example, transportation cost or travel time.(Dykstra 1984).When the problem is presented in a tableau form, the rows illustrate the sources (supply centres), and the columns are the destinations (demand centres).
A usual way to present the problem is to formulate it as a linear programming problem (see also Dykstra 1984): where M is the number of supply centres, N is the number of demand centres, C ij is the relevant link distance, and x ij is the number of units to be transported from the supply centre i to the demand centre j.Total cost of transporting (z) is minimised.It follows that in feasible solutions for this kind of a problem: Thus, in the transportation model, all constraints are equalities, and the sum of supplies and demands must be equal (Dykstra 1984).For solving this kind of problems, an algorithm that takes advantage of the special structure of the transportation model has been developed.The algorithm, used in this study, has been described in detail with an example case by Dykstra (1984), (see also Taha 1992).
To generate the initial solution for solving the transportation problem, feasible methods, such as the least cost method (LCM) and Vogel's approximation method (VAM), exist (see Taha 1992).In the least cost method, as much weight as possible is iteratively allocated to a variable (tableau cell) with the smallest unit cost in the entire tableau.After each allocation, supply and demand are adjusted accordingly and the satisfied row, or column, is crossed out.VAM is a heuristic method and usually provides a better starting solution than the least cost method and generally yields an optimum, or close to an optimum, starting solution (Taha 1992).For this study, VAM serves an interesting method for map generation in the case that it could produce a suitable map presentation even without solving the transportation problem.
The VAM method is an iterative method, based on computing the penalties -of not choosing the cheapest alternative -for each row and column in a transportation tableau (see Taha 1992).The initial solution is generated by choosing iteratively the row, or column, with the largest penalty and, then, allocating as much weight as possible to the variable (cell) having the least cost.

Applying the Methods in Raster Map Generation
Four different combinations of methods, to produce the raster map of mean volume, were tested and they are summarised in Table 3.For each method, the pixel level accuracy is computed for test sample plot pixels where a true value is known.Test sample plot data are used only in the evaluation part of the analysis.Thus, the test sample plot pixels have not been used as demand centres in the transportation problem either.
Local averaging in a 3 × 3 window was performed for each generated raster map as a postprocessing phase to smooth the map.Test plot results were calculated both for the unfiltered raster map and the filtered raster map.Only the pixels belonging to the forestry land stratum, in the digital map data, were used in the local averaging.
The first area weight estimation was computed using the RSP method.An inventory area A that contains n A pixels in the satellite image is assumed.For each of these n A pixels, i = 1,…,n A , a pixel level estimate is computed using Eq.3; the weight W kA which a sample plot k obtains in the area level estimate for area A is the sum of the weights it receives in the pixel level estimates in the inventory area The weights W kA , k = 1,…,n, for the area level estimate for area A, sum up to n A , and are the area weights.If inventory area A includes only a single pixel, the area level estimate is equal to a single pixel level RSP estimate.Thus, the estimates for an area A, containing several pixels, would be computed as weighted averages of the plot pixel field measurements, similarly as in Eq. 5. Instead of computing RSP area level estimates, the weights, from this nearest neighbour method, were stored and further used as initial weights for one of the SAE-based methods (Table 3).
In both of the methods, C1 and C2, for the area weight estimates by plots, the calibration estimator was applied, (see Lappi 2001, Kangas andMaltamo 2000).The methods differ only in the mode the initial weights are produced (Table 3).In the C1 method, the initial weights d k in Eq. 4 were computed using an approach of a 12-km geographical inclusion zone around the inventory area.Inside the inventory area, the prior weight was set to a constant value, and the weight decreases linearly with respect to the geographical distance to the inventory area.In the C2 method, the sums of weights by plots from the RSP estimation by pixels, were used as prior weights.Thus, the third (C2) area weight estimation was computed using the weights W kA , k = 1,…,n from the RSP as prior weights d k in Eq. 4. A sample plot is used in the estimation if it has been given a positive prior weight, in the pixel level estima- tion, for the pixels in the inventory area and, at area level, the prior weight d k for sample plot k (in Eq. 4) is, then, the sum of these weights.
In the methods C1 and C2, quite a large geographical search radius -60 km from each target pixel, where the spectrally nearest neighbours are sought -was used so as not to limit the usage of the sample plot data, in the surrounding zone of the area, in creating the raster map.The raster map generation was performed using the transportation problem model.For the analysis of this study (see Eqs. 8 and 9), the supply centres are the M pixels of the area and each pixel i has a supply S i equal to 1.The demand centres are the N sample plot pixels in the data, each plot pixel j having a demand D j equal to the area weight (noted with w k for plot k in Eq. 5) it has received in the area weight estimation.A transportation cost C ij , between source i and destination j, is a distance function value (Eq. 1) which was calculated from the satellite data, the number of neighbours being 1.As a result, total spectral distance is minimised when transporting the area weight of the pixels to the plots up to a given demand.In the RSP estimation, the raster map is generated as a by-product and no separate computation is needed for this part, as the inventory area is processed pixel by pixel, the number of neighbours being 5.
In the C3 method, VAM was applied in the raster map generation.The spectral distance function value, with 1 nearest neighbour, was used as a transportation cost (Table 3).

Pixel Level Accuracy
Pixel level accuracy was analysed using a cross validation technique.When defining the coefficients for the distance function, the mean volume for each sample plot pixel was estimated using the other plot pixels as a reference.The pixel level accuracy of the estimation was evaluated, first, with where n = number of sample plot pixels in the estimation; and, second, (see Tokola et al. 1996, Tokola 2000), with In R 2 , the original variation, in the test data under analysis, is taken into account, not only the RMSE of the estimation.
The possible bias of the estimation, at pixel level, was evaluated using bias Relative RMSE and bias were computed by dividing the RMSE and bias by the mean value of the y variable in the data set.The statistical significance of bias, at pixel level, was tested by computing where s = standard deviation of the differences ( ˆ) y y i i − .
In the evaluation phase, using the test sample plot pixels, the measured value and the estimate from the method under analysis, were drawn for the pixels, and the Eqs.11-13 were applied.

Results
The coefficients of the distance function, resulting in a minimum RMSE and an insignificant bias (t bias < 1.96) at pixel level, were sought by performing a grid search.The grid search was made separately for the cases by using 1 nearest neighbour -to be applied in the raster map generation using a transportation model -and 5 nearest neighbours -to be applied in the RSP estimation generating the raster map pixel by pixel.
It can be noted that using one nearest neighbour (m = 1 in Eq. 3), in the estimation, produces an RMSE value which is greater than the original standard deviation (Table 4).This means that the overall mean could give a better estimate than using the one nearest neighbour does.Thus, if the number of nearest neighbours can be chosen for the estimation, as in the RSP, a value greater than 1 would be selected.In the raster map generation using RSP, 5 nearest neighbours were used to produce a more proper basis for the comparisons of the usability of the methods.
Pooled pixel level results, in the test areas TPA and TPB, from the raster map generation techniques which were applied, show that the problem of using 1 nearest neighbour in the allocation of weights causes a poor pixel level quality (Table 5).RSP estimation produced the smallest RMSE in the pooled test data (Table 5).Local averaging with a 3 × 3 filter decreased the pixel level error -and the bias -in all cases, and the differences between methods are smaller (Table 5).Bias was not found to be significant (t bias < 1.96).
In the test area TPA, the test plot results show a larger relative pixel level error than in TPB, for the tested methods.Local averaging caused a decrease of RMSE in both areas, and all tested methods.For RSP, local averaging decreased bias in both areas, but for the other techniques, the effect of the local averaging on bias is not so clear in TPA and TPB (Fig. 1).After local averaging, C2 has the smallest bias in the pooled data (Table 5) -as a result of opposite bias values in TPA and TPB test data, but C1 is more stable (Fig. 1).The computationally simple method, C3, has largest bias, after local averaging, in the pooled data, resulting from the bias in TPB test area (Table 5, Fig. 1).
As an example of the raster map presentation of the estimated mean volume, generated by using the C1 and RSP methods, Fig. 2 illustrates the classified output of the applied methods in a selected 1.5 km × 1.5 km sample in the TPA area.The effect of local averaging can be seen, as expected, as a smoother output raster map appearance.

Discussion
In this study, necessary tools, concerning raster map generation, for the approach based on the  The allocation of the area weight of each pixel to sample plots was formulated as a transportation problem, using a spectral distance measure as a transportation cost, and solved using the transportation simplex algorithm (Dykstra 1984, Taha 1992).
The raster map generation is necessary for presenting results in a map form in the case of using the method presented by Lappi (2001) -that is, a method based on the calibration estimator -which does not directly produce a map.Another simple way to deal with the sub-area level reporting might be to compute the area weights for each sub-area separately, without a raster map generation, but it may suffer from the possibility of having quite a small number of sample plots available, even if including the plots from a larger geographical inclusion zone.It is clear, however, that the small number of neighbours is problem-atic in the sense of having accurate pixel level results (Tokola et al. 1996, Tomppo et al. 1998, Katila and Tomppo 2001).
In the raster map generation, the transportation cost was expressed by means of a spectral distance measure.To find a most suitable demand centre (sample plot) for each pixel, the spectral distance was minimised.This leads to an approach similar to a nearest neighbour estimation with the number of neighbours being 1 (C1, C2, C3) or 5 (RSP) and with no constraints for the weights allocated to each plot.Here, the weight to be allocated to a sample plot is constrained to have a maximum value which equals the area weight of the sample plot.This plot level area weighting is obtained from the calibration estimation procedure, where the calibration equations have to be fulfilled.
The selection of the channel weights, in the distance functions, was made by using a grid search method to have a minimum RMSE and an insignif- icant bias for the mean volume estimation.Other methods, for channel weight search, would be the use of multiobjective optimisation (Haara 2002), or canonical analysis (Moeur and Stage 1995).In these methods, the channel weights could be calculated by taking several forest variables into account.These forest variables could include both age and volume, for example.Another approach, for taking into account several forest variables at area level -and simultaneously also at pixel level -could be incorporating pixel level forest vari-able estimates into the estimation as calibration variables, as suggested by Lappi (2001).
In the selection of the channel weights, the whole 1999 data were used.However, the models calibrated in the large area, produced bias in the smaller test areas.A larger number of sample plots might be needed, in the smaller test areas, than used in this study; as in the study of Lappi (2001), the zone size could be computed to have a fixed number of sample plots.In this study, also using a part of the plots in the test areas as test plots decreased the number of the plots available in the estimation.Channel weights, for the distance function used in the transportation problem, produced large RMSE also in the whole 1999-year data set.The number of nearest neighbours, for this distance function, had to be this small, because the weight of a pixel is transported to one sample in most of the cases.This is the price, in raster map generation, of using the SAE approach -that is, calibration estimation, as used in this study -in computing the area weights.The precaution, reported by McRoberts et al. (2002), that a small number of nearest neighbours may result in an RMSE value greater than the standard deviation of the observations, came true in this study.It is known that with a smaller number of nearest neighbours RMSE increases (Tokola et al. 1996, McRoberts et al. 2002).Increasing the number of neighbours, to a greater value than 5, was not carried out, in order to preserve the estimation variance.When increasing the size of the neighbourhood, the estimation variance decreases and the bias increases (Altman 1992).Hence, the pixel level error of the methods C1-C3, in the study, without averaging was considerable, as could be expected.Lappi (2001) has suggested local averaging to improve the raster map generation by smoothing the map.This was performed in the present study to see whether this kind of postprocessing would be necessary to improve the pixel level accuracy of the results in the generated raster map.The filtering appeared to be necessary, especially for the methods based on the calibration estimation (C1 and C2), giving a better pixel level accuracy and a smaller bias to the results calculated over the test plots.
The results of RSP, in the test areas TPA and TPB together, show a bias of -1.7 m 3 /ha (-1.45%) also after local averaging (Table 5).The results from RSP, obtained for test plots in TPA and TPB, are biased (ca.5% in TPA and -7% in TPB), although the level of bias in the grid search of the channel weights is not so high (Table 4).The RMSE level of RSP is a result of similar nature as reported when using the Euclidean distance function and RSP (Muinonen andTokola 1990, Tokola et al. 1996).The estimation may suffer from too few observations in the spectral space, and the number of sample plots, in the analyses for test areas, was too small, in the study, to have a satisfactory pixel level accuracy.As a result, the calibration of the RSP technique was not good in the test areas, that is, at local level.Further analyses concerning, for instance, the zone size, weights of the covariates, number of neighbours, are not in the scope of this study.A larger inclusion zone could have been a usable choice, especially for RSP estimation.
The transportation simplex algorithm appeared to work as planned in the problem solving.For the purpose, a computer programme was developed due to the size of the problem being large.The problem solving took approximately 5 days with a computer equipped with a 2.4 GHz processor and 992 MB of RAM.Thus, it is clear that for larger areas, there is a need to compress the size of the problem and improve the problem solving, for example, by allowing only a limited set of possible target plots for each pixel.The large number of rows (pixels) is a problem, which also needs to be solved so that the approach could be used in larger scale and in practice.For these reasons, the nonparametric approach has a straightforward capability to produce the raster map output, because the inventory area is processed pixel by pixel.Furthermore, it has to be taken into account that the pixel level estimates of the methods C1-C3 are dependent on the properties of the inventory area.If a different area is selected over a pixel, this pixel can obtain a different estimate, as the area level transportation problem has changed.
In C1 and C2, calibration estimation, after having necessary initial weights, and the transportation model were applied.Geographical initial weighting was the basis of C1, and the initial weighting from the nearest neighbour analysis was the basis of the C2 method.The differences between C1 and C2 are quite small, but local averaging was necessary to lower the pixel level RMSE.(Table 5).The bias level, in the test areas TPA and TPB, is lower in C1 than in C2.Actually, the level of bias calculated over test plots is lowest in C1 than in any other methods tested.For C2, the pooled test results show a small bias, but C2 did not perform well in TPA and TPB test areas, and the average filter increased the bias of C2 (Fig. 1).
The usability of VAM in the map generation can also be discussed by comparing the pixel level results from the VAM solution and transportation model.Satisfactory results from using VAM, especially for larger applications, would offer a way to avoid the computation in solving the transportation For this reason, a raster map was generated and the estimates for the test sample plot pixels were calculated by using the VAM approach, too.It was computationally easier, as the time demand, in the test areas, was only approximately 2 hours for each test area.However, the C3 method had largest bias calculated over test plots.The results, however, indicate that it is still an interesting alternative, when the inventory area increases.It has to be noted that the results, concerning pixel level accuracy and bias, can be the outcome of the non-optimal calibration of the nonparametric estimation.More detailed search of the distance function parameters and variables, might improve the results of each method.
To perform calculations over a larger area, as a municipality, by using the approach of solving the transportation problem, the computational efficiency of the method has to be improved before it is an available tool for raster map generation.This kind of work was not possible in the present research project.For nearest neighbour methods, the area size is not such a problem, as the inventory area is processed pixel by pixel.As a result, the raster map generation approach applied in the study, can be seen as an optional part -followed possibly by the classification of the pixel level results -of the whole computation task, especially when it is based on the calibration estimation, also taking into consideration the properties of the method, and capabilities of satellite data in the estimation of forest variables.

Table 1 .
Mean volume of the growing stock (m 3 /ha) in the sample plot data: a) 1999 data (test data not included); b) test area TPA; and c) test area TPB.(n: number of plots; TPAtest, TPBtest: test plots only; TPA&zone, TPB&zone: plots in the test area and in the 12-km zone).

Table 2 .
Land use class distribution (number of plots) among the field plots belonging to the map stratum forestry land over 1999 data and the test areas.(Land use class Forestry includes forest land, other wooded land, waste land and other forestry land).a) 1999 data (without test data; n = 2760); b) test area TPA; and c) test area TPB.(TPAtest, TPBtest: test plots only; TPA&zone, TPB&zone: plots in the test area and in the 12-km zone).

Table 3 .
Techniques applied in the raster map generation.Local averaging was tested in postprocessing the generated raster map after each method.

Table 4 .
Coefficients for the distance function, produced by the grid search, and the pixel level accuracy of the nearest neighbour estimates for mean volume (m 3 /ha).Data: sample plots of 1999, test data excluded; n = 2760; Mean = 130.9 m 3 /ha; S.D. = 108.1 m 3 /ha; geographical reference zone = 60 km; m is the number of nearest neighbours; TM1…TM5 and TM7 are the satellite image channels 1…5 and 7.

Table 5 .
Pixel level estimation results in the pooled TPA and TPB test data for mean volume; a) without local averaging, and, b) with 3 × 3 local average filter.Bias is computed in the test data.(Mean 116.89 m 3 /ha; S.D. 107.60 m 3 /ha; n = 103).