Evaluation of DTM Generation in Surfer 8.0


Majid Hamrah
Assistant Professor
Islamic Azad University Estahban Branch
Estahban, Iran.
Hamrah@kntu.ac.ir

Davod Shojaee
Photogrammetry MSc Students
Dshojaee@gmail.com

Azam Mosavi
GIS MSc Student
Dept. of GIS Eng. K.N. Toosi University of Technology
Vali_asr St,
Tehran, Iran,
P.C. 19697
A_mousavi_m@yahoo.com


Abstract
Different data sources for the determination of a Digital Terrain Model (DTM) are available. New measurement techniques like Laser Scanning (LS) allow a high degree of automation and dense sensing. They pushed the development of new methods for data modelling. There are many methods for surface generation. Surfer 8.0 is powerful software for DTM generation, contouring and visualization. In Surfer there are 12 methods that each having specific functions for interpolation. In this paper accuracy of DTM generation for 3 subsets of data were investigated.

1-Introduction
The relief of the terrain surface is a continuous variation. It provides difficulties in representation because it is a continuous feature and three dimensional. Representation of relief is probably the most difficult problem in mapping and stems from the basic need to present the three dimensions of the terrain in two dimensional. Methods used to overcome this difficulty can be considered under three main headings: Land Surveying, Photogrammetry and Map Compilation. Height information may be provided by precise survey methods of leveling in which the exact height of a set of points is determined. Much mapping today is carried out from aerial photography. The advantage of a view of the three-dimensional scene in the overlapping pair of photographs tends to construct of "space model" in which the contour lines can be plotted by steering along them the floating point in the plotting machine. More advanced methods of mapping from air photographs enable contours to be plotted automatically. In the past centuries cartographers tried to implement different qualitative techniques to represent the relief.

With the invention of computer itself and development of digital methods and its potentials in handling huge data, the representation of continuous surface of terrain was possible, so the map users could perceive the terrain topography more effectively and also easier. The extraction of topographic data in digital form is considered under the topic of "DTM". The DTM is a computer representation of the earth surface and such provides a base data from which topographic parameters can be digitally generated. The DTM as its name implies, is a model of the elevation surface. In fact, The DTM can be defined as any digital representation of the continuous variation of relief over space. If the elevation data is sampled at point locations or digitized from the contour lines, a primary DTM generation concern is the interpolation method. A search for optimal DTM interpolation method has been of interest for quite some time.

2-Methods of Interpolation
Surfer uses 12 interpolation methods for create DEM. Some of these interpolators are: Inverse Distances to a power, Kriging, Minimum Curvature, etc. In the following sections all methods were reviewed and evaluated.

2-1-Inverse Distance to a Power
The Inverse Distance to a Power gridding method is a weighted average interpolator, and can be either an exact or a smoothing interpolator. With Inverse Distance to a Power, data are weighted during interpolation such that the influence of one point relative to another declines with distance from the grid node. Weighting is assigned to data through the use of a weighting power that controls how the weighting factors drop off as distance from a grid node increases. The greater the weighting power, the less effect points far from the grid node has during interpolation. As the power increases, the grid node value approaches the value of the nearest point. For a smaller power, the weights are more evenly distributed among the neighboring data points. Normally, Inverse Distance to a Power behaves as an exact interpolator. When calculating a grid node, the weights assigned to the data points are fractions, and the sum of all the weights is equal to 1.0. When a particular observation is coincident with a grid node, the distance between that observation and the grid node is 0.0, and that observation is given a weight of 1.0, while all other observations are given weights of 0.0. Thus, the grid node is assigned the value of the coincident observation. The Smoothing parameter is a mechanism for buffering this behavior. Inverse Distance to a Power is a very fast method for girding.

2-2-Kriging
Kriging is a geostatistical girding method that has proven useful and popular in many fields. This method produces visually appealing maps from irregularly spaced data. Kriging attempts to express trends suggested in data, so that, for example, high points might be connected along a ridge rather than isolated by bull's-eye type contours. Kriging is a very flexible gridding method [2]. By accepting the Kriging defaults, an accurate grid of the data can be obtained. Within Surfer, Kriging can be either an exact or a smoothing interpolator depending on the user-specified parameters. It incorporates anisotropy and underlying trends in an efficient and natural manner.

2-2-1-Variogram Overview
Surfer includes an extensive variogram modeling subsystem. This capability was added to Surfer as an integrated data analysis tool. Variogram modeling may also be used to quantitatively assess the spatial continuity of data even when the Kriging algorithm is not applied.Variogram modeling is not an easy or straightforward task. The development of an appropriate variogram model for a dataset requires the understanding and application of advanced statistical concepts and tools. In addition, the development of an appropriate variogram model for a dataset requires knowledge of the tricks, traps, pitfalls, and approximations inherent in fitting a theoretical model to real world data. Skill with the science and the art are both necessary for success. The development of an appropriate variogram model requires numerous correct decisions. These decisions can only be properly addressed with an intimate knowledge of the data at hand, and a competent understanding of the data genesis (i.e. the underlying processes from which the data are drawn).

