## Article info

**Volume ** 3

**Pages ** 128 - 139

**DOI **10.5027/jnrd.v3i0.13

**Published **27/12/2013

**Keywords **Detain topography, Digital elevation model, Interpolation, Tillage systems

Under a Creative Commons license

## Authors

**Muhammad Anggri Setiawan ^{a}, Martin Rutzinger ^{b}, Volker Wichmann ^{c}, Johann Stoetter ^{b, c}, Junun Sartohadi ^{a}**

^{a} Department of Environmental Geography, Universitas Gadjah Mada, Indonesia

^{b} Institute of Geography, University of Innsbruck, Austria

^{c} AlpS GmbH – Center for Climate Change Adaptation and Technologies

^{*}Corresponding author: panyidiksiti@gmail.com [stag_icon icon="envelope-o" url="" size="15px" new_window="no"]

## Abstract

There are very little attempts of DEM evaluation in such a disturbed or discontinuous surface (e.g., in tillage area). Present study aims to evaluate common interpolation methods (triangulation, nearest neighbor, natural neighbor, minimum curvature, multiquadratic radial basis function (MRBF), ordinary kriging, and inverse distance weight) in representing the detail topography of two different tillage types, namely bench terrace and furrow. Evaluation procedure was conducted through a stepwise analysis by using combination between the accuracy level (coefficient of determination (R2), mean error (ME) and standard deviation error (S)) and the shape similarity analysis. This study also shows the application of break-line function during the interpolation process in order to optimize some interpolation methods and the usage of drainage sink area as another step in evaluating DEM quality. To achieve the aim of this study, two field-size of dry-land agriculture (tegalan) were observed by using a set of total station Nikon DTM 322 with 3” angle accuracy. These plots, namely Tieng (1652 m²) and Buntu (673 m²), are situated in the upper part of Wonosobo District, Central Java Province, Indonesia. Tieng plot represents the bench terrace system embedded with stones on its terrace risers and showing relatively smooth ground surface. On the other side, Buntu plot shows the ridges and furrows system that lays perpendicularly to the contour lines. In terms of R², ME and S, there were slight differences in results between each method, except the multiquadratic radial basis function which was failed to generate terrace form in Tieng. The final result shows that triangulation is the best fit method followed by natural neighbour at representing the bench terraces in Tieng plot. In the case of furrow in Buntu plot, natural neighbour is the most accurate method. Despite its superiority at representing the bench terrace, triangulation has larger sink drainage area compared to natural neighbour. This study has confirmed the robustness of a stepwise analysis between quantitative and qualitative assessment techniques for DEM accuracy. A fine value of quantitative parameter does not necessarily mean that it will fairly possess a good spatial accuracy.

## Introduction

Representing tillage systems by Digital Elevation Model (DEM) is an important task in environmental modelling, particularly for soil erosion and hydrological processes within agricultural areas. It is common to represent the characteristics of tillage system by an index instead of using the real dimension of it. For instance, Wischmeir and Smith (1978) introduced the tillage system index (P) for the Universal Soil Loss Equation model (USLE). Morgan (2005) used the soil surface roughness (RFR) to represent the effect of different tillage systems on soil erosion process. In a detailed analysis, it is however necessary to maintain the real dimension of the tillage systems in DEMs by using appropriate interpolation algorithms.

Tillage activity changes the surface continuity, which can create specific morphologic structures such as ridges, terraces, and furrows. Some of those characteristics appear as regular patterns, others occur as abrupt changes i.e. discontinuities. Thus, it requires an adequate point sampling technique in order to maintain the surface structure of interest in the interpolated DEM (Aguilar et al., 2005). Surveying such an area could be efficiently done by means of topographic LiDAR technology (Fröhlich & Mettenleiter, 2004; Hack et al., 2004).

Nevertheless, traditional survey techniques acquiring single points such as theodolite or differential GPS are still common in practical work, especially if the acquisition of topographic LiDAR data is not affordable. The limited number of sampling points from such a traditional survey requires a comprehensive investigation of the pros and cons of different interpolation methods in order to be able to produce a reliable and accurate DEM.

