The effect of the distribution of measurement points around the node on the accuracy of interpolation of the digital terrain model

Contemporary land information systems allow the generation of models of numerical surfaces by applying a number of interpolation algorithms. The appropriate selection of measurement points included in the calculations is one of the key factors influencing the quality and accuracy of interpolation. Properly selected points contribute to improving the accuracy of the generated surface models and to shortening the computation process. This paper analyzes the effect of the location and density of measurement points on the accuracy of interpolation surfaces in view of the morphological differentiation of the generated models.


Introduction
Spatial information systems (SIT, GIS) are frequently used for processing, analyzing and maintaining large volumes of information originating from a variety of sources (Harmon and Anderson 2003;Heywood et al. 2002;Longley et al. 2005). Generally, data originating from the digital terrain model (DTM) provide a basis for spatial organization in such systems (Carlisle 2002;Li et al. 2005). Source data forming the DTM are most frequently obtained from direct measurements or indirectly from the digitalization of analog materials (Maune 2001;Wechsler 2003;Weng 2002). Contemporary measurement systems, such as aerial laser scanning or multiple beam sonic sounders, can obtain large volumes of information within a relatively short time (Aerts et al. 2003;Aruga et al. 2005;Cobby et al. 2001;Wack and Wimmer 2002). With a large volume of measurement data available, it is possible to generate a very accurate terrain model. The direct use of such sets, however, is difficult because of the volume and non-standardized nature of information recording (Heo 2003;Takefusa 2001;Wechsler 2003). Due to the high pace of development and the particular role played by information, the quality and reliability of data and dynamic data processing are equally important considerations. It is also desirable to conduct comprehensive comparative analyses for the same elements during different measurement periods. As a consequence, attempts are underway to organize the spatial structure and limit the number of points generating the digital terrain model while maintaining the required accuracy. The presentation of information describing the surface in the form of a regular GRID-type node structure (a regular grid of nodes-the points forming the DTM are located in the corners of squares with a specified side length) is one of the methods applied to organize the above information (Hutchinson and Gallant 2000;Raaflaub and Collins 2006;Zhou and Liu 2004). It enables a reduction in the amount of data, limits redundancy and significantly accelerates information processing. During the processing of data originating from different sources, particular attention should be paid to its distribution in the measurement spectrum, which is determined by the acquisition method. The quality and accuracy of data depends on the acquisition method. Irrespective of data type, in each case of GRID interpolation, we should aim at achieving the maximum accuracy of the DTM and reducing the number of points forming the model. Appropriate selection of interpolation algorithm parameters is an important consideration in the process of generating the grid squares (Maune 2001;Schabenberger and Gotway 2005). These algorithms support the interpolation of the values at individual nodes based on measurement points (Scherzinger et al. 2001;Shi et al. 2004;Wilson and Gallant 2000). Proper selection of measurement points used in the calculations is of key importance for precise determination of the values at the nodes.
The subject matter covered by the article presents the methodology and the process of analysis conducted on selected surface models and selected interpolation algorithms. Altogether, 18 different cases of surface have been analyzed. Based on each surface, the accuracy of 11 different interpolation algorithms has been examined in seven different configurations of search sectors, for 4 different values of measurement points density. The study was conducted for two cases, where entire surfaces were taken into account (together with the edge nodal points) and selected subareas on the surfaces (without edge nodes). Altogether, 11,088 calculated values were obtained, which were subsequently analyzed. Due to the limited size of this article, it is not possible to present all of the analyzed cases. For this reason, presentation of the study results was limited to 6 surfaces, 8 algorithms of interpolation, 7 sector systems and 4 densities for the 2 mentioned cases. Since the rule of procedure and analysis is the same for all cases, a detailed description of methodology has been presented for one selected surface with a different configuration of the measurement points.
From among all analyses, examples representative of the majority of the analyzed cases were selected. These examples illustrate the procedure of determining the accuracy of interpolation algorithms and estimating the effect of the location of measurement points on the calculation results. The results of the study provide an example of an analytical procedure and a presentation of the possibility of analyses performed in order to improve the quality of DTMs.

