INVERSE DISTANCE WEIGHT SPATIAL INTERPOLATION FOR TOPOGRAPHIC SURFACE 3D MODELLING

AbstrackTopographic modelling is an important aspect in engineering work and planning such as dam design, road planning, volumetric calculation, water flow analysis and so forth. To get topographic data, usually a land survey or topographic measurement is taking for a region or an area under study. A number of points that represent an area are measured to get a height dataset. The dataset which is consist of height points will be used to model topographic condition of the area with generating contour map and also 3D model. There are some methods can be used to generate topographic surface in 3D model like linear interpolation in Triangular Irregular Network (TIN), Kriging and Inverse Distance Weight (IDW). This research implemented IDW spatial interpolation algorithm to model earth surface or topographic model in 3D. The IDW method was implemented in Python programming language using Numpy library for computation and Plotly graphic library to visualize the 3D Model. Using 2500 and 10000 interpolation points with 100 random sampling points that extracted from Digital Surface Model (DSM), IDW was successfully estimated the height at unsampled locations. The results show, more higher interpolation point number will produce more detail surface texture.


Introduction
Topography could be described as a portrait of earth surface that shows all features on earth surface both natural and man-made features including the height of the earth surface. Commonly, all features are shown in a topographic map. Therefore topographic map is like a mirror of earth surface in a plain media or two dimension (2D). The height of earth surface on a topographic map is represented by contour lines. Contour line is a height isoline which has the same value along the line. Using contour lines ones could imagine the shape, pattern or height variation of the earth surface for a region or area. Along with computer and software advancement, currently the earth surface can be modelled and visualize in three dimension (3D). Thus ones not just can imagine how the earth surface looks like in 3D, but moreover it can be seen directly through a display device. Topographic surface in 3D is more vivid and realistic compare to 2D visualization. According to Stout and Blunt [1], there are some benefits can be obtained with 3D topographic visualization such as: it can represent natural characteristics of a surface, have powerful function and flexibility as required by user, the information obtained from 3D is significantly greater than can be achieved by 2D and more reliable and representative to represent earth surface. To model earth surface a height dataset is required. A height dataset consist of representative height points that cover a region to be mapped or modelled. The height of a point is measured using a measurement device like total station or theodolite. The measurement activity in order to get the height of a point is commonly called topographic mapping. The quality of 3D model that derived from height point measurement depends on the number of points that measured over an area and point's location. In general thumb, a more dense point that cover an area to be mapped, will give a better result that close to actual surface condition. But more points to be measured will require more cost and time. Point's location is another factor that need to be considered in a topographic mapping. Poorly placing points will cause 3D surface deviates with actual condition. In short, placing points in right location where gradient occurs and adequate number of points that cover an area to be mapped are keys to get a good topographic 3D model. The height dataset that taken from field measurement will be used and processed to get a contour map in 2D and also to generate 3D model. To generate contour or 3D model, a number of interpolation points at unsampled location must be populated to get a dense height point through the whole study area. There are some methods that can be used to estimate interpolation point namely: linear interpolation, IDW and Kriging. Linear interpolation is the simplest interpolation method to estimate an unsampled point along a line between two known points. In the implementation step, three points are connected to each other with three lines and formed a triangle which is called Triangular Irregular Network (TIN). The interpolation points will take place on all the lines of TIN. This method is simple and fast to implement, but the drawback in 3D topographic modelling could produce a coarse surface visualization, because interpolation points only take place along the TIN lines. Kriging is another interpolation method which is more complex and requires a longer time in computation. Kriging is a stochastic interpolation approach which considers inter-point correlation in the weight determination. Therefore the weight in Kriging is not only determined from the closeness of surrounding sampling points that involve in calculation, but also the variance-covariance value. Sampling point which has high covariance will get more weight than a lower one. IDW is a spatial interpolation approach that is used commonly to estimate an unsampled or unmeasured variable at any location in a study area. IDW is a deterministic interpolation approach which considers the distance of an unsampled point towards a set of surrounding sampling points in weights determination stage. In contrast with stochastic interpolation approach like Kriging, which uses inter-points correlation in weight determination, IDW is more simple and fast in computation. In this paper the algorithm of IDW was implemented in Python programming language using Numpy numeric library for computation and Plotly graphic library to visualize 3D topographic surface model.