Several interpolation methods have been developed and improved in order to provide high quality DEMs. Most of those methods were developed for certain purposes and therefore have their advantages and insufficiencies (Mitas & Mitasova, 1999; Moore et al., 1991; Webster & Oliver, 2007). Some studies have evaluated different interpolation methods by investigating (i) the influence of the sampling patterns of field measurements (Heritage et al., 2009; Merwade, 2009; Zimmerman et al., 1999), (ii) the density of sampling points (Aguilar, et al., 2005; Chaplot et al., 2006; Heritage, et al., 2009), and (iii) the morphological characteristics of the surface (Aguilar, et al., 2005; Chaplot, et al., 2006; Zimmerman, et al., 1999). Most evaluation studies on DEM interpolation are carried out on continuous surfaces. Thus, there is a demand on comprehensive investigations on point interpolation methods applied to surfaces with regular patterns and step edges such as in tillage areas.

Our study aims at finding the best fit interpolation method from the already published methods in order to represent two different tillage types, namely bench terrace and furrow. To improve the interpolation, some methods allow the integration of breaklines during the interpolation process. The resulted DEMs are investigated by a stepwise analysis using a combination between quantitative and shape similarity analysis. In terms of the quantitative analysis, parameters such as mean error (ME), error standard deviation (S), variables of linear regression (i.e. weighted coefficient of determination (wR²) and intercept value (a)), and sink drainage area are used in this study. After the calculation of quantitative parameters, the interpolated DEMs are visualized in 3D views and profile to analyze the shape similarity. Finally, rank classification is performed to decide on the best quality interpolation methods.

## Material and Methods

**Study site**

Two field plot areas in dry-land agriculture (tegalan) were observed. These plots, namely Tieng (1652 m²) and Buntu (673 m²), are situated in the upper part of Serayu Watershed in Wonosobo District, Central Java Province, Indonesia (Fig. 1). These plot areas were initially selected for measuring and modelling the soil erosion rate under both different land use and soil conservation methods. In this area, the limited availability of arable land area has forced the local farmers to modify the hillslope surface into agricultural field. Regardless of the slope steepness, bench terraces and furrows have widely been adopted as the common tillage systems within this area.

The Tieng plot has vertical terrace risers (ca. 2 m) and nearly plane terrace beds (Fig. 1a). An andesitic boulder from the early Holocene period is located on top of this plot. In addition, numbers of rocks are embedded along the terrace risers to reduce the soil erosion effect. The Buntu plot area (Fig. 1b) is identified by a set of ridges and furrows interrupted with a number of ditches, which are perpendicularly oriented to the main slope direction. Buntu plot has a relatively thin soil layer (30-40 cm) nearly approaching the bedrock.

**Data acquisition and pre-processing**

There is always a trade-off between the number of sampling points to be taken and the desired resolution and target scale. In order to carry out an effective survey, this study followed the principal work of (Aguilar, et al., 2005) and (Heritage, et al., 2009) who described the crucial role of the sampling strategy according to the morphological characteristics of the surface. In our study, point sampling was done at every major change within the plot’s area. Demanding a high resolution DEM (10x10 cm for Tieng and 5x5 cm for Buntu) a total station (Nikon DTM 322 with 3” angle precision) was employed.

At Tieng plot, the sampling points were taken along the edges of the terrace risers, on both the upper and lower part (Fig. 1a). At the spot with the andesitic boulder, the sampling points were surveyed right on its surface changes. As a consequence, each interpolation algorithm had to show its capability to interpolate both continuous (the boulder spot) and terraced surfaces (the bench terraces). The terrace bed was considered as plane surface with low surface roughness (2 to 5 cm differences in elevation). Hence, for practical reasons the terrace bed was not particularly observed since roughness did not exceed 10 cm.

Figure 2 shows the point data set of Tieng. The field survey campaign yielded 277 sampling points. For validation purposes, we chose randomly 29 points from the total sampling points. To optimize the interpolation processes, the remaining data (248 points) were then linearly densified along the terrace’s edges (0.5 m interval). The densification process resulted in 2254 points which were used in the interpolation process. In addition, to derive the breakline features, the points observed along the terrace edges (excepting the validation points) were simply connected and converted into blanking file format by using the module of Export Surfer Blanking File in SAGA-GIS (System for Automated Geoscientific Analyses) (Conrad, 2007).