Test models
Special theoretical test models were used in order to analyze the results of interpolation conducted on different surfaces by different algorithms, with different arrangements of search sectors and different densities of measurement points. The model construction principle was described for one selected surface. The function of two variables (1) necessary to generate the mathematical surface presented in Fig. 1a was used for model development. Based on that function, within the specified range of coordinates x, y (1), pseudo-measurement points (further referred to as measurement points) were generated at random.
They formed a standard surface, which was then rescaled and shifted to obtain a DTM of 500 m by 700 m with an area of 350,000 m 2 , containing exclusively measurement points with positive values of x, y and z coordinates (Fig. 1b). The measurement points generated in this way supported the calculation of practical values at nodal points using a variety of interpolation algorithms. Theoretical nodal points were determined based on function (1). The base square of the grid was determined in such a way that the morphology of the DTM generated using function (1) should not allow a precise interpolation of values at all nodes. The applied procedure makes it possible to obtain different accuracies during the interpolation of nodes in morphologically different parts of the model. Therefore, it enables an analysis of the entire spectrum of interpolation accuracy, which is lower at morphologically different parts of the model and higher at morphologically similar sites. Each of the tested interpolation algorithms responds differently to different surface morphology. Such an approach supports an accuracy analysis of a given algorithm for different cases. The spacing S between nodes was 10 m for each tested case (Fig. 1c). Next, based on the same function (1), theoretical nodular points were generated with the same positioning parameters (10 m 9 10 m) as the interpolated nodes. On the basis of the test model developed in this way, analyses were carried out to compare the practical value calculated by interpolation with the theoretical value calculated using the function for all nodular points of the GRID structure. All points generated based on a mathematical function (both measurement points and nodes) are theoretical points whose values have been rounded off. In the discussed example, the values were calculated to eight decimal places and are presented accurate to two decimal places (RMS). The error in the determination of a theoretical value (based on the function) at a given node is several orders of magnitude smaller than the interpolation error at the same node, which is why it has no significant effect on the presented results.
Based on different functions of two variables, there were 18 surface models generated using the methodology described above. Six surface models (S1-S6) (Fig. 2) were selected for presentation of analyses shown in this article. The selected models have different morphology, where changes in the surface morphology run in different directions. The selected subareas (S1-s-S6-s), which will be used for analyses described in chapter 5, were marked on all the surface models. The first algorithm (R2) assumed calculation of the value at the nodal point based on the sum of the values of n measurement points balanced by the reverse square of the distance from the node. The second algorithm (KR) assumed calculation of the value in the node based on the semi-variogram adjusted to the height distribution of measurement points around the node. The third algorithm (LP) calculated the value at the nodal point using the function of first-level local polynomial. The fourth algorithm (RF) for interpolation of the value at the nodal point used multiquadric radial functions. The fifth algorithm (PL) employed for interpolation of a plane which was approximated based on the nearest measurement points and a line running through a node. The sixth algorithm (TR) used triangulation with linear interpolation for calculations. In the seventh algorithm (MC), values in a node were calculated by the minimum curvature method. In the eighth algorithm (MA), values in a node were interpolated by the moving average method.
For comparative purposes, all algorithms relied on the same database of measurement points and calculations were performed for identically located nodal points. In order to determine the influence of measurement points' position around the node on the calculation results, seven cases of their positioning were analyzed (Fig. 3). Only in the first case presented in Fig. 3 were measurement points chosen from the entire space around the node without division into sectors.
In the other cases, the search space was divided into equal sectors. The angle of lines creating the first sector was fixed at 0°. The search radius around individual nodes was not limited and was the same in all directions.
In order to analyze the effect of different density of measurement points on the accuracy of interpolation by different algorithms, four submodels consisting of a different number of measurement points were generated for each analyzed surface. The RMS (2) coefficient was applied to compare the matching accuracy of the interpolated surface and the standard surface (Cressie and Pavlicova 2002;Schabenberger and Gotway 2005).
where • f(x, y)-value of function (1) in the theoretical nodal point with coordinates x y, • z-value calculated using a given interpolation algorithm on the basis of measurement points at the nodal point with coordinates x y, • n-number of nodal points.
A differential diagram was used to present the areas of deformation in the analyzed surface models. The diagrams present the absolute differences in the values calculated in all nodes between the practical interpolation model-DTM ''P'' (Fig. 5a)-and the theoretical function model-DTM ''T'' (Fig. 5b). Figure 5c presents the differential diagram representing the ranges of values of the differences determined based on the values calculated for each node. The RMS coefficient in [m], calculated jointly for the entire differential surface, is shown above the diagram. The legend under the diagram shows the assumed ranges of differential values in (m).
The use of differential models enables to locate differences in interpolation models as compared to the theoretical model. The application of the RMS coefficient allows, on the other hand, to determine the accuracy of matching the entire surface by means of a single value. That, in turn, offers the possibility of comparing, on a global scale, different interpolation surfaces by showing on the graphs the dependencies between the location of measurement points and the values Further analysis was divided into two steps. In the first one, the situation where all nodal points of the test model were taken into consideration was analyzed. In the second one, internal points of models selected using a rectangle were analyzed to eliminate the effect of edge nodes on the final results.