2-2-2-The Variogram
The variogram is a measure of how quickly things change on the average. The underlying principle is that, on the average, two observations closer together are more similar than two observations farther apart. Because the underlying processes of the data often have preferred orientations, values may change more quickly in one direction than another. As such, the variogram is a function of direction.

The variogram is a three dimensional function. There are two independent variables (the direction q, the separation distance h) and one dependent variable (the variogram value g (q, h)). When the variogram is specified for Kriging, by setting the sill, range and nugget parameters the anisotropy information can be specified. The variogram grid is the way this information is organized inside the program.

The variogram (XY plot) is a radial slice (like a piece of pie) from the variogram grid, which can be thought of as a "funnel shaped" surface. This is necessary because it is difficult to draw the three-dimensional surface, let alone try to fit a three dimensional function (model) to it. By taking slices, it is possible to draw and work with the directional experimental variogram in a familiar form a XY plot.

Remember that a particular directional experimental variogram is associated with a direction. The ultimate variogram model must be applicable to all directions. When fitting the model, the user starts with numerous slices, but must ultimately mentally integrate the slices into a final 3D model.

2-3-Minimum Curvature
Minimum Curvature is widely used in the earth sciences. The interpolated surface generated by Minimum Curvature is analogous to a thin, linearly elastic plate passing through each of the data values with a minimum amount of bending. For DTM generation the Minimum Curvature, however is not an exact interpolator, but can bring the most possible smooth surface. Minimum Curvature produces a grid by repeatedly applying an equation over the grid in an attempt to smooth the grid. Each pass over the grid is counted as one iteration. The grid node values are recalculated until successive changes in the values are less than the Maximum Residuals value, or the maximum number of iterations is reached (Maximum Iteration field).

2-4-Modified Shepard's Method
Modified Shepard's Method uses an inverse distance weighted least squares method. As such, Modified Shepard's Method is similar to the Inverse Distance to a Power interpolator, but the use of local least squares eliminates or reduces the "bull's-eye" appearance of the generated contours. Modified Shepard's Method can be either an exact or a smoothing interpolator.

2-4-1-Quadratic Neighbors
The Modified Shepard's Method starts by computing a local least squares fit of a quadratic surface around each observation. The Quadratic Neighbors parameter specifies the size of the local neighborhood by specifying the number of local neighbors. The local neighborhood is a circle of sufficient radius to include exactly this many neighbors.

2-4-2-Weighting Neighbors
The interpolated values are generated using a distance-weighted average of the previously computed quadratic fits associated with neighboring observations. The Weighting Neighbors parameter specifies the size of the local neighborhood by specifying the number of local neighbors. The local neighborhood is a circle of sufficient radius to include exactly this many neighbors.

2-5-Natural Neighbor
The Natural Neighbor gridding method is quite popular in some fields. Consider a set of Thiessen polygons (the dual of a Delaunay triangulation). If a new point (target) were added to the dataset, these Thiessen polygons would be modified. In fact, some of the polygons would shrink in size, while none would increase in size. The area associated with the target's Thiessen polygon that was taken from an existing polygon is called the "borrowed area." The Natural Neighbor interpolation algorithm uses a weighted average of the neighboring observations, where the weights are proportional to the "borrowed area. “The Natural Neighbor method does not extrapolate contours beyond the convex hull of the data locations (i.e. the outline of the Thiessen polygons).

2-6-Nearest Neighbor
The Nearest Neighbor gridding method assigns the value of the nearest point to each grid node. This method is useful when data are already evenly spaced, but need to be converted to a Surfer grid file. Alternatively, in cases where the data are nearly on a grid with only a few missing values, this method is effective for filling in the holes in the data.

2-7-Polynomial Regression
Polynomial Regression is used to define large-scale trends and patterns in the data. Polynomial Regression is not really an interpolator because it does not attempt to predict unknown Z values. There are several options to define the type of trend surface.

2-8-Radial Basis Function
Radial Basis Function interpolation is a diverse group of data interpolation methods. In terms of the ability to fit the data and to produce a smooth surface, the multiquadric method is considered by many to be the best. All of the Radial Basis Function methods are exact interpolators, so they attempt to honor the data. A smoothing factor can be introduced to all of the methods in order to produce a smoother surface.

