Horizontal accuracy assessment of a novel algorithm for approximate a surface to a DEM

This study evaluates the horizontal positional accuracy of a new algorithm that defines a surface that approximates DEM data by means of a spline function. This algorithm allows evaluating the surface at any point in its definition domain and allows analytically estimating other parameters of interest, such as slopes, orientations, etc. To evaluate the accuracy achieved with the algorithm, we use a reference DEM 2 m × 2 m (DEMref) from which the derived DEMs are obtained at 4 m × 4 m, 8 m × 8 m and 16 m × 16 m (DEMder). For each DEMder its spline approximant is calculated, which is evaluated at the same points occupied by the DEMref cells, getting a resampled DEM 2x2m (DEMrem). The horizontal accuracy is obtained by computing the area amongs the homologous contour lines derived from DEMref and DEMrem, respectively. It has been observed that the planimetric errors of the proposed algorithm are very small, even in flat areas, where you could expect major differences. Therefore, this algorithm could be used when an evaluation of the horizontal positional accuracy of a DEM product at lower resolution (DEMpro) and a different producing source than the higher resolution DEMref is wanted.


Introduction
Having a mathematical function that represents the terrain throughout its continuous definition domain has different advantages, among others the following: a) It is possible to sample regular meshes to generate digital elevation models (DEM) of both higher and lower resolution than the data from which the mathematical function was obtained; and this is achieved thanks to its definition domain is continuous, b) morphological variables of interest can be obtained from the corresponding mathematical formulas of the surfaces, such as slope, orientation, curvature and normal direction, c) You could intersect two surfaces corresponding to homologous DEMs from different dates and calculate the increase or decrease in the terrain volume. The provision of functions of this type has allowed resampling through bilinear (Maune, 2007) and bicubic (Keys, 1981) interpolations that have been used in different applications both to obtain DEMs of higher and lower resolution. In the first case, obtaining higher resolution has been used to, for example, improve urban flood zones in the absence of denser models (Shen and Tan, 2020); in the second case, its use has been frequent when it was intended to compare the altimetric accuracy of a lower resolution product model with another higher resolution reference model (Gao, 1998, Mukherjee et al. 2013, Wang et al. 2015, although, in most cases, the error introduced by the resampling from a higher resolution to a lower one was left unanalysed, as indicated by Mesa and Ariza, 2020. Other studies have addressed the influence of resampling techniques on products derived from DEMs such as streamflows. (Leong et al., 2015). Although procedures are available to extract information directly from a DEM, we propose a new algorithm to build a surface with low computational cost that adjusts the elevations in order to have an explicit expression (function) from which to find the elements of interest. As the terrain has many irregularities, a surface should not be constructed too regular. C1 continuity is sufficient. The regular structure of the DEM allows to define a piecewise surface defined on a quadrangular partition of the terrain. More precisely, we define a piecewise bicubic surface by providing simple rules that give the Bézier ordinates (cf. The reason for the lower number of investigations devoted to the horizontal component is due to the difficulty of finding a satisfactory method. In this work we will study the horizontal accuracy achieved by the new algorithm that we have proposed. We will use the automatic algorithm based on homologous contour lines introduced in Reinoso, 2010 and rigorously demonstrated in Reinoso, 2011 for evaluating the horizontal accuracy.

Material and Methodology
Our study was carried out on a 2x2 m resolution DEM ref produced by the Instituto Geográfico Nacional of Spain for the Navarra region. The following phases have been carried out:

Material
The DEM ref has a cell size of 2x2m and y occupies an area of 4.8x4.8 Km. Figure 1 shows the geographical characteristics of the environment, as well as the DEM ref that contains flat areas along with other steep slopes. The coordinates are referred in the ETRS89 system 30N UTM zone.

The approximation algorithm
We propose to construct a spline surface by means of a tensor product of 1D spline approximants, defining the surface patches directly in the Bernstein basis. Suppose that for a real function the values ( ), ∈ ℤ, are known, where = ℎ, with ℎ > 0 the size of the partition ≔ { : ∈ ℤ}. The 1D approximating spline reduces on each interval ≔ [ , +1 ] to a cubic polynomial, whose control polygon is formed by four control points with Bézier abscissae { , + where the masks = ( −1 , 0 , 1 ) ∈ ℝ 3 and = ( −1 , 0 , 1 ) ∈ ℝ 3 are determined to achieve 1 continuity as well as the reproduction of the quadratic polynomials. , 1, 1 6 ). Furthermore, the uniform norm of the corresponding operator is equal to 4/3.
From the exactness of , the following result regarding the approximation error holds.
Proposition 2 There exist constans , = 0,1, independent of and ℎ and , such that Now, given a 2D function ( , ) a bi-cubic piecewise surface is defined as a tensor product approximant: the operator is applied to as a function depending on (or ), and then is again applied to de resulting function, i.e. ( , ) = ( , ).
On each square × this function is a bi-cubic Bézier surface, so that it can be represented in Bernstein-Bézier form. It is a linear combination of functions  Also for Nielson test function (Nielson 1978) it provides good results for the same step length, shown in Fig.3.

Contours-based algorithm to measure the horizontal displacement
The the horizontal displacement computation of the fases with respect to the DEM ref is carried out in the following phases: 1. The contours of both DEMs are calculated (Fig. 4  a and b respectively), and their homologous curves are automatically identified, e.g. curves C4 a and C4 b in Fig. 4 a and b. 2. After superimposing the homologous contours ( Fig. 4 c), the areas enclosed between them are calculated (gray area in Fig. 4

Results and discussion
To calculate the DEMremXxX a 10 m interval between contour lines has been used, that in our DEM produces a total of 22 levels, specifically their heights ranging from 450 to 660 m. In Fig. 5 you can see the homologous contour drawn on a shadow map, as well as a detail where the area enclosed between those homologous contours are highlighted on green color.    Table 1 shows that our new algorithm produces better results than the traditional bicubic algorithm regardless of the cell size used as DEM der . While our algorithm seems to decrease the error due to horizontal displacement at a rate of ¼ as the resolution increases at a rate of 2, in the bicubic algorithm the rate of decrease is much lower. However, no large differences are observed in the values of the standard deviations if both algorithms are compared for each level of resolution. However, additional tests should be carried out with a greater number of DEMs, covering all types of terrain (flat, undulating and mountainous), in order to statistically validate the apparently better results or our algorithm respect the traditional bicubic one. Another advantage of the new algorithm with respect to the traditional bicubic one is that it has a lower computational cost, making it a candidate to be implemented in cartographic production software packages. We believe that this new algorithm can be used when you want to know the positional accuracy (horizontal and vertical) of a lower resolution DEM coming from a source other than the reference one or that has been created with a different method from the reference one. On the other hand, it would also be interesting to have an algorithm that not only reported the horizontal displacement with a scalar value, but also included information about direction in each of the cells, which would be possible by adapting the contours method by Reinoso 2011.

Conclusions
In this work a new algorithm (A apx ) is presented to approximate a DEM by means of a piecewise defined surface (S der ). It presents some advantages linked to its definition type, such as being able to obtain the altitude of a point in the entire definition domain of that surface, as well as morphological variables that characterize the terrain surface: slope, orientation, curvature or normal direction in an analytical way. An immediate application would be the possibility of resampling S der to obtain DEMs (DEM remXxX ) of higher or lower resolution than those used to create S der . One consequence of A apx resampling capabilities is being able to assess the accuracy of a product DEM (DEMpro) against a higher accuracy DEM ref . This assessment could be carried out both in the vertical component and in the horizontal component, which is the one studied in this work. DEM pro can come from both a source or a method other than the source or method used to create the DEM ref .
A apx has shown a lower horizontal displacement than the traditional bicubic interpolation algorithm, which can be interpreted as a lower error when resampling DEMs of lower resolution to others of higher resolution; These processes are necessary when trying to compare the accuracies of a DEM pro against a DEM ref , and whenever possible it will be necessary to choose those algorithms that produce the least error (horizontal displacement). The A apx computational cost is lower than other conceptually similar such as the traditional bicubic one. Finally, in the future an experiment will have to be designed with a sufficiently large number of DEMs on which to test A apx so that the results that appear in this first A apx study can be verified.