Analysis of the accuracy of DTMs for all nodal points
Using the differential diagrams, the influence of measurement point density on DTMs accuracy was analyzed. The diagrams presented in Fig. 6 show the location of deformations at different measurement point densities for two selected interpolation algorithms (R2 and KR). In all interpolation and density cases, 16 measurement points closest to the node were used. The measurement points were searched across the entire surface area around the node, without a division into sectors. For a later comparison of results, both algorithms relied on the same measurement points during interpolation. Interpolation using the R2 algorithm at the density of 0.3 p/100 m 2 (Fig. 6a) caused relatively large deformations across the entire surface (RMS = 2.51 m).
Increasing the measurement point density to 14.3 p/100 m 2 (Fig. 6b) significantly improved the matching of the model (RMS = 1.06 m). A similar situation occurred when the KR algorithm was used for interpolation. At a measurement point density of 0.3 p/100 m 2 (Fig. 6c), interpolation yielded clearly worse results, compared with a density of 14.3 p/100 m 2 (Fig. 6d). A comparison of the diagrams and the RMS coefficient for both algorithms revealed a clear advantage of the KR algorithm over the R2 algorithm for both a lower and higher density of measurement points. As in the discussed cases, the selection of points was carried out without a division into sectors and the influence of edge nodes on the accuracy of interpolation models was not observed. Similar dependences were obtained for the remaining algorithms. For all tested algorithms, interpolation accuracy increased with an increase in the density of measurement points. The diagrams shown in Fig. 7 show the influence of the number of search sectors on DTM accuracy. Further examples present deformations after interpolation with the RF algorithm. The density of measurement points in all cases was 0.3 p/100 m 2 . In all cases, 16 measurement points were used and were located in the following way: without sectors (Fig. 7a), in 4 sectors (Fig. 7b), in 8 sectors (Fig. 7c) and in 16 sectors (Fig. 7d).
The number of nodes rejected as a consequence of inaccuracies in calculations caused by an insufficient number of measurement points in individual sectors is given above the diagrams, next to the number of sectors. In consecutive diagrams, an increasing share of rejected edge nodes with an increasing number of sectors can be seen. On the other hand, in the area inside the diagrams, the accuracy increases with an increase in the number of sectors. Figure 7 shows an example representative of the entire group of tested algorithms-all algorithms display similar trends.
In order to analyze the effect of the number of search sectors and density of measurement points on the accuracy of interpolation by different algorithms, the RMS coefficient was used, whose values are given in Fig. 8. The diagrams shown in Fig. 8 were developed for individual surfaces (S1-S6) shown in Fig. 2   interpolation algorithms with symbols given in the corner (the symbols are explained in chapter 3). The vertical axis shows the values of the RMS coefficient, whereas the horizontal axis shows the number of sectors, and the density of measurement points is represented by the series lines. In order to make a uniform comparison of the accuracy for all the analyzed DTMs, all the diagrams present situations when one measurement point situated closest to the node is found in every sector.

and subsequent
An analysis of the graphs shows a clear improvement in the accuracy of all models with an increase in measurement points density (the RMS coefficient decreases). In all graphs, the influence of the number of sectors on RMS coefficient values decreases with an increase in the density of measurement points, indicating that the dependence of DTMs quality on the number of sectors decreases. This is a consequence of the fact that the number of points closer to the node increases with the density of the measurement points.
For the density of 0.3 p/100 m 2 and 1.4 p/100 m 2 for most diagrams, there is a clear effect visible of the number of sectors on the DTM quality. This applies especially to R2, KR and RF algorithms on all the presented surfaces (S1-S6). The RMS coefficient in those cases decreases gradually with an increasing number of search sectors. When the density of measurement points is higher (5.7 p/100 m 2 and 14.3 p/100 m 2 ), distinct stabilization of the model accuracy takes place on the level of 4 sectors. This applies especially to R2, KR, RF and MA algorithms. The higher stability of the accuracy of all the DTMs, regardless of the number of sectors, can be observed for TR and MC algorithms, where the RMS coefficient remains at the same level. For most cases, regardless of the surface morphology (from S1 to S6), it is typical that the accuracy is improved during interpolation with the use of 4 search sectors. When the number of search sectors is further increased, the value of the RMS coefficient increases, that is, the accuracy of DTM decreases. This applies mainly to LP, PL and MA algorithms. More points may be found when the number of sectors increases, but the position of the points relative to the node is less advantageous (they are usually farther from it), which has a decisive effect on interpolation with those algorithms. For most cases, when the number of search sectors is increased above 4, the model accuracy is not significantly improved. Increasing the number of search sectors considerably extends the time of interpolation between individual interpolations. The analyzed examples show that it is the optimum solution to determine 4 search sectors for most examined interpolation algorithms, regardless of the surface morphology.