A thorough ground survey was also carried out at the Buntu plot. Visualizing every ridge and furrow as a straight line, the points were sampled both at the starting and the end point of each imagined line. This resulted in two sample points for each ridge and furrow. The point data set of Buntu is shown in Figure 3. In order to validate this sampling strategy some points were also measured in the middle part of the ridges and furrows located on the western part of the plot. In addition, sampling points were also taken on the positions of newly planted trees and along the middle site of the ditches. In total, the field survey yielded 2016 sample points. For validation purposes, we excluded 235 points from the original dataset. From those validation points, 171 points were derived from the middle site of furrows and ridges of the Buntu’s western part and 64 points were taken from the trees’ positions. Applying the same technique used at Tieng, linear densification with a 0.5 m interval was also carried out to the remaining original data set resulting in 6508 points for interpolation.

**Interpolation process**

Interpolation methods were selected in this study based on two considerations. First, is their availability in some readily GIS software packages and, secondly, is according to their complexities algorithm from the simplest principle to the most compound one. The selected interpolation methods for this study were nearest neighbour, triangulation, natural neighbour, inverse distance weighting (IDW), minimum curvature, multiquadratic radial basis function (MRBF) and ordinary kriging.

To process those interpolation routines, GIS software of Surfer Version 8.0 (www.goldensoftware.com) was used in this study due to its complete functions in surface interpolation process. It is possible to incorporate the so-called break-line function, which can be used to integrate discontinuity features in the interpolation process. The SAGA-GIS (Conrad, 2007) was also deployed in proportion of validating and visualizing the DEMs. In the following a brief description to each interpolation method is provided. Some examples from studies using these surface interpolation methods in various fields are mentioned as well.

**Assessment of the DEM accuracy**

DEM accuracies are often confirmed quantitatively with the root mean square error (RMSE). However, RMSE merely provides a standard deviation value of the elevation error calculated from sampled points and higher accuracy data. It assumes two conditions: First, the values of sampled points are in normal distribution and secondly, the mean error is equal to zero which might not always be the case (Fisher & Tate, 2006). The mean error can be vary because of inaccurate interpolation result. Furthermore, a DEM with a good RMSE does not always represent similarity in comparison to the real surface (Declercq, 1996; Yang & Hodler, 2000) due to a limited sampling data set. To overcome this shortcoming, (Fisher & Tate, 2006) recommended to apply the error standard deviation (S) in which the real value of mean error is considered.

Another alternative technique in assessing the DEM quality can be represented through measuring the total area of sink drainages. Such sink drainages can occur due errors in the input data or imperfect interpolation (Wang & Liu, 2006). The size of those sinks can range from single cells up to groups of connected cells which do not have any down-slope path on its surrounding cells (Wang & Liu, 2006). As a consequence drainage paths end in such sinks instead of reaching the main outlet. A critical question, however, always arises on how to distinguish between real and spurious sinks, in particular, when working on overview scales (e.g. catchment area). It may be easier to recognize the spurious sinks as noise on a field-plot scale where sampling points are conducted at every major surface change. The basic consideration is that the more spurious sinks appear, the more spatial errors occur in the DEM.

Other authors (Desmet, 1997; Fisher & Tate, 2006) have agreed that evaluating the DEM accuracy should consider both the quantitative and the shape similarity analysis. It is, in fact, a matter of practical reason whether to merely use the quantitative level or to apply the shape similarity analysis. On the one hand, the quantitative level is easier to compare but it lacks any measure on shape similarity (Yang & Hodler, 2000). On the other hand, shape analysis is more like a descriptive or qualitative assessment but has the ability to describe the similarity between interpolated and reference surface (Wood & Fisher, 1993). It remains a challenging task to combine these two assessment methods in order to get a more reliable measure of the spatial accuracy of DEMs.

As mentioned in the previous section, it is not sufficient to merely use the statistical test for assessing the DEM accuracy. Thus, current study followed the work of (Desmet, 1997; Wood & Fisher, 1993; Yang & Hodler, 2000) that combined the quantitative and qualitative parameters in order to assess the DEM accuracy. The investigation phase of this study is depicted in Figure 4.

