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

Horizontal accuracy assessment of a novel algorithm for approximate a surface to a DEM Barrera, D. a,b, Ibáñez, M.J. a, Eddargani, S. a, c, Romero, R. d, Ariza-López, F.J.e, Reinoso-Gordo, J.F.f,* a Department of Applied Mathematics, University of Granada, Spain – [email protected], [email protected], [email protected] b IMAG -- Institute of Mathematics, Granada, Spain c University Hassan First, MISI Laboratory, Morocco d Department of Computer Science and Artificial Intelligence, University of Granada, Spain – [email protected] e Department of Cartographic, Geodesic and Photogrammetry Engineering – [email protected] f Department of Architectural and Engineering Graphic Expression – [email protected] * [email protected] Abstract: 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. Keywords: Spline, Bernstein basis, control points, Bézier ordinates, tensor product, resample, horizontal accuracy, DEM 1. 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. Farin 2001) of the surface patches relative to the squares that make up the partition. Most studies on DEM accuracy of DEM THT compare a product with a reference have analyzed the altitude Proceedings of the International Cartographic Association, 4, 2021. 30th International Cartographic Conference (ICC 2021), 14–18 December 2021, Florence, Italy. This contribution underwent single-blind peer review based on submitted abstracts. https://doi.org/10.5194/ica-proc-4-11-2021 | © Author(s) 2021. CC BY 4.0 License. 2 of 5 component, leaving the horizontal component untouched. 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. 2. Material and Methodology Our study was carried out on a 2x2 m resolution DEMref produced by the Instituto Geográfico Nacional of Spain for the Navarra region. The following phases have been carried out: 1. Starting from the 2x2 m resolution DEMref, a resampling has been performed using the nearest neighbor method at resolutions of 16x16, 8x8 and 4x4 m (DEMder) which we respectively denote DEM4x4, DEM8x8, DEM16x16. In this way, the DEMder have the same values as the DEMref at the points interpolated by the nearest neighbor method. 2. To evaluate the capacity of our approximation algorithm (Aapx) to estimate the DEMref values, the surface of each DEMref se ha calculado la superficie de cada DEMder (Sder) has been calculated using Aapx. Sder is assessed to obtain DEMs with the same resolution as DEMref which we call DEMremXxX (DEMrem4x4, DEMrem8x8, DEMrem16x16). 3. The horizontal displacement between DEMref and each one of the DEMremXxX is calculated using the Reinoso 2011 contour method. 4. The results obtained with Aapx are compared with the traditional bicubic resampling algorithm. 2.1 Material The DEMref 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 DEMref that contains flat areas along with other steep slopes. The coordinates are referred in the ETRS89 system 30N UTM zone. Figure 1: DEM used as reference and its geographical environment 2.2 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 {𝑣𝑖,𝑣𝑖+ℎ 3,𝑣𝑖+1 −ℎ 3,𝑣𝑖+1}. If 𝐷3 is the union (without repetitions) of these Bézier abscissae and 𝑢𝑖=𝑣𝑖−ℎ 3y 𝑤𝑖= 𝑣𝑖+ℎ 3, then 𝐷3=⋃{𝑢𝑖,𝑣𝑖,𝑤𝑖} 𝑖∈ℤ . For 𝑥 ∈𝐼𝑖, it holds 𝐴𝑓(𝑥)=𝑉𝑖𝐵0(𝑡)+𝑊𝑖𝐵1(𝑡)+𝑈𝑖+1𝐵2(𝑡)+𝑉𝑖+1𝐵3(𝑡), with 𝑡 = 𝑥 ℎ−𝑖 and 𝐵𝑘(𝑡)=(3 𝑘)𝑡𝑘(1−𝑡)3−𝑘, 0≤𝑘 ≤3. The Bézier ordinates 𝑉𝑖,𝑊𝑖,𝑈𝑖+1 and 𝑉𝑖+1 are defined as linear combinations of point values 𝑓(𝑣𝑖): 𝑉𝑖=𝑓(𝑣𝑖), 𝑈𝑖=∑𝛼𝑟𝑓(𝑣𝑖+𝑟) 1 𝑟=−1 , 𝑊𝑖=∑𝛽𝑟𝑓(𝑣𝑖+𝑟) 1 𝑟=−1 , 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. Proceedings of the International Cartographic Association, 4, 2021. 30th International Cartographic Conference (ICC 2021), 14–18 December 2021, Florence, Italy. This contribution underwent single-blind peer review based on submitted abstracts. https://doi.org/10.5194/ica-proc-4-11-2021 | © Author(s) 2021. CC BY 4.0 License. 3 of 5 Proposition 1 The unique masks yielding 𝐶1 continuity and the reproduction of quadratic polynomials are 𝛼=(1 6,1,−1 6) and 𝛽=(−1 6,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 ‖(𝑓−𝐴𝑓)(𝑙)‖∞,𝐼𝑖≤𝐾𝑙ℎ3−𝑙‖𝑓(3)‖∞,𝛺, where 𝛺 =[𝑣𝑖−1,𝑣𝑖+2]. 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 𝐵𝑚(𝑥 ℎ−𝑖)𝐵𝑛(𝑥 ℎ−𝑗),0 ≤𝑚,𝑛 ≤3, whose coefficients are expressed in terms of values 𝑓(𝑣𝑖+𝑘,𝑣𝑗+𝑙),−1≤𝑘,𝑙 ≤1. For Franke function (see Franke 1982) (top), the results provided by the method proposed for ℎ =2−7 are shown in Fig. 2. Figure 2: Results provided by our Aapx for h=2-7 according to the Franke funcion Also for Nielson test function (Nielson 1978) it provides good results for the same step length, shown in Fig.3. Figure 3: Results from our Aapx algorithm according to the Nielson test function 2.3 Contours-based algorithm to measure the horizontal displacement The the horizontal displacement computation of the fases with respect to the DEMref 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 C4a and C4b 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 d). The horizontal displacement (Hdi) computed by the ith pair of homologous contours (Cia, Cib) is formulated as the area enclosed by both curves (Ai) divided by the mean length of those contours (𝐿𝑚𝑖= 𝐿𝑖𝑎+𝐿𝑖𝑏 2) 𝐻𝑑𝑖=𝐴𝑖 𝐿𝑚𝑖 3. And the mean displacement of the DEMremXxX respect to DEMref (𝐻𝑑DEMremXxX) is computed as the weighted average of the displacements of all the homologous contours, the weighting factor being the average length of those contours, divided by the total length of the average contours being the total length 𝐿𝑇𝑜𝑡 =∑𝐿𝑚𝑖 𝑛 𝑖=1 : 𝐻𝑑DEMremXxX =1 𝐿𝑇𝑜𝑡 ∑𝐴𝑖∗ 𝑛 𝑖=1 𝐿𝑚𝑖 Proceedings of the International Cartographic Association, 4, 2021. 30th International Cartographic Conference (ICC 2021), 14–18 December 2021, Florence, Italy. This contribution underwent single-blind peer review based on submitted abstracts. https://doi.org/10.5194/ica-proc-4-11-2021 | © Author(s) 2021. CC BY 4.0 License. 4 of 5 Figure 4: Hologous contours and area between them 3. 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. Figure 5: Area between homologous contours corresponding to the DEMref and the DEMrem16x16 Table 1 shows the mean horizontal displacements as well as their standard deviations calculated for 4x4, 8x8 and 16x16m resolutions using our new Aapx algorithm and the traditional bicubic resampling method. DEMremXxX 𝐻𝑑DEMremXxX Aapx (m) Bicubic (m) DEMrem4x4 Mean 0.07 0.99 Std 0.03 0.06 DEM rem8x8 Mean 0.28 1.05 Std 0.13 0.06 DEM rem16x16 Mean 0.81 1.38 Std 0.38 0.27 Figure 6: Horizontal displacement from the 𝐻𝑑𝐷𝐸𝑀𝑟𝑒𝑚𝑋𝑥𝑋 respect to the DEMref, comparing our algorithm and the traditional bicubic. Table 1 shows that our new algorithm produces better results than the traditional bicubic algorithm regardless of the cell size used as DEMder. 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 Proceedings of the International Cartographic Association, 4, 2021. 30th International Cartographic Conference (ICC 2021), 14–18 December 2021, Florence, Italy. This contribution underwent single-blind peer review based on submitted abstracts. https://doi.org/10.5194/ica-proc-4-11-2021 | © Author(s) 2021. CC BY 4.0 License. 5 of 5 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. 4. Conclusions In this work a new algorithm (Aapx) is presented to approximate a DEM by means of a piecewise defined surface (Sder). 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 Sder to obtain DEMs (DEMremXxX) of higher or lower resolution than those used to create Sder. One consequence of Aapx resampling capabilities is being able to assess the accuracy of a product DEM (DEMpro) against a higher accuracy DEMref. This assessment could be carried out both in the vertical component and in the horizontal component, which is the one studied in this work. DEMpro can come from both a source or a method other than the source or method used to create the DEMref. Aapx 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 DEMpro against a DEMref, and whenever possible it will be necessary to choose those algorithms that produce the least error (horizontal displacement). The Aapx 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 Aapx so that the results that appear in this first Aapx study can be verified. 5. References Farin, G. (2001). Curves and Surfaces for CAGD: A practical guide, 5th Edition. Morgan Kaufmann, San Francisco, USA. Franke, R. (1982) Scattered data interpolation: tests of some methods, Math. Comp. 157, pp. 181-200. Gao, J. (1998). Impact of sampling intervals on the reliability of topographic variables mapped from grid DEMs at a micro-scale. Int. J. Geogr. Inf. Sci., 12, 875– 890. Keys, R. (1981). "Cubic convolution interpolation for digital image processing". IEEE Transactions on Acoustics, Speech, and Signal Processing. 29 (6): 1153– 1160 Leong Tan, M.; Ficklin, D.L.; Dixon, B.; Ibrahim, A.L.; Yusop, Z.; Chaplot, V. (2015). Impacts of DEM resolution, source, and resampling technique on SWATsimulated streamflow. Appl. Geogr., 63, 357–368 Maune D (2007) Digital elevation model technologies and applications: the DEM user’s manual, 2nd edn.Bethesda, USA Mesa-Mingorance, J.L.; Ariza-López, F.J. (2020) Accuracy Assessment of Digital Elevation Models (DEMs): A Critical Review of Practices of the Past Three Decades. Remote Sens. , 12, 2630. Mukherjee, S.; Joshi, P.K.; Mukherjee, S.; Ghosh, A.; Garg, R.D.; Mukhopadhyay, A (2013). IU accuracy of open source Digital Elevation Model (DEM). Int. J. Appl. Earth Obs. Geoinf., 21, 205–217. Nielson G. M. (1978) A first order blending method for triangles based upon cubic interpolation, Int. J. Numer. Meth. Engrg. 15, pp. 308-318. Reinoso, J. F. (2010) A priori horizontal displacement (HD) estimation of hydrological features when versioned DEMs are used, Journal of Hydrology,384: 1–2, pp. 130141. Reinoso, J. F. (2011). An algorithm for automatically computing the horizontal shift between homologous contours from DTMs, ISPRS Journal of Photogrammetry and Remote Sensing, 66: 3, pp. 272-286. Shen, J., Tan, F. (2020). Effects of DEM resolution and resampling technique on building treatment for urban inundation modeling: a case study for the 2016 flooding of the HUST campus in Wuhan. Nat Hazards 104, 927– 957. Wang, B.; Shi, W.; Liu, E. (2015). Robust methods for assessing the accuracy of linear interpolated DEM. Int. J. Appl. Earth Obs. Geoinf., 34, 198–206. Proceedings of the International Cartographic Association, 4, 2021. 30th International Cartographic Conference (ICC 2021), 14–18 December 2021, Florence, Italy. This contribution underwent single-blind peer review based on submitted abstracts. https://doi.org/10.5194/ica-proc-4-11-2021 | © Author(s) 2021. CC BY 4.0 License.