Analysis of DTM accuracy for selected nodal points
Previous analyses accounted for the influence of edge nodes on the calculation results. Since the number of rejected edge nodes has a significant effect on the values of the RMS coefficient, a situation should be analyzed in which the number of measurement points around the nodes is always sufficient for their determination. The number of sectors affects the number of rejected edge nodes (Fig. 7). An increase in the number of sectors is accompanied by an increase in the number of nodes that cannot be determined due to an insufficient number of measurement  The examples discussed in this section support the analysis of a situation when the effect of errors due to edge nodes has been eliminated and the size of GRID has been adjusted to terrain morphology.
In the second stage of the study, to eliminate the influence of edge points, an area selected using a 250 m 9 290 m rectangle (i.e., a surface area of 72,500 m 2 ) was analyzed (Fig. 9a, b). Additionally, the area for analysis was selected in a way to minimize the influence of a mismatch between the size of the base square of the node grid and the morphology of the DTM. Figure 9c presents a differential diagram showing model deformation after interpolation using the R2 algorithm for 6 search sectors. Above the diagram, the RMS coefficient is shown for the entire interpolation area after the elimination of edge points. Figure 9d shows the differential diagram generated for selected 672 (24 9 28) nodal points calculated in the same way as above. The lower the RMS coefficient, the better the model matching. The entire area (Fig. 9c) contains more differences; therefore, the RMS is higher.
The selected fragment is free from the effect of edge errors, and the size of its GRID is better adjusted to terrain morphology at the analyzed site; thus, the RMS is lower. The RMS coefficient presented above the diagram indicates a much more accurate match of the area selected using the rectangle to the standard surface. The coefficient value is lower because the selected area had a base square size better matching the DTM morphology and the influence of errors developing at edge points was eliminated. A similar situation also occurs in the cases presented in Fig. 10. Further differential diagrams present the results of analyses for different interpolation algorithms: R2 (Fig. 10a), KR (Fig. 10b), LP (Fig. 10c) and RF (Fig. 10d).
In all cases, 16 measurement points, 4 search sectors and a density of 14.3 p/100 m 2 were used for interpolation of each node. The worst matching of the DTM Fig. 9 Selection of measurement points in the DTM and in the differential diagram Interpolation of the digital terrain model 525 was achieved for the LP algorithm, while the best was offered by KR and RF algorithms. Similar to the preceding case, the dependence of the RMS coefficient on the number of sectors was analyzed for selected areas. To this end, subsequent diagrams ( Fig. 11) of the relationships between RMS and the number of sectors were formed with a different density of measurement points, for 6 selected subareas (S1-s-S6-s), shown in Figs. 2, and 8 interpolation algorithms with symbols given in the corner (the symbols are explained in chapter 3). As in the previous case, in order to perform a uniform comparison of the accuracy for individual subareas, all the diagrams show situations when one point located closest to the node is found in each sector. The vertical axis of the diagrams shows the values of the RMS coefficient, while the horizontal axis shows the number of sectors, and the density of measurement points is represented by the series lines. Due to area selection, interpolation results for all algorithms were corrected. Since the quality of DTMs (lower RMS) increased considerably, the scale of the vertical axis was changed to improve legibility.
Also, in this case, an improvement in the model quality can be observed in all graphs with an increase in the number of measurement points. However, unlike in the previous case, a significant increase in the DTM's quality was not observed with increasing density of the measurement points. Similar characteristics of the RMS values can be observed in all the diagrams (Fig. 11) for individual algorithms, regardless of the area (S1-s-S6-s) and in all the surfaces. Stabilization of the model accuracy (similar values of RMS) was visible, especially for PL, TR and MA algorithms. This group of algorithms is especially sensitive to edge node errors and eliminating them allowed similar results to be achieved, regardless of the number of search sectors. The most significant effect on the change of the number of sectors was observed in R2, KR, LP and RF algorithms. In those cases, stabilization of the model accuracy takes place with 4 search sectors. A further increase in the number of sectors does not improve the quality of the TMs. For LP and MC algorithms, increasing the number of search sectors above 4 resulted in decreasing the accuracy of the surface models. Those algorithms are the most sensitive to location (particularly distance) of the points around a node.  The analyzed examples show that after eliminating the effect of edge nodes, establishing 4 search sectors for most examined interpolation algorithms still remains the optimum solution, regardless of the subarea surface morphology. Only for PL, TR and MA algorithms is it possible to claim that after eliminating the edge nodes and fitting the GRID to morphology, the number of sectors used for finding measurement points does not significantly affect the accuracy of the calculations or the quality of interpolation models. The values of RMS for those algorithms are similar for different numbers of sectors.