Instead of the RMSE, we used ME and S for the statistical test routine, as recommended by (Fisher, 1998; Li, 1994). The value of ME was kept without absolute value ((Fisher & Tate, 2006) to define whether the DEM is underestimate (negative) or overestimate (positive). Those equations are described as follows:

Incorpo

rated with ME and S, further assessment was also carried out through the weighted coefficient of determination (wR²) and intercept value (a) based on the linear regression (Krause et al., 2005). The wR² is obtained from the calculation of R² and the gradient b through the following equation:

In terms of R² and its function, earlier authors (Caruso & Quarta, 1998; Desmet, 1997; Heritage, et al., 2009) have also encompassed them to identify DEM accuracy generated from different sampling strategies and interpolation methods. Finally, the last quantitative assessment parameter used was the total area of sinks drainage (Wang & Liu, 2006). It is considered that a larger area (m²) of sink drainage occurs if the DEM is less accurate.

After all quantitative parameters were identified; visual analyses were then conducted to select the most similar DEM compared to the original shape of the terraces and furrow shape. Two simple visualization techniques were used in namely cross-section profiling and combination of shaded relief maps with 3D views. To produce such a representative profile for the test sites Tieng and Buntu, two perpendicular cross-section lines were drawn in each of the plot area (Figure 2 and Figure 3). For Tieng plot, line A-B represents the shape of the boulder and the terrace risers while line C-D depicts the terrace bed. The line E-F at Buntu plot represents the longitudinal shape of furrow and ditches, whilst the G-H illustrates the regular repetition of the ridges and the furrows.

A DEM, which has relatively extreme quantitative value or dissimilar shape, was deliberately rejected from the final analysis. Rank classification was then performed among the potentially realiable DEMs based on every quantitative parameter. The best DEM on certain quantitative parameter will gain the smallest score. Rank classification has already been addressed by the work of (Chaplot, et al., 2006), which merely focused on RMSE value. To emphasize the reliability effect of each parameter to the DEM accuracy, weighting by factor of 2 was embedded to the parameter of S.

## Result

In order to derive optimum DEM result, crucial parameters for each interpolation method were properly adjusted beforehand, particularly for minimum curvature, IDW, MRBF and ordinary kriging. For minimum curvature, the maximum number of iterations was set between one to two times from the total grid nodes, i.e. 350,000 for Tieng and 400,000 for Buntu. The maximum residual value was set by default 0.022 for Tieng and 0.017 for Buntu yielded from 0.001 x (Zmax - Zmin). In the case of IDW, some combination values of distance power, search radius and maximum numbers were simulated to get the most appropriate value both of Tieng and Buntu. It was confirmed that combination value of 3, 5, and 2 gave the best performance for Tieng plot, while Buntu plot was 5, 1, and 4 for distance power, search radius and maximum numbers, respectively. For the MRBF, we used smoothing factor (R²) of 0.12 for both Tieng and Buntu. In the case of ordinary kriging, linier (Fig. 5a) and gaussian variogram (Fig. 5b) were likely more fit to Tieng’s data trend whilst power variogram (Fig. 6) was fit for Buntu area. In addition, to optimize the interpolation result in Tieng plot, breakline function was set into the minimum curvature, MRBF and ordinary kriging.

Both of Tieng and Buntu’s data set were described in Table 1. The coefficients of variation (CV) for both areas were significantly low due to the regular repetition of the tillage forms. However, Tieng possessed higher CV (0.34%) than Buntu (0.08%) due to the occurrence of terrace riser that has distinct height to the terrace bed.

Following the stepwise analysis in evaluating the DEM accuracy (Figure 4), the quantitative validation data from all of the interpolation results were initially calculated. Table 2 and 3 compile all of those data for both Tieng and Buntu plot. In order to observe the effectiveness of breakline function in Tieng plot, the DEMs resulted by minimum curvature, MRBF and ordinary kriging without breakline function were also included in this quantitative analysis.

Overall, validation result by means of linear regression parameter in the Tieng plot exhibited fine results. The R² values of all DEMs were nearly perfect – close to 1. Likewise, the b values were also fine except the DEM resulted by MRBF-breakline method. It only produced 0.6, which then resulted in low value of wR². Among others, its intercept value was also extremely large and similar to the elevation average provided in Table 1. In contrast, the MRBF showed a better value of linear regression parameter when the breakline function is not included during the interpolation process. Meanwhile, the best value for the linear regression parameter was derived by the triangulation method with 0.99 of wR² and -7.06 of a value.

The next quantitative analysis for Tieng plot was carried out by considering the S and ME values. The lower value of S and ME is gained, the more accurate DEM will be resulted. Focusing on the S value, both triangulation and natural neighbour showed the best value (0.29). However, the ME value of natural neighbour was slightly lower rather than that of triangulation. Although gaining lower accuracy than natural neighbour and triangulation, the S value of ordinary kriging linear variogram with and without breakline function were better compared to the remaining methods. Moreover, among the DEMs resulted by the group of ordinary kriging, implementation of breakline function could provide better S value. The S values between IDW, ordinary kriging gaussian variogram-breakline, and nearest neighbour were slightly similar showing from 0.54, 0.58, and 0.53, respectively. Meanwhile, the ordinary kriging gaussian variogram tended to give an underestimate DEM result (-0.16 of ME) with high S value (0.89). Likewise, the MRBF also showed the most unsatisfactory result of S value (1.09).

Analysing the sink area was the final step for the quantitative evaluation between the DEMs in Tieng plot. Among the DEMs, both small sink area was resulted by natural neighbour (4.08 m²) and triangulation (5.16 m²). Although both DEMs resulted by MRBF showed low value of S, they had smaller total sink area (19.58 m² for MRBF with breakline and 23.26 m² for MRBF without breakline function) compared to nearest neighbour (50.04 m²) and IDW (120.41 m²). Such small sink areas in MRBF could be affected by the smoothing factor during its interpolation process which is neither considered by nearest neighbour nor IDW. The exceptionally large area of sink drainage in IDW and nearest neighbour gave a reason to reject this method at representing terrace feature in Tieng plot. Likewise, the ordinary kriging gaussian variogram had also a potential reason to be rejected due to its large sink drainage (100.5 m²).

After acquiring the result from quantitative analysis in Tieng plot, the qualitative analysis was then conducted by using the 3D view and cross profile. All of the 3D views of DEMs are depicted in Figure 7. In order to create a contrast view over the terrace features, the hill-shading effect was also incorporated into each of DEMs. Those 3D views provided a general overview of each DEM in representing the step form of bench terrace. Meanwhile, Figure 8 illustrates the cross profile of Tieng plot in order to observe the detailed feature of the terrace riser and terrace bed in every DEM.

First observation was given to the DEMs resulted from the variation of ordinary kriging method. As illustrated in Figure 5, the Tieng data set was fit to the linear and gaussian variogram. Based on the 3D view in Figure 7, both variogram obviously failed in representing the step feature of bench terrace. Although the DEM from linear variogram showed better view rather than that of gaussian variogram, it created such undulating surface instead of step form. In contrast, implementation of breakline function could result in better feature of bench terrace in the ordinary kriging with linear variogram but not for the gaussian variogram. For this reason, the ordinary kriging linear variogram without breakline and both DEMs from gaussian variogram were rejected for representing the Tieng plot.

Other obvious failures in representing the bench terrace were shown in the MRBF and minimum curvature method. It was evident that the implementation of breakline function did not give any better result for both methods. All of them produced undulating surface instead of step form. In this sense, the DEMs resulted from both methods were potentially rejected.

Beside the ordinary kriging linear variogram-breakline, there were others method that could depict the step form of the bench terrace in Tieng plot i.e. triangulation, IDW, nearest neighbour, and natural neighbour. Although they posed similar result, the gradation colour through the hill-shading technique showed different effect. The darker colour of hill-shading shows the steeper feature. Between the triangulation and natural neighbour methods had similar pattern showing dark colour on the terrace riser and smooth graduation of grey colour on the terrace bed. For the IDW, scattered black and white spot were found along the edge of terrace riser and on the terrace bed. In the 3D view of nearest neighbour, there were two distinct features embedded in its DEM. First, the stripping shadow on the terrace riser and secondly, the narrow black colour spot along the terrace bed.

Some discrepancies showed in the 3D view were then thoroughly observed by cross profile illustrated in Figure 8. In order to simplify the analysis, the DEMs of MRBF, minimum curvature, ordinary kriging linear variogram and both ordinary kriging gaussian variogram were not included in the cross profile analysis due to their incapability based on the previous step analysis. For a visual comparison, we inserted some observed points into the cross-section rofile.Through the A-B profile, it was confirmed that both minimum curvature-breakline and MRBF-breakline failed to represent the plain shape of terrace bed despite implementing the break-line function. In the case of nearest neighbour, its profile exhibited such stepwise form through the boulder shape, even though it showed better form on the terrace form. Over the cross profile, IDW tended to produce small micro-topography right on the terrace’s edge. There were numbers of artefacts along the line indicating more sink occurrence compared to natural neighbour and triangulation.

The C-D profile depicted a better view of inaccuracy of some DEMs. All of the profiles originated from minimum curvature-breakline, MRBF-breakline and ordinary kriging linear variogram-breakline were overestimating the observed points. Through this profile, the IDW also produced small topography along the plane terrace bed. In contrast, both triangulation and natural neighbour were pretty reliable in representing the terrace bed.

In the case of Buntu plot, the wR² and intercept value were initially compared as the first step of the quantitative analysis. Among the tested interpolation methods have the same value of wR² (0.98) except for the nearest neighbour that only resulted in 0.93. Additionally, the nearest neighbour method exhibited the worst intercept value (-27.11). This was the first indication of inaccurate DEM resulted from nearest neighbour method.

After analysing the value of ME and S in Buntu plot, it was more evident that the nearest neighbour method was less accurate compared to other methods. It exhibited 0.29 of S and 0.05 of ME. That was another reason to reject the nearest neighbour in representing Buntu plot area. In contrast, other methods resulted in 0.1 of S. Between them the IDW was slightly better on the ME value (0.01). These results confirmed that based on the S value the IDW method was superior compared to other methods.

The last quantitative analysis for Buntu plot was the sink drainage area. According to the values mentioned in Table 3, IDW showed the largest area of sinks (46.7 m²) while nearest neighbour resulted in smallest area (1.97 m²). Those values were contrary with earlier mentioned result of S, ME and wR². Thus, this required further analysis through the qualitative method. Meanwhile, between the ordinary kriging, minimum curvature and MRBF produced relatively similar sink area, which were actually larger than that of natural neighbour (13.06 m²) and triangulation method (10.24 m²).

The first impression through the E-F and G-H profiles comparison was shown by the stepwise form yielded by nearest neighbour (Figure 10). Sharp tip were also identified at the ridge and furrow’s edges (figure 10a) which was in fact not found in Buntu plot. Figure 9a enhances those sharp edges depicted as distinct shadow areas right at the upper edge. Those sharp edges were exhibited by the triangulation method. Due to their contrived form, DEMs interpolated by means of nearest neighbour and triangulation were rejected for further analysis in Buntu plot.

Another spurious rough surface was also found in the E-F profile of the IDW method. This rough surface could promote the occurrence of sinks (Table 3). Through the Figure 9f, group of depressions occurred along the furrows. Likewise, number of stripping shadow feature were also found on top of the ridges. Those features were actually the common effect of local extrema resulted by IDW method. Hence, IDW was also rejected to represent the furrow form in Buntu plot.

Based on the E-F cross profile, three methods were likely to exhibit similar form, namely minimum curvature, MRBF and ordinary kriging. Unexpectedly, they exhibited undulating surface at the end of ridge’s edge and lower site of furrow before reaching the ditches area. In reality, this occurrence was not found along the furrow area in Buntu plot. Due to that dissimilarity, those three irregularly rough DEMs were also rejected for further analysis.

Combination of quantitative and qualitative assessment of DEM accuracy have filtered some potential DEM and rejected spurious DEM at representing the terrace and furrow form. The terrace form in Tieng plot was well represented by means of triangulation and natural neighbour. The scoring procedure was then conducted between those DEMs and showing that the triangulation is better compared to natural neighbour (Table 4). In case of Buntu plot, there was only one method left, namely natural neighbour, which was regarded as the best fit DEM at representing the furrow. It fairly means that scoring and rank classification was not carried out in case of Buntu plot.

## Discussion

Two issues should be addressed first before describing the result in detail. First, this study convinces the work of Heritage et al. (2009) that morphological changes can effectively be used as the basic sampling framework during the survey. The point sampling taken along the terrace’s edges (Tieng plot) and at the tip of ridge and furrow (Buntu plot) can be an efficient and effective technique for creating proper DEM data of bench terrace and furrow tillage. Secondly, the point density is of great importance to the DEM result. Without linear densification through the observed points, any tested interpolation method in this study can not produce optimum DEM. In this sense, our study is in accordance to Aguilar et al. (2005) who mentioned the principle key in generating the DEM data, i.e. sampling technique and number of point sampling.

Each of the interpolation method tries to accurately represent the regular repetition of bench terraces and furrow forms. Any interpolation method, which is overdoing the smoothing process, will yield less accurate DEM either in Tieng or Buntu plot. Even for a more sophisticated interpolation method such as ordinary kriging can not provide optimum result to the discontinuous surface of bench terrace. That result corresponds to Mitas and Mitawova (1999) who underlined the incapability of kriging method in representing the local geometry – in this case, the step form of bench terrace. Natural neighbour and triangulation, in contrast, show better performance in representing the Tieng plot despite its simple algorithm.

The break-line function is implemented to the ordinary kriging, minimum curvature, and MRBF to optimize the interpolation process in Tieng plot. This function constraints the grid calculation nearby the edges of terrace risers to have value as close as to the value along the break-line. Among them, ordinary kriging linear variogram exhibit relatively close to the bench terrace form, even though it is not as accurate as both triangulation and natural neighbour. Meanwhile, Minimum curvature and MRBF always attempts to give the smoothest surface to its DEM along the break-line area. As a result, there are irregularity shapes in the minimum curvature and MRBF.

Through this study, we can confirm that assessing the DEM accuracy by mere statistical value (ME, S, and wR²) does not always provide reliable analysis. For instance, in the case of Buntu plot, all of the tested interpolation methods – excluding the nearest neighbour – gain same value of S and wR². In fact, we found different final result after conducting a thorough qualitative technique by using the 3D view and cross profile analysis. Thus, this finding follows the earlier studies by Wood and Fisher (1993), Declercq (1996), Desmet (1997), Yang and Hodler (2000), and Fisher and Tate (2006) that stated the importance of combination between the quantitative and qualitative technique in assessing the DEM accuracy.

Distinguished with other quantitative parameters which depend on the point’s comparison, the sink drainage area is promoted in this study as an indication of spatial error. Indeed, it has strong correlation to the shape similarity. The more artefacts occur at the DEM surface, the larger the sink area will be. It applies to all methods except the nearest neighbour at the furrow system. Despite its stepwise form on furrow, it yields the smallest sink area compared to others. It indicates that a more flat area occurred rather than a sink due to those stepwise forms. Those flat features are also considered as noise within the DEM because they can introduce a discontinuity to the surface flow (Garbrecht & Martz, 1997).

The final result shows that triangulation is the best fit method followed by natural neighbour at representing the bench terraces in Tieng plot. In the case of furrow in Buntu plot, natural neighbour is the most accurate method. Despite its superiority at representing the bench terrace, triangulation has larger sink drainage area compared to natural neighbour. It means that triangulation will require more pre-processing to eliminate such sinks rather than natural neighbour. Eventually, there are found such a combination tillage form between bench terrace and furrow in the study area. In that features, natural neighbour might perform as the best fit method with minimum efforts.

## Conclusions

This study has confirmed the robustness of a stepwise analysis between quantitative and qualitative assessment techniques for DEM accuracy. A fine value of quantitative parameter does not necessarily mean that it will fairly possess a good spatial accuracy. Thus, this study suggests to not independently apply the quantitative parameters without crosschecking it by the shape similarity analysis. For further applications, it is always important to compare the accuracy of the DEM yielded from the readily possible methods instead of depending to one favourite method. In addition, the applicability of breakline function for optimizing such more sophisticated interpolation method such as ordinary kriging requires additional studies. However, with a proper point sampling technique, a simple interpolation can result in more reliable DEM. In this study, we conclude that natural neighbour performs as the best fit method at the furrow and triangulation generates the most accurate DEM at the bench terrace area.

## References

Aguilar, F. J., Agüera, F., Aguilar, M. A., & Carvajal, a. F. (2005). Effects of Terrain Morphology, Sampling Density, and Interpolation Methods on Grid DEM Accuracy. Photogrammetric Engineering & Remote Sensing, 71(7), 805-816.

Caruso, C., & Quarta, F. (1998). Interpolation methods comparison. Computers & Mathematics with Applications, 35(12), 109-126. doi: 10.1016/s0898-1221(98)00101-1

Chaplot, V., Darboux, F., Bourennane, H., Leguédois, S., Silvera, N., & Phachomphon, K. (2006). Accuracy of interpolation techniques for the derivation of digital elevation models in relation to landform types and data density. Geomorphology, 77(1-2), 126-141.

Conrad, O. (2007): SAGA - Entwurf, Funktionsumfang und Anwendung eines Systems für Automatisierte Geowissenschaftliche Analysen. electronic doctoral dissertation, University of Göttingen

Declercq, F. A. N. (1996). Interpolation Methods for Scattered Sample Data: Accuracy, Spatial Patterns, Processing Time. Cartography and Geographic Information Science, 23, 128-144.

Desmet, P. J. J. (1997). Effects of Interpolation Errors on the Analysis of DEMs. Earth Surface Processes and Landforms, 22, 563-580. doi: 10.1002/(sici)1096-9837(199706)22:6<563::aid-esp713>3.0.co;2-3

Fisher, P. (1998). Improved Modeling of Elevation Error with Geostatistics. GeoInformatica, 2(3), 215-233. doi: 10.1023/a:1009717704255

Fisher, P. F., & Tate, N. J. (2006). Causes and consequences of error in digital elevation models. Progress in Physical Geography, 30(4), 467-489. doi: 10.1191/0309133306pp492ra

Fröhlich, C., & Mettenleiter, M. (2004). Terrestrial laser scanning - new perspective in 3D surveying. International archives of photogrammetry, remote sensing and spatial information sciences, XXXVI(8/W2).

Garbrecht, J., & Martz, L. W. (1997). The assignment of drainage direction over flat surfaces in raster digital elevation models. Journal of Hydrology, 193(1-4), 204-213. doi: 10.1016/s0022-1694(96)03138-1

Hack, R., Azzam, R., Charlier, R., & Slob, S. (2004). 3D Terrestrial Laser Scanning as a New Field Measurement and Monitoring Technique Engineering Geology for Infrastructure Planning in Europe (Vol. 104, pp. 179-189): Springer Berlin / Heidelberg.

Heritage, G. L., Milan, D. J., Large, A. R. G., & Fuller, I. C. (2009). Influence of survey strategy and interpolation model on DEM quality. Geomorphology, 112(3-4), 334-344.

Krause, P., Boyle, D. P., & Bäse, F. (2005). Comparison of different efficiency criteria for hydrological model assessment. Advances in Geosciences, 5, 89-97.

Li, Z. (1994). A comparative study of the accuracy of digital terrain models (DTMs) based on various data models. ISPRS Journal of Photogrammetry and Remote Sensing, 49(1), 2-11.

Merwade, V. (2009). Effect of spatial trends on interpolation of river bathymetry. Journal of Hydrology, 371(1-4), 169-181.

Mitas, L., & Mitasova, H. (1999). Spatial interpolation. In P. Longley, K. F. Goodchild, D. J. Maguire & D. W. Rhind (Eds.), Geographical Information Systems: Principles, Techniques, Management and Applications (pp. 481-492). New York: Wiley.

Moore, I. D., Grayson, R. B., & Ladson, A. R. (1991). Digital terrain modelling: A review of hydrological, geomorphological, and biological applications. Hydrological Processes, 5(1), 3-30.

Morgan, R. P. C. (2005). Soil erosion and conservation (3rd ed.). Oxford, UK: Blackwell Publishing.

Wang, & Liu. (2006). An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. International Journal of Geographical Information Science.

Webster, R., & Oliver, M. A. (2007). Geostatistic for Environmental Scientist (2 ed.): John Wiley & Sons, Ltd.

Wischmeier, W. H., & Smith, D. D. (1978). Predicting rainfall erosion losses-a guide to conservation planning US Department Agriculture Handbook (Vol. No. 537, pp. 57). Washington DC: USDA.

Wood, J. D., & Fisher, P. F. (1993). Assessing Interpolation Accuracy in Elevation Models. IEEE Computer Graphics and Applications, 13, 48-56.

Yanalak, M. (2003). Effect of gridding method on DTM profile based on scattered data. Journal of Computing in Civil Engineering, 17(1), 58-67.

Yang, X., & Hodler, T. (2000). Visual and Statistical Comparisons of Surface Modeling Techniques for Point-based Environmental Data. Cartography and Geographic Information Science, 27, 165-176.

Zimmerman, D., Pavlik, C., Ruggles, A., & Armstrong, M. (1999). An Experimental Comparison of Ordinary and Universal Kriging and Inverse Distance Weighting. Mathematical Geology, 31(4), 375-390.