2-8-1-Function Types
The basis kernel functions are analogous to variograms in Kriging. The basis kernel functions define the optimal set of weights to apply to the data points when interpolating a grid node.

2-9-Triangulation with Linear Interpolation
The Triangulation with Linear Interpolation method in Surfer uses the optimal Delaunay triangulation. The algorithm creates triangles by drawing lines between data points. The original points are connected in such a way that no triangle edges are intersected by other triangles. The result is a patchwork of triangular faces over the extent of the grid. This method is an exact interpolator.

Each triangle defines a plane over the grid nodes lying within the triangle, with the tilt and elevation of the triangle determined by the three original data points defining the triangle. All grid nodes within a given triangle are defined by the triangular surface. Because the original data are used to define the triangles, the data are honored very closely. Triangulation with Linear Interpolation works best when the data are evenly distributed over the grid area. Datasets that contain sparse areas result in distinct triangular facets on the map.

2-10-Local Polynomial
The Local Polynomial gridding method assigns values to grid nodes by using a weighted least squares fit with data within the grid node's search ellipse.

2-11-Data Metrics
The collection of data metrics gridding methods creates grids of information about the data on a node-by-node basis. The data metrics gridding methods are not, in general, weighted average interpolators of the Z-values.

Data metrics use the local dataset including break lines, for a specific grid node for the selected data metric. The local dataset is defined by the search parameters. These search parameters are applied to each grid node to determine the local dataset.

2-12-Moving Average
The Moving Average gridding method assigns values to grid nodes by averaging the data within the grid node's search ellipse.

2-12-1-Search Ellipse
To use Moving Average, define a search ellipse and specify the minimum number of data to use. For each grid node, the neighboring data are identified by centering the search ellipse on the node. The output grid node value is set equal to the arithmetic average of the identified neighboring data. If there is fewer than the specified minimum number of data within the neighborhood, the grid node is blanked.

3-Results
For evaluating of DTM generation three types of region applied. The datasets collected from respectively Laser Scan, photogrammetry and ground surveying. The surfaces were made from these datasets. The first dataset is related to a rather flat region and second dataset is belonged to a semi-mountainous region and third dataset is related to a region with various topography. Figure 3-1 shows the surfaces of three datasets.


Figure 3-1. Surfaces of dataset


At first, each of the twelve methods on each of these data was tested and then the RMSEs were calculated. The results of these methods are displayed in table 3-1.

Table 3-1. Results of datasets from 12 methods


Table 3-1 shows that result of first dataset of Kriging method has a better RMSE with respect to other methods. In second dataset, Radial Basis has the lowest RMSE. So In third dataset this method has a rather quite consequence than the other methods. And Minimum Curvature in third dataset has a better result with respect to other methods.

In Kriging method with using different variogram functions for first dataset, following results were obtained. Table 3-2 shows these results.

Table 3-2. Results of different variogram functions for first dataset


This table shows that Linear and Power variogram functions Kriging has better results.
For second dataset, Radial Basis method has better results. This method has 5 Basis Function (Inverse Multi quadric, Multilog, Multi quadric, Natural cubic spline and thin plate spline). For second dataset different types of basis function was evaluated. Table 3-3 shows these results.

Table 3-3. Results of different Basis functions for second dataset


This table shows that Thin Plate Spline for second dataset has a better RMSE with respect to other Basis function.
For analyzing other Radial Basis Functions the third dataset was applied. The Result of third dataset shows that Multi Quadratic has better consequences in this dataset. Table 3-4 shows these results.

Table 3-4. Results of different Basis functions for third dataset


4-conclusion
On the basis of RMSE minimizing, a test has been carried out to find out an appropriate interpolation method for region of various topography. Three datasets, belonging to the three regions were examined in different interpolation functions. The results showed that radial basis function and thin plate spline had minimum RMSE, although the linear variogram of kriging gave nearest result.

For semi-mountainous region the radial basis function and minimum curvature function presented the suitable result. The power and linear variogram methods of Kriging had nearly the best RMSE.

It can be concluded that in flat region the linear variogram and power with drift type linear of kriging method give an acceptable result. In general, radial basis function is suitable for various topography regions.

References

  1. Surfer 8.0, Help Documents, www.goldensoftware.com
  2. The Discontinuous Nature of Kriging Interpolation for Digital Terrain Modeling, Thomas H. Meyer, Department of Natural Resources Management, and Engineering, University of Connecticut, Storrs, Connecticut, W.B
  3. Franke R. Scattered Data Interpolation: Test of Some Methods [J]. Mathematics of Computations, 1982,33(157):181.
  4. Franke R, Nielson G. Smooth Interpolation of Large Sets of Scattered Data [J]. International Journal for Numerical Methods in Engineering, 1980,15(2):1691.