Comparison of the accuracy of interpolation algorithms
For the purpose of comparing the accuracy of interpolation algorithms, DTMs generated for all (Fig. 12) and selected (Fig. 13) nodal points were analyzed. The interpolation was carried out for measurement points positioned in 4 search sectors for each nodal point. The vertical axis shows the RMS coefficient value, and the density was placed along the horizontal axis while the series consisted of consecutive interpolation algorithms. The surface to which each diagram applies (as per Fig. 2) is provided in the corners.
The graphs illustrate an increase in DTM accuracy with increasing density of measurement points. A clear improvement in accuracy for all algorithms was noted during interpolation performed for selected points (free from edge errors and with an adjusted size of the base square) (Fig. 13).
This list allows a comparison of the accuracy of each algorithm with the same interpolation parameters. Comparison of the RMS values shows that the most accurate interpolation, both for the entire surfaces and for the selected subareas, for an entire surface was performed with the TR algorithm. Accurate models were developed also by means of KR, RF and TR algorithms. MA algorithm proved to be equally accurate, but only for higher densities of measurement points. The lowest accuracy of DTMs, both for the entire surfaces and for the selected subareas, was achieved with R2 and MC algorithms. The accuracy of the interpolation models achieved with the other algorithms is comparable. In general, it can be stated that the RMS ratio decreases (hence increases the accuracy of the model, Figs. 12 and 13) when edge effects are eliminated; however, the relative differences between the algorithms remain constant. Eliminating the edge effects has little effect on the relative differences in accuracy between the interpolation algorithms, however, significantly affects the accuracy of the models in all cases.
For the purpose of their practical validation, the obtained results were used to develop real DTMs. The models were generated using all sector variants and methods of measurement point determination. Attention was paid to the duration of calculations and the quality of the generated DTMs. The models shown in Fig. 14 were based on the measurement points provided by a multibeam echo-sounder while probing the sea bottom. Interpolation was performed with 4 search sectors, using RF algorithm, which provided an accurate interpolation model for the selected subareas. Edge nodes were eliminated from the models. When the number of sectors was lower than 4, the quality of the DTMs deteriorated. With an increase in the number Comparison of accuracy of interpolation algorithms (with 4 sectors) for 6 selected subareas with no edge nodes of sectors, the calculation process was more time-consuming. When the number of search sectors was limited to 4, the calculation process was shorter and the high quality of DTMs was maintained. The elimination of edge nodes reduced deformations along the edges of the models. This method improved the quality of work and optimized the time needed to develop individual digital terrain models.

Conclusions
The methodology applied in the study enabled an interpretation of the results and the choice of the optimum interpolation parameters for individual algorithms while using morphologically diverse surface models. The application of differential diagrams enabled the location of deformations of interpolation models. Such diagrams allow the presentation of the scope of deformations within the determined accuracy ranges. The application of the RMS coefficient, on the other hand, allows the accuracy of matching of the entire surface to be determined using a single value. That, in turn, enables comparative analyses to be performed of the influence of the distribution of measurement points on the accuracy of the generated DTMs. During generation of digital terrain models using the GRID structure, a significant role is played by correct selection of measurement points used for calculations. Such selection is possible using sectors. Eighteen different cases of surface morphology and 11 different interpolation algorithms were analyzed. Examples representative of a larger group of analyzed cases were selected. The determination of the optimum number of search sectors was important for the majority of analyzed cases. As the analyzed examples show, the optimum solution involves establishing 4 search sectors for most of the examined interpolation algorithms, regardless of the surface morphology. A similar situation exists after eliminating the effect of edge nodes and fitting the size of the base network of node location to the morphology of a surface by the choice of a subarea. During studies aimed at determining interpolation accuracy, one should consider the influence of the edge effects. Each of the tested algorithms has allowed to increase the accuracy of the interpolation after eliminating the influence of the edge effects. The adjustment of the base square to the morphology of the studied surface also has a major influence on the results. The selection of a less morphologically diversified area using a rectangle enabled better adjustment of the base square of the grid (10 m 9 10 m) to surface morphology. As a result, the accuracy of the interpolation model increased for all applied algorithms. Interpolation accuracy using the optimum number of search sectors was confirmed in practice during interpolation on a real digital terrain model. Due to the fact that the optimum number of 4 search sectors was used for each node, it was possible to achieve the desired quality of DTMs with simultaneous optimization of the required development time.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.