Methodology
To estimate an unsampled point using IDW method, a number of sampling points are required. Based on the distance between an unsampled point to each sampling point that involves in the calculation process, the weight will be determined with Shepard's rule [2] as in equation 1, where wi is the weight of i-th sample point, dij is the Euclidian distance from the unsampled to the sampling point which can be calculated using equation 2. Lastly p stands for power. In Shepard's rule the default p value is 2. (1) After getting all weights for each sampling point, the weight will be used to estimate the value at the unsampled location using the equation 3. To ensure an unbiased interpolator the sum of weights needs equal to 1 [3].
As mentioned above, some points are required to perform IDW interpolation. There are two approaches can be chosen in selecting sampling point around an estimated point location. First approach can be done with a defined radius and the second one just define a minimum number of points to be selected. This research will implement the second approach. Figure 1 shows the second approach implementation flowchart [4]. The process starts with initiating 3 input variables, there are initial radius, p values, minimum number of selected points and location of unsampled point to be estimated. Based on initial radius then a filter block will be defined using minimum and maximum coordinate of the block, where minimum coordinate can be obtained with subtracting abscissa and ordinate of estimation point by radius value, and vice versa the maximum coordinate can be obtained by adding respected abscissa and ordinate with initial radius value. After the block is defined, then the algorithm will test if a point is within the block. A sampling point will be selected if it lies between the minimum and maximum block coordinate. In each selection step, the number of selected points will be counted. If the number of selected points less than defined input variable, then the size of filtered block will be enlarged by adding radius value with a distance. The process will iterate in a loop until the number of selected sampling points reach the minimum input value. After selecting the sampling points, then the algorithm will continue to calculate the weights. In this step, firstly the Euclidian distance between the unsampled point to each selected sampling point will be calculated. If the distance equal to 0, it means the unsampled point's location is the same with the current sampling point, therefore the sampling point will be excluded. Otherwise the algorithm will proceed to calculate the weight for the sampling point. When weight calculation step is finished, then the algorithm will continue to estimate the unsampled value using equation 3. The IDW method that described above is implemented in Python 3.7.4 programming language. All the steps in the flowchart will be put in a function. For each unsampled point to be interpolated or estimated the function will be called. The unsampled/estimation points are generated through an area to be modelled with a certain interval. More closer the interval will generate more interpolation points and vice versa. In Figure 2 can be seen 5 m interval interpolation points in an area to be modelled.

Result and Discussion
The algorithm that had been implemented in Python was used to model a 3D surface of an area. The dataset which used was 100 random height point that extracted from Digital Surface Model (DSM) of garbage disposal dumping site at Gampong Jawa, Banda Aceh municipality. The DSM was generated from drone survey on September 9 th , 2018. The area of disposal dumping site to be modelled is around 9 hectares, with minimum-maximum coordinate 756405E, 617023N-756724E, 617289N in Universal Transverse Mercator (UTM) coordinate system zone 46N. Figure 3 shows the dumping site area and sampling points in blue dot.
To model the site area in 3D, 50 and 100 interpolation points were tested for both x axis and y axis. Thus the total of interpolation points will be 50x50= 2500 points, with x and y interval respectively 6.4 m and 5.3 m. Next for the another one, the total number points are 100x100=10000 with 3.2 and 2.7 m for x and y interval . In figure 4 can be seen the result of IDW interpolation for both scenario. From the figure can be observed that the first case with 2500 interpolation points produced a coarser output compare to the second one, it happened because the interval for first case is bigger compare to another one. 10000 interpolation points Figure 5 shows the 3D surface topographic modelling from both cases. The 3D surface model was generated using minimum 5 closest sampling points with power value (p) 2. The surface function from Plotly graphic library was used to visualize the 3D model. From figure  5 can be seen that the Plotly's surface function produced smooth surface, but figure 5 (b) gives more detail surface texture

Conclusion
IDW spatial interpolation has been successfully estimated a number of unsampled height points that used in topographic surface modelling.
As a deterministic interpolation approach that used distance of estimated location to each neighborhood sampling points in weight determination stage, IDW also proved as a robust interpolation approach which can bridge the drawback of linier interpolation and the complexity of stochastic Kriging interpolation. Using the default Shepard's power value 2, the IDW interpolation method was implemented in Python programming language by creating a function. For each unsampled point to be estimated the function will be called. The output of estimation values are stored in a list and visualized by Plotly's surface function library which produces a smooth surface. More interpolation points will produce more detail surface texture.

Recommendation
The default Shepard's power value 2 was used in this paper in weight determination. For next research it is suggested to use an optimized power value by finding the minimum error between interpolation result and sampling value. More over comparing with other interpolation method like Kriging is also suggested to get more insights how some interpolation methods work in 3D topographic surface modelling.