Computational modelling of terrains

Computational modelling of terrains Hugo Ledoux Ken Arroyo Ohori Ravi Peters Maarten Pronk v2025.0 ©2025 Hugo Ledoux, Ken Arroyo Ohori, Ravi Peters, and Maarten Pronk c b This work is available under a Creative Commons Attribution 4.0 International Licence. For licence details, see http://creativecommons.org/licenses/by/4.0/ v2025.0 [2025-10-27] v2024.0 [2024-10-29] v2023.0 [2023-11-12] v0.9 [2022-11-14] v0.8 [2021-11-08] v0.7 [2020-11-09] v0.6 [2019-11-11] Download latest version The latest version of this book can be downloaded in PDF at https://github.com/tudelft3d/terrainbook/releases Extra material Most chapters have a short YouTube video explaining the key concepts, and some chapters have extra material. Available at https://tudelft3d.github.io/terrainbook/videos Source code The source code of the book, in L A TEX, is available at https://github.com/tudelft3d/terrainbook Errors? Feedback? Please report errors or potential improvements at https://github.com/tudelft3d/terrainbook/issues Colophon This book was typeset with L A T E X using the kaobook class. The figures were created using Ipe, OmniGraffle, Affinity Designer, or Blender. The front cover images was generated by DALL·E 3, and we made some custom colour changes. Preface This open book presents an overview of algorithms and methodologies to reconstruct terrains, to manipulate them, and to extract information from them. It covers different representations of terrains (eg TINs, rasters, point clouds, contour lines), discusses different applications (visibility analysis, runoff modelling, etc.), presents techniques to handle large datasets, and discusses related topics such as global elevation models and bathymetric datasets. The book presents the theory and gives examples of libraries and software to perform certain tasks, but no code examples are provided. We wanted the book to be language-agnostic. The course website for which this book was developed ( https://3d.bk.tudelft.nl/courses/geo1015 ) provides assignments where students use Python and/or C++. Open material. This book was primarily developed for the course Digital terrain modelling (GEO1015) in the MSc Geomatics programme at Delft University of Technology in the Netherlands. The course is tailored for MSc students who have already followed introductory courses in GIS, programming, and acquisition of geographical datasets. Each chapter corresponds to a lesson in the course, whose content is also openly available: https://3d.bk.tudelft.nl/courses/geo1015 Accompanying videos. Most of the chapters have a short video explaining the key concepts, and those are freely available online: https://tudelft3d.github.io/terrainbook/videos Who is this book for? The book is written for MSc Geomatics students, but it can also be used at the BSc level. Prerequisites include: knowledge of GIS, background in linear algebra, and an introductory programming course. Acknowledgements. We thank the many students of the GEO1015 course over the years who have helped us by pointing—and often fixing with a pull request—the errors, typos, and weird sentences of this book. A special thanks to the students of the year 2018–2019 who had to deal with the first version of this book, and to Balázs Dukai for thoroughly proofreading the first drafts of this book. Contents 1 What is a terrain? 1 1.1 Dimensionality of DTMs ...................................... 2 1.2 2.5D terrain == field ......................................... 3 1.3 Data models for terrains ....................................... 4 1.4 TIN versus raster for modelling terrains .............................. 8 1.5 Notes and comments ........................................ 9 1.6 Exercises ............................................... 9 2 Acquisition of elevation measurements 11 2.1 Principles of lidar .......................................... 11 2.2 Other acquisition techniques .................................... 15 2.3 Artefacts ............................................... 17 2.4 Notes and comments ........................................ 22 2.5 Exercises ............................................... 22 3 Global digital elevation models 23 3.1 Acquisition of global elevation data ................................ 24 3.2 Most common products available ................................. 26 3.3 Specific characteristics of gDEMs .................................. 28 3.4 Notes and comments ........................................ 31 3.5 Exercises ............................................... 32 4 Delaunay triangulations & Voronoi diagrams 33 4.1 Voronoi diagram ........................................... 33 4.2 Delaunay triangulation ....................................... 34 4.3 Duality between the DT and the VD ................................ 37 4.4 Incremental construction of the DT ................................ 38 4.5 Data structures for storing a DT .................................. 42 4.6 Constrained and Conforming Delaunay Triangulations ..................... 44 4.7 Notes and comments ........................................ 46 4.8 Exercises ............................................... 46 5 Spatial interpolation: deterministic methods 49 5.1 What is a good interpolation method for terrains? ........................ 49 5.2 Fitting polynomials ......................................... 50 5.3 Weighted-average methods ..................................... 51 5.4 Assessing the interpolation results ................................. 58 5.5 Overview of all methods ...................................... 59 5.6 Notes and comments ........................................ 61 5.7 Exercises ............................................... 62 6 Spatial interpolation: kriging 63 6.1 Statistical background ........................................ 63 6.2 Geostatistics and the standard geostatistical model ....................... 65 6.3 Covariance, dissimilarity and the semivariogram ........................ 66 6.4 Simple kriging ............................................ 70 6.5 Ordinary kriging ........................................... 72 6.6 Other types of kriging ........................................ 73 6.7 Implementation details ....................................... 74 6.8 Notes and comments ........................................ 74 6.9 Exercises ............................................... 75 7 Conversions between terrain representations 77 7.1 Conversion of PC/TIN to raster .................................. 77 7.2 Conversion to isolines ........................................ 78 7.3 Simplification of a TIN ....................................... 81 7.4 Wedding cake effect ......................................... 84 7.5 Notes and comments ........................................ 85 7.6 Exercises ............................................... 85 8 Topographic properties and features 87 8.1 Topographic properties ....................................... 87 8.2 Topographic features ........................................ 92 8.3 Properties and features used in practice .............................. 94 8.4 Notes and comments ........................................ 96 8.5 Exercises ............................................... 96 9 Visibility queries on terrains 99 9.1 Rendering + ray casting ....................................... 99 9.2 For 2.5D terrains, the problem is simpler ............................. 100 9.3 Notes and comments ........................................ 102 9.4 Exercises ............................................... 103 10 Runoff modelling 105 10.1 Computing the flow direction ................................... 106 10.2 Computing the flow accumulation ................................. 107 10.3 Solving issues with sinks ...................................... 108 10.4 Assigning flow direction in flats .................................. 109 10.5 Drainage networks and basins ................................... 110 10.6 Notes and comments ........................................ 111 10.7 Exercises ............................................... 111 11 Processing point clouds 113 11.1 Thinning ............................................... 113 11.2 Outlier detection ........................................... 115 11.3 Ground filtering ........................................... 116 11.4 Shape detection ........................................... 122 11.5 Notes and comments ........................................ 127 11.6 Exercises ............................................... 127 12 Handling and processing massive terrains 129 12.1 Raster pyramids ........................................... 129 12.2 Indexing points in 3D space with the kd-tree ........................... 132 12.3 Streaming paradigm to construct massive TINs/grids ...................... 135 12.4 Notes and comments ........................................ 140 12.5 Exercises ............................................... 140 13 Spatial extent of a set of points 141 13.1 Properties of the region ....................................... 142 13.2 Convex hull .............................................. 143 13.3 Moving arm ............................................. 143 13.4 𝜒-shape ................................................ 145 13.5 𝛼-shape ................................................ 145 13.6 Clustering algorithms ........................................ 146 13.7 Notes and comments ........................................ 147 13.8 Exercises ............................................... 147 14 Processing bathymetric data to produce hydrographic charts 149 14.1 How are depth-contours produced in practice? .......................... 149 14.2 Common methods used in practice are not satisfactory ..................... 151 14.3 A Voronoi-based surface approach ................................. 153 14.4 Some examples of results with real-world datasets ........................ 157 14.5 Notes and comments ........................................ 159 14.6 Exercises ............................................... 159 Appendices 161 A Point cloud file formats 163 A.1 ASCII formats ............................................ 163 A.2 PLY format .............................................. 164 A.3 LAS format .............................................. 165 A.4 LAZ format .............................................. 167 B Extra information about the AHN datasets 169 C Estimating the normals in a point cloud 171 D Some useful equations 173 D.1 Centre of a circle defined by 3 points ............................... 173 Bibliography 175 Index 181 6 1 What is a terrain? TIN: triangulated irregular network Figure 1.4: A TIN is obtained by lifting the vertices to their elevation. All the triangles are usually Delaunay, ie their circumcircle(green) is emptyofany other points in the plane. 1.3.3 Irregular tessellations The cells of an irregular tessellation can be of any shape and size, and they usually ‘follow’—or are constrained by—the samples points that were collected, albeit this is not a requirement. Subdividing the space based on the samples has the main advantage of producing a tessellation that is adaptive to the distribution of the samples. The subdivision is potentially better than that obtained with regular tessellations (which subdivide arbitrarily the space without any considerations for the samples). The most known examples of the use of irregular tessellations in terrain modelling is the triangulated irregular network, or TIN. As shown in Figure 1.4, a TIN refers to an irregular tessellation of the 𝑥𝑦 -plane into non-overlapping triangles (whose vertices are formed by three sample points), and to the use of a linear interpolation function for each triangle. One way to explain the 2.5D properties of a TIN is as follows: if we project vertically to the 𝑥𝑦 -plane the triangles in 3D space forming the TIN, then no two triangles will intersect. While not a requirement, the triangulation is usually a Delaunay triangulation (more about this in Chapter 4). The main reason is that Delaunay triangles are as “fat” as possible (long and skinny triangles are avoided), and thus they behave better for interpolation. As can be seen in Figure 1.5, the estimated value can be significantly different, and in this case the right one would make more sense since sample points that are closer to the interpolation location are used (in the TIN on the left, the value of 95m is not used). Every point (which becomes a vertex in the triangulation) is lifted to its elevation to create a surface, embedded in three dimensions, approximating the morphology of the terrain. The value of elevation at an unsampled location 𝑝 is obtained by linearly interpolating on the plane passing through the three vertices of the triangle containing 𝑝 . TINs are the most popular alternatives to 2D grids for modelling elevation; both representations have advantages and disadvantages. A TIN in which a linear interpolation function is used yields a 𝐶0 piecewise representation, ie it is a continuous function but at the edges of the triangles the first derivative is not possible. It is possible to use higher-order functions in each triangle of a TIN, to construct a 𝐶1 or 𝐶2 field, ie where the first and second derivative of the surface can be obtained. Chapter 5gives more details about spatial interpolation and continuity. Figure 1.5: Two TINs (left is nonDelaunay; right is Delaunay) and the result of estimating with linear interpolation in the TIN. estimated height = 38m estimated height = 99m 110m 29m 95m 41m 110m 95m 41m 29m 1.3 Data models for terrains 7 quadtree contour lines raster TIN point cloud Figure 1.6: Four most common data models for terrains. 1.3.4 Hierarchical tessellations Hierarchical tessellations attempt to reduce the number of cells in a tessellation by merging the neighbouring cells having the same value (thus yielding cells of different sizes). While both regular and irregular tessellations can be hierarchical, in the context of the representation of terrains, the former is more relevant and is sometimes used in practice. A commonly used hierarchical structure in two dimensions is the quadtree, which is a generic term for a family of tessellations that recursively subdivide the plane into four quadrants. As is the case for grids, quadtrees are relatively easily implemented in a computer because they are trees in which each node has exactly four children, if any. The shortcomings of regular hierarchical tessellations are similar to those of regular tessellations: the rotation and scaling operations are difficult to handle. The main advantage of using them—saving memory space—is present only when there is spatial coherence between cells having the same attribute value, ie when they are clustered together. Indeed, the size of a quadtree is not dependent on the number of cells, but on their distribution. The quadtree of a 2D grid having no two adjacent cells with the same value (eg a checkers board) contains the same number of cells as the grid, and its size would most likely be worse because of the overhead to manage the tree. Another disadvantage is that the notion of neighbours, which is straightforward in regular tessellations, is less trivial. 1.3.5 Other common terrain representations used in GIS In the GIS literature, besides the ones above, different representations for terrains are often listed, the two most relevant being: 1. irregularly spaced sample points, such a point cloud; 2. contour lines. It should be noticed that these two are however incomplete: the set of rules to reconstruct the surface at unsampled locations is not explicitly given, they are not continuous surfaces. Conceptually speaking, these should therefore not be considered valid representations of a terrain. While this might seems odd, this is in line with the consensus among practitioners today, where a point cloud or contour lines would typically be used as an input to a process to generate a terrain. In Chapter 11 we present and discuss several algorithms and techniques to process raw point clouds, so that the points can be used to construct terrains. We will nevertheless consider these in the course; the four representations we will use are shown in Figure 1.6. Contour lines. Given a bivariate field 𝑓(𝑥, 𝑦)=𝑧 , an isoline (commonly named contour line) is the set of points in space where 𝑓(𝑥, 𝑦)=𝑧0 , where 𝑧0 is a constant. Isolines have been traditionally used to represent the elevation in topographic maps and the depth in bathymetric maps for navigation at sea. One particular property of an isoline is that its direction is always perpendicular to the direction of the steepest slope of the terrain. Another property that follows from the 2.5𝐷 property of the field is that contours neither intersect themselves nor each other. 8 1 What is a terrain? Figure 1.7: An example of a terrain (top) and its contour lines at every 100m (bottom). 100m 200m 300m 400m 500m 600m 600 500 400 600300 200 200 300 400 100 Figure 1.8: Cross-sectionofaterrain(left), and the 200m isoline extracted from a TIN representation of it (right). 200m The purpose of isolines on a map is to reveal the shape of the underlying terrain. By observing the shape and interrelation of neighbouring contours, the presence and significance of surface features becomes apparent; see Figure 1.7 for an example. It should be noticed that data between contours is absent in the contour map. Yet, in case of good contours the reader will still be able to deduct the general morphology of the field. It is even so that the use of contours will speed up the map reading process, as it conveys just that relevant bit of data to the map reader rather than ‘flooding’ the reader with information which essentially makes the user do his own cartographic selection. Contouring is a form of discretizing the field that makes it easier to use a map. Naturally, this comes at a price. The level of approximation of the field can (dramatically) differ between contours, the biggest error would be midway in between contour lines. But, depending on the relation between the spacing between contours (the contour interval) and the map scale, which in turn is dependent on the map application, this effect may be neglected. In practice, isolines are only approximated from the computer representation of a field. They are usually extracted directly from a TIN or a regular grid. As shown in Figure 1.8, the idea is to compute the intersection between the level value (eg 200m ) and the terrain, represented for instance with a TIN. Each triangle is scanned and segment lines are extracted to form an approximation of an isoline. Chapter 7gives more details. 1.4 TIN versus raster for modelling terrains There is an ongoing debate about whether TINs or rasters are the better data model to model terrains. Practitioners and scientists are probably split 50/50 on the issue. A data model will be judged more efficient than another if it represents a surface more accurately within the same amount of storage space, 1.5 Notes and comments 9 measured in bytes. This of course depends on the data structure used to store that data model. It should be said that both TIN and raster have advantages and disadvantages (as we will see during this course), and in practice one should choose the most appropriate model for the task at hand. This means converting between the two data models when it is necessary (topic of Chapter 7). 1.5 Notes and comments Kumler (1994) carried out a 4-year comparison between TINs and rasters. He states that the common belief that a TIN is more space-efficient than raster is handicapped by the fact that a TIN must have at least 3 times less points to be of equal space. His conclusions are also that rasters can estimate elevation more accurately than comparably-sized TINs. However, he still finishes with by stating: “Yeah, well. . . TINs still look better.” Fisher (1997) discusses the disadvantages of rasters, in a GIS and remote sensing context. Frank (1992) and Goodchild (1992) discuss at length the issue of data model, data structure and representation of reality. Tse and Gold (2004) coined the term ‘2.75D GIS’ and show an example of a where a triangulation is used to represent the surface of the Earth, with holes (for tunnels), cliffs and caves. The same idea is also referred to as a ‘2.8D GIS’ by Gröger and Plümer (2005). While more an academic exercise then something used in practice, multi-resolution triangulation have been described and used for some application by De Floriani and Magillo (2002). Akima (1978) shows the advantages of using higher-order functions in each region of a TIN, to construct a 𝐶1or 𝐶2field. Dakowicz and Gold (2003) demonstrate that using simple rules (nearestneighbour for instance) yields fields that are not realistic and have bad slope, which is in practice problematic for several applications. Obtaining good slope from contour lines is possible, but is in practice a complex process. 1.6 Exercises 1. Explain in your own words why a point cloud (eg collected with airborne lidar) is not considered a complete representation of a terrain. 2. What is a bivariate function? 3. Assume you have a 2.75D terrain of an area. Is it possible to extract the isolines from it? What properties will these have? Will they intersect? While ‘lidar’ is often treated as the acronym of light detection and ranging, it actually originated as a portmanteau of ‘light’ and ‘radar’. (from Wikipedia) Backscattering is the natural phenomenon of the reflection of (electromagnetic) waves or signals back to the direction they came from. Figure 2.1: Lidar range measurement Acquisition of elevation measurements 2 uyoutu.be/_BSuNu3Ahw0 2.1 Principles of lidar . . . . . . 11 2.2 Other techniques . . . . . . 15 2.3 Artefacts . . . . . . . . . . . 17 2.4 Notes and comments . . . . 22 2.5 Exercises . . . . . . . . . . . 22 The very first step in the process of terrain modelling is the acquisition of elevation measurements. Nowadays, these measurements are usually collected in large quantities using some form of remote sensing, ie sensors that measure—in our case—the distance to the Earth’s surface from an airborne or even a spaceborne platform. In raw form, elevation measurements are typically stored as a point cloud, ie a collection of georeferenced 3D points with each point representing one elevation measurement on the Earth’s surface. There are a number of remote sensing techniques that are used to measure elevation on Earth or other planets. Typically, these techniques measure 1. the distance to the target surface; 2. their own position and orientation with respect to some global reference system. By combining these, we can compute the 3D coordinates of the measured location on the target surface. In this chapter we will focus primarily on lidar, the most common acquisition technique for large scale terrain models with centimetre level accuracy. But we also give an overview of other acquisition techniques, for example photogrammetry, InSAR, and sonar. And to conclude we will look at typical artefacts that you might encounter while working with elevation data. This is because, as with any kind of real-world measurements, there are various uncertainties and restrictions in the acquisition process that lead to distortions—the artefacts—in the acquired data. These artefacts need to be taken into account when further processing and using the elevation data. 2.1 Principles of lidar A lidar system measures the distance to a target by illuminating it with pulsed laser light and measuring the reflected or backscattered signal with a sensor (see Figure 2.1). By measuring the time-of-flight, ie the difference in time between emitting a pulse and detecting its return or echo, the distance to the target that reflected the pulse can be found using a simple formula. To be exact, the time-of-flight 𝑇is equal to 𝑇=2𝑅 𝑐(2.1) where 𝑐 is the speed of light (approximately 300,000 km/s), and 𝑅 is the distance or range between the lidar scanner and the target object that reflects the laser pulse. Therefore the range 𝑅 can be found from the measured time-of-flight 𝑇using 𝑅=1 2𝑇𝑐. A typical lidar systems performs hundreds of thousands of such range measurements per second. 12 2 Acquisition of elevation measurements Figure 2.2: An airborne lidar system. Figure from Dowman (2004). NASA has used space lidar on Earth,on the Moon, and on Mars. inertial navigation system (INS) differential GPS inertial measurement unit (IMU) Lidar scanners exist in various forms. They can be mounted on a static tripod (terrestrial lidar) for detailed local scans, or on a moving platform such as a car (mobile lidar) or an aircraft (airborne lidar) for rapid scanning of larger areas. Nowadays, also hand-held lidar systems exist, and even some of the latest smartphones have a lidar sensor. Furthermore, lidar can also be used from a satellite in space. However, in the remainder of this text we will focus on airborne lidar. 2.1.1 Georeferencing the range measurements Apart from the laser scanner itself, a lidar system uses a GPS receiver and an inertial navigation system (INS), see Figure 2.2. These devices, which respectively provide the global position and orientation of the laser scanner, are needed for georeferencing, ie to convert the range measurements of the laser scanner to 3D point measurements in a global coordinate system such as WGS84. To obtain an accurate global position, differential GPS (DGPS) is employed. DGPS is a technique to enhance the accuracy of GPS by using GPS stations on the ground (one is visible in Figure 2.2). These DGPS stations have a known position and they broadcast the difference between that known position and the position at the station as indicated by GPS. This difference is essentially a correction for errors in the GPS signal. The aircraft receives these differences from nearby DGPS stations and uses them to correct the GPS position of the aircraft. Using DGPS the accuracy of the GPS position on the aircraft can be improved from around 15 meters to several centimetres. To obtain the accurate orientation of the laser scanner, the INS of the aircraft is used. The INS accurately measures the orientation, ie the yaw, pitch and roll angles of the aircraft, by means of an inertial measurement unit (IMU). Only when we accurately know the orientation of the laser scanner, can we know the direction (in a global coordinate system) in which a laser pulse is emitted from the aircraft. By combining the global position and the global orientation of the laser scanner with the range measurement from the laser scanner, the 2.1 Principles of lidar 13 Emited pulse t t echo threshold Detected ehoesBackscattered waveform P P (a) (b) Figure 2.3: The emitted laser pulse, (a) the returned signal, and (b) the recorded echoes. Figure adapted from Bailly et al. (n.d.). return georeferenced 3D position of the point on the target object that reflected the lase pulse can be computed. 2.1.2 Echo detection A lidar system performs ranging measurements using the time-of-flight principle that allows us to compute range from a time measurement using the known speed of light in the air. The time measurement starts when the laser pulse is emitted and is completed when a backscattered echo of that signal is detected. In practice one emitted pulse can even lead to multiple echoes in the case when an object reflects part of the laser pulse, but also allows part of the pulse to continue past the object. Notice that lidar pulses are typically emitted in a slightly divergent manner. As a result the footprint of the pules at ground level is several centimetres in diameter, which increases the likelihood of multiple echoes. Figure 2.3 illustrates what the backscattered signal looks like when it hits a target object in the shape of a tree. A tree is particularly interesting because it often causes multiple echoes (one or more on its branches and one on the ground below). The lidar sensor observes a waveform that represents the received signal power ( 𝑃 ) as a function of time ( 𝑡 ). With the direct detection lidar systems that we focus on in this book, the echoes are derived from the backscattered waveform by using a thresholding technique. This essentially means that an echo is recorded whenever the power of the waveform exceeds a fixed threshold (see Figure 2.3b). An echo can also be referred to as a return. For each return the return count is recorded, eg the first return is the first echo received from an emitted laser pules and the last return is the last received echo (see Figure 2.3). The return count can in some cases be used to determine if an echo was reflected on vegetation or ground (ground should then be the last return). 2.1.3 Anatomy of a lidar system A lidar system consists of an optical and an electronic part. As shown in Figure 2.4, each part consists of several components. 14 2 Acquisition of elevation measurements Figure 2.4: Conventional architecture of a direct detection lidar system. Figure from Chazette et al. (2016). atmospheric scattering In the optical part, a pulse of a particular wavelength (typically nearinfrared) is generated by the laser source for each lidar measurement. It then passes through a set of optics (lenses and mirrors) so that it leaves the scanner in an appropriate direction. After the pulse interacts with the scattering medium, it is reflected back into the scanning optics which then directs the signal into a telescope. The telescope converges the signal through a field diaphragm (essentially a tiny hole around the point of convergence). The field diaphragm blocks stray light rays (eg sunlight reflected into the optics from any angle) from proceeding in the optical pipeline. Next, the light signal is recollimated so that it again consists only of parallel light rays. The final step of the optical part is the inference filter which blocks all wavelengths except for the wavelength of the laser source. This is again needed to block stray light rays from distorting the measurement. The electronic part consists of a photodetector, which first transforms the light signal into an electrical current, which is then converted to a digital signal using the analogue-to-digital converter. Once the digital signal is available, further electronics can be used to interpret and record the signal. 2.1.4 Laser wavelength Choosing the optimal laser wavelength is a compromise of several different factors. One needs to consider atmospheric scattering, ie how much of the signal is lost simply by travelling through the atmosphere, and the absorption capacity of vegetation, ie how much of the signal is lost because it is absorbed by vegetation. In addition, there is the stray signal due to direct and scattered contributions of sunlight. While it is possible to filter such stray signals in the lidar system to some degree, it remains wise to choose a wavelength that is only minimally affected by it. Finally there are regulations that limit the laser radiance values permissible to the eye. This means that the power of emitted signal needs to be carefully controlled, and/or a wavelength must be chosen that is not absorbed by the eye so much. As a result, most lidar systems use a wavelength in the near-infrared spectrum, usually between 600 and 1000 nm. A notable exception is made for bathymetric purposes, in which case a green (532 nm) laser is used because that has a greater penetration ability in water. 2.1.5 Scanning patterns In order to improve the capacity to quickly scan large areas, a number of rotating optical elements are typically present in a lidar system. Using 2.2 Other acquisition techniques 15 Rotating mirror or prism Oscillating mirror Rotating polygonal mirror 2axis galvanometer NON-UNIFORM UNIFORM Figure 2.5: Different configurations of rotating mirrors and the associated scanning patterns from a moving platform. Arrows indicate the direction of the emitted laser signal. Figure from Chazette et al. (2016). To read or to watch This YouTube video explains the principles of an aerial lidar system: https://youtu.be/EYbhNSUnIdU Figure 2.6: Photogrammetry these optical elements, ie mirrors or prisms, the emitted laser pulse is guided in a cross-track direction (ie perpendicular to the along-track direction in which the aircraft moves, see Figure 2.1), thereby greatly increasing the scanned ground area per travelled meter of the aircraft. Figure 2.5 depicts a number of possible configurations of rotating optics and shows the resulting scanning patterns. It is clear that density of points on the ground is affected by the scanning pattern. The top example for example, yields much higher densities on edges of the scanned area. In practice more uniform patterns, such as the bottom two examples are often preferred. 2.2 Other acquisition techniques Apart from lidar there are also other sensor techniques that can be used to acquire elevation data. Some of these are active sensors just like lidar (a signal is generated and emitted from the sensor), whereas others are passive (using the sun as light source). And like lidar, these sensors themselves only do range measurements, and need additional hardware such as a GPS receiver and an IMU to georeference the measurements. What follows is a brief description of the three other important acquisition techniques used in practice. 2.2.1 Photogrammetry Photogrammetry allows us to measure the distance from overlapping photographs taken from different positions. If a ground point, called afeature, is identifiable in two or more images, its 3D coordinates can be computed in two steps. First, a viewing ray for that feature must be reconstructed for each image. A viewing ray can be defined as the line from the feature, passing through the projective centre of the camera, to the corresponding pixel in the image sensor (see Figure 2.6). Second, considering that we know the orientation and position of the camera, 22 2 Acquisition of elevation measurements To read or to watch This is a paper that compares lidar and photogrammetry derived point clouds for the generation of a DEM. It shows that even when artefacts seem to be under control, both techniques may measure different elevations C. Ressl et al. (2016). Dense image matching vs. airborne laser scanning — Comparison of two methods for deriving terrain models. Photogrammetrie - Fernerkundung - Geoinformation 2016.2, pp. 57–73 PDF: https://3d.bk.tudelft.nl/courses/geo1015/data/others /Ressl16.pdf 2.4 Notes and comments If you would like to learn more about how a lidar scanner works, the chapter from Chazette et al. (2016) is recommended. More details on InSAR can be found in the manual from Ferretti et al. (2007). Reuter et al. (2009) give an elaborate overview of the processing that needs to be done to derive a high quality (raster) DTM from raw elevation measurements. 2.5 Exercises 1. Name three differences between point cloud acquisition with lidar and with photogrammetry. 2. Explain what the time-of-flight principle entails. 3. How can you minimise occlusion effects in a point cloud during acquisition? 4. Why does positioning, using for instance GPS, play such an important role in acquisition? https://opentopography.org Global digital elevation models 3 3.1 Acquisition of global data 24 3.2 Most common products . . 26 3.3 Specific characteristics . . . 28 3.4 Notes and comments . . . . 31 3.5 Exercises . . . . . . . . . . . 32 We define as “global digital elevation models” (or global terrains; we refer to them as “gDEMs” in this book) the datasets that cover (most of) the Earth. Those datasets require different acquisition methods from local datasets, since flying an airplane or performing local surveys at the scale of the Earth has not been done yet, and it would be prohibitively expensive. The acquisition instruments for a global coverage must be space-borne, ie mounted on a satellite for instance. Notice that the orbit of some satellites makes them technically non-global, but that we still refer to their datasets as global since they have a wide coverage, just restricted to certain latitudes. It cannot be understated that before the introduction of gDEMs (with SRTM v1 in 2000, see below for more information), we had no way of knowing the elevation for a given location on the Earth. Indeed, looking at the datasets listed and/or hosted at OpenTopography (Figure 3.1), we can observe that, even in 2022, local elevation datasets are mostly limited to developed countries. Global DEMs enable us to perform global environmental studies, such as geological studies, hydrological modelling, ecosystems dynamics, the understanding of volcanic processes, and flow simulations (see Chapter 10). While gDEMs are elevation models like local ones (and can be modelled with essentially the same formats and tools as local ones), they have several properties and characteristics that apply only to them, and we report in this chapter on the main ones. Figure 3.1: OpenTopography coverage with some European datasets added in pink; there are in fact more European datasets but there is no global registry for them. 24 3 Global digital elevation models https://icesat-2.gsfc.nasa.gov/ https://gedi.umd.edu/ 3.1 Acquisition of global elevation data The acquisition of gDEMs requires the use of a sensor mounted on a satellite. The three most used instruments to collect elevation information are: 1. Photogrammetry from optical satellite images 2. Interferometric synthetic-aperture radar (InSAR) 3. Lidar 3.1.1 Photogrammetry from high-resolution satellite images See Section 2.2.1. 3.1.2 InSAR The first gDEM (SRTM v1) was collected with InSAR in 2000, aboard the Space Shuttle Endeavour. As InSAR requires two images of the same area for stereoscopy to derive a DEM, the Shuttle was equipped with two antennas, one in the payload bay, and one on a mast protruding 60m from the Shuttle. Similarly, the TerraSAR-X satellite was joined by TanDEM-X in 2010, a twin satellite orbiting only a few hundred meters (!) away, to generate InSAR data in a single pass. See Section 2.2.2 for more information. 3.1.3 Spaceborne lidar (ICESat-2 + GEDI) Spaceborne lidar is a relatively new technique, and has not yet been used to produce a gDEM. We still include it here because it enables terrain elevation measurements globally (vegetation can be filtered out), and we believe it will help us create gDEMs in the near future. Lidar was first used in space on the Apollo missions, and with further technological developments, it has been used extensively from the 1990s onwards. For example, Mercury, Mars, near-Earth asteroids, and lately again the Moon have been scanned using lidar. Earth surface elevation lidar measurements have also been developed, often flown on the Space Shuttles. ICESat was the first Earth-based lidar satellite, launched in 2003, with the primary goal of ice sheet monitoring. It had an elevation accuracy of several cm and was operational for five years. Most recently, NASA launched in 2018 two missions measure the elevation of the Earth globally with lidar instruments: ▶ ICESat-2 (Ice, Cloud, and Land Elevation Satellite-2) is in a low Earth and polar orbit to investigate ice sheets, it covers the Earth between −88◦ and 88◦ latitude. Its instrument to measure altimetry is called Advanced Topographic Laser Altimeter System (ATLAS). Apart from terrain retrieval, ICESat-2 also measures the surface, such as canopy height, and has many other applications such as measuring bathymetry and estimating biomass. 3.1 Acquisition of global elevation data 25 ICESat-2(92°inclination) GEDI(51.6°inclination) Figure 3.2: Ground tracks for three successive orbits of ICESat-2 and GEDI. The satellite is a represented by a triangle and past orbits fade out. Note the increased density of ground tracks at the latitude of inclination, as well as the lack of coverage beyond 51.6◦latitude for GEDI. Table 3.1: Key characteristics of GEDI and ICESat-2 missions in comparison with a typical airborne lidar mission. ICESat-2 GEDI airborne lidar type discrete photon full waveform either main objective cryosphere monitoring ecosystems — duration 2018–2024(ongoing) 2019–2023, 2024-(ongoing) single flight(s) orbit inclination 92◦51.6◦NA laser pulse power 120𝜇J/30𝜇J 15 000𝜇J/4500 𝜇JNA altitude ∼480km ∼420 km 0.5 km beam footprint 11m 23 m 0.05 m along track spacing 0.7m 70 m 0.1 m across track spacing 3 km/90 m between pair 0.6km 0.1m swath width 6.6km 4.2km 1 km beam frequency 532nm (green) 1064 nm (near-infrared) either # laser(s) 1 3 1 # beams 6 4 1 # ground tracks 6 (in 3 strong/weak pairs) 8 (4 strong, 4 weak) 1 ▶ GEDI (Global Ecosystem Dynamics Investigation) is attached to the international space station (ISS) and its primary goal is to investigate global ecosystems. It does not have global coverage since it collects measurements only between 51.6◦ N and 51.6◦ S. GEDI has been combined with TanDEM-X data to produce biomass estimates and with Landsat imagery to produce a global canopy height map. The characteristics of both missions are summarised in Table 3.1 ICESat-2 laser split into six beams, divided into three pairs, each pair 90m apart and the pairs 3.3km apart, for a total swath width of 6.6km . Along-track, it can measure each 0.7m , while its beam footprint is ∼11m , so each measurement overlaps. GEDI instrument has three lasers, forming 4 beams and eight tracks, each 600m apart, for a total swath width of 4.2km . GEDI measures a point every 70m along-track, with a beam footprint of 23m . Furthermore, whereas GEDI employs a full-waveform laser at a typical near-infrared wavelength of 1064nm , ICESat-2 employs a single-photon lidar at a bathymetric “green” wavelength of 532nm. Product levels. The data from the ICESat-2 and GEDI missions is made publicly available in several data products, categorised in 3 levels (Level 1, 26 3 Global digital elevation models Figure 3.3: Filtered ICESat-2 and GEDI points from a single granule each at the 47th latitude, demonstrating the beam patterns. NotethatICESat-2hasasmaller beam footprint and a much higher pulse repetition, but a more uneven spatial coverage than GEDI. The gaps between data here will decrease by using multiple granules, but will never disappear completely. ICESat-2ATL03 strongbeam weakbeam GEDIL2A strongbeam weakbeam SRTM 2, and 3 data products), where a higher-level is derived from a lower-level product. 1. Level 1 products contains the raw telemetry; 2. Level 2 products contain directly usable geolocated data to which several corrections—such as accounting for atmospheric effects— are applied. 3. Level 3 data are aggregated versions of Level 2 products, which are smaller in filesize and easier to process. ICESat-2 differentiates between a Level 3A, which are aggregated Level 2 data products per granule, and a Level 3B, which are gridded versions of the aggregated Level 3A data products. GEDI’s Level 3 data product are gridded versions of Level 2 data products, like ICESat-2’s Level 3B. GEDI also has Level 4 data products, which are model outputs—like carbon estimates—based on Level 2 data. Comparison to typical airbone lidar. These space borne lasers also differ considerably from airborne lasers, most notably so in their platform, resulting in significant differences in beam footprint and ground coverage. The altitude increase results in a wider beam footprint, from ∼0.5m (at 500m ) for airborne platforms to ∼15m for space platforms. Although much wider, it is a small increase compared to the increase in altitude, going from 0.5km to 500km . A comparison is given in Table 3.1. Airborne lidar often focuses on maximizing coverage ( points/m2 ) of smaller areas, whereas the coverage for space lasers is the ground track of the satellite. While both ICESat-2 and GEDI employ instruments with multiple (split) laser beams, including the ability to point the laser away from the ground track, all to maximize coverage, this still results in very sparse and uneven coverage as shown in Figure 3.3. 3.2 Most common products available An overview of the most common gDEMS are given in Figure 3.4. Note that all these products differ considerably in terms of coverage, resolution, accuracy and licensing. Even the same product can have different versions, with different resolutions and licenses. For example, SRTM is freely available, including its derived NASADEM, and was first introduced at 90m , with subsequent versions at 30m . Tandem-X, and its derived WorldDEM, is a commercial product, with a 3.2 Most common products available 27 SRTM ASTER ALOS NASADEM CopernicusDEM FABDEMMERIT WorldDEM Tandem-X 2000 2010 2020 radar optical derived DTM* not-open Tandem-X 2020 BIOMASS FathomDEM GEDTM30 future Figure 3.4: An overview of current gDEMs Table 3.2: Overview of global DEMS, see Figure 3.5 for their lineage. Year released By Sensor Type License Resolution SRTM 2001 NASA InSAR DSM Open 30–90m ASTER 2009 NASA optical DSM Open 30m Tandem-X 2014 DLR InSAR DSM Closed 12m WorldDEM 2014 Airbus InSAR DSM/DTM′Closed 5–12m ALOS 2016 JAXA optical DSM Open 30m MERIT 2017 Yamazaki et al. (2017) InSAR DTM′Open 90m NASADEM 2019 NASA InSAR DSM Open 30m CopernicusDEM 2020 ESA InSAR DSM Open 30–90m Tandem-X 2020 2022 Wessel et al. (2022) InSAR DSM Closed 30m FABDEM 2022 Hawker et al. (2022) InSAR DTM′Closed 30m FathomDEM 2025 Uhe et al. (2025) InSAR DTM′Closed 30m GEDTM30 2025 Ho et al. (2025) InSAR DTM′Open 30m resolution of ±12m . WorldDEM-NEO, a newer version of WorldDEM with more Tandem-X data (as used in Tandem-X 2020), even has a resolution of ±5m . CopernicusDEM is a resampled WorldDEM—bought with your taxpayer money by ESA and freely distributed—at 30m resolution. The pseudo DTMs FABDEM (Forest And Building removed) and FathomDEM, while derived from the freely available CopernicusDEM, are only free for research purposes. Similarly, while ALOS World3D is freely available at 30m , it also comes in a commercial version at 5m resolution. There is even a 0.5m commercial version, based on multiple optical satellites, available on request. 28 3 Global digital elevation models (a) NASADEM (b) CopernicusDEM Figure 3.5: NASADEM and CopernicusDEM for the Indus delta in Pakistan. Note the striped noise in NASADEM, and how CopernicusDEM has more detail. There is ∼twelve years between these images. GeoTIFF 3.3 Specific characteristics of gDEMs 3.3.1 Global often means near-global The data collected depend on the orbit of the satellite. Some satellites have a polar orbit and can therefore completely measure the Earth (ICESat-2 is one example, see 3.2), while some will cover only certain latitudes (eg GEDI, see 3.2). Similarly, while the Tandem-X mission produced a DEM with a global coverage, the SRTM DEM was measured from the Space Shuttle and only covers up to 60◦latitude. 3.3.2 Format The formats used to store and exchange DEMs have historically been defined by the military, which are still used by many government agencies. One example is the Digital Terrain Elevation Data (DTED) format, developed in the 1970s, which stores elevation in integers (which tells us a lot about the accuracy and precision possible 50 years ago). It specifies several possible levels in terms of resolution (in arcseconds), from level 0 at ±1km to level 2 at 30m . More recently, in 2016, the Defence Gridded Elevation Data (DGED) has been defined, specifying more levels to higher resolutions and allowing GeoTIFFs to used (which removes the integer constraint). DGED also defines the structure and the specific tiling of the data at higher latitudes, as resolutions in arcseconds become smaller near the poles. Most of the gDEMS are tiled in a similar way. CopernicusDEM—adhering tothe DGED level 3 standard—hastilesof3601 × 3601pixelson the equator, but 3600 × 2400 (height, width) pixels at 50◦ latitude, 3600 × 1800 at 60◦ 3.3 Specific characteristics of gDEMs 29 2Downloading gDEMs Most of the gDEMS are available in a standard and easily accessible format, such as GeoTIFF. However, for broad compatibility, recent advances such as new compression techniques and tiling strategies are not yet widely used. One such advance is Cloud Optimized GeoTIFF (COG: https://www.cogeo.org ), which is a normal GeoTIFF with a specific structure that allows it to be read piecewise from the cloud. Without such a structure, a GeoTIFF has to be downloaded in its entirety—even if one is only interested in a small part of it—before it can be read. latitude, 3600 × 1200 at 70◦ latitude, becoming as small as 3600 × 360 for the last 5 degrees of latitude. This tiling scheme results in pixels being as square as possible, but makes it hard to work with tiles from different latitudes. In this context it becomes clear that the resolution should not be discussed in terms of meters, but terms of degrees (or divisions of a degree). As the Earth has a circumference of ±40000 km (measured around the Equator), 1◦ of latitude is ±111km and 1◦ of longitude is 111cos 𝜙km at latitude 𝜙 . Degrees are further divided into 60 arcminutes (‘), which themselves are divided into 60 arcseconds (“). In practice, the highest resolution for SRTM ( 30m ) is actually 1′′ , and its 90m product has a resolution of 3′′ . DTED level 0 thus has a resolution of 30′′ , while level 2 has a resolution of 1′′. Datasets such as SRTM and NASADEM can be provided as .hgt (height) files, which are not even a format, but are simply a binary file with the elevation values listed in a given order, with the geographic extent to be derived from the filename. 3.3.3 Accuracy TheaccuracyofgDEMsisoftenbrokendownintoseveralcomponentsand related metrics, which are not always well-defined. The DTED and DGED specifications differentiate between horizontal and vertical accuracy, and within each defines both relative and absolute accuracy. Relative accuracy describes the consistency of the measurements, specified as the random error component of the uncertainty between two DEM pixels. Absolute accuracy describes the total error of a measurement compared to a reference. The DGED standard specifies a relative vertical accuracy of less than 12m for level 2 (resolution of 1′′ , or ±30m ) and an absolute accuracy (goal) of 18m . CopernicusDEM reports a mean error of less than 2m for 65% of its tiles, and another 19% with an error of less than 5m in terms of absolute accuracy. 3.3.4 Errors As with any measurements taken, gDEMs contain errors and outliers. For example, SRTM contains a lot of noise, and has (diagonal) striping artefacts, visible in Figure 3.5a for the derived NASADEM. 30 3 Global digital elevation models CopernicusDEM[m+MSL] Figure 3.6: Low outliers in CopernicusDEM, with orthophoto on the right for context. Electricity poles, visible by their shadows, are the cause for these errors here. voids in DEMs Earth Gravitational Model (EGM) mean sea level (MSL) Another example is CopernicusDEM that suffers from multipath errors in urban areas, which lead to small pits in the DEM, see Figure 3.6. Most gDEMS suffer from voids in steep terrain, as peaks can occlude valleys below (called the shadow effect). These voids are often filled not only by interpolation, but with the help of other gDEMS, as they measured the same area from a different angle. Furthermore, gDEM products are often accompanied by a quality layer, which indicates where data has been void-filled. Similarly, error masks with the calculated instrument error and water masks are often provided alongside the elevation data. 3.3.5 Vertical datums Global DEMs are vertically referenced to a specific geoid, specifically EGM96 (EPSG:5171) for SRTM and EGM2008 (EPSG:3855) for CopernicusDEM. Be aware that there is a small difference—generally less than 0.5m—between these versions of the EGM geoid. Note that the geoid is not the same as the mean sea level (MSL), as it does not take into account the dynamic effects of temperature and currents. Depending on the location, the difference between the geoid and MSL can be up to 1.5m. 3.4 Notes and comments 31 DTM′=pseudo DTM 3.3.6 gDEMs are DSM (more than DTM) In contrast to local DEMs, which are often provided as either a classified point cloud, or as separate raster DSMs and DTMs, gDEMs should be classified as DSMs. The current measurement techniques will measure the top of canopy and buildings, and not the ground below. This is the largest source of error in gDEMS, and can considerably limit the applicability of the data. Several attempts have been made to correct gDEMS for vegetation and buildings, leading to what we denote as a pseudo DTM (DTM ′ ). While these methods improve the accuracy of gDEMs considerably, they are not perfect, and resulting terrain still contain vegetation and/or buildings. Recent work has also suggested that while these datasets have improved vertical accuracy, the accuracy of derived geomorphometric parameters such as slope and curvature suffers. 3.4 Notes and comments Space lidar is promising as a source to reconstruct gDEMs because lidar penetrates vegetation, and thus obtaining a DTM is an easier process. Recently, several global coastal pseudo DTMs—covering only areas near or below sea-level—have been produced using space lidar (to some extent): CoastalDEM (Kulp and Strauss, 2024), DiluviumDEM (Dusseau et al., 2023) and DeltaDTM (Pronk et al., 2024). Likewise, the new ESA Biomass mission (Quegan et al., 2019)—launched in 2025—will use P-band radar to penetrate vegetation and measure the ground below, potentially enabling a truly global DTM in the future. Yang et al. (2011) provide a detailed list of applications where gDEMs (SRTM, but when the paper was written (2011) SRTM was still the main product available globally) are necessary as input. Schumann and Bates (2018) make a case for the need for high-accuracy open-access DEMs, demonstrating SRTM is not good enough for many applications. Hancock et al. (2021) investigates the requirements for a global lidar DEM. Furtherreadingabout DEM terminologycan be foundinGuthetal. (2021), which is one of the products of the Digital Elevation Model Intercomparison eXperiment (DEMIX) group. They also published a paper comparing the vertical accuracy and derived parameters of several gDEMs (Guth et al., 2024). Arguably the best place to download DEMs (gDEMS, local ones, lidar datasets, etc.) is OpenTopography ( https://opentopography.org ). Otherwise each gDEM has its own download portal, with its own registration systems and its specific ways of searching and downloading the data. In case of the most recent gDEMS, such as CopernicusDEM, the data must be downloaded as .tar (archives) for 1◦×1◦ tiles via FTP, in folders for each continent and then country. Data is duplicated for the border areas, totalling 2TB. To learn more about the DGED format, read https://dgiwg.org/do cuments/dgiwg-standards . Similarly, reading a user guide on any of the gDEMS is a good idea, such as the one for CopernicusDEM Product Handbook. 38 4 Delaunay triangulations & Voronoi diagrams p Figure 4.12: (top) The DT before and (bottom) after a point 𝑝 hasbeen inserted. Noticethat the DTis updatedonly locally (only the yellow triangles are affected). 4.4 Incremental construction of the DT Since the VD and the DT are dual structures, the knowledge of one implies the knowledge of the other one. In other words, if one has only one structure, she can always extract the other one. Because it is easier, from an algorithmic and data structure point of view, to manage triangles over arbitrary polygons (they have a constant number of vertices and neighbours), constructing and manipulating a VD by working only on its dual structure is simpler and usually preferred. When the VD is needed, it is extracted from the DT. This has the additional advantage of speeding up algorithms because when the VD is used directly intermediate Voronoi vertices—that will not necessarily exist in the final diagram—need to be computed and stored. While there exists different strategies to construct at DT, we focus in this book on the incremental method since it is easier to understand and implement. An incremental algorithm is one where the structure is built incrementally; in our case this means that each point is inserted one at a time in a valid DT and the triangulation is updated, with respect to the Delaunay criterion (empty circumcircle), after each insertion. Observe that the insertion of a single point 𝑝 in a DT modifies only locally the DT, ie only the triangles whose circumcircle contains 𝑝 need to be deleted and replaced by new ones respecting the Delaunay criterion (see Figure 4.12 for an example). In sharp contrast to this, other strategies to construct a DT (eg divideand-conquer and plane sweep algorithms, see Section 4.7), build a DT in one operation (this is a batch operation), and if another point needs to be inserted after this, the whole construction operation must be done again from scratch. That hinders their use for some applications where new data coming from a sensor would have to be added, or where we want to delete points because they are outliers. The incremental insertion algorithm, and the other well-known algorithms, can all construct the DT of 𝑛 points randomly distributed in the Euclidean plane in O(𝑛log 𝑛). Figure 4.13 illustrates the steps of the algorithm, and Algorithm 1its pseudo-code. In a nutshell, for the insertion of a new point 𝑝 in a DT( 𝑆 ), the triangle 𝜏 containing 𝑝 is identified and then split into three new triangles by joining 𝑝 to every vertex of 𝜏 . Second, each new triangle is tested—according to the Delaunay criterion—against its opposite neighbour (with respect to 𝑝 ); if it is not a Delaunay triangle then the Algorithm 1: Algorithm to insert one point in a DT 1Input: A DT(𝑆)T, and a new point 𝑝to insert Output: T 𝑝=T∪{𝑝}// the DT with point 𝑝 2find triangle 𝜏containing 𝑝 3insert 𝑝in 𝜏by splitting it in to 3 new triangles (flip13) 4push 3 new triangles on a stack 5while stack is non-empty do 6𝜏={𝑝, 𝑎, 𝑏} ←pop from stack 7𝜏𝑎={𝑎, 𝑏, 𝑐} ←get adjacent triangle of 𝜏having the edge 𝑎𝑏 8if 𝑐is inside circumcircle of 𝜏then 9flip22 𝜏and 𝜏𝑎 10 push 2 new triangles on stack 4.4 Incremental construction of the DT 39 Figure 4.13: Step-by-step insertion, with flips, of a single point in a DT in two dimensions. o1 o2o3 DT(S) Figure 4.14: The set 𝑆 of points is contained by a big triangle formed by the vertices 𝑜1 , 𝑜2 and 𝑜3 . Many triangles outside conv(𝑆) are created. ∞ Figure 4.15: The infinite vertex ( ∞ ) is used to ensure that the triangles in DT( 𝑆 ) are always adjacent to exactly 3 triangles. This DT contains 7 finite triangles and 5 infinite triangles. edge shared by the two triangles is flipped (a flip is an operation to modify adjacent triangles, see below) and the two new triangles will also have to be tested later. This process stops when every triangle having 𝑝 as one of its vertices respects the Delaunay criterion. 4.4.1 Initialisation: the big triangle or the infinite vertex The DT of a set 𝑆 of points subdivides conv( 𝑆 ), which means in practice that the triangles on the boundary of conv( 𝑆 ) will not be adjacent to exactly 3 neighbouring triangles. Because it is convenient to store and manipulate triangles having exactly 3 neighbours, in practice most DT construction algorithms will use one of these two “tricks”: ▶ Big triangle: 𝑆 is entirely contained in a big triangle 𝜏𝑏𝑖𝑔 several times larger than the spatial extent of 𝑆 ; conv( 𝑆 ) therefore becomes 𝜏𝑏𝑖𝑔 . Figure 4.14 illustrates this. The construction of DT( 𝑆 ) is for example always initialised by first constructing 𝜏𝑏𝑖𝑔 , and then the points in 𝑆are inserted one by one. ▶ Infinite vertex: a fictitious vertex is inserted at the “infinity”, and therefore the edges on the boundary of conv( 𝑆 ) are incident to “infinite triangles” formed by a convex hull edge and the infinite vertex, see Figure 4.15. This can be conceptually seen as embedding 𝑆 on a sphere, and adding the infinite vertex on the other side of the sphere. The infinite vertex is conceptually the same as the big triangle but is numerically more stable since the size of the big triangle does not need to be defined. Observe however that since the infinite vertex has no coordinates, the predicates Orientation and InCircle used to construct a DT (see Section 4.4.5) cannot be used with the infinite vertex and infinite triangles, instead one should handle those with specific cases. Using a big triangle or an infinite vertex has many practical advantages. First, since an edge is always guaranteed to be shared by two triangles, point location algorithms never “fall off” the convex hull. Second, when a single point 𝑝 needs to be inserted in DT( 𝑆 ), this guarantees that 40 4 Delaunay triangulations & Voronoi diagrams 2How are DT created in practice? Several implementations of the DT use a big triangle or the infinite vertex, CGAL ( https://www.cgal.org/ ) and startinpy ( https: //github.com/hugoledoux/startinpy ) are two examples. Those will refer in their API to “finite” and “infinite” vertices, edges, and triangles. It is therefore essential to understand the mechanism to use those librairies, even if one is not constructing the DT herself. Figure 4.16: The Walk algorithm for a DT in two dimensions. The query point is 𝑝 . starting triangle p 𝑝 is always inside an existing triangle; we thus do not have to deal explicitly with vertices added outside the convex hull. Third, identifying the vertices that bounds conv( 𝑆 ) is easy: they have one incident triangle that has one or more of the big triangle vertices (or it contains the infinite vertex). Fourth, the Voronoi cells of the points that bounds conv( 𝑆 ) will be bounded, since the only unbounded cells will be the ones of the 3 points of 𝜏𝑏𝑖𝑔 . This can help for some of the spatial analysis operations, for instance interpolation based on the VD (see Chapter 5). The main disadvantage is that more triangles than needed are constructed. For example in Figure 4.14 only the shaded triangles would be part of DT( 𝑆 ). The extra triangles can nevertheless be easily marked as they are the only ones containing at least one of the 3 points forming 𝜏𝑏𝑖𝑔. 4.4.2 Point location with walking To find the triangle containing the newly inserted point 𝑝 , we can use the point-in-polygon test for every triangle (the standard GIS operation), but that brute-force operation would be very slow (complexity would be O(𝑛)or a single point location since each triangle must be checked). A better alternative is to use the adjacency relationships between the triangles, and use a series of Orientation tests, as described in Section 4.4.5, to navigate from one triangle to the other. The idea, called “walking”, is shown in Figure 4.16 and details are given in the Algorithm 2. The idea is as follows: in a DT( 𝑆 ), starting from a triangle 𝜏 (it can be any), we move to one of the adjacent triangle of 𝜏 ( 𝜏 has three neighbours, we choose one neighbour 𝜏𝑖 such that the query point 𝑝 and 𝜏 are on each side of the edge shared by 𝜏 and 𝜏𝑖 ) until there is no such neighbour, then the simplex containing 𝑝 is the current triangle 𝜏 . Notice that this algorithm is not affected by degenerate cases, and that if an Orientation test returns 0 (collinearity), then it is simply considered a positive result. This will ensure that if the query point 𝑝 is located exactly at the same position as one point in 𝑆, then one triangle incident to 𝑝will be returned. 4.4 Incremental construction of the DT 41 Algorithm 2: Walk(T,𝜏,𝑝) 1Input: A DT(𝑆)T, a starting triangle 𝜏, and a query point 𝑝 Output: 𝜏𝑟: the triangle in Tcontaining 𝑝 2𝜏𝑟= None 3while 𝜏𝑟== None do 4visitededges = 0 5for 𝑖←0to 2do 6𝜎𝑖←get edge opposite to vertex 𝑖in 𝜏 7if Orientation (𝜎𝑖, 𝑝)<0then 8𝜏←get neighbouring triangle of 𝜏incident to 𝜎𝑖 9break 10 visitededges += 1 11 if 𝑣𝑖𝑠𝑖𝑡𝑒𝑑𝑒𝑑𝑔𝑒𝑠 == 3then // all the edges of 𝜏have been tested 12 𝜏𝑟=𝜏 13 Return(𝜏𝑟) flip22 flip22 b c d a b c d a Figure 4.17: Aflip22. It should be mentioned that while it appears straightforward, the point location step is the biggest computational bottleneck for a DT implementation. For a large dataset (eg a lidar point cloud, see Chapter 12 for some massive examples), if several thousands/millions of triangles must be tested to find the one containing a give point, then it will be very slow; the insertion itself with a series of flips is generally since around 4 flips will be performed for a normal distribution of points. In practice, because most real-world datasets will have a high spatial coherence (in simple terms, two consecutive points in the dataset are close in reality; see Section 12.3.4), the time spent on walking will be minimised since most library will start the walk from the previously inserted point. 4.4.3 Flips Flips are operations that modify locally the triangulation. There are 3 flip operations (the numbers refer to the number of triangles before and after the flip): ▶ aflip22 modifies the configuration of two adjacent triangles. Consider the set 𝑆={𝑎, 𝑏, 𝑐, 𝑑} of points in the plane forming a quadrilateral, as shown in Figure 4.17. There exist exactly two ways to triangulate 𝑆 : the first one contains the triangles 𝑎𝑏𝑐 and 𝑏𝑐𝑑 ; and the second one contains the triangles 𝑎𝑏𝑑 and 𝑎𝑐𝑑 . Only the first triangulation of 𝑆 is Delaunay because 𝑑 is outside the circumcircle of 𝑎𝑏𝑐 . A flip22 is the operation that transforms the first triangulation into the second, or vice-versa. It is performed in constant time O(1). ▶ aflip13 is the operation of inserting a vertex inside a triangle, and splitting it into three triangles (see Figure 4.18). ▶ aflip31 is the inverse operation that deletes a vertex (see Figure 4.18). 42 4 Delaunay triangulations & Voronoi diagrams flip13 flip31 b c a b c a d Figure 4.18: Aflip13 and its inverse operation flip31. 4.4.4 Controlling the flips To control which triangles have to be checked and potentially flipped, we use a stack ‗ . When the stack is empty, then there are no more triangles to be tested, and we are guaranteed that all the triangles in the triangulation have an empty circumcircle. 4.4.5 Predicates Constructing a DT and manipulating it essentially require two basic geometric tests (called predicates): Orientation determines if a point 𝑝 is left, right or lies on the line segment defined by two points 𝑎 and 𝑏 ; and InCircle determines if a point 𝑝 is inside, outside or lies on a circle defined by three points 𝑎 , 𝑏 and 𝑐 . Both tests can be reduced to the computation of the determinant of a matrix: Orientation(𝑎, 𝑏, 𝑝)= 𝑎𝑥𝑎𝑦1 𝑏𝑥𝑏𝑦1 𝑝𝑥𝑝𝑦1(4.2) InCircle(𝑎, 𝑏, 𝑐, 𝑝)= 𝑎𝑥𝑎𝑦𝑎2 𝑥+𝑎2 𝑦1 𝑏𝑥𝑏𝑦𝑏2 𝑥+𝑏2 𝑦1 𝑐𝑥𝑐𝑦𝑐2 𝑥+𝑐2 𝑦1 𝑝𝑥𝑝𝑦𝑝2 𝑥+𝑝2 𝑦1 (4.3) 4.5 Data structures for storing a DT A triangulation is simply a subdivision of the plane into polygons, and thus any data structure used in GIS can be used to store a triangulation. Simple Features: while many use this (PostGIS and any triangulation you see in Shapefiles), this is not smart: (1) the topological relationships between the triangles are not stored; (2) the vertices are repeated for each triangle (and we know that for a Poisson distribution of points in the plane a given point has exactly 6 incident triangles). Edge-based structures: all the edge-based topological data structure used for storing planar graphs (eg DCEL, half-edge, winged-edge, etc) can be used. These usually lead to large storage space. Observe that in practice, if only the DT is wanted (and not the constrained one, see below), practitioners will often simply store the sample points and reconstruct on-the-fly the DT, since it is unique (if we omit points not in general position that is). However,becauseitissimpler tomanage triangles over arbitrary polygons (they always have exactly 3 vertices and 3 neighbours), data structures specific for triangulations have been developed and are usually used. The simplest data structure, as shown in Figure 4.19, considers the triangle as being its atom and stores each triangle with 3 pointers to its vertices and 3 pointers to its adjacent triangles. Observe that the order in which the vertices and adjacent triangles stored correspond to each other. This ‗ The first-in-last-out data structure: https://en.wikipedia.org/wiki/Stack_(abstra ct_data_type) 4.5 Data structures for storing a DT 43 τ τa τb τc a b c d triangle v1v2v3adj1adj2adj3 τa b c τaτbτc τab d c τ... τ τ... ... ... ... ... ... ... ... Figure 4.19: The triangle-based data structure to store efficiently a triangulation (and the adjacency relationships between the triangles). is an important property that allows an efficient retrieval of triangles in the Walk algorithm (Algorithm 2) for instance. 44 4 Delaunay triangulations & Voronoi diagrams Figure 4.20: (top) A set 𝑆 of points and straight-line segments. (middle) Constrained DT of 𝑆 .(bottom) Conforming DT of 𝑆 ; the Steiner points added are in red. Steiner point 4.6 Constrained and Conforming Delaunay Triangulations Given as input a set 𝑆 of points and straight-line segments in the plane, different triangulations of 𝑆 (so that the segments are respected) can be constructed. We are mostly interested in the constrained Delaunay triangulation (ConsDT) and the conforming Delaunay triangulation (ConfDT), see Figure 4.20 for one example. Constrained DT (ConsDT). Given a set 𝑆 of points and straight-line segments in ℝ2 , the ConsDT permits us to decompose the convex hull of 𝑆 into non-overlapping triangles, and every segment of 𝑆 appears as an edge in ConsDT( 𝑆 ). ConsDT is similar to the Delaunay triangulation, but the triangles in ConsDT are not necessarily Delaunay (ie their circumcircle might contain other points from 𝑆 ). The empty circumcircle for a ConsDT is less strict: a triangle is Delaunay if its circumcircle contains no other points in 𝑆 that are visible from the triangle. The constrained segments in 𝑆act as visibility blockers. Figure 4.21 shows one example. Without going into details about one potential algorithm, one way to construct a ConsDT(𝑆) is (see Figure 4.22): 1. construct DT( 𝑆𝑝 ), where 𝑆𝑝 is the set containing all the points in 𝑆 and the end points of the line segments (Figure 4.22b) 2. insert each line segment, each insertion will remove edges from DT(𝑆𝑝). In Figure 4.22c 3 edges are removed. 3. this creates 2 polygons that need to be retriangulated, in Figure 4.22d there is a blue and a green one. 4. retriangulate each separately, the Delaunay criterion needs to be verified only for the vertices incident to the triangles incident to the hole/polygon. Observe that the ConsDT can be used to triangulate polygons with holes (see Figure 4.23), it suffices to remove the triangle outside the exterior boundary, but inside the convex hull. Conforming DT (ConfDT). A ConfDT adds new points to the input 𝑆 (called Steiner points) to ensure that the input segments are present in the triangulation. As Figures 4.20 and 4.23 show, the input straight-line segments will be potentially split into several collinear segments. The 3 buildings DT of the vertices of the buildings Constrained DT Figure 4.21: The ConsDT of a set of segments. On the right, the triangle whose circumcircle is green is a Delaunay (no other points in its interior) and so is the triangle whose circumcircle is in purple (there is one point in its interior, but it cannot be seen because of the constrained segment). 4.6 Constrained and Conforming Delaunay Triangulations 45 (a) (b) (c) (d) (e) Figure 4.22: Steps to construct a ConsDT. Steiner points have to be carefully chosen (where to put them is beyond the scope of this course). Observe that each triangle in a ConfDT respects the Delaunay criterion, but that more triangles are present. If 2 segments are nearly parallel, many points could be necessary (for 𝑚 segments, up to 𝑚2 could be necessary). (a) (b) (c) Figure 4.23: (a) One polygon with 4 holes (interior rings). (b) its ConsDT. (c) its ConfDT (the Steiner point added is in red). 46 4 Delaunay triangulations & Voronoi diagrams 4.7 Notes and comments The DT and the VD have been discovered, rediscovered and studied many times and in many different fields, see Okabe et al. (2000) for a complete history. The VD can be traced back to 1644, when Descartes used Voronoi-like structures in Part III of his Principia Philosophiæ. The VD was used by Dirichlet (1850) to study quadratic forms—this is why the VD is sometimes referred to as Dirichlet tessellation—but was formalised and defined by Voronoi (1908). The first use of the VD in a geographical context is due to Thiessen (1911), who used it in climatology to better estimate the precipitation average around observations sites; the DT was formalised by Delaunay (1934). For the construction of the DT, the incremental algorithm was first described by Lawson (1972). Fortune (1987) describes a sweep-line one, and Guibas and Stolfi (1985) a divide-and-conquer algorithm. The local optimality of a DT, which implies globally optimality in the case of the DT, was proven by Delaunay (1934) himself. The max-min angle optimality of the DT was firstly observed by Sibson (1978). This parabolic lifting was first observed by Brown (1979) (who used a spherical transformation), further described by Seidel (1982) and Edelsbrunner and Seidel (1986). Liu and Snoeyink (2005) explains the details of the infinite vertex. The walking algorithm described in this chapter, with a few modifications, can perform point location in O(𝑛1/3 ). However, it is in theory not the fastest solution: Mücke et al. (1999) and Devillers et al. (2002) discuss alternatives that are optimal (ie O(log 𝑛) ). However, they both note that optimal algorithms do not necessarily mean better results in practice because of the amount of preprocessing involved, the extra storage needed, and also because the optimal algorithms do not always consider the dynamic case, where points in the DT could be deleted. Several criteria for constructing data-dependent triangulations are discussed in Dyn et al. (1990). While these can be used, in practice it was proven that the Delaunay triangulation is still the triangulation that minimises the roughness of a surface (Wang et al., 2001; Rippa, 1990) Shewchuk (1997) shows that while the triangle-based data structure requires twice as much code as with the quad-edge (to store and construct a ConsDT), the result is that the code runs twice as fast and the memory requirement as about 2X less. CGAL ( https://www.cgal.org/ ), among many others, uses the triangle-based data structure. Since a DT can be locally modified by adding one point (and not reconstructing the whole structure from scratch, see Figure 4.12), it is also possibletodelete/remove one vertexfrom aDTwith only local operations. Mostafavi et al. (2003) and Devillers (2009) describe algorithms. 4.8 Exercises 1. A DT has 32 triangles and we insert a new point 𝑝 that falls inside one of the triangles. If we insert and update the triangulation (for Delaunay criterion), what is the number of triangles? 4.8 Exercises 47 2. Given the input formed of elevation points and breaklines below (both projected to the 𝑥𝑦 -plane), draw both the constrained and conforming Delaunay triangulation (an approximation is fine). Constained DT Conforming DT 3. If a given vertex 𝑣 in a DT has 7 incident triangles, how many vertices will its dual polygon contain? 4. Identify the 5 infinite triangles in Figure 4.15. 5. A DT has 6 vertices, and 3 of these are forming the convex hull. How many triangles does the DT have? 6. Assume you have 8 points located on a circle. Draw the DT and the VD of these 8 points. 7. When inserting points in a DT (Algorithm 1), what happens if a new point is inserted directly on an edge? Line 2 states that the triangle is split into 3 new triangles, does it still hold? 54 5 Spatial interpolation: deterministic methods Figure 5.7: A TIN is obtained by lifting the vertices to their elevation. All the triangles are usually Delaunay, ie their circumcircle(green) is emptyofany other points in the plane. p0 p2 p1 x A0 A1 A2 wi(x) = Ai Figure 5.8: Barycentric coordinates. 𝐴𝑖 defines the area of a triangle. to their elevation value). The number of samples used in the interpolation is therefore always 3, and their weight is based on the barycentric value (see below). To obtain satisfactory results, this method is usually used in 2D with a Delaunay triangulation because, among all the possible triangulations of a set of points in the plane, it maximizes the minimum angle of each triangle. Themethodisexact,continuous,local,adaptative,efficient,andautomatic. Only the property #3 is not fulfilled (at the edges of the triangles). If the point location strategy is used to identify the triangle containing 𝑥 (Section 4.4.2), then O(𝑛1/3 ) on average is used. The interpolation itself is performed in constant time. Data-dependent triangulations. It was shown in Chapter 4and in Figure 1.5 that, for terrain modelling, the Delaunay triangulation is preferred over other triangulations because it favours triangles that are as equilateral as possible. However, it should be noticed that the elevation of the vertices are not taken into account to obtain the DT, ie if we changed the elevation of the samples we would always get the same triangulation. One might therefore wonder whether the DT best approximates the morphology of a terrain. A triangulation that considers the elevation (or any 𝑧 coordinate) is called adata-dependent triangulation. The idea is to define a set of criteria (instead of the empty circumcircle). One example is trying to minimise the change in normals for the two incident triangles of an edge. While such methods will yield longer and skinnier triangles, these might better approximate the shape of the terrain for some specific cases. One drawback of these methods is that different criteria will be required for different cases, and that computing such triangulation can be computationally expensive. In practice, one would need to first compute the DT, and then take each edge (and the two incident triangles), and perform a local flip based on the elevation values; the final triangulation is obtained by optimisation the wished criterion. Barycentric coordinates. The linear interpolation in a triangle can be efficiently implemented by using barycentric coordinates, which are local coordinates defined within a triangle. Referring to Figure 5.8, any point 𝑥 inside a triangle 𝑝0𝑝1𝑝2 can be represented as a linear combination of the 3 vertices: 𝑥=𝑤0𝑝0+𝑤1𝑝1+𝑤2𝑝2(5.3) and 𝑤0+𝑤1+𝑤2=1(5.4) The coefficients 𝑤𝑖 are the barycentric coordinates of the point 𝑥 with respect to the triangle 𝑝0𝑝1𝑝2 . Finding the coefficients 𝑤0 , 𝑤1 , and 𝑤2 can be done by solving a system of linear equations. If we subtract 𝑝2 from 𝑥 , and we use 𝑤2=1−𝑤0−𝑤1, we obtain 𝑥−𝑝2=𝑤0(𝑝0−𝑝2)+𝑤1(𝑝1−𝑝2)(5.5) We obtain 2 vectors ( 𝑝0−𝑝2 and 𝑝1−𝑝2 ), which represent 2 edges of the triangle. This equation can be solved and we find that the 3 coefficients are equal of the area of the 3 triangle subdividing the original triangle (as shown in Figure 5.8). 5.3 Weighted-average methods 55 x Figure 5.9: The VD of a set of points with an interpolation location 𝑥. x p1 p2 p3 p4 p5 p6 w6(x) wi(x) = stolenarea Figure 5.10: Natural neighbour coordinates in 2D for 𝑥 . The shaded polygon is V + 𝑥. Higher-order function in each triangle (TIN-c1). Itispossibletomodify the linear function inside each triangle by a higher-order function. As is the case for splines, there are several ways to achieve this, and the details of these is out of scope for this course. These methods are usually used more for finite element analysis where the flow of a certain fluid (eg wind) around or through a mechanical piece is studied. Most methods would define a cubic Bézier polynomial inside each triangle (which is 𝐶1 ), and then ensure that the function is 𝐶1 along the edges and at the 3 vertices of the triangles. To achieve this the normals of each vertex is calculated by averaging the normals of the incident triangles, and the normal along an edge is computed similarly with the 2 incident triangles. 5.3.4 Natural Neighbour Interpolation (NNI) This is a method based on the Voronoi diagram for both selecting the data points involved in the process, and assigning them a weight. It is also called Sibson’s interpolation, after the name of its inventor. It uses two VDs: one for the set 𝑆 of data points (Figure 5.9), and another one where a point 𝑥 is inserted at the estimation location (Figure 5.10). The insertion of 𝑥 modifies locally a VD( 𝑆 ): the Voronoi cell V 𝑥 of 𝑥 ‘steals’ some parts of some Voronoi cells of VD(𝑆). This idea forms the basis of natural neighbour coordinates, which define quantitatively the amount V 𝑥 steals from each of its natural neighbours (Figure 5.10). Let D be the VD( 𝑆 ), and D+=D∪{𝑥} . The Voronoi cell of a point 𝑝 in D is defined by V 𝑝 , and V + 𝑝 is its cell in D+ . The natural neighbour coordinate of 𝑥with respect to a point 𝑝𝑖is 𝑤𝑖(𝑥)=𝐴𝑟𝑒𝑎(V 𝑝𝑖∩V + 𝑥) 𝐴𝑟𝑒𝑎(V + 𝑥)(5.6) where 𝐴𝑟𝑒𝑎(V 𝑝𝑖) represents the area of V 𝑝𝑖 . For any 𝑥 , the value of 𝑤𝑖(𝑥) will always be between 0 and 1: 0 when 𝑝𝑖 is not a natural neighbour of 𝑥 , and 1 when 𝑥 is exactly at the same location as 𝑝𝑖 . A further important consideration is that the sum of the areas stolen from each of the 𝑘 natural neighbours is equal to 𝐴𝑟𝑒𝑎(𝑉+ 𝑥), in other words: 𝑘 X 𝑖=1 𝑤𝑖(𝑥)=1.(5.7) Therefore, the higher the value of 𝑤𝑖(𝑥) is, the stronger is the ‘influence’ of 𝑝𝑖 on 𝑥 . The natural neighbour coordinates are influenced by both the distance from 𝑥to 𝑝𝑖and the spatial distribution of the 𝑝𝑖around 𝑥. Natural neighbour interpolation is based on the natural neighbour coordinates. The points used to estimate the value of an attribute at location 𝑥 are the natural neighbours of 𝑥 , and the weight of each neighbour is equal to the natural neighbour coordinate of 𝑥 with respect to this neighbour. The natural neighbour interpolant possesses all the wished properties from above, except that the first derivative is undefined at the data points. Its main disadvantage is that its implementation is rather complex, and obtaining an efficient one is not simple and involves complex manipulation of the VD. From Section 4.4 we know that one insertion of a single 56 5 Spatial interpolation: deterministic methods Figure 5.11: Top: The NNI interpolant in 1D is equivalent to a linear interpolation. Bottom: If the gradient at each sample points are calculated/estimated, then it is possible to modify the weights so that a𝐶1interpolant is obtained. point 𝑝 in a DT can be done in O(log 𝑛) , but the deletion of a point is a more complex operation (outside the scope of this book). Higher-order function (NNI-c1). The NNI method can be thought of performing linear interpolation, in the 1D case (where we have one independent variable) then it is equivalent to a linear interpolant (see Figure 5.11). It can be modified so that the first derivative is possible everywhere, including at the data points. This is achieved by modifying the weights so that they are not linear anymore. The gradient of the surface at each sample point is taken into account, ie for each data point we can estimate the slope (with a linear function, a plane) and modify the weights; how this is done is out of scope for this book. The resulting interpolant is 𝐶1 , and Figure 5.12 shows an example. 5.3.5 Laplace interpolant The Laplace interpolant, or non-Sibsonian interpolation, is a computationally faster variant of the natural neighbour interpolation method. It is faster because no (stolen) areas need to be computed, instead the lengths of the Delaunay and the Voronoi edges are used. Figure 5.12: Notice how the NNI interpolant creates “inverted cups” around each sample point, and how NNI-c1 results in a more rounded surface. (a) NNI (b) NNI (𝐶1) 5.3 Weighted-average methods 57 x p1 p2 p3 p4 p5 p6 w6(x) = |e| |xp6| e Figure 5.13: The weight for the Laplace interpolant for one neighbour (𝑝6.) input grid output grid Figure 5.14: Resampling of an input grid, theoutputgrid hasadifferentorientation and a different resolution. For a given interpolation location 𝑥 , the natural neighbours 𝑝𝑖 of 𝑥 are used for the Laplace interpolant. The weight 𝑤𝑖 of a 𝑝𝑖 is obtained, as shown in the Figure 5.13, by: 𝑤𝑖(𝑥)=|𝑒𝑑𝑔𝑒𝑖(V + 𝑥)| |𝑥𝑝𝑖|(5.8) where |𝑒𝑑𝑔𝑒𝑖(V + 𝑥)| represents the length of the Voronoi edge dual to the Delaunay edge 𝑥𝑝𝑖 (the orange edge in Figure 5.13 for one neighbour); and |𝑥𝑝𝑖|the Euclidean distance (in 2D) between 𝑥and 𝑝𝑖. If we consider that each data point in 𝑆 has an attribute 𝑎𝑖 (its elevation), the interpolation function value at 𝑥is: 𝑓(𝑥)=P𝑘 𝑖=1𝑤𝑖(𝑥)𝑎𝑖 P𝑘 𝑖=1𝑤𝑖(𝑥)(5.9) Note that the fraction becomes indeterminate when 𝑥 equals one of the sample points 𝑝𝑖 . In this case the Laplace interpolant therefore simply defines that 𝑓(𝑥)=𝑎𝑖. Firstly the Laplace interpolant is exact: the interpolation method returns the exact value, rather than some estimate, of a sample point when it is queried at that precise location. Secondly, it is continuous and continuously differentiable ( 𝐶1 ) everywhere except at sites where finitely many Voronoi circles intersect. Thirdly, it is local, ie it uses only a local subset of data for the interpolation of a point. This limits the computational cost and supports efficient addition or removal of new data points. Finally, like the VD itself, it is adaptive to the spatial configuration of sample points. Unlike other methods such as IDW interpolation, the Laplace interpolant requires no user-defined parameters. 5.3.6 Bilinear Interpolation When one wants to know the value of the elevation at a location 𝑝 , one can simply look at the value of the pixel (which is equivalent to using nearest neighbour interpolation), but this method has many drawbacks, for example when one needs to resample a grid. Resampling means transforming an input grid so that the resolution and/or the orientation are different, see Figure 5.14. Bilinear interpolation has been shown to give better results than using the value of the pixel. The method, which can be seen as an ‘extension’ of linear interpolation for raster data, performs linear interpolation in one dimension (say along the 𝑥 axis), and then in the other dimension ( 𝑦 ). Here one has to be careful about the meaning of a grid: does the value of a pixel represent the value of the whole pixel? or was the grid constructed by sampling the values at the middle of each pixel? In most cases, unless metadata are available, it is not known. But in the context of terrain modelling, we can assume that the value of a pixel represents the value at the centre of the pixel. Suppose we have 4 adjacent pixels, each having an elevation, as in Figure 5.15. Bilinear interpolation uses the 4 centres to perform the interpolation at location 𝑝=(𝑝𝑥, 𝑝𝑦) ; it is thus a weighted-average method because the 4 samples are used, and their weight is based on the linear interpolation, as explained below. We need to linearly interpolate the values at locations 𝑞 and 𝑟 with linear interpolation, and then linearly interpolate along the 58 5 Spatial interpolation: deterministic methods Figure 5.15: Bilinear interpolation. 10 12 15 18 p n4 = 10 n3 = 12 n1 = 15 n2 = 18 p q r cross-validation root-mean-square error 𝑦 axis with these values. Also, notice that the result is independent of the order of interpolation: we could start with interpolating along the 𝑦 axis and then the 𝑥 axis and we would get the same result. For the case in Figure 5.15, the calculation would go as follows: 𝑞𝑧=𝑝𝑥−𝑛4𝑥 𝑛3𝑥−𝑛4𝑥×(𝑛3𝑧−𝑛4𝑧)+𝑛4𝑧 𝑟𝑧=𝑝𝑥−𝑛1𝑥 𝑛2𝑥−𝑛1𝑥×(𝑛2𝑧−𝑛1𝑧)+𝑛1𝑧 𝑝𝑧=𝑝𝑦−𝑟𝑦 𝑞𝑦−𝑟𝑦∗(𝑞𝑧−𝑟𝑧)+𝑟𝑧 (5.10) 5.4 Assessing the results of an interpolation method and/or fine-tuning the parameters Finding the “best” interpolation method for a given dataset, and the most suitable parameters (if any are needed), can be a rather tricky task in practice because we most often do not have extra control points. One simple technique, which is also very easy to implement, is called jackknife, or cross-validation. It is a simple statistics resampling technique to estimate the bias and the variance of an estimation. Imagine you have a dataset 𝑆 consisting of 𝑛 sample points. The main idea is to remove/omit from 𝑆 one sample point 𝑝 and calculate the estimation ˆ 𝑎𝑝 obtained for the elevation at the location ( 𝑥, 𝑦 ) of 𝑝 , and to compare this value with the real value 𝑎𝑝 . And then to repeat this for each of the 𝑛 points in 𝑆 ; each estimation is thus obtained with 𝑛−1 points. One method (with given parameters) for a given dataset can be characterised by computing the root-mean-square error: 𝑅𝑀𝑆𝐸 =rP𝑛 𝑖=1(ˆ 𝑧𝑖−𝑧𝑖)2 𝑛(5.11) And it is a good idea to plot the results to observe where the largest differences between the estimation and the real values are obtained, this can help in identifying which parameters should be fine-tuned. See for instance one example in Figure 5.16. It can be seen in Figure 5.16c that the largest differences between the observed and estimated values are (mostly) concentrated around the two peaks of the terrain, which is not surprising. The differences in the lower areas (which is water) are smaller since these areas have a flatter morphology. Figure 5.16d shows the same absolute differences but in a scattered plot of the observed values versus the estimated ones. 5.5 Overview of all methods 59 (a) 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 x (km) 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 y (km) 0 200 400 600 800 (b) 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 x (km) 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 y (km) 25 50 75 100 125 150 175 (c) 0 200 400 600 800 Real value 0 100 200 300 400 500 600 700 800 Estimated value 25 50 75 100 125 150 175 (d) Figure 5.16: (a) A terrain of a given area containing 2 hills. (b) A sample of 1000 points of this terrain. (c) A plot of the errors (absolute values) obtained from the jackknife (with IDW and a given search radius and power). (d) A plot of the absolute elevation versus the estimated ones . 5.5 Overview of all methods Figure 5.17 shows the result of 8 different interpolants for the same (real-world) sample points, and Table 5.1 give an overview of their properties. 60 5 Spatial interpolation: deterministic methods (a) Nearest neighour (b) IDW (radius=1500m; pow=2) (c) IDW (radius=1500m; pow=4) (d) TIN (linear) (e) TIN (𝐶1)(f) Natural neighbours (g) Natural neighbours (𝐶1)(h) Laplace Figure 5.17: Results of a few interpolation methods for the same dataset; the samples are shown on the surface (red dots). 5.6 Notes and comments 61 Download the terrains Download the 8 resulting terrains and explore them in a GIS software, eg QGIS. https://tudelft3d.github.io/terrainbook/extra/interpol/. Table 5.1: Overview of the interpolation methods discussed in this chapter, with their properties (as described in Section 5.1). exact continuous local adaptable efficient automatic global function ×𝐶2+×– – × splines ×𝐶2+depends 0 - × nearest neigh. ✓×✓+ ++ ✓ IDW ✓×✓- 0 × TIN ✓𝐶0✓+ ++ ✓ NNI ✓𝐶0✓++ 0 ✓ NNI-c1 ✓𝐶1✓++ - ✓ Laplace ✓𝐶0✓++ + ✓ bilinear ✓𝐶0✓++ ++ ✓ 5.6 Notes and comments Watson (1992), in his authoritative book, lists the essential properties of an ‘ideal’ interpolation method for bivariate geoscientific datasets; we have added computationally efficient and automatic to the list. Mitasova and Mitas (1993) gives a full description of the regularised splines with tension (RST) interpolation method. This method has also been implemented in the open-source GIS GRASS. For a discussion about influence of the power in IDW on the resulting surface, please see Watson (1992). The description of the barycentric coordinates is mostly taken from Eberly (2018). The natural neighbour interpolation method is also called Sibson’s interpolation, after the name of the inventor (Sibson, 1981). An excellent summary of the methods to modify Equation 5.6 to obtain a continuous function is found in Flötotto (2003). The Laplace interpolant was discovered independently by Belikov et al. (1997) and Hiyoshi and Sugihara (1999) and Hiyoshi and Sugihara (2000). The regularised spline with tension (RST) is available in the open-source GIS GRASS. The construction of a polynomial inside each triangle of a TIN can be done with several methods. The simplest method is the Clough-Tocher method (Clough and Tocher, 1965; Farin, 1985). It splits each triangle into 3 sub-triangles (by inserting a temporary point at the centroid of the triangle) and a cubic function is built over each. Dyn et al. (1990) shows how to obtain a data-dependent triangulation. Rippa (1990) proves that the DT is the triangulation that minimizes the roughness of the resulting terrain, no matter what the actual elevation of the data is. Here, roughness is defined as the integral of the square of the 𝐿2 -norm of the gradient of the terrain. Gudmundsson et al. (2002) shows that a variation of the DT (one where 𝑘 vertices can be inside the 62 5 Spatial interpolation: deterministic methods circumcircle of a given triangle) can yield fewer local minima; whether it yields a “better’ terrain is an open question. 5.7 Exercises 1. Given a triangle 𝜏 with coordinates (20.0, 72.0, 21.0), (116.0, 104.0, 32.0), and (84.0, 144.0, 26.0), estimate the elevation at 𝑥 = (92.0, 112.0) with linear interpolation in the triangle (both by finding the equation of the plane and with barycentric coordinates). 2. What happens when the search distance is very large for inverse distance weighting interpolation (IDW)? 3. For grids, can IDW or others be used instead of bilinear? If yes, how does that work? 4. The 15 elevation samples below have been collected. You want to interpolate at two locations: a) at location ( 7,6 ) with IDW (radius=3; power=2); the purple circle. b) at location ( 15,6 ) with linear interpolation in TIN; the orange cross. What are the respective answers? 10 155 5 10 0 15 12 9 5 2 3 5 5 67 18 16 14 10 8 5 random process stochastic process probability distribution random variable Spatial interpolation: kriging 6 uyoutu.be/ZaadzBvgE2s 6.1 Statistical background . . . 63 6.2 Geostatistical model . . . . 65 6.3 Covariance, dissimilarity and the semivariogram . . 66 6.4 Simple kriging . . . . . . . 70 6.5 Ordinary kriging . . . . . . 72 6.6 Other types of kriging . . . 73 6.7 Implementation details . . 74 6.8 Notes and comments . . . . 74 6.9 Exercises . . . . . . . . . . . 75 Kriging is a spatial interpolation method that was developed mostly by Georges Matheron based on the earlier work of Danie Krige, who created it to estimate the yield of gold mines in South Africa. In contrast to other spatial interpolation methods, it involves creating a custom model that is fine-tuned using the statistical properties of each dataset. In this way, kriging can take into account the specific characteristics of a dataset, often yielding better results than other interpolation methods. Like other techniques based on geostatistical models, kriging relies on the fact that when one moves across space, values such as the gold content in rock or the elevation in a terrain have both a general spatial trend (eg a flat mean value, a fitted plane or a more complex polynomial defining a surface) and a certain spatially correlated randomness (ie closer points tend to have more similar values). Both of these elements are modelled in kriging. More than a single method, kriging comprises a family of related methods. Within this chapter, we will look at two related types of kriging: simple kriging and ordinary kriging. These treat the spatially correlated randomness in a similar way, but they make different assumptions about the trend in a dataset. 6.1 Statistical background The physical processes that shape the world can be considered to be at least partly deterministic. In the case of a terrain, the elevation is determined by processes that we can model (more or less accurately), such as plate tectonics, volcanic activity, and erosion. However, these processes are too complex and not understood well enough to use them to obtain accurate elevation values. Imagine, for instance, how difficult it would be to get an accurate elevation map of the world using only the shape of the tectonic plates and some other parameters (eg their direction and speed of movement). Because of this complexity, the value of complex properties, such as the elevation of a terrain, are usually treated in geostatistics as the result of what is known as a random or stochastic process. In this context, randomness can be understood as the fact that the value of a property at an unsampled location is not known exactly, and so we cannot assign it an exact number. Instead, we can make an educated guess of the value at that location by creating a statistical model of its possible values using a probability distribution, which we can associate with a function (ie a probability distribution function) or with a set of standard statistical measures, such as the mean and variance. This situation is phrased in mathematical terms by saying that the value of the elevation property 𝑍 at a location 𝑥 isarandom variable 𝑍(𝑥) . For the sake of simplicity, we will usually omit the location and denote it just as 𝑍 ; or when working with multiple locations (eg 𝑥𝑖 and 𝑥𝑗 ), we will shorten their respective random variables ( 𝑍(𝑥𝑖) and 𝑍(𝑥𝑗) ) using subscripts (𝑍𝑖and 𝑍𝑗). 70 6 Spatial interpolation: kriging minimisation of the variance Beforemoving on to applythese tokriging, there are a couple of important points. First, note that these theoretical functions are often only applied when |ℎ|>0 , since setting 𝛾(0)=0 helps to ensure that kriging passes exactly through the sample points (the exact property as explained in Section 5.1). Second, all of the semivariogram-related functions seen in this section can be converted to covariance functions as well, taking into account that 𝛾(ℎ)=sill −𝐶(ℎ) . Note that this means that the covariance is high when |ℎ|is small and it decreases as |ℎ|increases. 6.4 Simple kriging Simple kriging is similar to other spatial interpolation methods that use a weighted average. It starts from the assumption of second-order stationarity. Moreover, the expectation is also known, and so the general procedure to perform it is to: (i) subtract it from the sample points to obtain residuals, (ii) use the residuals to define a function that estimates the value of the residual term at any location, and (iii) interpolate at the desired locations using the function added to the expectation. Thus, simple kriging defines a function ˆ 𝑅0 that estimates the value of the residual 𝑅 of the random variable 𝑍 at a location 𝑥0 as a weighted average of its residuals at the 𝑛 sample points 𝑥𝑖 that we will use for the interpolation (where 1≤𝑖≤𝑛). We denote this as: ˆ 𝑅0=ˆ 𝑍0−𝐸[𝑍0]= 𝑛 X 𝑖=1 𝑤𝑖(𝑍𝑖−𝐸[𝑍𝑖] | {z } 𝑅𝑖 ).(6.19) Simple kriging is unbiased, and therefore the expected value of the estimation at a location 𝑥0 is equal to the expected value at that location. In mathematical terms, we can formulate this as: 𝐸ˆ 𝑍0−𝑍0=0or 𝐸[𝑍0]=𝐸ˆ 𝑍0.(6.20) Then, in order to derive the equations used in simple kriging, we start from the fact that it minimises the variance of the estimation error, which in this case is given by var ˆ 𝑅0−𝑅0 . If we use the definition of the variance from Equation 6.1, this can be instead put in terms of an expectation: var ˆ 𝑅0−𝑅0=𝐸ˆ 𝑅0−𝑅0−𝐸ˆ 𝑅0−𝑅02 However, we know from the unbiased criterion from Equation 6.20 that 𝐸ˆ 𝑅0−𝑅0=0, and so we can simplify the previous equation as: var ˆ 𝑅0−𝑅0=𝐸ˆ 𝑅0−𝑅02. If this is expanded, it results in: 6.4 Simple kriging 71 simple kriging system var ˆ 𝑅0−𝑅0=𝐸ˆ 𝑅2 0−2ˆ 𝑅0𝑅0+𝑅02 =𝐸ˆ 𝑅2 0−2𝐸ˆ 𝑅0𝑅0+𝐸𝑅02 =𝐸"𝑛 X 𝑖=1 𝑛 X 𝑗=1 𝑤𝑖𝑤𝑗𝑅𝑖𝑅𝑗#−2𝐸"𝑛 X 𝑖=1 𝑤𝑖𝑅𝑖𝑅0#+𝐸𝑅02 = 𝑛 X 𝑖=1 𝑛 X 𝑗=1 𝑤𝑖𝑤𝑗𝐸𝑅𝑖𝑅𝑗−2 𝑛 X 𝑖=1 𝑤𝑖𝐸[𝑅𝑖𝑅0]+𝐸𝑅02. Here, we can use the definitions of the variance based on residuals from Equations 6.7 and 6.8 together with our covariance formula from Equation 6.9, which yields: var ˆ 𝑅0−𝑅0= 𝑛 X 𝑖=1 𝑛 X 𝑗=1 𝑤𝑖𝑤𝑗cov(𝑅𝑖, 𝑅𝑗)−2 𝑛 X 𝑖=1 𝑤𝑖cov(𝑅𝑖, 𝑅0)+cov(𝑅0, 𝑅0) (6.21) = 𝑛 X 𝑖=1 𝑛 X 𝑗=1 𝑤𝑖𝑤𝑗𝐶(𝑥𝑖−𝑥𝑗)−2 𝑛 X 𝑖=1 𝑤𝑖𝐶(𝑥𝑖−𝑥0)+𝐶(𝑥0−𝑥0). (6.22) In order to minimise this equation, we can find where its first derivative is zero. This is: 𝜕var ˆ 𝑅0−𝑅0 𝜕𝑤𝑖 =2 𝑛 X 𝑗=1 𝑤𝑗𝐶(𝑥𝑖−𝑥𝑗)−2𝐶(𝑥𝑖−𝑥0)=0for all 1≤𝑖≤𝑛, which yields the set of 𝑛simple kriging equations: 𝑛 X 𝑗=1 𝑤𝑗𝐶(𝑥𝑖−𝑥𝑗)=𝐶(𝑥𝑖−𝑥0).(6.23) While these equations can be used to perform simple kriging, it is often easier to deal with these in matrix form: ©« 𝑤1 . . . 𝑤𝑛ª®®¬ | {z } 𝑤 ©« 𝐶(𝑥1−𝑥1) ··· 𝐶(𝑥1−𝑥𝑛) . . ..... . . 𝐶(𝑥𝑛−𝑥1) ··· 𝐶(𝑥𝑛−𝑥𝑛)ª®®¬ | {z } 𝐴 =©« 𝐶(𝑥1−𝑥0) . . . 𝐶(𝑥𝑛−𝑥0)ª®®¬ | {z } 𝑑 (6.24) which is known as the simple kriging system. Finally, if we invert the matrix 𝐴, the weights are given by: 𝑤=𝐴−1𝑑. (6.25) 72 6 Spatial interpolation: kriging These weights can be applied to interpolate the value of the residual term at any location as a weighted average of the sample points, where the correlation between the points is given by a covariance function, which can be obtained from the semivariogram. However, simple kriging does not tell us what value we should use for the expectation 𝐸[𝑍] , since we start from the assumption that it is known, which is often not the case. Even a seemingly reasonable value, such as average of all points, can be very inaccurate if the sample points are unevenly distributed across the domain. 6.5 Ordinary kriging Ordinary kriging is similar to simple kriging in that it estimates values using a weighted average function with weights computed from the semivariogram/covariance based on the distance between the points. However, it estimates a value at a location using only sample points in the neighbourhood of the location and relies only on local second-order stationarity, ie a constant mean within a moving neighbourhood, with a mean that is computed within the method itself. Rather than relying on the residuals, it defines a function ˆ 𝑍0 that directly estimates the value of the random variable 𝑍 at a location 𝑥0 as a weighted average of its value at the 𝑛neighbouring sample points 𝑥𝑖that we will use for the interpolation (where 1≤𝑖≤𝑛). We denote this as: ˆ 𝑍0= 𝑛 X 𝑖=1 𝑤𝑖𝑍𝑖.(6.26) Like simple kriging, ordinary kriging is unbiased, which is achieved by making sure that the interpolation weights add up to one, ie P𝑛 𝑖=1𝑤𝑖=1 . Ordinary kriging also minimises the variance of the estimation error, which is given by var ˆ 𝑍0−𝑍0 . For this, we can use the same derivation as for simple kriging up to Equation 6.21 but using the variogram for the final step. This is: var ˆ 𝑅0−𝑅0= 𝑛 X 𝑖=1 𝑛 X 𝑗=1 𝑤𝑖𝑤𝑗cov(𝑅𝑖, 𝑅𝑗)−2 𝑛 X 𝑖=1 𝑤𝑖cov(𝑅𝑖, 𝑅0)+cov(𝑅0, 𝑅0) =− 𝑛 X 𝑖=1 𝑛 X 𝑗=1 𝑤𝑖𝑤𝑗𝛾(𝑥𝑖−𝑥𝑗)+2 𝑛 X 𝑖=1 𝑤𝑖𝛾(𝑥𝑖−𝑥0)−𝛾(𝑥0−𝑥0). (6.27) Using the previous equation and the unbiased criterion from Equation 6.20, we can apply the minimisation method known as Lagrange multipliers‗and arrive at the set of 𝑛+1ordinary kriging equations: ‗https://en.wikipedia.org/wiki/Lagrange_multiplier 6.6 Other types of kriging 73 ordinary kriging system 𝑛 X 𝑗=1 𝑤𝑖𝛾(𝑥𝑖−𝑥𝑗)+𝜇(𝑥0)=𝛾(𝑥𝑖−𝑥0)for all 1≤𝑖≤𝑛 𝑛 X 𝑗=1 𝑤𝑖=1(6.28) where 𝜇(𝑥0) is a Lagrange parameter that was used in the minimisation process. Like with simple kriging, these equations can be used to perform ordinary kriging, but it is often easier to deal with these in matrix form: ©« 𝛾(𝑥1−𝑥1) ··· 𝛾(𝑥1−𝑥𝑛)1 . . ..... . .1 𝛾(𝑥𝑛−𝑥1) ··· 𝛾(𝑥𝑛−𝑥𝑛)1 1··· 1 0 ª®®®®¬ | {z } 𝐴 ©« 𝑤1 . . . 𝑤𝑛 𝜇(𝑥0)ª®®®®¬ | {z } 𝑤 =©« 𝛾(𝑥1−𝑥0) . . . 𝛾(𝑥𝑛−𝑥0) 1ª®®®®¬ | {z } 𝑑(6.29) which is known as the ordinary kriging system. Finally, if we invert the matrix 𝐴 , the weights and the Lagrange multipliers are given by: 𝑤=𝐴−1𝑑(6.30) 6.6 Other types of kriging Directional kriging is useful when the similarity of points differs according to different directions, eg north-south versus east-west. It involves creating variograms for different directions. Block kriging attempts to obtain a local trend for an area around a point (rather than just at a point). It can be used to obtain smoother results. Cokriging applies kriging to multiple correlated variables. It is particularly useful when you want to interpolate one variable with limited sample points but there are more sample points for another correlated variable. Indicator kriging applies thresholds to obtain discrete values, such as presence/absence of a variable or a set of distinct classification classes. Poisson kriging applies kriging to data involving cumulative counts or rates. It is often applied together with polygonal datasets. Universal kriging fits the trend using a predefined deterministic function. It is tricky to use in practice because it imposes conditions on the underlying variogram. 74 6 Spatial interpolation: kriging 6.7 Implementation details There are a few important details with respect to the implementation of kriging methods in practice. First of all, within this chapter, we have assumed that you always use all sample points to interpolate any point on the plane. While this is optimal in theory, if a large number of sample points are used to interpolate every point, kriging can be very slow in practice. The reason for this is because matrix 𝐴 will be very large, and inverting a matrix is a computationally expensive process. Since the weights of far-away sample points are usually very small, the usual solution is to limit the number of sample points used, either by using a search radius, or by selecting only a given number of its closest sample points. However, when the weights of far-away points are not negligible, this will cause artefacts in the final result. Another related issue is that kriging is often said to be exact in theory, ie it passes exactly through the sample points. However, this is only true if the nugget of the theoretical variogram function is set to zero. Some authors and implementations get around this by hard-coding 𝛾(0)=0 , but this can create a discontinuity between the value at the sample point and its immediate neighbourhood. Finally, it is worth noting that kriging can be directly applied to any point on the plane, yielding a result such as the one in Figure 6.5. However, much like other interpolation methods, kriging is only reliable in the domain (ie roughly the convex hull of the points). It can extrapolate (often by using negative weights), but that does not mean that the results outside the domain are accurate. Figure 6.5: The result of using ordinary kriging to interpolate on a grid of points using the sample dataset using only the sample points within 15 units of each interpolated point. 0 50 100 150 200 250 x 0 50 100 150 200 250 y 6.8 Notes and comments Krige (1951) is the original publication by Danie Krige, which was later formalised by Georges Matheron (Matheron, 1962; Matheron, 1965). How this came to be is best explained in Cressie (1993). 6.9 Exercises 75 https://pyinterpolate.readthedoc s.io/ If you have trouble following the derivations of the kriging equations or want to know more about them, Lichtenstern (2013) explains this well. If you feel like your statistics background is a bit weak, you first might want to have a look at Fewster (2014), particularly Chapter 3. A relatively simple explanation of kriging with agricultural examples is given by Oliver and Webster (2015). A standard reference textbook that is good but not so easy to follow is Wackernagel (2003). The mathematics covered in this chapter is partly based on the latter. Pyinterpolate is a good Python library to perform kriging and is used to generate some of the example figures from this chapter. Two other good Youtube videos that explain kriging: ▶https://www.youtube.com/watch?v=CVkmuwF8cJ8 ▶https://www.youtube.com/watch?v=98zz25kTteQ 6.9 Exercises 1. Why can using a search radius create artefacts in the interpolated terrain? 2. If kriging generally provides better results than other interpolation methods, why would you use something else (eg IDW)? 3. What does a nugget of zero say about a dataset? What about a large nugget? 4. What kind of dataset would yield a flat variogram (ie a horizontal line)? Conversions between terrain representations 7 uyoutu.be/Nqfp94KpcUA 7.1 Conversion to raster . . . . 77 7.2 Conversion to isolines . . . 78 7.3 Simplification of a TIN . . 81 7.4 Wedding cake effect . . . . 84 7.5 Notes and comments . . . . 85 7.6 Exercises . . . . . . . . . . . 85 We consider in this chapter the following four terrain representations and discuss the conversions between them: raster TIN isolines point cloud (PC) from/to PC raster TIN isolines PC — interpolate at middle points of cells (§7.1) create DT using 2D projection of points (ie using 𝑥 and 𝑦 only) convert to TIN + extract from triangles (§7.2.2) + structure output (§7.2.3) raster keep middle points only — create TIN using middle points of cells + TIN simplification (§7.3) extract from grid cells (§7.2.1) + structure output (§7.2.3) TIN keep only vertices interpolate at middle points of cells (§7.1) — extract from triangles (§7.2.2) + structure output (§7.2.3) isolines keep only vertices —warning: ‘wedding cake’ effect (§7.4) convert lines to points + interpolate (§7.1)—warning: ‘wedding cake’ effect (§7.4) create DT using points — warning: ‘wedding cake’ effect (§7.4) — 7.1 Conversion of PC/TIN to raster As shown in Figure 7.1, this step is trivial: one needs to interpolate at the locations of the centre points of the raster cells. The interpolation method can be any of the ones described in Chapters 5and 6. (a) (b) (c) (d) (e) Figure 7.1: (a) input sample points. (b) size/location of output raster. (c) 9 interpolations must be performed (at locations marked with ◦ ): at the middle of each cell. (d) the convex hull of the sample points show that 2 estimations are outside, thus no interpolation. (e) the resulting raster. 78 7 Conversions between terrain representations isoline components 7.2 Conversion to isolines Reading a contour map requires some skill, however it is considerably easier to learn to interpret a contour map than to manually draw one from a limited set of sample points. Yet this was exactly the task of many cartographers in the past couple of centuries: it was intuitively done by imagining a local triangulation of sample points. Isolines are usually directly extracted from either a TIN or a grid representation of a terrain. The basic idea, as shown in Figure 7.2, is to compute the intersection between a level value (eg 200m ) and each cell of the terrain (triangle or grid cell in our case). Notice that the cells are ‘lifted’ to their elevation. Each cell of the terrain is thus visited, one after the other, and for each cell if there is an intersection (which forms a line segment) then it is extracted. The resulting set of segment lines forms an approximation of the isoline. This process is then repeated for every level value. Notice that an isoline can have several components, for instance when the terrain has more than one peak. Therefore the number and size of the line segments in the resulting isoline are dependent on the resolution of the data representation. The basic algorithm for extracting one isoline is shown in Algorithm 3. Note that since the algorithm visits every grid cell or triangle individually and requires only local information, it is very easy to parallelise. It is thus a scalable algorithm. Its time complexity is O(𝑐) , where 𝑐 is the number of cells. Recall from Chapter 4that for 𝑛 points a DT contains about 2𝑛 triangles. The same idea can be used to extract all the isolines: for each triangle/cell and each level value, extract all the necessary line segments. 7.2.1 Conversion of raster to isolines Observe that, for a raster, the dual of the raster must be constructed (see Section 4.3), that is we consider the centre each pixel as a sample point, and we join with an edge the centres of two adjacent pixels (in Figure 7.3, the four values are centres of 4 adjacent raster cells). Intersections are computed by linearly interpolating the elevations of the vertex pairs along the edges of this grid. Figure 7.3 illustrates the different possible configurations. The top-left case indicates the case for which there are no intersections: all vertices are either higher or lower than 𝑧0. Observe that when two vertices are exactly at 𝑧0 , then the extraction of these is in theory not necessary because the neighbouring cell could also extract them. However, we do not want to obtain an output with duplicate line segments, and thus a simple solution to this is to only Figure 7.2: Vertical cross-section of a terrain (left), and a 2D projection of the terrain TIN with the extracted 200m isoline (right). 200m the 200m isoline has 2 components 7.2 Conversion to isolines 79 Algorithm 3: Simple extraction of one isoline Input: a planar partition 𝐸formed of cells (either rectangular or triangular cells); the elevation value 𝑧0 Output: a list of unstructured line segments representing the contour lines at 𝑧0 1𝑠𝑒𝑔𝑚𝑒𝑛𝑡𝐿𝑖𝑠𝑡 ← [ ]; 2for 𝑒∈𝐸do 3if 𝑧0intersects 𝑒then /*See Figures 7.3 and 7.4 */ 4extract intersection 𝜒of 𝑧0with 𝑒; 5add 𝜒to 𝑠𝑒𝑔𝑚𝑒𝑛𝑡𝐿𝑖𝑠𝑡; 15 17 19 21 10 10 10 10 10 17 19 10 12 8 9 15 12 8 6 8 10 17 7 8 10 8 10 15 12 8 14 7 12 8 14 7 12 8 14 7 Figure 7.3: Different cases when extracting an isoline at elevation 10m (in blue) for a regular grid. The blue lines are the ones extracted for that cell. extract such line segments if they are for instance the lower and/or left segments of a given cell. The most interesting case is the bottom-left one in Figure 7.3, it occurs when the two pairs of opposing points are respectively higher and lower than 𝑧0 . This forms a saddle point. The ambiguity arises here since there are two ways to extract a valid pair of contour line segments (only one of the 2 options must be extracted). This can be resolved by simply picking a random option or consistently choose one geometric orientation. 7.2.2 Conversion TIN to isolines Since a triangle has one fewer vertices than a square grid cell, there are less possible intersection cases and, more importantly, there is no ambiguous case. When one or more vertices of the triangle are at the same elevation as 𝑧0 , then one must be careful. As shown in Figure 7.4, if only one vertex is at 𝑧0 then nothing should be extract; if two vertices are at 𝑧0 then the edge between these can be extracted; if all three vertices are at 𝑧0 then the triangle is flat/horizontal and nothing should be extracted (because adjacent triangles will have edges extracted). To avoid extracting twice the same line segment when two vertices are at 𝑧0 (case on the right in Figure 7.4), then we can simply look at the normal 86 7 Conversions between terrain representations 4. Given a raster terrain (GeoTiff format) that contains several cells with no_data values, describe the methodology you would use to extract contour lines from it. As a reminder, contours lines should be closed curves, except at the boundary of the dataset. 5. Assume you have the small terrain formed of 3 triangles below, draw the isoline in this TIN for an elevation of 10 m. 8m 12m 9m 9m 15m 6m tinput x toutput 3×3 filter Figure 8.1: Example of a 3 × 3 filter. The new value of the cell 𝑥 of the input (in blue) is calculated by using its 8 neighbours (eg averaging the values) and the output terrain 𝑡𝑜𝑢𝑡𝑝𝑢𝑡 contains that value for its cell 𝑥 . This operation is usually performed for all cells in the input 𝑡𝑖𝑛𝑝𝑢𝑡 . Topographic properties and features 8 uyoutu.be/iF64Qeb_isw 8.1 Topographic properties . . 87 8.2 Topographic features . . . . 92 8.3 Properties used in practice 94 8.4 Notes and comments . . . . 96 8.5 Exercises . . . . . . . . . . . 96 While a terrain is a (2.5D) surface, it can also be conceptualised as an aggregation of many topographic features that are inter-related. Common examples of features are peaks, ridges, valleys, lakes, cliffs, etc., but one can think of application-specific ones such as the navigational channels in bathymetry, buildings in city modelling, or dikes for flood modelling. Identifying the different features forming a terrain enhances our understanding of the raw dataset. To help us extract and identify features, some properties (or characteristics) need to be extracted from terrains. We describe in this chapter the main properties of terrains, explain how they can be extracted, and how they are used in practice (for specific use-cases in different fields, and to identify features). 8.1 Topographic properties We describe in this section the main topographic properties that are commonly extracted from terrains: ▶slope ▶curvature ▶roughness Since these differ significantly from the data model used (TINs and grids), we give examples for both. For grids. The extraction of these properties is more common with grids, and in most cases kernel filters are used, Figure 8.1 shows one example. A filter is conceptually the same as a focal operation in map algebra (in GIS), or as a convolutional filter in computer vision. Observe that in Figure 8.1 a 3 × 3 window is shown, but a 5 × 5or larger window could also be used; the size to use depends on the scale at which one wants to extract a property. The main advantage of grids is that the value of a property can be easily calculated and the output is a new grid having the same resolution and orientation, but the values are the property (eg the gradient at this location) instead of the elevation. For TINs. A given property is in most cases as easy as for a grid to calculate for a given location ( 𝑥, 𝑦 ), but the issue is how to store the results: in a grid? only a set of points? This makes the use of TINs more cumbersome in practice. 88 8 Topographic properties and features Figure 8.2: The slope at a given location 𝑝𝑖 is defined by the tangent plane 𝐻𝑖 to the surface. Here are 3 examples for a profile view of a terrain. z=0 H1 H2 H3 p1 p2 p3 α α α Figure 8.3: One DTM with contour lines, and the gradient and aspect for a given location (blue cross). ∆z ∆xy α gradient aspect N E W S θ slope = gradient + aspect cartographical azimuth 8.1.1 Slope The slope at a given location 𝑝 on a terrain is defined by the plane 𝐻 that is tangent at 𝑝 to the surface representing the terrain (see Figure 8.2). What we casually refer to as ‘slope’ has actually two components (see Figure 8.3). 1. gradient 2. aspect Gradient. The gradient at a given point 𝑝 is the maximum rate of change in elevation. It is obtained by the angle 𝛼 between 𝐻 and the horizontal plane (Figure 8.3). From a mathematical point-of-view, the gradient is the maximum value of the derivative at a point on the surface of the terrain (maximised over the direction). The gradient will most often be expressed in degrees or in percentage. If a percentage is used, the following is used (see Figure 8.3): 𝑝𝑒𝑟𝑐𝑒𝑛𝑡𝑎𝑔𝑒 =Δ𝑧 Δ𝑥𝑦 (8.1) which means that a hill with a gradient of 100% is equal to 𝛼=45◦. Notice that if we calculate the gradient at every location for a terrain, then we obtain a new field since the gradient is a continuous phenomena (values from 0◦ to 30◦ for instance). This means in practice that for a given terrain in raster, calculating its gradient will create a new raster file that can be further processed. Aspect. At a given point 𝑝 on the terrain the gradient can be in any direction, the aspect is this direction projected to the 𝑥𝑦 -plane. It is basically a 2D vector telling us the direction of the steepest slope at a given point; at a given location the aspect will always be perpendicular to the contour line. Observe that for the parts of the terrain that are horizontal (eg a lake) the value of the aspect is unknown. The aspect is usually expressed as a cartographical azimuth, which is 8.1 Topographic properties 89 ci,j ci−1,j ci−1,j+1 ci,j+1 ci+1,j+1 ci+1,j ci+1,j−1 ci,j−1 ci−1,j−1 113 112 113 108 109 110 105 104 108 #2 #3 r= 5m Figure 8.4: (top) Given a cell 𝑐𝑖,𝑗 , the 3 × 3 kernel and its 8 neighbours. (bottom) A hypothetical case with some elevations; orange = aspect for method #2 below, purple = aspect for method #3 below. expressed in degrees, from the North, clockwise: 0◦ means North, 90◦ East, 180◦South, and 270◦West. Slope in TINs. Calculating the slope in a TIN is fairly straightforward: for a point 𝑝=(𝑥, 𝑦) find the triangle 𝜏 containing this point, and compute the normal vector ® 𝑛 of 𝜏 (pointing outwards). The projection of ® 𝑛 on the 𝑥𝑦 -plane is the aspect (this is done by simply ignore the 𝑧 - component of the vector). And the gradient is obtained by calculating the angle 𝛾 between ® 𝑛 and the horizontal plane, and taking the complement of 𝛾. If 𝑝 is directly on a edge of the TIN then the solution cannot be obtained directly; it is common practice to calculate the normal vector of the 2 incident triangle and average them to obtain one ® 𝑛 . The same is applied if 𝑝 is directly on a vertex 𝑣 of the TIN: the average of all the normal vectors of all the incident triangle to 𝑣is used. Slope in grids. If the terrain is represented as a regular grid (say of resolution 𝑟 ), then there exist several algorithm to obtain the slope at a given cell 𝑐𝑖,𝑗 . We list here a few common ones. It should be noticed that most algorithms use a 3 × 3 kernel, ie the value for the gradient/aspect at cell 𝑐𝑖,𝑗 is computed by using (a subset of) the 8 neighbours. 1. Local triangulation + TIN method. It is possible to locally triangulate the 9 points, calculate the normal of the 8 triangles, and then use the method above for TINs. 2. Maximum height difference. This method simply picks the maximum height difference between 𝑐𝑖,𝑗 and each of its 8 neighbours, the maximum absolute value is the direction of the aspect and the gradient can be trivially calculated. Notice that this means that there are only 8 possibilities for the slope (at 45°intervals). For the case in Figure 8.4, the aspect would be facing south (180°) and the gradient would be 45°. 3. Finite difference. With this method, the height differences in the 𝑥 -direction (west-east) and in the 𝑦 -direction (south-north) are calculated separately, and then the 2 differences are combined to obtain the slope. This means that only the direct 4-neighbours of 𝑐𝑖,𝑗 are used. 𝜕𝑧 𝜕𝑥=𝑧𝑖−1,𝑗 −𝑧𝑖+1,𝑗 2𝑟,𝜕𝑧 𝜕𝑦=𝑧𝑖,𝑗−1−𝑧𝑖,𝑗+1 2𝑟(8.2) The gradient is defined as: tan 𝛼=s(𝜕𝑧 𝜕𝑥)2+(𝜕𝑧 𝜕𝑦)2(8.3) and the aspect as: tan 𝜃= 𝜕𝑧 𝜕𝑦 𝜕𝑧 𝜕𝑥 (8.4) 90 8 Topographic properties and features 2How does it work in practice? The GDAL utility gdaldem ( https://www.gdal.org/gdaldem.html ) does not have the best documentation and does not explicitly mention which method is used. After some searching, we can conclude that the method “4. Local polynomialfitting” is usedby default forslope/aspect,and specifically the Horn’s method is used (Horn, 1981). This uses a 3 × 3 window, and fits a polynomial; the centre pixel value is not used. If the option -alg ZevenbergenThorne is used, then the algorithm of Zevenbergen and Thorne (1987) is used. This uses only the 4 neighbours, and is a variation of the method “3. Finite difference” above. The documentation of gdaldem states that: “literature suggests Zevenbergen & Thorne to be more suited to smooth landscapes, whereas Horn’s formula to perform better on rougher terrain.” 2nd derivative of the surface The value of 𝜃 should be resolved for the correct trigonometric quadrant, and if 𝜕𝑧 𝜕𝑥=0 then it means the aspect should be handled differently (considering only the variation in the south-north direction). For the case in Figure 8.4, the gradient would be 39.5°and the aspect would be 194.0°. 4. Local polynomial fitting. Based on the9 elevation points, it is possible to fit a polynomial (as explained in Chapter 5) that approximate the surface locally; notice that the polynomial might not pass through the point if a low-degree function is used. A quadratic polynomial could for instance be defined: 𝑓(𝑥, 𝑦)=𝑎𝑥2+𝑏𝑦2+𝑐𝑥𝑦 +𝑑𝑥 +𝑒𝑦 +𝑑(8.5) , and thus: 𝜕𝑓 𝜕𝑥=2𝑎𝑥 +𝑐𝑦 +𝑑(8.6) 𝜕𝑓 𝜕𝑦=2𝑏𝑦 +𝑐𝑥 +𝑒(8.7) and if a local coordinate system centered at 𝑐𝑖,𝑗 is used, then 𝑥=𝑦=0 , and thus 𝜕𝑓 𝜕𝑥=𝑑and 𝜕𝑓 𝜕𝑦=𝑒. 8.1.2 Curvature The curvature is the 2nd derivative of the surface representing the terrain, it represents the rate of change of the gradient. We are often not interested in the value of the curvature itself ( ° 𝑚 ) but whether the curvature is: convex, concave, or flat. The curvature at a point 𝑝is often decomposed into types: 1. profile curvature: the curvature ofthevertical cross-sectionthrough 𝑝 perpendicular to the contour line passing through 𝑝 (or of the vertical plane along the 2D vector of the aspect at 𝑝) 8.1 Topographic properties 91 Figure 8.5: Nine curvatures (Figure adapted from Kreveld (1997)). 2. plan curvature: the curvature along the contour line passing through 𝑝 (or along the line segment perpendicular to the 2D vector aspect and passing through 𝑝) Because there are 2 types of curvatures and each have 3 potential values, there are 9 possible options (as Figure 8.5 shows). Computing for grids. Computing the curvature is a complex operation and we will not describe one specific method. The idea is to reconstruct locally the surface (eg with the polynomial fitting from Section 8.1.1 above, or with a TIN), and then verify whether the 2 curvature types are convex/concave/flat. Observe that the curvature, as it is the case for the slope, is heavily influenced by the scale of the terrain (its resolution) and thus having a 3 × 3 kernel might be influenced by the noise in the data, or by small features. Computing for TINs. For a TIN, it is possible to define for each vertex 𝑣 the profile and the plan curvatures by using the triangles that are incident to 𝑣 and extract the contour line for the elevation of 𝑣 (as is shown in Figure 8.9). The idea is to classify each vertex into one of the 9 possibilities in Figure 8.5. If there is no contour segment, then 𝑣 is either a peak or a pit. A peak will be profile and plan convex; a pit will be profile and plan concave. If there are 2 segments, then we can use these to estimate the direction of the aspect, it will be perpendicular (thus the bisector between the 2 segments is a good estimate) in the direction of lower elevations. If we simply look at the elevations higher and lower than 𝑣 along this direction, then we can easily verify whether 𝑣 is profile convex or concave. For the plan curvature, we can simply walk along one of the 2 edges so that higher elevations are on our left, 𝑣 is plan convex if the contour line makes a left turn at 𝑣 , if it makes a right turn it is concave, and if it is straight then it is plan flat. If there are >2 segments, then 𝑣 is a saddle point and thus no curvatures can be defined. 92 8 Topographic properties and features (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (f,f) (f,f) (v,c) (v,v) (v,v) (f,f) (v,v) (v,v) (v,v) (v,f) (v,f) (v,f) (v,f) (v,f) (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (v,c) (f,f) (f,f) (v,c) (v,v) (v,v) (f,f) (v,v) (v,v) (v,v) (v,f) (v,f) (v,f) (v,f) (v,f) (v,c) (v,v) (v,f) (v,f) (v,v) Figure 8.6: (top) Points from a TIN classified according to their curvatures (convex, concave, flat). (middle) The VD of the points. (bottom) The Voronoi edges between the cells having the same label are removed, to create polygons. Figure 8.7: The green profile of a terrain has a lower roughness than the orange one (normals locally deviate less). When each point has been assigned a curvature—a pair (𝑝𝑟𝑜 𝑓 𝑖𝑙𝑒, 𝑝𝑙𝑎𝑛) — we can use for instance the Voronoi diagram, as shown in Figure 8.6. It suffices to remove the Voronoi edges incident to cells having the same label, and polygonal zones are obtained. 8.1.3 Roughness & ruggedness The terms ‘roughness’ and ‘ruggedness’ are often used interchangeably and have slightly different definitions depending on the software and/or the documentation. We can however claim that they both refer to how “undulating” or “regular” a (part of a) terrain is. A terrain with a high roughness will have small local deviations, while one with low roughness will be “smoother”. In other words, the normals of the surface of the terrain will deviate from each others greatly for a high roughness, and less for low roughness. Another way to measure roughness, is to think of it as the ratio between the surface area and its projection into a plane. For gridded terrains, the roughness is most often calculated by simply looking at the differences in elevations for the cells inside a 3 × 3 filter (or larger filter). The roughness value for one location can be one of these variations (more exist): ▶ the standard deviation of the 9 values (or 25 if a 5 × 5 filter is used) in the filter; ▶ the largest difference in elevation between the value in the centre of the filter and one its neighbouring cell in the filter; ▶ the difference between the elevation of the central pixel and the mean of its surrounding cells; It should be observed that these methods are highly influenced by the scale, ie the resolution of the grid and the size of the filter will yield potentially very different results. Notice also that since the differences in the elevations are used, a terrain that would be a constant slope (eg a talus) would get a roughness that is not zero (in comparison to a perfectly flat terrain). A solution to this would be to fit a plane with least-square adjustment to all the points involved in a filter, and then compare the differences of the elevations to plane. For TINs, the same three variations above can be used for a single location ( 𝑥, 𝑦 ), if for instance we pick the natural neighbours or if all the points within a certain distance threshold are used. However, as mentioned above, how to create a new field of roughness is not as trivial as for a grid. One could recreate a TIN with the values of the vertices having the roughness, or create a grid. 8.2 Topographic features 8.2.1 Peak A point 𝑝 whose surrounding is formed only of points that are of lower elevation is a peak. The size and shape of the surrounding is dependent on the application and on the data model used to represent the terrain. If a grid is used, this surrounding could be the 8 neighbours; if a TIN is used they could be the vertices that of the triangles incident to 𝑝 . Observe that 8.2 Topographic features 93 Figure 8.8: (a) Peaks and pits. (b) A saddle point (Figure from https://www.armystudyguide.com) peaks and pits are local and influenced by the scale of the data 8 9 14 17 5 5 6 15 18 12 10 Figure 8.9: A saddle point at elevation 10m , and its surrounding points. The triangulation of the area is created and used to extract the contour line segments at 10m (red lines). a peak can be local, that is one point that happens to be a few centimetres higher than all its neighbours would be classified as a peak (the small terrain in Figure 8.8 contains several peaks), while if we consider a hill we would surely consider only the top as the peak. A peak is therefore on the scale of the data. The contour line through the 𝑝does not exist. 8.2.2 Pit A point 𝑝 whose surrounding is formed only of points that are of higher elevation is a pit. The same remarks as for peak apply here. The contour line through the 𝑝does not exist. 8.2.3 Saddle point As shown in Figure 8.8b, a saddle point, also called a pass, is a point whose neighbourhood is composed of higher elevations on two opposite directions, and 2 lower elevations in the other two directions. From a mathematics point-of-view, it is a point for which the derivatives in orthogonal directions are 0, but the point is not maximum (peak) or a minimum (pit). If we consider the contour line of a saddle point 𝑝 , then there are 4 or more contour line segments meeting at 𝑝 ; for most point in a terrain this will be 2, except for peaks/pits where this is 0. Figure 8.9 shows an example for a point with an elevation of 10m , the contour lines at 10m is drawn by linearly interpolating along the edges of the TIN of the surrounding (see Chapter 7). 8.2.4 Valleys & ridges Valleys and ridges are 1-dimensional features. If a terrain is represented as a TIN, we can extract the edges of the triangles that form a valley or a ridge. An edge 𝑒 , incident to 2 triangles, is considered a valley-edge if the 94 8 Topographic properties and features Figure 8.10: Edges in a TIN can be classified as valley, ridge, or neither projection of the 2 normals of the triangles, projected to the 𝑥𝑦 -plane, point to 𝑒 . If the 2 normals projected point in the opposite direction, then 𝑒is a ridge. If they point in different directions, then 𝑒is neither. 8.3 Properties and features used in practice 8.3.1 Slope The slope (gradient + aspect) are a cornerstone of runoff modelling (see Chapter 10), the prediction of the flow and accumulation of water on a terrain. The slope is used to calculate the flow direction at a given location, which is the direction with the steepest descent at that location. The slope can also be used to predict the irradiation (from the sun) that a given location at a given day/time would receive. This is often the input of (local) meteorological models, can be used to optimise the location of solar panels or to predict land surface temperature. 8.3.2 Curvature While curvature is used implicitly to calculate the flow direction in runoff modelling, there are use-cases where the value is useful. One of them is for the predicting of where snow covering will be. The values of the curvature can help a practitioners understand and characterise the drainage basins, once extracted from a terrain (see Section 10.5). 8.3.3 Roughness The roughness can be used directly as a predictor for the habitats of different species. The variations in roughness in a terrain can be used to delineate the terrain into geomorphological and geological areas. 8.3 Properties and features used in practice 95 Figure 8.11: Left: a DTM visualised with height as a shade of blue. Right: when hillshading is applied. 2Why does the sunlight come from the North-West? The source of the light for hillshading is usually set at the NorthWest, but in reality the sun is never located there (in the northern hemisphere). Why is this a common practice then? The main reason is because the human brain usually assume that the light comes from above when looking at picture. Doing so reduces the chances of relief inversion, ie when mountains are perceived as valleys, and vice-versa. This website  gives a clear example where a valley is interpreted as a mountain ridge by many if the sun is coming from the South. 8.3.4 Hillshading Hillshading is a technique used to help visualise the relief of a gridded terrain (see Figure 8.11 for an example). It involves creating an image that depicts the relative slopes and highlights features such as ridges and valleys; a hillshade does not depict absolute elevation. This image assumes that the source of light (the sun) is located at a given position (usually North-West). While it would be possible to use advanced computer graphics methods (see Chapter 9) to compute the shadows created by the terrain surface, in practice most GIS implements a simplified version of it which can be computed very quickly. Given a regular gridded terrain, hillshading means that each cell gets a value which depicts the variation in tone, from light to dark. The output of a hillshade operation is thus a regular gridded DTM, usually with the same extent and resolution as the original grid (for convenience). The values computed for each cell need as input the gradient and the aspect of the terrain. The formula to compute the hillshade of a given cell 𝑐𝑖𝑗 differs from software to software, and we present here one (it is used in QGIS and ArcGIS for example, and surely others). It assumes that the output hillshade value is an integer in the range [0,255] (8-bit pixel), and that the direction (azimuth) and the height (given as an angle) of the illumination source is known. Notice that the position of the sun is relative to the cell, its position thus changes for different cells of a terrain. As above and in Figure 8.12, for a cell 𝑐𝑖𝑗 , its gradient is 𝛼𝑖𝑗 , its aspect is 𝜃𝑖𝑗 , the azimuth of the sun is 𝜓 (angle clockwise from the north, like the aspect), and the height of the sun is 𝛾 (0 rad is the horizon, 𝜋 2 rad is the 102 9 Visibility queries on terrains v τ1 τ2 τ3 Figure 9.8: The 3 triangles 𝜏1 , 𝜏2 , and 𝜏3 form a cycle when viewed from the viewpoint 𝑣 , and it is not possible to sort them from furthest to closest (without decomposing them). τ3 τ1 τ2 τ4 τ5 τ6 v Figure 9.9: The triangles in a DT can be ordered in an in-front/behind manner when viewed from a viewpoint. One solution involves using so-called tangents. The current tangent 𝑡𝑐𝑢𝑟 is first initialised as a vector pointing downwards. Then, starting at 𝑣 , we walk along the ray 𝑣𝑞𝑖 , and for each cell intersected its elevation 𝑧 is compared to the elevation of 𝑡𝑐𝑢𝑟 at that location. If 𝑧 is lower, then the cell is invisible. If 𝑧 is higher, then the cell is visible and 𝑡𝑐𝑢𝑟 is updated with a new tangent using the current elevation. Viewsheds with several viewpoints 𝑣𝑖 are also very useful, think for instance of obtaining the viewshed along a road. This can be computed by sampling the road at every 50m and computing the viewsheds from each of the points. Each viewshed yields a binary grid, and it suffices to use a map algebra operator to combine the results into one grid (if one cell is visible from any viewpoint, then it is visible). 9.2.2 Visibility in TINs Using the depth-sort algorithm for arbitrary triangles would require using a BSP-tree for indexing and sorting the triangles, and some triangles would need to be decomposed, as explained above. Figure 9.8 shows one simple example. However, it has been proven that Delaunay triangulations are acyclic for any fixed viewpoint. In other words, the in-front/behind relationship for the triangles of a DT, with respect to a given viewpoint, is acyclic (see Figure 9.9). Therefore, to obtain the triangles intersecting a ray coming out of a viewpoint (ordered from the closest to farthest), it suffices to modify slight the point location algorithm from Section 4.4.2. This operation can be performed in 2D, by projecting the triangles of the TIN to the 𝑥𝑦-plane. This means that visibility queries in TINs—like in grids—are greatly simplified compared to the general case where the ordering of objects is the main difficulty (and handling overlapping objects like in Figure 9.5). 9.3 Notes and comments The ‘tangent algorithm’ to compute viewsheds was first described by Blelloch (1990). The description here is inspired by that of De Floriani and Magillo (1999). Newell et al. (1972) first proposed the depth-sorting algorithm and the decomposition necessary when polygons in the scene cannot be sorted from furthest to closest. Edelsbrunner (1990) proved that Delaunay triangulations, in any dimensions, are acyclic. 9.4 Exercises 103 9.4 Exercises 1. Explain why the spacing in Figure 9.6c along the profile has points that are not equally spaced. 2. You are given a 2.75D terrain of an area, it is composed of triangles, and your aim is to perform line-of-sight queries between some locations. Describe the algorithm that you will implement to perform the queries. flow direction flow accumulation Runoff modelling 10 uyoutu.be/0NzZoJATFjc 10.1 Flow direction . . . . . . 106 10.2 Flow accumulation . . . 107 10.3 Solving issues with sinks 108 10.4 Flow direction in flats . . 109 10.5 Drainage networks . . . 110 10.6 Notes and comments . . 111 10.7 Exercises . . . . . . . . . . 111 Many interesting DTM operations are based on runoff modelling, ie the computation of the flow and accumulation of water on a terrain. Examples include: knowing where streams will form in the case of heavy rainfall, finding the areas that will be affected by a waterborne pollutant, tracing the areas that could become submerged by floodwater, or calculating the rate of erosion or sedimentation in a given area. In hydrology, runoff modelling can be very complex (Figure 10.1). Hydrological models usually consider different precipitation scenarios, model various types of overland and subsurface flows, and take into account many locationand time-dependent factors, such as the depth of the water table and the permeability of the soil. Such models can be quite accurate, but they require high-resolution data that is often not available, they are difficult to create without specialised knowledge, and they involve substantial manual work. By contrast, the simpler GIS models of runoff can be performed automatically in large areas with only a DTM. These models mostly use gridded raster terrains, and so we will generally refer to these in this chapter, but the methods described here mostly work just as well with other representations (eg a TIN). In order for the GIS models of runoff to achieve their results, two big assumptions are usually made: 1. that all water flow is overland, thus ignoring all subsurface flows and dismissing factors such as evaporation and infiltration; and 2. that a good estimate for the total flow at any point is the drainage area upstream from it, ie the area above the point which drains through/to it, which is roughly equivalent to rain that is falling evenly all over a terrain. Basedontheseassumptions,runoff modellingissimplifiedbyconsidering only two values, which are computed for every cell in a DTM: Flow direction Given a DTM cell, towards which nearby cells and in which proportions does water flow from it? Flow accumulation Given a DTM cell, what is the total water flow that passes through it? We look at a few different methods to compute these values in the next two sections. 106 10 Runoff modelling Figure 10.1: Different types of water flows as modelled in hydrology. Figure adapted from Beven (2012). single flow direction (SFD) D8 flow direction rho8 (𝜌8) 10.1 Computing the flow direction Theoretically, the flow direction of a point is the direction with the steepest descent at that location, which does correspond to the direction towards which water would naturally flow. However, the discretisation of a terrain into DTM cells means that some kind of an approximation needs to be made. There are two broad approaches that can be followed to do this: computing a single flow direction, which assumes that all the water in a DTM cell flows to one other cell, or multiple flow directions, which assumes that the water in a DTM cell can flow towards multiple other cells. 10.1.1 Single flow direction The earliest and simplest method to compute the flow direction of a cell is to compute the slope between the centre of the cell and the centre of all its neighbouring cells (using the distance between the centres and the difference in elevation), then assign the flow direction towards the neighbour with the steepest descent. The method is known as the single flow direction (SFD) approach, and when applied to a raster grid, it usually considers that there are eight neighbours to each pixel (left, right, up, down and the diagonals). For this reason, it is also known as the eight flow directions (D8) approach. On one hand, the method is very fast and easy to implement, and it avoids dispersing the water flow between multiple cells. On the other hand, it can have significant errors in the flow direction, and it does not allow for divergent flows. For instance, in a square grid, the errors can be of up to 22.5◦ (because the method is forced to choose a neighbouring cell in increments of 45◦ ). This method can therefore easily create artefacts in certain geometric configurations (Figure 10.2). Many of these artefacts can be eliminated by using the rho8 ( 𝜌8 ) method, which modifies D8 to assign the flow direction to one of its lower neighbours randomly with probability proportional to the slope. However, it produces non-deterministic results, which is often a sufficient reason not to use it. 10.2 Computing the flow accumulation 107 312 TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS A. Theoretical. C. MS. B. D8. •" ....... .?'•.'• ::• • •i:•:•;• ........ .... •:.•.......:.• .......................... .• ...................... .• .... .... . ............ • ...... ,•• •,, ........... :::::::::::::::::::::: ..... '" '•;- •3' • .... '•-2{ '"?'":":":' '":"'•: D. Lea's method E. DEMON F. Doo Figure 4. Top half of an outward draining circular cone. Elevation was defined as 200 minus the radius from the center on a 16 x 16 grid with grid spacing 10 units. Specific catchment area is theoretically radius/20 ranging from 0 at the center to 53 on the corners. Contours (10-unit interval) depict elevation. Gray scale depicts contributing area (1 white to 60 black). (a) Theoretical values, (b) single direction (D8) procedure, (c) Quinn et al.'s [1991] procedure (MS), (d) Lea's [1992] method, (e) DEMON [Costa-Cabral and Burges, 1994], and (f) new procedure (D•). and e i and d• are elevations and distances between pixels as labeled in Figure 3. The slope direction and magnitude are r = tan -• (s2/s0 S = (S12 -'1 c S22) 1/2 (3) If r is not in the range (0, tan -• (d2/d•)), then r needs to be set as the direction along the appropriate edge and s assigned as the slope along that edge. if r < 0, r = 0, s = s• (4) if r > tan -• (d2/di), r = tan -• (d2/di), s = (eoe2)/(d• 2 + d22) 1/2 (5) Next recognize that each of the eight facets depicted in Figure 2 can be mapped by appropriate selection of corner elevations and rotation/transformation onto the facet in Figure 3. Table 1 gives the node elevations corresponding to the corners of each of the triangular facets used to calculate slopes and angles in (1)-(5). These are arranged such that e o is the center point, e • is the point to the side, and e 2 is the diagonal point. The local angle associated with the largest downwards slope from the eight facets (r' = r from facet with maximum s) is then adjusted to reflect an angle counterclockwise from east (Figure 2) to obtain the flow direction angle. ra = air' + a cVr/2 (6) The multiplier af and constant a c depend on the facet selected and are listed in Table 1. The procedure that searches for the facet with the largest slope proceeds in the order of facets 1 to 8 shown in Figure 1 and in the case of ties (facets with equal slope) picks the first one. In nature ties are extremely rare so the bias introduced by this is deemed negligible. In the case where no slope vectors are positive (downslope) a flow direction angle of -1 is used to flag the pixel as "unresolved," that is, a flat area or pit. Unresolved flow directions are resolved iteratively by making them flow toward a neighbor of equal elevation that has a flow direction resolved. This is the same approach for resolving pits and flats as used in the D8 method [e.g., Mark, 1988; Jenson and Domingue, 1988]. I therefore use the calculation of D8 flow directions as a preprocessor to raise the elevation of all pixels in a pit to the level of the overflow. Then where pixels are flagged as "unresolved," the flow angle returned by the D8 procedure is used. This ensures that flat pixels drain to a neighbor that ultimately drains to a lower elevation, eliminating the possibility of inconsistencies such as loops in the flow direction angles. This method of representing flow directions based on triangular facets is designated Dc• (an infinite number of possible single flow directions). Figure 10.2: The D8 method creates artefacts when water is draining from a circular cone. Figure from Tarborton (1997). multiple flow directions (MFD) Figure 10.3: The flowwidth 𝐿 canbe computed using the geometry of the DTM cells. In the case of a square grid with spacing 𝑑 , it is √2 4𝑑 for the diagonals ( 𝐿2 ) and 1 2𝑑 for the adjacencies ( 𝐿1 ), where 𝑑 is the grid spacing. Based on Quinn et al. (1991). Despite its age and limitations, the SFD method is still widely used and available in many GIS tools. 10.1.2 Multiple flow directions In an attempt to overcome the limitations of the SFD method, a variety of methods assign the flow direction of a DTM cell fractionally to some or all of its lower neighbouring cells according to some criteria. These methods are collectively known as multiple flow directions (MFD), and they usually use a variation of this equation: 𝐹𝑖=(𝐿𝑖tan 𝛼𝑖)𝑥 P𝑛 𝑗=1𝐿𝑗tan 𝛼𝑗𝑥(10.1) where 𝐹𝑖 is the flow towards the i-th neighbouring cell, 𝐿𝑖 is the flow width (Figure 10.3), 𝛼𝑖 is the gradient towards the i-th neighbouring cell (and so tan(𝛼𝑖) is the slope), 𝑥 is an exponent that controls the dispersion, and 𝑛is the number of neighbours of the cell. As shown in Figure 10.4, MFD methods show characteristically wider flows compared to SFD methods. D8 does not disperse the flow, but the path is constrained to the 8 possible grid directions. By contrast, Quinn et al. (1991) (an MFD method) can model the flow direction in a way that matches the topography better, but it also introduces substantial dispersion. More modern approaches try to combine some of the advantages of both approaches. TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS 313 A. Single direction procedure, D8 B Quinn et al. (1991) procedure, MS C. Lea's (1992) method D. DEMON ......... E. New Procedure, Doo Figure 5. Circular cone influence maps for the circled pixels. Gray scale ranges from white (0 or no influence) to black (1 or 100% influence). Calculation of Upslope Areas Upslope area is calculated using a rccursivc procedure that is an extension of the very efficient rccursivc algorithm for single directions [Mark, 1988]. The upslope area of each pixel is taken as its own area (one) plus the area of upslope neighbors that have some fraction draining to the pixel in question. The flow from each cell either all drains to one neighbor (if the angle falls along a cardinal or diagonal direction) or is on an angle falling between the direct angle to two adjacent neighbors. In the latter case the flow is proportioned between these two neighbor pixels according to how close the flow direction angle is to the direct angle to those pixels. The following pseudocode gives the logic of this algorithm: Procedure DPAREA(i, j) if AREA(i, j) is known then no action else AREA(/, j) = 1 (the area of a single pixel) for each neighbor (location in, jn) p = proportion of neighbor (in, jn) that drains to pixel (i, j) based on angle if (p > 0) then call DPAREA(in, jn) (this is the recursive call to calculate area for the neighbor) AREA(i, j) = AREA(i, j) + p x AREA(in, jn) Return The calculation is initiated by calling this function for the outlet pixel. It then recursively calls itself for all pixels that contribute to the upslope area at the outlet. The recursion stops when it reaches a pixel that has no pixels upslope. Illustrative Examples This section gives examples of results from this method, D% compared to the single direction approach, D8; Quinn et al.'s [1991] multidirection algorithm, MS; Lea's [1992] method; and DEMON [Costa-Cabral and Burges, 1994]. In these examples we use the notion of influence and dependence maps. The influence function I(x, Xo) is defined as the upslope area at each pixel x from a specific pixel Xo. It maps where flow from pixelxo goes and how it is dispersed. It is computed by running a modified version of the procedure for calculating upslope area that uses an area contribution of one from pixel x o but zero for all other pixels. The dependence function D(x, Xo) is the opposite of the influence function, defined as D(x, Xo) = I(xo, x). The upslope area at pixelx o is composed of the sum of the area of upslope pixels that have some proportion of their flow go through pixelxo. D(x, Xo) maps the contribution from pixel x to the calculation of upslope area at Xo. It is calculated through repeated evaluation of the influence function. Figure 4 shows the upslope area by each approach for a circular cone. Figure 5 shows the influence maps from each of the five algorithms (D8, MS, Lea's [1992] method, DEMON, and Din) applied to the circular cone. D8 results in no spreading, but flow paths (which are what the influence map plots) Figure 10.4: Flows in a circular cone: SFD (D8) vs. MFD (Quinn et al. (1991)). Figure from Tarborton (1997). 10.2 Computing the flow accumulation After the flow directions in all the cells of a DTM have been computed, the usual next step is to use this information to compute the flow 108 10 Runoff modelling sink outlet accumulation in all of them. As stated in the assumptions we make for GIS models of runoff, the flow accumulation at a given DTM cell can be estimated by the area that drains to it. Note that in the case of a square grid, it is simply the number of cells that drain to it. In practical terms, the flow accumulation is defined based on a recursive operation: 𝐴0=𝑎0+ 𝑛 X 𝑖=1 𝑝𝑖𝐴𝑖(10.2) where 𝐴0 is the accumulated flow for a cell, 𝑎0 is the area of the cell, 𝑝 is the proportion of the i-th neighbour that drains to the cell, 𝐴𝑖 is the accumulated flow for the i-th neighbour, 𝑛 is the total number of the neighbouring cells. Note that this calculation can be sped up substantially by: (i) storing the accumulated flows that have already been computed, and (ii) not following the recursion when 𝑝𝑖=0. 10.3 Solving issues with sinks Sinks, which are also known as depressions or pits, are areas in a DTM that are completely surrounded by higher terrain. Some of these are natural features that are present in the terrain (eg lakes and dry lakebeds) and where water would flow towards (and stagnate) in reality, and are thus not a problem for runoff modelling. However, they can also be artefacts of the DTM (eg noise and areas without vegetation can create depressions), or they can be very small areas that easily filled (ie flooded), after which water would flow out of them. In the latter case, we need to implement a mechanism to route water flows out of these depressions, since otherwise our runoff model could have very large water flows stopping at even tiny depressions. We will look at two common options to solve this problem: modifying a DTM by filling in (certain) sinks, and implementing a flow routing algorithm that allows water to flow out of sinks. 10.3.1 Filling in sinks The aim of the algorithms to fill in sinks is to increase the elevation of certain DTM cells in a way that ensures that all the cells in the DTM can drain to a cell on its boundary (Figure 10.5), or possibly to a set of cells that are known to be valid outlets, eg lakes and oceans. At the same time, the elevation increases should be minimised in order to preserve the original DTM as much as possible. In the best case scenario, we can imagine that the resulting DTM is one that resembles follows the terrain elevation of terrain where there is no water and the top of natural water bodies (but no artificial features as in a DEM). One efficient method to fill in sinks is the priority-flood algorithm (Barnes et al., 2014b). It works by keeping: (i) a list of DTM cells that are known to drain, which is kept sorted by elevation; and (ii) a raster marking whether each cell of the DTM is known to drain yet. The list is initialised with the cells on the boundary of the DTM (which are assumed to be able to drain outwards), as well as other specially marked cells (eg those forming part of a lake or a large river). Then, it iteratively: (i) removes the lowest cell from the sorted list of cells that are known to drain, (ii) increases the 10.4 Assigning flow direction in flats 109 been resolved. If an unresolved neighbor nis at a lower elevation than c(such as cells D;E;Fin Fig. 1), it is raised to the elevation of c before it is placed on the queue. In this way, an ϵ-descending path is constructed for each cell, leading to the fulfillment of the second depression-filling criterion. Because nis always brought up only to the elevation of the lowest cell which still drains to the edge of the DEM, the third depression-filling criterion is fulfilled. The algorithm terminates when the priority queue is empty (Fig. 1k). Since all cells are given elevations greater or equal to their original values in the DEM, the first property of the depressionfilling criteria is fulfilled. Because all three criteria are fulfilled at each step, it follows that the result of the Priority-Flood Algorithm is a solution to the depression-filling problem. Fig. 1. A conceptual cross-section of a DEM demonstrating the Priority-Flood Algorithm. The algorithm is initialized in (a), with all the edge cells having been pushed into a priority queue from which the cell with the lowest elevation is popped. In (b)–(e), cells are popped from the priority queue, marked as processed (dark lines on the landscape), and their neighbors pushed. In (f), cell Cis the lowest cell which is guaranteed to drain, therefore its neighbor Dmust be in a depression. To resolve the depression, cell Dis raised to the elevation of Cbefore being pushed into the priority queue. D's old elevation is no longer needed, but is shown with a dotted line for clarity. In (g), Emust also, by the same logic, be in a depression; therefore, it is raised to the new elevation of D. This continues through (i). In (j), the final cells are popped, though there is nothing to do. (k) shows the final elevations of the cells. Since the cells are raised as the algorithm progresses, no extra work is necessary between (j) and (k); the algorithm is complete when the priority queue is empty. R. Barnes et al. / Computers & Geosciences 62 (2014) 117–127120 Figure 10.5: A vertical cross-section of a DTM with filled sinks. The dashed line represents the original terrain, whereas the thicker solid line represents the filled terrain. From Barnes et al. (2014b). least-cost paths (LCP) flat elevation of its neighbours that are not yet known to drain to the level of the cell, (iii) adds the neighbours that are not yet known to drain to the list. Note that implicit in this last step is the fact that the neighbour cells are deemed to be able to drain through the current (lowest) cell. 10.3.2 Least-cost (drainage) paths An alternative to modifying a DTM to eliminate sinks is to implement a more complex water routing algorithm that allows water to flow out of sinks. For this, the usual approach is to implement a variation of the 𝐴∗ search algorithm, which in this context is known as the least-cost paths (LCP) algorithm (Metz et al., 2011). The LCP algorithm is similar to priority-flood in that it keeps a sorted list of DTM cells that are known to drain, which is also initialised to the boundary pixels (and possibly other cells). Then, it iteratively: (i) removes the lowest cell from this list, (ii) sets the drainage direction of its neighbours that are not yet known to drain towards itself, (iii) adds the neighbours that are not yet known to drain to the list. 10.4 Assigning flow direction in flats Flats are areas in a DTM that have the same elevation. They therefore do not have a well-defined flow direction, which causes problems for many water routing algorithms. Flats can sometimes occur naturally, but they are more often the result of precision limits, noise removal, or sink filling algorithms. It is thus often necessary to apply a method that assigns a flow direction to flats, either by: (i) modifying the DTM to eliminate them, and then assigning them a flow direction in the usual way, or (ii) assigning them a flow direction directly. After all flats in a DTM have been identified and their extent is known, algorithms usually work by (i) assigning an artificial gradient away from higher terrain (Figure 10.6), ie terrain in a flat is assumed to become lower as we move farther away from its neighbouring higher terrain; and/or (ii) assigning an artificial gradient towards lower terrain (Figure 10.7), ie terrain in a flat is assumed to become lower as we move closer to its neighbouring lower terrain. Barnes et al. (2014a) is a good example of an efficient method that combines both of these approaches, resulting in more natural flow directions and better results than would be possible with either approach individually. 110 10 Runoff modelling Figure 10.6: In a flat surrounded by higher terrain (dark grey) with a single lower-elevation outlet (light grey), we can use a gradient away from the higher terrain to route water out of the flat and towards the outlet. For this, we can iteratively assign (tiny or symbolic) elevation decreases in the flat starting from the higher terrain until all non-draining cells in the flat have been covered. Note that in this case, a sink is produced by the procedure. Figure from Barnes et al. (2014a). Since every cell surrounding a flat must be either a low edge or a high edge cell, if none of either are found, then the DEM has no flats. If there are no low edge cells, but there are some high edge cells, then the DEM has flats but they cannot be drained. In both such instances, the algorithm indicates this and proceeds no further. Otherwise, the algorithm then labels each unique flat (see Algorithm 4). Let the data structure holding these labels be called Labels, and let it be initialized to a value N O L ABEL . To apply the labels, each cell cin the LowEdges queue is used as the seed to a flood fill algorithm. If chas not yet been labeled then a new label l is given to it, its elevation eis noted, and all cells which can be reached from cwhile traversing only cells of elevation eare given the same label l. A FIFO queue could be used for this, or a D8 variant of Heckbert (1990), or one of many other solutions to the “connected components”problem (Grana et al., 2010 provide a good overview). Recursive solutions should be avoided in implementations using languages without proper tail recursion, as stack overflows are likely when processing large DEMs. In the end each flat will have a unique label and all member cells of a flat will bear its label. The LowEdges queue will not contain any non-draining flats; however, such flats must be eliminated from the HighEdges queue. Fig. 1. Step 2: gradient away from higher terrain (see Section 2.3 and Algorithm 5). The number of increments applied to each cell is shown. (a) Numbers of increments after the first iteration. This first set of incremented cells are the edge cells adjacent to higher (but not to lower) terrain; they were found in Step 1 and stored in the queue HighEdges. After being incremented, each of these cells will add its unincremented neighbors—those just beyond the thick black line—to the queue. In (b) the neighbors of the edge cells have been popped off of the queue and incremented, advancing the black line. This process continues through (c), which shows the final number of increments, though the gradient is the inverse of what is desired. By noting that at most 3 increments were applied to this flat it is possible to obtain the desired gradient (d) by subtracting each cell of (c) from 3. This is done in Step 3 (Fig. 3). Note that the cells (F2, F3, F4) adjacent to lower terrain are ignored. (e) shows the resulting drainage field. The flow directions of cells E3–E5 will be uniquely determined by a gradient towards lower terrain in Step 3 (see Fig. 3d). Fig. 2. Step 3: gradient towards lower terrain (see Section 2.4 and Algorithm 6). The number of increments applied to each cell is shown. (a) Numbers of increments after the first iteration. This first set of incremented cells are the edge cells adjacent to lower terrain; they were found in Step 1 and stored in the queue LowEdges. After being incremented, each of these cells will add its unincremented neighbors—those just beyond the thick black line—to the queue. In (b) the neighbors of the edge cells have been popped off of the queue and incremented, advancing the black line. This process continues through (c), (d), and (e). (e) shows the final number of increments. (f) shows the resulting drainage field. There is a second part to Step 3, shown below, which combines the gradients. R. Barnes et al. / Computers & Geosciences ∎(∎∎∎∎)∎∎∎–∎∎∎ 3 Please cite this article as: Barnes, R., et al., An efficient assignment of drainage direction over flat surfaces in raster digital elevation models. Computers & Geosciences (2013), http://dx.doi.org/10.1016/j.cageo.2013.01.009i Figure 10.7: In a flat surrounded by higher terrain (dark grey) with a single lower-elevation outlet (light grey), we can use a gradient towards from the outlet to route water out of the flat and towards the outlet. For this, we can iteratively assign (tiny or symbolic) elevation increases in the flat starting from the outlet until all non-draining cells in the flat have been covered. Figure from Barnes et al. (2014a). Since every cell surrounding a flat must be either a low edge or a high edge cell, if none of either are found, then the DEM has no flats. If there are no low edge cells, but there are some high edge cells, then the DEM has flats but they cannot be drained. In both such instances, the algorithm indicates this and proceeds no further. Otherwise, the algorithm then labels each unique flat (see Algorithm 4). Let the data structure holding these labels be called Labels, and let it be initialized to a value N O L ABEL . To apply the labels, each cell cin the LowEdges queue is used as the seed to a flood fill algorithm. If chas not yet been labeled then a new label l is given to it, its elevation eis noted, and all cells which can be reached from cwhile traversing only cells of elevation eare given the same label l. A FIFO queue could be used for this, or a D8 variant of Heckbert (1990), or one of many other solutions to the “connected components”problem (Grana et al., 2010 provide a good overview). Recursive solutions should be avoided in implementations using languages without proper tail recursion, as stack overflows are likely when processing large DEMs. In the end each flat will have a unique label and all member cells of a flat will bear its label. The LowEdges queue will not contain any non-draining flats; however, such flats must be eliminated from the HighEdges queue. Fig. 1. Step 2: gradient away from higher terrain (see Section 2.3 and Algorithm 5). The number of increments applied to each cell is shown. (a) Numbers of increments after the first iteration. This first set of incremented cells are the edge cells adjacent to higher (but not to lower) terrain; they were found in Step 1 and stored in the queue HighEdges. After being incremented, each of these cells will add its unincremented neighbors—those just beyond the thick black line—to the queue. In (b) the neighbors of the edge cells have been popped off of the queue and incremented, advancing the black line. This process continues through (c), which shows the final number of increments, though the gradient is the inverse of what is desired. By noting that at most 3 increments were applied to this flat it is possible to obtain the desired gradient (d) by subtracting each cell of (c) from 3. This is done in Step 3 (Fig. 3). Note that the cells (F2, F3, F4) adjacent to lower terrain are ignored. (e) shows the resulting drainage field. The flow directions of cells E3–E5 will be uniquely determined by a gradient towards lower terrain in Step 3 (see Fig. 3d). Fig. 2. Step 3: gradient towards lower terrain (see Section 2.4 and Algorithm 6). The number of increments applied to each cell is shown. (a) Numbers of increments after the first iteration. This first set of incremented cells are the edge cells adjacent to lower terrain; they were found in Step 1 and stored in the queue LowEdges. After being incremented, each of these cells will add its unincremented neighbors—those just beyond the thick black line—to the queue. In (b) the neighbors of the edge cells have been popped off of the queue and incremented, advancing the black line. This process continues through (c), (d), and (e). (e) shows the final number of increments. (f) shows the resulting drainage field. There is a second part to Step 3, shown below, which combines the gradients. R. Barnes et al. / Computers & Geosciences ∎(∎∎∎∎)∎∎∎–∎∎∎ 3 Please cite this article as: Barnes, R., et al., An efficient assignment of drainage direction over flat surfaces in raster digital elevation models. Computers & Geosciences (2013), http://dx.doi.org/10.1016/j.cageo.2013.01.009i drainage network drainage basin drainage divide 10.5 Drainage networks and basins Interpreting DTM cells as nodes and the flow direction as directed edges connecting them yields the drainage network of a DTM. However, it is usually best to filter out the least important parts of the network using a flow accumulation threshold. A good rule of thumb for this threshold is the mean flow accumulation in the DTM, but an exact value is usually set by trial and error until the desired parts of the network are kept. Based on a computed drainage network, it is then possible to extract the drainage basins of a DTM by considering the areas that are drained by one or more nodes of the network (Figure 10.8). This operation can be performed in many different places, such as the end node of a river (yielding its river basin), the nodes just before junctions in the network (yielding the drainage basins of the tributaries of a river), or the end nodes of a selected part of the network (yielding the drainage basin of a sea or ocean). The lines that separate adjacent drainage basins are drainage divides, which form topographical ridges. 10.6 Notes and comments 111 Figure 10.8: The areas that drain to all the oceans can be computed by selecting the DTM cells on the coastline of these oceans and finding the areas that drain through them. Note the endorheic basins that drain to none of these cells. These actually form sinks in the DTM. From Wikimedia Commons. 10.6 Notes and comments Beven (2012) is a good reference book on hydrology. It covers how to make much more complex runoff models than the ones described here. O’Callaghan and Mark (1984) was the original paper to describe the D8 method. Fairfield and Leymarie (1991) modify D8 into the stochastic rho8 method. Quinn et al. (1991) describes the original MFD method. Tarborton (1997) describes the alternative (D ∞ ) MFD method and contains nice figures comparing multiple methods. Barnes et al. (2014b) describes how to fill in sinks, while Metz et al. (2011) describes how to use a variation of 𝐴∗ search algorithm to route water out of them. Barnes et al. (2014a) describes how to assign the drainage direction over flats. 10.7 Exercises 1. Given a raster map of precipitation values, how would you be able to improve the flow accumulation estimates? 2. Why is the flow width important? 3. You have a cycle in your drainage network. How can that happen? How would you solve it? 4. How can you detect endorheic basins without finding all other basins first? 5. Come up with an algorithm to identify flats in a DTM. 118 11 Processing point clouds invert Figure 11.6: Basic idea behind the CSF algorithm for ground filtering of a point cloud: inverting the data and letting a cloth fall. t0 t1 t2 t3 Figure 11.7 process, we aim to control the tension (or rigidity) of the cloth, so that areas where there is no sample point (eg where there are large buildings or water) can be filled realistically. The CSF algorithm is a simplification of an algorithm in computer graphics to simulate a piece of cloth falling on an object. The cloth is modelled as a surface formed of particles (vertices) that are regularly distributed on a grid, these particles have a mass and they are connected to their neighbours (4-neighbours in this case). For terrains (2.5D objects), the particles are constrained to only move vertically. Two factors influence the 𝑧 -value of a particle during the cloth falling process: 1. external forces: in this case this is the gravity pulling down a particle; 2. internal forces: the tension in the cloth, which is modelled by the interactions between a particle and its neighbours. As particles fall down, some will reach the ground and become unmovable. These will potentially be neighbours to movable ones, whose elevation will be controlled by how we define the rigidity of the cloth. As shown in Figure 11.7, the process is iterative. We first define a cloth formed of particles, and then for each iteration we calculate the next 𝑧 -value of each particle based on the vector of displacement from the external and internal forces at the previous step. If a particle is movable (ie it has not reached the ground yet), then the gravity force is applied (a vectorpointingdownwards;itsmagnitudewilldependonthemomentum of the particle) and afterwards the internal forces are applied. Notice that in Figure 11.7, the particle in red at 𝑡3 was moved downwards because of the gravity, but its internal forces are a vector pointing upwards since its 2 neighbours (it would be 4 for a 2D case) have higher 𝑧-values. The algorithm is detailed in Algorithm 6. Initialisation of the cloth. The cloth is first initialised at an arbitrary height above the highest points are the point cloud. The cloth is formed of particles regularly distributed according to a user-defined parameter. We assume that all particles have the same mass, and we define arbitrarily a first displacement vector due to the gravity (pointing downwards) (line 6). For each particle 𝑝 , we need to define the lowest elevation it can move, once it reaches it it is labelled as unmovable (line 5). The lowest elevation of one particle 𝑝 is defined as the original elevation of the closest sample point 𝑠after projecting both to the 2D-plane. 11.3 Ground filtering 119 Algorithm 6: CSF algorithm 1Input: A set 𝑆of sample points from a point cloud; resolution 𝑟of the cloth grid; tolerance 𝜖𝑧𝑚𝑎𝑥 to stop the iterations; tolerance 𝜖𝑔𝑟𝑜𝑢𝑛𝑑 to classify points in 𝑆 Output: The points in 𝑆are classified as ground/non-ground 2invert 𝑆 3initialise the cloth 𝐶at an elevation 𝑧0higher than the highest elevation 4for all 𝑝∈𝐶do 5𝑝𝑧𝑚𝑖𝑛 =lowest possible elevation based on 𝑆 6𝑝𝑧𝑝𝑟𝑒𝑣 =𝑧0+𝑑𝑖𝑠𝑝𝑙𝑎𝑐𝑒𝑚𝑒𝑛𝑡 7𝑝𝑧𝑐𝑢𝑟 =𝑧0 8while Δ𝑧> 𝜖𝑧𝑚𝑎𝑥 do /*external forces */ 9for all 𝑝∈𝐶do 10 if 𝑝is movable then 11 tmp = 𝑝𝑧𝑐𝑢𝑟 12 𝑝𝑧𝑐𝑢𝑟 =(𝑝𝑧𝑐𝑢𝑟 −𝑝𝑧𝑝𝑟𝑒𝑣)+𝑝𝑧𝑐𝑢𝑟 13 𝑝𝑧𝑝𝑟𝑒𝑣 =𝑡𝑚𝑝 /*internal forces, process once each pair 𝑒of adjacent particles */ 14 for all 𝑒∈𝐶do 15 𝑝0=𝑒𝑠𝑡𝑎𝑟𝑡 16 𝑝1=𝑒𝑒𝑛𝑑 17 update 𝑝0𝑧𝑐𝑢𝑟 and 𝑝1𝑧𝑐𝑢𝑟 if they are movable /*calculate the max Δ𝑧*/ 18 for all 𝑝∈𝐶do 19 if (𝑝𝑧𝑐𝑢𝑟 −𝑝𝑧𝑝𝑟𝑒𝑣)>Δ𝑧then 20 Δ𝑧=(𝑝𝑧𝑐𝑢𝑟 −𝑝𝑧𝑝𝑟𝑒𝑣) 10m 2m v = (0,0,4) 10m 2m v = (0,0,4) v = (0,0,−4) (a) niunmovable (b) nimovable ni p ni p Figure 11.8: Internal forces in the CSF algorithm: 2 cases are possible. Internal forces. Theinternalforcesareappliedonlytomovableparticles; once a particle has been labelled as unmovable it cannot be moved again. Given a movable particle 𝑝 , we apply the internal forces by individually looking at its 4 neighbours 𝑛𝑖 . For each neighbour 𝑛𝑖 , there are two cases (see Figure 11.8): 1. 𝑛𝑖is unmovable: only 𝑝is moved, towards 𝑛𝑖. 2. 𝑛𝑖 is movable: the idea is that both 𝑝 and 𝑛𝑖 will try to move towards each other to the same height. The vector applied to each will thus be in opposite direction. Controlling the tension/rigidity. Notice that in Figure 11.8b, both particles are moved to the same elevation, but that it is also possible to scale the internal forces displacement vector, eg to 0.8 of its length (and thus decrease the tension in the cloth). Lower internal force displacement means that the particles will move more during a single iteration and so the tension in the cloth is effectively reduced. The same idea applies to Figure 11.8a, the displacement vector can be controlled by scaling the displacement vector. In Figure 11.8a, it is 0.5 of the difference in elevation, butif lesstension iswanted, then thescale could be forinstance0.4 (so that 𝑝 has an internal displacement ® 𝑣=(0,0,3.2) , because 0.4∗(10 −2)=3.2 ) or 0.3. 120 11 Processing point clouds 2CSF is implemented in several open-source libraries The description of the CSF algorithm in this book is a simplification of the original algorithm (see Zhang et al. (2016)), it omits the post-processing to take into account steep slopes. The complete algorithm is implemented in the open-source software CloudCompare ( https://cloudcompare.org ) and in the open-source library PDAL (https://pdal.io). How the process ends. This iterative process is repeated until the maximum displacement of all particles is less than a user-defined parameters (𝜖𝑧𝑚𝑎𝑥); or until a certain number of iterations has been performed. Two possible outputs. When the process is completed, the surface of the cloth can be used to obtain two outputs (see Figure 11.9): 1. classification of points into ground/non-ground 2. a surface representing the ground (the cloth) If the surface of the cloth is used, it is for instance possible to triangle it or to create a grid from it. If a segmentation/classification of the input points is wanted, then the distance between a sample point of the original point cloud and the cloth can be used (this is the parameter 𝜖𝑔𝑟𝑜𝑢𝑛𝑑 in Algorithm 6). If this distance is less than a given user-defined threshold, then the sample point is a ground point. 11.3 Ground filtering 121 (a) Original point cloud (b) Output #1: the ground points (with the ground surface shown in grey) (c) Output #2: the ground surface Figure 11.9: Two outputs of the CSF algorithm for a given area. 122 11 Processing point clouds spheres, cones, cylinders, planes, etc. can be detected minimal set 11.4 Shape detection Shape detection is used to automatically detect simple shapes—such as planes—in a point cloud. See for example Figure 11.10 where the points are randomly coloured according to the corresponding planar surfaces. Shape detection is an important step in the extraction and reconstruction of more complex objects, eg man-made structures such as buildings are often composed of planar surfaces. In this section, three shape detection methods will be introduced: 1. RANSAC 2. region growing 3. Hough transform First, some common terminology. Let 𝑃 denote a point cloud, if we perform shape detection on 𝑃 we aim to find a subset of points 𝑆⊂𝑃 that fit with a particular shape. Most shape detection methods focus on shapes that can be easily parametrised, such as a line, a plane, or a sphere. If we specify values for the parameters of such a parametrised shape, we define an instance of that shape. For example, a line in the plane can be parametrised using the equation 𝑦=𝑚𝑥 +𝑏 , in this case 𝑚 and 𝑏 are the parameters. We can create an instance of a line by specifying values for its parameters 𝑚 and 𝑏 , respectively fixing the slope and the position of the line. In the following, the methods are described in a general way, ie without specialisations for one particular shape. Only for illustrative purposes specific shapes such as a line or a plane are used to (visually) explain the basic concept of each shape detection method, but the same could be done with spheres, cones, or other shapes. 11.4.1 RANSAC RANSAC is short for RANdom SAmpling Consensus and, as its name implies, works by randomly sampling the input points. In fact it starts by picking a random set of points 𝑀⊂𝑃 . This set 𝑀 is called the minimal set and contains exactly the minimum number of points that is needed to uniquely construct the shape that we are looking for, eg 2 for a line and 3 for a plane. From the minimal set 𝑀the (unique) shape instance I is constructed (see Figures 11.11b and 11.11c). The algorithm then checks for each point 𝑝∈ {𝑃\𝑀} if it fits with I . This is usually Figure 11.10: Planar regions in the AHN3 point cloud. Each region was assigned a random colour. 11.4 Shape detection 123 (a) Input points  score=0 (b) 1st minimal set  score=11 (c) 2nd minimal set (d) Detected line instance Figure 11.11: RANSAC for line detection (𝑘=2iterations) inlier consensus set score done by computing the distance 𝑑 from 𝑝 to I and comparing 𝑑 against a user-defined threshold 𝜖 . If 𝑑< 𝜖 we say that 𝑝 is an inlier, otherwise 𝑝 is an outlier. The complete set of inliers is called the consensus set, and its size is referred to as the score. The whole process from picking a minimal set to computing the consensus set and its score, as shown in Algorithm 7, is repeated a fixed number of times, after which the shape instance with the highest score is outputted (Figure 11.11d). Algorithm 7: The RANSAC algorithm Input: An input point cloud 𝑃, the error threshold 𝜖, the minimal number of points needed to uniquely construct the shape of interest 𝑛, and the number of iterations 𝑘 Output: the detected shape instance I 𝑏𝑒𝑠𝑡 1𝑠𝑏𝑒𝑠𝑡 ←0; 2I 𝑏𝑒𝑠𝑡 ←nil; 3for 𝑖←0. . . 𝑘 do 4𝑀←𝑛randomly selected points from 𝑃; 5I←shape instance constructed from 𝑀; 6𝐶← ∅ ; 7for all 𝑝∈𝑃\𝑀do 8𝑑←distance(𝑝, I); 9if 𝑑< 𝜖 then 10 add 𝑝to 𝐶; 11 𝑠←score(𝐶); 12 if 𝑠>𝑠𝑏𝑒𝑠𝑡 then 13 𝑠𝑏𝑒𝑠𝑡 ←𝑠; 14 I 𝑏𝑒𝑠𝑡 ←I; The most touted benefit of RANSAC is its robustness, ie its performance in the presence of many outliers (up to 50%). Other algorithms to identify planes, eg fitting a plane with least-square adjustment, are usually more sensitive to the presence of noise and outliers (which are always present in real-world datasets). The probability that a shape instance is detected with RANSAC depends mainly on two criteria: 1. the number of inliers in 𝑃, and 2. the number of iterations 𝑘. Naturally, it will be easier to detect a shape instance in a dataset with a relatively low number of outliers. And it is more likely that a shape instance is found if more minimal sets are evaluated. Picking a sufficiently 124 11 Processing point clouds (a) Input points with normals and three seed points 0◦ 0◦ 0◦ (b) Start growing. Add neighbours if the normal angle is small. 0◦ 0◦ 0◦ (c) Continue growing from new region point 90◦ 0◦ 0◦ (d) Stop growing where the normal angle is too great (e) Finalregions fromallthree seedpoints Figure 11.12: Region growing for plane detection based on the angle between neighbouring point normals seed point 1: seeAppendix Cto estimatethenormal in a point cloud high 𝑘 is therefore important for the success of the algorithm, although a higher 𝑘also increases the computation time. Because of the random nature of RANSAC, the minimal sets that it will evaluate will be different every time you run the algorithm, even if the input data is the same. The detected shape instance can therefore also be different every time you run the algorithm; RANSAC is therefore said to beanon-deterministic algorithm. This could be a disadvantage. Time complexity. The time complexity of RANSAC is O(𝑘𝑛) , where 𝑛 is the size of 𝑃. 11.4.2 Region growing Region growing works by gradually growing sets of points called regions that fit a particular shape instance. A region 𝑅 starts from a seed point, ie a point that is suspected to fit a shape instance. More points are added to 𝑅 by inspecting candidate points, ie points in the neighbourhood of the members of 𝑅 . To check if a candidate point 𝑐 should be added to 𝑅 , a test is performed. In the case of region growing for plane detection (see Figure 11.12) this test entails computing the angle between the normal vector of 𝑐 and the normal vector 1 of its neighbour in 𝑅 . If this angle is small it is assumed that 𝑐 lies in the plane instance that corresponds to 𝑅 , and that it can therefore be added to 𝑅 . Otherwise 𝑐 is ignored (Figure 11.12d). This process of growing 𝑅 continues until no more candidates can be found that are compatible with 𝑅 . When this happens, the algorithm proceeds to the next seed point to grow a new region. Algorithm 8gives the pseudo-code for the region growing algorithm. Noticethat the set 𝑆 isused tokeeptrackofthe points in the current region whose neighbours still need to be checked. Also notice that candidate points that are already assigned to a region are skipped. The seed points can be generated by assessing the local neighbourhood of each input point. For example in case of plane detection one could fit a 11.4 Shape detection 125 Algorithm 8: The Region growing algorithm Input: An input point cloud 𝑃, a list of seed points 𝐿𝑆, a function to find the neighbours of a point 𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑢𝑟𝑠() Output: A list with detected regions 𝐿𝑅 1𝐿𝑅← []; 2for each 𝑠in 𝐿𝑆do 3𝑆← {𝑠}; 4𝑅← ∅; 5while 𝑆is not empty do 6𝑝←pop(𝑆); 7for each candidate point 𝑐∈neighbours(𝑝)do 8if 𝑐was not previously assigned to any region then 9if 𝑐fits with 𝑅then 10 add 𝑐to 𝑆; 11 add 𝑐to 𝑅; 12 append 𝑅to 𝐿𝑅; plane through each point neighbourhood and subsequently sort all points on the fitting error. Points with a low plane fitting error are probably part of a planar region so we can expect them to be good seeds. To compute the point neighbourhoods a k-nearest neighbour search or a fixed radius search can be used, which can both be implemented efficiently using a 𝑘d-tree (see Section 12.2). Notice that region growing is based on the idea that we can always find a path of neighbouring points between any pair of points within the same region. This does mean that two groups of points that fit the same shape instance but are not connected through point neighbourhoods will end up in different regions. Other shape detection methods described in this chapter do not need point neighbourhood information. Time complexity. If we assume that 1. the number of seeds in 𝐿𝑆is linear with 𝑛, ie the size of 𝑃, 2. the size of 𝑆is at most 𝑛, and that 3. a 𝑘 nn search takes O(𝑘log 𝑛) , where 𝑘 is the number of neighbours, we come to a worst-case time complexity of O(𝑛2𝑘log 𝑛) . In practice it should be better since 𝑆 is not likely to be 𝑛 large, and it will get smaller the more regions have been found. 11.4.3 Hough transform The Hough transform uses a voting mechanism to detect shapes. It lets every point 𝑝∈𝑃 vote on each shape instance that could possibly contain 𝑝 . Possible shape instances thus accumulate votes from the input points. The detected shape instances are the ones that receive the highest number of votes. To find the possible shape instances for 𝑝 , the algorithm simply checks all possible parameter combinations that give a shape instance that fits with 𝑝. It is thus important to choose a good parametrisation of the shape that is to be detected. For instance when detecting lines one could use the slope-intercept form, ie 𝑦=𝑚𝑥 +𝑏 . However, this particular 126 11 Processing point clouds O y x r φ (a) Line parametrisation O 1 2 3 4 1 2 3 4 y x (b) Input points O 1 2 3 4 1 2 3 4 y x (c) Line instances for each point 𝑟𝜙 0◦90◦ 0.0 0 0 0.511 1.0 0 0 1.513← 2.03←0 2.5 0 1 3.0 0 0 3.510 4.012 4.511 (d) Accumulator contains the number of votes for each line instance. O 1 2 3 4 1 2 3 4 y x (e) Detected line instances with a minimal vote count of 3. Figure 11.13: Hough transform for line detection with a 10 ×2 accumulator. The (𝜙, 𝑟) line parametrisation is chosen because this form can represent vertical lines (unlike the 𝑦=𝑚𝑥 +𝑏form for example). Hesse normal form accumulator quantisation parametrisation can not easily represent vertical lines, because 𝑚 would need become infinite which is computationally difficult to manage. A better line parametrisation is the Hesse normal form which is defined as 𝑟=𝑥cos 𝜙+𝑦sin 𝜙. As illustrated in Figure 11.13a, (𝑟, 𝜙) are the polar coordinates of the point on the line that is closest to the origin, ie 𝑟 is the distance from the origin to the closest point on the line, and 𝜙∈ [0◦,180◦] is the angle between the positive 𝑥 -axis and the line from the origin to that closest point on the line. This parametrisation has no problems with vertical lines (ie 𝜙=90◦ ). Similarly, for plane detection we can use the parametrisation 𝑟=𝑥cos 𝜃sin 𝜙+𝑦sin 𝜙sin 𝜃+𝑧cos 𝜙. Where (𝑟, 𝜃,𝜙) are the spherical coordinates of the point on the plane that is closest to the origin. Figure 11.13 shows an example for line detection with the Hough transform and Algorithm 9gives the full pseudo-code. The votes are saved in an accumulator which is essentially a matrix with an axis for each parameter of the shape, eg for detecting lines we would need two axes (See Figure 11.13d). Notice that each element in the accumulator represents one possible shape instance. Because each axis only has a limited number of elements, each parameter is quantised. This means that each parameter is restricted in the possible values it can have. The chosen quantisation determines the sensitivity of the accumulator. The accumu- 11.5 Notes and comments 127 Algorithm 9: The Hough transform algorithm Input: An input point cloud 𝑃 , an accumulator matrix 𝐴 , a detection threshold 𝛼 Output: A list with detected shape instances 𝐿𝐼 1for each 𝑝in 𝑃do 2for each instance 𝑖from 𝐴that fits with 𝑝do 3increment 𝐴[𝑖]; 4𝐿𝐼←all shape instances from 𝐴with a more than 𝛼votes; 2: https://en.wikipedia.org/wiki/ Random_sample_consensus#Paramete rs lator of Figure 11.13d for example, can only detect horizontal and vertical lines, because the 𝜙 parameter is quantised in only two possible values. Notice that the accumulator can be made more sensitive by choosing a finer quantisation, effectively increasing the size of the accumulator (although that will also make the algorithm run slower). Time complexity. The time complexity of the Hough transform algorithm as discussed here is O(𝑛𝑚) , where 𝑚 is the number of elements in the accumulator. 11.5 Notes and comments Arge et al. (2010) introduced the outlier detection method for echosounding datasets by cutting long edges in a TIN. Axelsson (2000) originally proposed the greedy TIN densification algorithm forground filtering. Healsodescribes howtohandlediscontinuities in the terrain such as cliffs. It should be said that his paper is scarce on details, and many variations of the algorithms have been proposed so that small/low objects are filtered out and so that it performs well in densely forested areas. See for instance Lin and Zhang (2014). The cloth simulation filter (CSF) algorithm is from Zhang et al. (2016). The original paper has a somewhat complex definition that has been simplified and modified for this book. Also, the original has a postprocessing step for steep slope that is omitted in this book. A comparison with several other ground filtering methods can be found in the work of Meng et al. (2010). Fischler and Bolles (1981) originally introduced the RANSAC algorithm and applied to cartography in that same paper. On Wikipedia you can read how you can compute the required number of RANSAC iterations to achieve a certain probability of success given that you know how many outliers there are in your dataset2. Limberger and Oliveira (2015) describes how to efficiently do plane detection in large point clouds using a variant of the Hough transform. 11.6 Exercises 1. The LAS standard gives a global point offset in the header. What is the benefit of using such a global offset? 134 12 Handling and processing massive terrains 5 10 10 5 0 (1, 3) (2, 7)(2, 7) (4, 8) (5, 6) (8, 9) (9, 5) (6, 4) (9, 1) (5, 6) (2, 7) (1, 3) (4, 8) (9, 5) (6, 4) (9, 1) (8, 9) x x y y q= (4.5, 4) (4.5, 4) far (a) 5 10 10 5 0 (1, 3) (2, 7)(2, 7) (4, 8) (5, 6) (8, 9) (9, 5) (6, 4) (9, 1) (5, 6) (2, 7) (1, 3) (4, 8) (9, 5) (6, 4) (9, 1) (8, 9) x x y y q= (4.5, 4) (b) 5 10 10 5 0 (1, 3) (2, 7)(2, 7) (4, 8) (5, 6) (8, 9) (9, 5) (6, 4) (9, 1) (5, 6) (2, 7) (1, 3) (4, 8) (9, 5) (6, 4) (9, 1) (8, 9) x x y y ? q= (4.5, 4) dist (c) 5 10 10 5 0 (1, 3) (2, 7)(2, 7) (4, 8) (5, 6) (8, 9) (9, 5) (6, 4) (9, 1) (5, 6) (2, 7) (1, 3) (4, 8) (9, 5) (6, 4) (9, 1) (8, 9) x x y y ? q= (4.5, 4) (d) Figure 12.6: Several states for the nearest neighbour query based on a 𝑘 d-tree, 𝑞=(4.5,4.0) is the query point and 𝑐=(6,4) is the nearest point. 12.3 Streaming paradigm to construct massive TINs/grids 135 swapping external memory algorithms 𝑑𝑖𝑠𝑡(𝑞, 𝑐𝑡𝑒𝑚𝑝) , 𝑑𝑖𝑠𝑡() being the Euclidean distance between 2 points. If that distance is shorter, then it is possible that one point in the subtree is closer than 𝑐𝑡𝑒𝑚𝑝 , and thus that subtree must be visited. If not, then the whole subtree can be skipped, and the algorithm continues. Figure 12.6c shows this idea after ( 1,3 ) has been visited. 𝑐𝑡𝑒𝑚𝑝 is ( 5,6 ), and we must decide whether the subtree right of ( 2,7 ) must be visited. In this case it must not be visited because the bounding box (light blue region) is 3.0 units from 𝑞 , and 𝑑𝑖𝑠𝑡(𝑞, 𝑐𝑡𝑒𝑚𝑝) is around 2.07; it is thus impossible that one point inside the subtree be closer than (5,6). The next step is verifying whether the subtree right of the root could contain a point closer than 𝑐𝑡𝑒𝑚𝑝 . In the Figure 12.6d, this is possible since the bounding box is only 0.5 unit from 𝑞 , and thus the subtree must be visited. The algorithm continues until all subtrees have either been visited or eliminated. At the end, 𝑐is (6,4). Time complexity. To insert a new point, and to search for a nearest neighbour, the time complexity on average is O(log 𝑛) ; this is assuming the tree is balanced, if not it could be O(𝑛) in the worst-case. The tree stores one node per point, thus the space complexity is O(𝑛). 𝑚 -closest neighbours. The algorithm can be extended in several ways by simple modifications. It can provide the 𝑚 nearest neighbours to a point by maintaining 𝑚 current closest points instead of just one. A branch is only eliminated when 𝑚 points have been found and the branch cannot have points closer than any of the 𝑚 current bests. This can help improve significantly the running time of several operations described in this book: IDW with a fixed number of neighbours (Section 5.3), extracting shapes from point clouds (Section 11.4), estimating normals in point clouds (Appendix C), calculating the spatial extent (Chapter 13), are only but a few examples. 12.3 Streaming paradigm to construct massive TINs and grids from point clouds The incremental construction algorithm for the Delaunay triangulation (DT), presented in Chapter 4, will not work if the size of the input dataset is larger than the main memory. Or if it works, it will be very slow. The reason for this is that the data structure for the DT (to store the points coordinates, and also the data structure for the triangles and their topological relationships) will not fully fit in the main memory. Therefore, part of it will be in the main memory (say 16GB of RAM) and the other part will be on the harddrive. The operating system controls which parts are in memory and which parts are on the harddrive, and we call swapping the operations to transfer between them. One solution to this problem is to design external memory algorithms. These basically do not rely on the operating system to decide which parts of the data structure are stored on the disk, but improve the process by explicitly storing temporarily files and having explicit rules for the swapping of data between the disk and the memory. The main drawbacks of this approach are that the design of such algorithms is rather complex, 136 12 Handling and processing massive terrains streaming data 4: https://en.wikipedia.org/wiki/ Stream_(computing) batch processing Tobler (1970) that for different problems different solutions have to be designed, and that for problems like the DT construction a lot of large temporary files need to be created. We discuss in this section an alternative approach to dealing with massive datasets: streaming.Astream is a sequence of data—in theory it can be infinite!—that is available over a period of time, and “can be thought of as items on a conveyor belt being processed one at a time rather than in large batches” 4 . One concrete example is YouTube: to watch a given video a user does not need to first download the whole file, she can simply start watching the video as soon as the first KB are downloaded. The content of the video is downloaded as she watches the video, and if she fast-forwards to, for instance, 5:32s then only the KB of content from where the cursor is need to be downloaded to watch the video. At no moment is the full video downloaded to the user’s device. Batch algorithms, like the incremental insertion algorithm described in Section 4.4, require to have all the points in memory to work. By contrast, a streaming algorithm operates only locally and can thus, in theory, process infinitely large datasets. The streaming paradigm can be used to process geometries (points, meshes, polygons, etc.) but it is slightly more involved than for a simple video. Since the First Law of Geography of Tobler stipulates that “everything is related to everything else, but near things are more related than distant things”, if we wanted to calculate the slope at one location in a point cloud we would need to retrieve all the neighbouring points and potentially calculate locally the DT. The question is: is it possible to do this without reading the whole file and only process one part of it? We focus in the following on the creation of a DT. The main issue that we are facing is that a triangle is respecting the Delaunay criterion if its circumcircle is empty of any point, therefore while constructing the DT we need to test all the other points in the dataset to ensure that a given triangle is Delaunay (or not). Streaming would mean here: can we assess that a given triangle is Delaunay without having to read/download the whole file? 12.3.1 Overview of streaming DT construction Figure 12.7 shows the overview of the processes involved for the construction of a DT with the streaming paradigm. Think of the stream as a 1D list of objects (points, triangles, tags, etc.) and the aim is to be able to perform an operation without having in memory the whole stream. 12.3.2 Finaliser: adding finalisation tags to the stream The key idea is to preprocess a set 𝑆 of points and insert finalisation tags informing that certain points/areas will not be needed again. For the DT construction, as shown in Figure 12.8, this can be realised by constructing aquadtree of 𝑆 ;a finalisationtag basically informs the following processes that a certain cell of the quadtree is empty, that all the points inside have been processed (have already appeared in the stream). In practice, this is performed by reading a LAS/LAZ file (or any format with points) twice from disk: 12.3 Streaming paradigm to construct massive TINs/grids 137 2Streaming is realised with Unix pipes The key to implementing streaming of geometries is to use Unix pipes (also called pipelines). Pipelines were designed by Douglas McIlroy at Bell Labs during the development of Unix, and they allow to chain several processes together. The output of a process becomes the input of the next one, and so on (the data flowing through the process is the stream). Given 2 processes, the 2nd one can usually start before the 1st one has finished processing all the data. In Unix, the pipe operator is the vertical line “ | ”, and several commands can be chained with it: “ cmd1 | cmd2 | cmd3 ”. A simple example would be “ls -l | grep json | wc -l” which would: 1. list all the files in the current directory (one file name per line); 2. send this to the operator grep which would discard all lines not having the keyword "json"; 3. send this to the operator “ wc -l ” which counts the number of line. finalisation of triangles 1. the first pass will count how many points are in each cell of the quadtree (and store the results) 2. the second pass will read again sequentially each point (and send it in the stream), and decrement the counter for each cell. When it is empty, a finalisation tag will be added to the stream. 12.3.3 Triangulator The input of the triangulator is the output of the finaliser: a set of points with finalisation tags. The triangulator will triangulate the points as described in Section 4.4, but will attempt to remove from memory the triangles that are final, those that we are sure will never be modified (since it is guaranteed that no new points will fall inside their circumcircle). This is performed with the finalisation tags and the following observation (see Figure 12.9): a triangle inside a finalised quadtree cell (ie where all the points in the streams inside that cell have been read) is final if its circumcircle does not encroach on an active quadtree cell. If its circumcircle overlaps with an active quadtree cell, then it is possible that later in the stream a new point will be added inside the circle, and thus the triangle will not be Delaunay. Final triangles can be removed from memory and written directly to disk; it is however possible to add another process to the pipeline and send the final triangles to them (eg to create a grid or to extract isolines). Notice also that the memory the triangle was using can be reused to store another new triangle created by new points arriving in the stream. 12.3.4 Spatial coherence The construction of a DT with the streaming paradigm will only succeed (in the sense that the memory footprint will stay relatively low) if the spatial coherence of the input dataset is high. It is defined by Isenburg 138 12 Handling and processing massive terrains input points finaliser triangulator output triangulation isoline extraction Finalisation tags One of these 2 Figure 12.7: Overview of the streaming pipeline to construct a DT (or extract isolines). 1. x y z 2. x y z 3. x y z 4. x y z 5. x y z 6. x y z 7. x y z 8. x y z 9. x y z 1. x y z 2. x y z 3. x y z finalise c3 4. x y z finalise c2 5. x y z 6. x y z 7. x y z 8. x y z finalise c0 9. x y z finalise c1 9 points and a ”normal” stream 9 points and a stream with finalisation tags 1 4 2 3 5 6 7 9 8 7 8 c0 c1 c2 c3 1 4 2 3 5 6 9 Figure 12.8: How the finaliser modifies to input and injects finalisation tags. Left: 9 points and the related streams (just a list of the points with coordinates). Right: if a quadtree of depth 1 (in orange) is used (with 4 cells), then the stream would be augmented with finalisation tags. Isenburg et al. (2006a) et al. as: “a correlation between the proximity in space of geometric entities and the proximity of their representations in [the file]”. They demonstrate that real-world point cloud datasets often have natural spatial coherence because the points are usually stored in the order they were collected. If we shuffled randomly the points in an input file, then the spatial coherence would be very low and the finalisation tags in the stream coming out of the finaliser would be located at the end (and not distributed in the stream). Itis possible tovisualisethespatial coherence ofa dataset by colouring,for an arbitrary grid, the positions of the first and last points; Figure 12.10 12.3 Streaming paradigm to construct massive TINs/grids 139 Figure 12.9: The DT at a given moment during the triangulation process. Blue quadtree cells are not finalised yet, white onesare;yellowtriangles arestillinmemory (their circumcircles (in red) encroach onunfinalised cells); white triangleshave been written to disk since their circumcircles do not encroach on an active cell (some green circles shown as example). position in file/stream first point position last point position low spatial coherence high spatial coherence Figure 12.10: The colour map used for the position of a point in the file, and 3 examples of cells. gives an example. The idea is to assign a colour map based on the position of the points in the file, and to colour the centre of the cells with the position of the first point inside that cell, and to colour the boundary of the cell with the position of the last point. Figure 12.11 illustrates the spatial coherence for 2 tiles of the AHN3 dataset in the Netherlands. Notice that the cells are generally of the same colour, which means that the spatial coherence is relatively high. It is interesting to notice that the two datasets have different patterns probably because they were compiled by different companies, who used different equipment and processing software to generate the datasets. 12.3.5 Streaming cannot solve all problems related to terrains The ideas behind streaming are very useful for certain local problems (eg interpolation, creation of grids, extraction of contour lines), but unfortunately they cannot be used directly (or it would be extremely challenging) for global problems such as visibility or flow modelling. 140 12 Handling and processing massive terrains positionin thefile 0 120M 240M 360M 480M AHN3 - 37EN1 (a) position in the file 0 100M 200M 300M 400M AHN3 - 07BZ2 (b) Figure 12.11: Spatial coherence of 2 AHN3 tiles. The inner cell colour indicates the position in the stream of first point in that cell, and the outer cell colour indicates the position in the stream of the last point in that cell. 12.4 Notes and comments The description of the 𝑘 d-tree and the nearest neighbour query is adapted from Wikipedia ( https://en.wikipedia.org/wiki/K-d_tree ) and the lecture notes entitled “kd-Trees—CMSC 420” from Carl Kingsford (available at https://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/ kdtrees.pdf). Vitter (2001) provides an overview of external algorithms. Agarwal et al. (2005) construct massive TINs by designing external algorithms, and Arge et al. (2006) and Agarwal et al. (2008) have implemented spatial analysis functions on TINs based on that paradigm. The streaming computation of the DT algorithm is a simplification of the algorithm described in Isenburg et al. (2006a), and some figures were basically redrawn. Isenburg et al. (2006b) explains in details how large rasters can be constructed with spatial interpolation by modifying the streaming pipeline of Figure 12.7. 12.5 Exercises 1. The tree in Figure 12.4 is balanced, but if ( 1,3 ) had been selected as the root, how would the tree look like? 2. Real-world point cloud datasets often have natural spatial coherence. Explain why that is the case for lidar datasets. 3. How to construct a 𝑘d-tree that is as balanced as possible? 4. “The ideas behind streaming are very useful for certain local problems, but unfortunately they cannot be used directly for global problems such as visibility or flow modelling”. Explain why that is with a concrete example. Spatial extent of a set of points 13 uyoutu.be/dqRb32iGZ_c 13.1 Properties of the region 142 13.2 Convex hull . . . . . . . . 143 13.3 Moving arm . . . . . . . . 143 13.4 𝜒-shape . . . . . . . . . . 145 13.5 𝛼-shape . . . . . . . . . . 145 13.6 Clustering algorithms .146 13.7 Notes and comments . . 147 13.8 Exercises . . . . . . . . . . 147 Given a point cloud, one operation that practitioners often need to perform is to define the spatial extent of the dataset. That is, they need to define the shape of the region that best abstracts or represents the set of points. As seen in Figure 13.1, this region is often in two dimensions, for example in the case of an aerial lidar datasets we may want to know where the ground is (after removing the points on the water), or in the case of the scanning of the façade of a building, we would like to obtain a polygon that represents where the wall is (omitting the windows). Another example is if we consider a subset of a point cloud that represents the points for a given building, we can interested in creating the polygon that represents its footprint (Figure 13.1c). Calculating the spatial extent is useful to calculate the area covered by a dataset, to convert it to other formats (eg raster), or to get an overview of several datasets it is faster to load a few polygons instead of billions of points, etc. The spatial extent is often called by different names, for instance: envelope, hull, concave hull, or footprints. It is important to notice that the spatial extent is not uniquely defined and that it is a vague concept. As Figure 13.2 shows, there are several potential regions for a rather simple set of points, and most of these could be considered ‘correct’ by a human. In this chapter we present methods that are used in practice to define the spatial extent of a set of points in ℝ2 , which implies that the points in a point cloud are first projected to a two-dimensional plane. (a) (b) (c) Figure 13.1: Three point cloud datasets for which we would like to find the spatial extent. (a) An aerial point cloud with several canals (dark colour). (b) A scan of a façade containing several windows. (c) For the subset formed of the building footprint (green points) we would like to obtain its footprint (line in red). 142 13 Spatial extent of a set of points Figure 13.2: Different methods to obtain the spatial extent of a given set of points in the plane. (a) A set of points in ℝ2(b) Its convex hull (c) A𝜒-shape (d) An 𝛼-shape 13.1 Properties of the region Let 𝑆 be a set of points in ℝ2 , and 𝑅(𝑆) the region that characterise the spatial extent of 𝑆 . The region is potentially formed by a set of polygons (if 𝑆 forms two distinct clusters for instance),and in practice most algorithms will compute a linear approximation of 𝑅(𝑆) , so the polygons will have straight edges as boundaries. To evaluate the different algorithms to create R( 𝑆 ), we list here different properties that one must consider when defining the spatial extent of a set of points. P1. Regular polygons? Are polygons allowed to have dangling parts (lines), such as the one in Figure 13.3b P2. All points part of the region? Can outliers be ignored? Or do they have to be part of the region? In Figure 13.3a and Figure 13.3b they are all part of the region, in Figure 13.3c and Figure 13.3d one outlier is not. P3. Region is one connected component? Or are more components allowed? In Figure 13.3a–Figure 13.3c there is one component, but Figure 13.3d has two. P4. Are holes allowed in a polygon? Polygons in Figure 13.3a–Figure 13.3c have only an exterior boundary, while in Figure 13.3d one polygon has an interior boundary too (a hole). (a) (b) (c) (d) Figure 13.3: Different properties for the spatial extent 13.2 Convex hull 143 a b c d e (a) (b) Figure 13.4: (a) First four steps of the gift wrapping algorithm to compute the convex hull. (b) The resulting convex hull. a b c d e l=10m Figure 13.5: First four steps of the moving arm algorithm (witha length 𝑙 )tocompute the spatial extent. a c b d e Figure 13.6: First four steps of the moving arm algorithm (with a knn where 𝑘=3 ) to compute the spatial extent. P5. Computational efficiency What is the time complexity of the algorithm, and does it require large and complex auxiliary data structures? 13.2 Convex hull As explained in Section 4.2.1, given 𝑆 , a set of points in ℝ2 , its convex hull, which we denote conv( 𝑆 ), is the minimal convex set containing 𝑆 . Two examples of convex hulls are in Figures Figure 13.2b and Figure 13.3a. For a given set of points, the convex hull is uniquely defined and does not require any parameters (unlike the other methods listed below). It is also relatively easy to compute: it can be extracted from the Delaunay triangulation, or computed directly using a specific algorithm. An example of the latter is the well-known gift wrapping algorithm, shown in Figure 13.4. It begins with a point that is guaranteed to be on conv( 𝑆 ) (we can take an ‘extreme’, such as 𝑎 in Figure 13.4, because it is the point with the lowest 𝑦 -coordinate), and then picks the point in 𝑆 (omitting the ones already on conv( 𝑆 )) for which the polar angle between the horizontal line and that point ( 𝑎 at this step) is the largest ( 𝑏 in this case), and adds it to conv( 𝑆 ). Then for 𝑏 , the polar angle is calculated from the line 𝑎𝑏 and the 𝑐 is chosen since it forms the largest angle. The algorithm continues this way until 𝑎is visited again. If 𝑆 has 𝑛 points and conv( 𝑆 ) is formed of ℎ points, then the gift wrapping algorithm has a time complexity of O(𝑛 ℎ) ; each of the ℎ points are tested against all 𝑛 points in 𝑆 . However, there exist more efficient algorithm that have a time complexity of O(𝑛log 𝑛). Properties convex hull: P1. The sole polygon is guaranteed to be regular (and convex). P2. All points are on or inside the region. P3. One component. P4. No holes in the region. P5. O(𝑛log 𝑛) 13.3 Moving arm Arm of length 𝑙 .Themovingarmisageneralisationofthegift wrapping algorithm (see Section 13.2) where the infinite line, used to calculate the polar angles, is replaced by a line segment of a given length 𝑙 (the “moving arm”). This means that, unlike the original gift wrapping algorithm, only a subset of the points in 𝑆 are considered at each step. This also means that potentially the result is a polygon that is non-convex. Figure 13.5 shows the first few steps for a given 𝑙 , and it can be observed that 1 point is not part of the final region. Observe also that if 𝑙 had been larger then conv(𝑆) could have been obtained.