UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Terrain drainage features and queries Yu, Sidi 1996

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1996-0501.pdf [ 2.21MB ]
JSON: 831-1.0051298.json
JSON-LD: 831-1.0051298-ld.json
RDF/XML (Pretty): 831-1.0051298-rdf.xml
RDF/JSON: 831-1.0051298-rdf.json
Turtle: 831-1.0051298-turtle.txt
N-Triples: 831-1.0051298-rdf-ntriples.txt
Original Record: 831-1.0051298-source.json
Full Text

Full Text

T E R R A I N D R A I N A G E F E A T U R E S A N D QUERIES by Sidi Yu B.S., Lock Haven University of Pennsylvania, 1993  A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E D E G R E E OF M a s t e r of Science in T H E F A C U L T Y O F G R A D U A T E STUDIES (Department of Computer Science)  We accept this thesis as conforming to the required standard  The University of British Columbia June 1996 © Sidi Yu, 1996  In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  Computer Science The University of British Columbia 2366 Main M a l l Vancouver, Canada V 6 T 1Z4  Date:  Abstract  Terrain drainage characteristics are of interest to a number of fields such as hydrology, hydraulic engineering, natural resources management, flood control, environmental sciences, geographic information systems, etc. The traditional way of obtaining terrain drainage characteristics information is to have human operators visually identify and delineate terrain features of interest from aerial photography or.contour maps. The disadvantages of this mode of operation are obviously the tediousness, the inefficiency of the process, and the inaccuracy of the information thus obtained. Over the past twenty years, considerable research effort has been devoted to automatic generation of terrain drainage information using digital computers. Significant progress has been observed in this field. Yet most of the research done in this area does not share common definitions of terms. Some publications have attempted to address this problem by proposing general definitions that are not tied to a particular data representation or algorithm. We perceive this as the direction to follow for all future research in this area. This thesis sets out to further develop the general definitions of terrain drainage characteristics.  Particular terrain features such as pits, peaks, water-  ii  courses, ridges are defined, and methods for identifying them are presented. Suggestions for data structures and algorithms are also given. The methods proposed not only solve the standard problems, but answer additional interesting queries as well. Proofs of computational complexity are given for the more important algorithms and problems.  iii  Contents Abstract  ii  Contents  iv  List of Tables  viii  List of Figures  ix  Acknowledgements  xi  Dedication 1  xiii  Introduction  1  1.1  Models and Approaches  2  1.2  What Problems Can Be Addressed Using Spatial Analysis  3  1.3  Digital Terrain Representations  4  1.3.1  6  1.4  1.5  Arguments for TINs  Characteristics of Our Approach  8  1.4.1. Implicit vs. Explicit Surface  8  1.4.2  9  Global vs. Local Processing  How Real is the Model  10 iv  2  3  Previous Research 2.1  In D E M s  12  2.2  In TINs  14  2.3  Topological Structure  3.2  15 16  Assumptions: How the Water Flows  16  3.1.1  Assumptions by Other Researchers  16  3.1.2  Our Assumptions  17  Land Forms  18  3.2.1  General Drainage Definitions  18  3.2.2  Local Drainage Features  18  3.2.3  Global Drainage Features  19  3.3  Terrain Partition  23  3.4  Graph Structure of Global Features  24  Extracting Terrain Features in TINs  26  4.1  Local Features  26  4.2  Definitions Particular to TINs  27  4.3  Global Drainage Features  29  4.3.1  Drainage Network in TINs  29  4.3.2  General Separating Lines in TINs  31  4.4  5  .  Definitions 3.1  4  12  Dubious Dualities  31  4.4.1  31  Dualities in a Terrain  Data Structure and Complexity of Drainage Processing 5.1  T I N Data Structure  33 34  v  6  5.2  Local Terrain Features  34  5.3  Global Terrain Features  35  5.3.1  Processing Complexity Is the Same as Output Complexity . .  35  5.3.2  Upper Bound on Drainage Network Complexity  36  5.3.3  A Worst-case Example . .  38  5.3.4  Average-case Complexity  39  5.3.5  Terrain Planar Graph  Drainage Queries 6.1  7  .  40 41  Simple Queries  41  6.1.1  41  Finding Basin Boundaries  6.2  Finer Partition  42  6.3  Querying the Watershed of a Point .  42  6.4  Querying the Size of the Watershed of a Point  44  6.4.1  Definitions  44  6.4.2  Lemma and Corollaries  45  Conclusions 7.1  48  Level of Abstraction  -48  7.1.1  Advantages  48  7.1.2  Disadvantages  49  7.2  Reflection of Reality  50  7.3  Future Research  50  7.3.1  Making the Model More Realistic  51  7.3.2  Data Structures  51  7.3.3  Handling Special Cases  . . .  vi  51  7.3.4  Adding Other Factors  7.3.5  Merging Data Sources  52 . .  52  Bibliography  53  Appendix A Terrain Function  58  Appendix B Horizontal Slicing Plane Method  60  vn  List of Tables Observed Drainage Network Complexity  viii  40  List of Figures 3.1  Drainage Forest  20  3.2  Three Types of General Separating Lines within a Basin  21  3.3  A Non-pass Point  22  3.4  Terrain Graph  25  4.1  Terrain Feature Around a Vertex  27  4.2  Drainage and Non-drainage Edges  28  4.3  Local Channel  28  4.4  Three Types of Edges  29  4.5  Two Types of Watercourse Merge Points  30  4.6  Dividing Line Coincides with Watercourse  30  4.7  Breaking Down of Dualities  32  5.1  Repeating a Triangle Repeats the Subsequent Sequence  37  5.2  Worst Case Example for Drainage Network Complexity: the Pyramid  38  6.1  New Trace-up Lines  43  6.2  Query Watershed of a Point  43  6.3  Strips  45  ix  6.4  Sweeping Trace-up Line  46  6.5  Interpolation Points  47  A. l  Sinusoidal Directional Derivative .  59  B. l  Horizontal Plane of a Pit  61  x  Acknowledgements I would like to express my special thanks to my thesis supervisor D r . Jack Snoeyink. I am very grateful for the fruitful discussions with him, the constructive guidance from him, and just the fun of working with him. I can always get my drafts back, marked with lots of comments, from him sooner than I dare to hope. He has offered materials directly to this thesis: Theorem 1, its proof, and the worst case example; the formal presentation of the concept "strip" in querying sizes of watersheds of drainage network points in Chapter 6; proof of terrain partition; Figure 5.1, Figure 5.2, Figure 6.3, Table 5.1, etc. Many thanks to D r . Dave Kirkpatrick for squeezing time to be the second reader of this thesis. I would also like to thank Dr. Kelly Booth for giving me much valuable advice and guidance in the Virtual Hand Project. Many friends and mentors have made my pursuit easier and possible. Here is my salute to them all for their generous help and support.  SIDI Y U  XI  The University of British Columb June 1996  This thesis is dedicated to my mom, who's the best.  xiii  Chapter 1  Introduction  Where does water flow on a terrain? If we know the underlying dry surface of the terrain, are we able to know where the streams and rivers lie? Terrain drainage characteristics include watercourse network, ridges, basins, lakes, peaks, hills, valleys, pits, watersheds, river volume, basin outlet (lowest point on the basin boundary), etc. They provide important information on water resources, possible flood areas, erosion and other natural processes. This information can be used in a number of fields such are hydrology, hydraulic engineering, natural resource management, etc.  1  1.1  Models and Approaches  Let's revisit our opening question: where does water flow on a terrain? The answer, in the physical world, depends on a variety of factors: gravity, viscosity of water, sediment concentration, undersurface water, permeability of the soil, vegetation, physical laws such as dynamics, and, of course, the terrain surface geometry. Thus, to answer the question "correctly" requires an all encompassing model of the entire physical process. First of all, that model has yet to be proposed; secondly, with the computational resources we have, we can only settle for a less ambitious approach, namely to ignore some of the factors, or to make assumptions about them and proceed with what we believe important. We have to choose the right kind and level of abstraction for the problem. Besides the mathematical model we use, we need to decide on the means of processing. Traditionally, human operators either visually inspect aerial photography to decide where rivers, ridges and basin boundaries lie; or they have a contour map at hand, and estimate terrain features from contour information; or they combine the two pieces of information to get better results. Manual quantification of terrain drainage characteristics is tedious and time consuming, and inevitably subject to inaccuracies. A lot of research has been dedicated to automatic generation of terrain drainage information using computers. Most algorithms use terrain surface geometry as the main data source in their processing. This is due to a number of factors: 1. Availability  of data: many government geographic survey agencies, such as  USGS (The United States Geological Survey), have been surveying terrain elevations at densely placed sampling locations (for instance, on a grid with 30-  2  meter inter-grid spacing). Such information has become available in the public domain in large volumes. Other types of information, such as vegetation, soil absorption, are harder to obtain. 2. Importance of geometric data: among all factors that determine terrain drainage, surface geometry is probably the most important. 3. Computational cost: a model that incorporates other factors, such as erosion, can be expensive computationally. The model itself will be complex, and the computational cost will be prohibitive. Considering all the above factors, we will stick with the geometric approach in this thesis. Our raw data will be the widely available digital elevation data. We call our approach spatial analysis of digital representations of topographic surfaces. By "spatial analysis" we mean the analysis is based solely on the geometric information contained in the digital terrain representations.  One of the goals of this  thesis is to outline certain sound choices in carrying out such spatial analysis which, we believe, help generate more consistent and accurate determination of drainage features.  1.2  What Problems Can Be Addressed Using Spatial Analysis  Without formally defining drainage network features (they are defined in later chapters), we claim that the following questions can be answered by spatial analysis: • Standard Problems: We can identify the drainage network (network of watercourses), network of ridges, pits, peaks, passes (saddle points). We can delin3  eate basins and hills (the duals of basins). With formal definitions, we can obtain a consistent graph representation of the overall structure of drainage features. • Non-standard Problems: With a richer representation, an additional set of problems can be solved correctly and efficiently. The watershed of a terrain surface point is the set of terrain surface points that contribute to the flow of that point. Watersheds of watercourses completely partition the terrain. By precomputing watersheds of certain critical points on the terrain, we can answer many queries efficiently, such as: what is the watershed of any point on the terrain; what is the watershed of a finite area; how much is the flux into a point on the terrain (or the size of that points watershed); where is the outlet of a basin, what is the volume of a basin, thus how much flood will come out from the basin if we know the level of precipitation.  1.3  Digital Terrain Representations  To perform spatial analysis correctly and accurately, it is essential to use an appropriate digital terrain representation, because the representation's properties have a bearing on what kind of information we can get and how structured the information is. Thus to describe the characteristics of our approach, we need to discuss the' different forms of digital terrain representations. "Digital" here refers mainly to the fact that the sample points on the terrain surface are discrete as opposed to continuous. Digital terrain representations contain only surface points. There are three main forms of digital terrain representations. In these representa-  4  tions, each surface sample point has an x-y location and an elevation z. • digital contour: represents the surface as a set of contours. Like traditional contour maps, each contour is of a particular elevation. There are sample points along the contour in the form of a sequence of x-y locations. There is an implied adjacency relation between neighboring x-y locations. • D E M (digital elevation model or regular grid): the sample points are regularly spaced in their x-y locations, like a regular grid. When storing them in a computer, the memory adjacency relations can be used to imply the spatial adjacency relations and orientation relations. Thus, we only need to explicitly store the elevations. • T I N (triangulated irregular network): uses a triangulation to represent the terrain surface. The x-y locations of the sample points are irregularly spaced. A triangulation of the x-y locations is performed, and that triangulation defines the adjacency relations between points. The most popular scheme for the triangulation is Delaunay Triangulation. Delaunay Triangulation maximizes the smallest angle of the triangles, and thus avoids long and skinny triangles. The resulting triangulation with triangles all having finite and local extent gives better approximations to the actual terrain surface. The triangulation surface with all the points elevated to their original heights becomes the T I N surface, expressed as a function z = f(x, y). It is continuous, but not different i a t e at the edges and vertices of the triangulation. To store a T I N , we need to store the x-y locations, elevations and adjacency relations. Due to processing efficiency and low complexity, the latter two representations are more popular. Virtually all current GIS systems support D E M s ; more and more are 5  supporting TINs and the conversion between the two. A critical study by M . P . Kumler [19] compares the relative strengths of the two representations. The author maintains that for the selected representative terrains in the United States, D E M s are more accurate, given the same amount of storage. TINs, because of its need for storing x-y locations and adjacency relations, are less space efficient. In general, for the same amount of storage, a D E M stores 10 times as many points as a T I N . The study uses a particular set of procedures in acquiring T I N vertices (there are many approaches to acquiring T I N vertices), thus the author contends that the study is not a conclusive statement on the relative merit of the two representations. The author also indicates further study is necessary to determine which of the two gives better estimates on the derivative statistics such as slope and aspect. These derivatives are important in the context of determining terrain drainage.  1.3.1  A r g u m e n t s for T I N s  • Data Acquisition From data acquisition point of view, since raw data generally are not collected at regular grid points, the grid point elevations in D E M s must be obtained by some kind of interpolation scheme. In TINs, the irregularly acquired points can be used as vertices directly. Furthermore, locations of certain important points such as pits, peaks and saddle points on the terrain can be represented exactly in TINs, but not in D E M s . . .. • Storage. Efficiency Traditionally, an argument for TINs. is that they adapt to the surface complexity better than D E M s : for instance, at flat areas, few triangles are needed in TINs, while, because of the uniformity, a large number of grid points are used in DEMs(the traditional fixed grid D E M s ) anyway. By the 6  same token, a rugged area, an area of high relief can be represented accurately by a large number of T I N vertices, which is not possible for D E M s . This is referred as variable resolution of the TINs. Nevertheless, the disadvantages of TINs lie mainly in the inefficiency of storage (of the x-y coordinates and the adjacency relations), and the consequent processing inefficiency (for example, having to calculate orientation between triangles in certain queries). • Richness of Information Kumler's paper[19] compares the two representations mainly by the criterion of efficiency defined as the elevation accuracy of the representation per storage space used. By accuracy he means the absolute deviations of the model surface and the actual measured elevations at many test points on the studied terrains. • Simplification Possibilities For D E M s , simplification entails adjusting the uniform grid size and resampling with some interpolation scheme. There is little we can do to exploit the local characteristics of the terrain in that process. For TINs, however, since we can control the x-y locations of the vertices, we can sparsely sample terrain regions where little information is lost in so doing, and densely sample regions where many vertices are necessary to represent the complexity of the surface. In essence, simplification amounts to choosing another set of parameters (in this case, vertices) to represent the surface. This is the classic resampling problem. TINs' superiority over D E M s in this respect arises from the non-uniformity of vertex selection, which affords them a more flexible, powerful and localized resampling scheme. The comparison between D E M s and TINs is similar to that between raster and vector representations in image processing and graphics. D E M s are simple yet at 7  times imprecise; TINs can be complex yet exact.  1.4  Characteristics of Our Approach  Our approach to the drainage problems has a number of characteristics.  1.4.1  Implicit vs. E x p l i c i t Surface  We call our approach an explicit surface approach, or defined surface approach as opposed to implicit surface or undefined surface approach. We use an explicitly defined continuous surface (the T I N triangulation surface) on which water flows. Many proposed methods [14, 33, 29, 7, 41] do not have an underlying surface, and that makes them mathematically unsound.  For example, a number of methods  based on D E M s treat each D E M grid point as a cell, and do not have the notion of a surface. In our method, with an explicit surface, many terrain features and processes become well defined and easily extracted. In a large number of previous approaches, a terrain feature is simply defined as the output of a particular algorithm. This makes it impossible for different researchers to reach the same result, and there is no objective measure of correctness. Several researchers have realized this difficulty, and tried to remedy it[10, 17, 46]. We follow the suggestions made by Frank et al. [10] that formal definitions should be used to define terrain-specific features so that properties of the structure can be established mathematically and contradictory definitions (or at least those in • disagreement with-other related-research) can be-avoided. rOur definitions, apply, to any.piecewise planar surface. This includes T I N surfaces and D E M surfaces if a piecewise planar interpolated surface is specified. This mathematical approach provides a common ground for discussions of terrain features. 8  Furthermore, by following a mathematical approach, we are able to generate consistent global structures for the terrain.  1.4.2  G l o b a l vs. L o c a l Processing  Many approaches, particularly the ones based on D E M s , generate drainage information by performing only local operations. They are similar to raster image processing, where you apply local operators on a image and then use line-connection or region-growing algorithms to generate global structures. The hope is that global structures will manifest themselves locally everywhere they lie, and turn into local structures, so they can be detected locally. Unfortunately, many global drainage features do not always reveal themselves as local structures; furthermore, the same type of local features may belong to different global structures. Thus, the local approach often results in broken or ambiguous structures.  •  We prefer to generate global structures, such as watercourse network, basin boundaries and watersheds, by performing global operations, and hope that the global operations are not too expensive computationally. Fortunately, with our definitions, it turns out that the operations to generate global structures do no extra work, because linear-complexity local operations give us all and only the necessary starts for global structures, and the tracing of the global structures is a linear operation to their own size. In natural terrains, these global features typically have linear complexity. This fact makes our approach both "correct" and efficient.  <j  9  1.5  How Real is the Model  Let's return to our discussion of level of abstraction and representativeness of the model. Our approach assumes that water flows on the terrain surface; certain water bodies, such as streams and rivers, have no depths or widths, thus they are onedimensional mathematical objects; streams merge, but never bifurcate. There are three assumptions in this model that are at variance with the physical world: 1. The sample points we get from the raw data are not on the dry surface underneath water bodies. If the sample points are on a river or a lake, they are sampled on the surface of these water bodies. Thus to say the given surface is the dry surface water runs on is not accurate. But this discrepancy is largely due to the limitation of our data source, and we have to live with it. 2. Real water bodies have widths and depths, so they are three-dimensional rather than one-dimensional. Our model currently treats rivers as one dimensional objects. A 3D model will be much more costly computationally, and hard to define. A n d the benefit of this complexity may not be relevant in the context of huge amount of data, which is the usual case for our purposes. 3. Streams and rivers, as we know them, in the real world, do bifurcate. We do not allow that in our model, because it will turn the river network from rooted forest, in graph theoretic sense, into network. That will make the structure and the algorithm more complex, and could potentially make the computational cost too high. Furthermore, such features are more prone to.inaccuracies and artifacts of the input data, thus may be better left out, at least for the scope of this thesis.  10  Artifacts in the data can cause problems in processing. For example, many small pits in the data are the result of inaccurate measurement of elevations. Given the precision of the survey data, such artifacts are inevitable. Our model assumes that the source data have already been put through some filtering operation, and most artifacts are removed. We will have more discussion on this in the next chapter.  11  Chapter 2  Previous Research We divide previous works in extracting terrain drainage features into two categories based on the digital representations they use, namely D E M s and TINs.  2.1  In DEMs  Methods based on D E M s do not have an explicit surface that water runs on, so typical models of how water flows bear little resemblance to the physical process. For example, an early approach proposed by Peucker [33] has no concept of surface. Many approaches [14, 15, 33, 29, 7, 36, 20, 41] based on D E M s extract features use techniques similar to those in raster image processing: local operators and region growing, etc.  Typically, local operators are used to detect possible pits, peaks,  channel cells (cells that have the potential to be channels), etc., then the channel cells are linked according to some heuristics to form the drainage network. Douglas [7] describes several methods for detecting channel (ridge) cells. Basin boundaries are identified in a similar fashion[14]. As an example, a channel is usually detected with a local operator similar to a 12  line detector in image processing. The idea is if a grid point and its 8 neighbor grid points form a V shaped local surface, we say the grid point is on a channel. One then attempts to use some heuristic line-connection algorithm to connect these channel grid points into a linked channel line. As with image processing, there are inherent ambiguities with this raster approach, which emphasizes local characteristics. The global nature of the terrain drainage is not adequately reflected, and consistent global structures are hard to obtain. Different heuristics result in different sets of extracted terrain features. Mark [22] proposed a drainage network detection method that does account for the global nature of drainage. Each grid point represents a cell, and the flow behavior of the grid point represents the flow behavior of the entire cell. The water on a cell goes to one of the eight neighbors: the one that forms the steepest down slope with the current grid point. This last assumption is not so unreasonable per se, it is similar to our assumption of no bifurcation, which keeps the complexity tractable. But the combination of the these assumptions is unrealistic. The algorithm goes as follows: sort all the cells in order of decreasing elevations. Each cell has water quantity one. For each cell in that order, decide to which of its eight neighboring cells the water will flow (the one that forms the steepest down slope with the current cell). A d d the current cell's water quantity to that cell. In the end, cells with water quantity over certain level is claimed to be on the water course. These cells are linked naturally to form the drainage network. This method takes into account the fact that drainage process has a global nature. Because of that, it is more costly: the sorting process is of order O ( n l o g n ) , which is higher than the local methods' O(n), where n is the number of grid points. These two types of D E M methods have the same complexity during their second  13  phase of processing: the line-connection phase. The complexity for the second phase is proportional to the size of the drainage network. In our approach the complexity is proportional to the resulting drainage network too. The difference lies in the fact that in these approaches, the drainage network complexity is bounded by the number of grid points, because they do not introduce new points. But in our approach, we introduce new points, thus the complexity can be potentially high. Fortunately, as we will see later, in natural terrains, the complexity is more likely to be linear or sub-linear in the input size. Takahashi, et al. [41] proposed a method that first extracts critical points (pits, peaks and passes) that satisfy Euler's relation, then traces channels and ridges from those critical points. This method takes advantage of the fact that there is a global structure in drainage features.  2.2  In TINs  There are fewer papers dealing with extracting terrain drainage features from TINs [30, 16]. In TINs, we do have an explicit surface. In [30], what later we define as local channels are extracted, then joined to form the drainage network, although how to do the joining is not explained. To partition the terrain into tributaries of channel segments, the terrain is subdivided into cells (not dissimilar to Mark's grid cells), all the water in a single cell flows into a single watercourse. In this paper, we will outline explicitly how to do the joining, and discard the notion . of cells; because in reality; water on. different parts of a cell can flow into different . watercourses. In [16], the drainage network is formed by tracing up local channels, which, as we will see in detail later, can not account for all the watercourses, because watercourses 14  include not only local channels. Tracing up local channels doesn't work when local channels do not connect. In the paper, when delineating tributaries, tracing up along the steepest accent is employed, but details are lacking. We will show that because water flows down, its more natural and complete to trace down instead of up. Tracing down gives us the complete set of watercourses. We will also explain explicitly how the trace-downs and trace-ups are performed.  2.3  Topological Structure  In [7], a new type of digital elevation model is proposed which extracts important ("information rich") terrain features such as ridges and channels, and use them to form a "richline" representation of the terrain. This is in the same spirit as the works on graphs of surface networks [34, 37, 41, 46]. Pfaltz [34] defines a graph structure with vertices at critical points (pits, peaks and passes) and edges at ravine and ridge lines. Takahashi et al. [37] give a different graph structure. We are interested in the quantitative aspects of terrain drainage, yet these graphs represent more of the qualitative ones. Nevertheless they provide an economical way of representing the overall terrain structure. The construction of these graphs will be addressed in this thesis.  15  Chapter 3  Definitions We have given justifications to some of the assumptions and choices we made in our approach. In this chapter, we present these assumptions and definitions in detail. These definitions are applicable to any kind of surfaces, not limited to T I N surfaces or piecewise planar surfaces. Their realization in TINs will be discussed in next chapter.  3.1 3.1.1  Assumptions: How the Water Flows Assumptions by Other Researchers  Defining how water flows is essential to computing terrain drainage. In raster D E M s , as proposed by Douglas [7] and Mark [22], a network of discrete point elements is used. In [22], each element represents an entire cell of the grid. Then it is decided to which of its four or eight neighbor elements the water on this element flows.' Thus there are two levels of discretization:' one in the use-of cells, the other in limiting the possible directions water can flow. Douglas discarded the  16  second discretization by allowing water to flow along the steepest descent instead of to the eight neighbor elements. To have an explicit surface for the steepest descent calculation, Douglas inserts a virtual grid point at the center of each square cornered at four neighboring grid points. The elevation of the center point is the average of the four corner points. Then four triangle surfaces each consisting of two points of an edge of the square and the center point are used in the steepest descent calculation. Nevertheless, Douglas's method still keeps the concept of cell.  3.1.2  O u r Assumptions  Assume we have a continuous surface which can be modeled as a function z = f(x, y) over the x-y plane, and that: 1. Water flows on the surface. 2. A t any point on the surface, the water follows the steepest descent at that point . 1  3. A t any point on the surface, the steepest descent is either unique or nonexistent . 2  Because the steepest descent is unique, we obtain a tree structure of watercourses. The uniqueness of steepest descent forbids bifurcation. As pointed out by Poiker in Chapter 39 of the N C G I A Core Curriculum, this assumption is also the basis of the traditional definition of watershed in hydrology. 'For the precise definition of "steepest descent", see the Appendix. We leave special cases such as flat surfaces out of our present discussion. They can be accounted for by assigning a steepest descent by some criteria. 2  17  To account for some of the special cases, in our processing, we may desire that local depressions be filled to form lakes and that their outlet points be identified, so that water can pass through the depression and flow down from its outlet, and find its final destination. Smoothing operations can be performed to remove small depressions. O'Callaghan and Mark [29] reported 90% removal rate for one pass of a simple binomial smoothing operator. For larger depressions, explicit methods have to be used to fill them.  3.2  Land Forms  In this section, we give definitions of important land forms.  3.2.1  General Drainage Definitions  • Trickle path of a point p: the path that a drop of water follows starting at p. The destination point is either a pit, which is a local minimum, or a point on the boundary of the terrain. • Watershed of a point p: the area formed by all the points whose trickle paths go through or end at p. A l l and only the water falling on the watershed of p flows through or into p. • Basin: the watershed of a pit.  3.2.2  L o c a l Drainage Features  • Pit: local minimum of the surface. • Peak: local maximum of the surface.  18  3.2.3  G l o b a l Drainage Features  Definition of Drainage Network(or Watercourses) Drainage network is the network of streams and rivers in a terrain. Formally, • Drainage network: the set of all points whose watersheds have nonzero projected area onto the x-y plane.  In other words, any point that has a  measurable amount of water flowing through it is on the drainage network. A watercourse is a monotonically descending branch of the drainage network. • Watercourse origin: the limit of drainage network points whose watershed sizes approach zero.. Because any point's trickle path is unique and trickle paths do not cross, the drainage network is a forest in graph theoretic sense, see Figure 3.1. Each tree in the forest is rooted at a pit. The nodes (vertices in graph theory) in the forest are watercourse origins (called leaves) and merge points of watercourses (called interior nodes). The portion of the watercourse connecting two nodes is a, link (edge in graph theory). Because the interior nodes are defined as merge points, they all have degrees greater than 2 (except for possibly pits). A parent node's watershed includes the watersheds of all its descendents' watersheds.  Compliment of Watercourses In a terrain, there are features that form the boundaries of watersheds. These are: dividing lines, crests, and watercourse trace-ups.  19  Figure 3.1: Drainage Forest Dividing Lines If an interior node A on the drainage network is not a pit, then one and only one of its adjacent nodes has elevation lower than A (it is the parent of A). Now we examine the projections (on the x-y plane) of all the incident links of A. If two consecutive (counting clockwise or counterclockwise) links lead to two nodes of higher elevations, ie. if two watercourses merge at A, we define a separating line between the two watercourses by tracing up along the maximum ascent between them starting from the node and ending at a peak. We call this separating line a dividing line (see Figure 3.2).  • Dividing line: a separating line between two merging watercourses. It starts at the merge point of the watercourses, traces up along the steepest ascent and ends at a peak.  20  A dividing line has the property that no measurable amount of water from one side flows to the other side by crossing any point on the dividing line.  watercourse trace-up dividing line  general separating line  watercourse  Figure 3.2: Three Types of General Separating Lines within a Basin  Crests A crest is what is generally called a ridge line between two peaks. It starts at a pass (saddle point), and ascends in two distinct directions following steepest ascents, and reaches two peaks, which need not be distinct. • Pass: any non-pit point p that has at least two maximum ascents A and B, 3  and there is no measurable amount of water that flows from one side of A-B to the other through p. It can be proved that the latter condition is satisfied if and only if the following is false: at p, there is an incoming watercourse on one side of A-B, and the direction of steepest descent is on the other side of A-B, see Figure3.3. A t a pass, it's also true that there must be at least one minimum descent on each side of A-B, otherwise the water falling on one or both sides of A-B at For the precise definition of maximum ascent, see Appendix.  21  Figure 3.3: A Non-pass Point the vicinity of p will have to flow across A-B. • Crest: a line consisting of two trace-ups that each originate from the same pass, ascend up along the steepest ascent, and end at a peak. No measurable amount of water from one side of the crest flows to the other side by crossing any point on the crest. A pass is the lowest point on a crest. We soon will see that crests are important division lines between basins.  Watercourse  Trace-ups  • Watercourse trace-up: is a line that starts at a watercourse origin, traces up along the steepest ascent, ends at a peak. Watercourses trace-ups are used to connect watercourse origins to crests and eventually to peaks. They complete the partition of the terrain. Because dividing lines, crests, and watercourse trace-ups all have the function of separating watercourses, they can be classified into the category of general separating line. • General separating line: a dividing line, a crest or a watercourse trace-up. 22  3.3  Terrain Partition  Watercourses and general separating lines form the boundaries of a complete partition of the terrain. • source area: a contiguous area of the terrain surface bounded by a portion of a watercourse and a portion of a general separating line which has no portions of other watercourses or general separating lines lying inside it. Lemma 1 Source areas partition the terrain. All and only the water falling on a source area flows into its bounding watercourse portion directly. Furthermore, the bounding watercourse portion descends monotonously from one end to the other. When we say the water flows "directly" into a watercourse, we mean the water makes its initial entry into the drainage network at that watercourse. A pair of source areas arise when two watercourses merge: the dividing line together with each of the two watercourses bounds one source area. The upper end of the watercourse is a watercourse origin. The lower end is a watercourse merge point. Therefore, each source area is related to a leaf node and an interior node of the drainage network forest. Each leaf is related to two source areas. Each interior node is related to two or more source areas. Each source area is part of the watershed of its related interior node. This fact is essential to using partition to solve drainage problems. It's possible that a small basin can lie completely within a source area spatially, but it is not part of the source area. While using the partition to solve drainage problems, we should not make the opposite assumption. We observe that when two dividing lines that bound the two sides of a watercourse trace up to a single peak, the boundaries for the two source areas are closed. But  23  when the two dividing lines end up at two distinct peaks, we do not have closed boundaries. This is where crests' role comes into play: there exists some crest that completes the boundaries. Conversely, if there is a crest, some watercourse will need it to complete its two source areas' boundaries. Therefore the work of finding crests are part of the work of finding the boundaries of source areas.  3.4  Graph Structure of Global Features  There is a planar graph structure coming out of the partition. The crests partition the terrain. A basin is bounded by crests. A crest connects two summits. Each summit is a vertex in graph theoretic sense. A crest is an edge that connects two vertices. A basin is a face in the graph. We can obtain a dual graph of the above described graph. Each basin has a pit in it, which becomes a vertex of the dual graph; watercourses are the dual graph's edges, they connect pits; hills are faces, which are bounded by watercourses. See Figure 3.4. Two basins are adjacent if they share a common crest. A watercourse and a crest meet at a pass. This graph is a compact representation of the overall structure of the terrain. Douglas's "richline" representation [7] is precisely an attempt on using this structure to represent the most important features of the terrain. This graph structure arises as a result of our definitions. Generating this graph takes little extra processing after we obtain all-the watercourses and general separating lines.  24  cresl  watercourse  o  peak  •  pit  Figure 3.4: Terrain Graph  25  Chapter 4  Extracting Terrain Features i n TINs We now apply the definitions we presented in the previous chapter to TINs. We start by outlining extraction methods for local features in TINs, then we will introduce new definitions that only apply to TINs. Most of these definitions and methods apply to any piecewise planar surface. We also assume that there are no degenerate cases such as flat areas or horizontal edges. These are special cases that can be handled by assigning steepest descent direction on the basis of other information such as flow direction, medial axes of a lake, etc.  4.1  Local Features  Local features such, as pits, peaks, passes can be easily extracted by examining vertices and their adjacent vertices, edges and faces. For instance, we can redefine the following: 26  • Pit: any vertex that has lower elevation than all its adjacent vertices. • Peak: any vertex that has higher elevation than all its adjacent vertices. In T I N , certain local features can only occur at certain types of triangulation objects (vertices, edges, faces). This can be decided by enumerating all T I N objects, and examining what features can lie on what objects. For instance, as we have seen, pits can only be on vertices. Another important observation is that according to the definition in the previous chapter, a "pass" in T I N can only lie on some vertex of the T I N . The identification of passes is a simple local test at T I N vertices. Figure 4.1 gives some idea about what types of objects can occur around a vertex.  upward channel  upward ridge  upward,slope(on a face)  downward watercourse  Figure 4.1: Terrain Feature Around a Vertex  4.2  Definitions Particular to TINs  We define the following terms in preparation for discussion on global feature extraction. In T I N ,  27  • Drainage edge of a face A: any incident edge of A that intersects the steepest descent line of some point on A. Some water on A drains out into that edge. See Figure 4.2. Non-drainage Edge  Drainage Edge  Figure 4.2: Drainage and Non-drainage Edges  • Non-drainage edge of a face A: any incident edge of A that intersects the steepest ascent line of some point on A. A n edge can then be classified into one of the following three types (see Figure 4.4): 1. Local channel: A n edge E that is a drainage edge of both its incident faces. Some water on both those faces drains into E, and flows down along E to its lower vertex.. Local channels collect water by themselves. See Figure 4.3 . 1  local channel  other sides of the triangles  Figure 4.3: Local Channel h is higher than /, ma and na are slope lines down to hi. All water in the shaded area flows through a  28  2. Local ridge: A n edge that is a non-drainage edge of both its incident faces. Water never drains into it. It is the dual of a local channel. 3. Plain edge: A n edge that is a drainage edge of one incident face, and a non-drainage edge of the other.  Local Channel  Local Ridge  Plain Edge  direction of steepest descent  Figure 4.4: Three Types of Edges  4.3  Global Drainage Features  Now we are ready to present algorithms for global feature extractions.  4.3.1  Drainage Network in T I N s  According to our definition of drainage network, in a T I N , every local channel is on the drainage network, because they collect water by themselves. For a point p on a triangular face to be on the network, since p does not collect water by itself, p has to lie on the trickle path of some local channel. This fact makes it easy to identify watercourses: they can start only at local channels, whose identification can be achieved by a local test. Specifically, the origin  29  of a watercourse is a local channel C whose upper vertex v has no incoming local channels or v doesn't drain into C (ie. C is not the steepest descent for v). By tracing down the steepest descents from watercourse origins, we obtain the complete drainage network.  local channels merge at a vertex  watercourse on a face merges with a local channel  Figure 4.5: Two Types of Watercourse Merge Points A n interior node (upon which two or more watercourses merge) of the drainage network can lie on either a vertex or a local channel in T I N (see Figure 4.5). If it's on a local channel, the dividing line that separates the two merging watercourses will coincide with one of them, but we can consider the dividing line as-being shifted up along the face by an infinitesimal offset (see Figure 4.6).  local channel  watercourse along a face  dividing line  Figure 4.6: Dividing Line Coincides with Watercourse  30  4.3.2  General Separating Lines in T I N s  Since crests start at passes, which can be identified with local operations, we can trace up along the steepest ascent from passes to obtain all the crests. Watercourse Trace-ups can be traced starting at watercourse origins, which can be identified with local operations. Watercourse dividing lines start at watercourse merge points, which are a by-product of identifying drainage networks.  '  We can identify the starting points of these general separating lines easily. We then perform tracing-ups from these points to obtain all general separating lines.  4.4 4.4.1  Dubious Dualities Dualities in a Terrain  There are many dual pairs of features in a terrain. For instance, pits and peaks, basins and hills, channels and ridges. A strict duality exists if by flipping the terrain up-side-down (by, say, negating all elevations), the definition we use for one feature can now identify the dual of that feature. In reality, if we insist on maintaining strict duality, we end up with incomplete features or infinite regress[44, 48]. For instance, in T I N , if we keep the definition of drainage network, and modify the definition of general separating lines to meet the strict duality requirement, then drainage network and general separating lines will not interlock to form a complete partition of the terrain. Instead, it can be shown, some needed separating lines will be missing, and some will be superfluous for the purpose of separating catchments. If, on the other hand, we keep our definition of general separating lines, and modify  31  the definition of drainage network by insisting that there must be a watercourse between two merging general separating lines (to form a dual of dividing lines, which separate two merging watercourses), we can get into troubles within even a single triangle. Say we have a pyramid triangle that has two local ridges and a local channel at the base, as in Figure 4.7. To separate the ridges, we have to introduce a watercourse, which, unfortunately inevitably, merges with the local channel. Now we need a general separating line to divide the two merging watercourses, ad infinitem.  Figure 4.7: Breaking Down of Dualities If we use differentiable surfaces instead of piecewise planar ones as an attempt to preserve duality, we end up facing a bigger problem: we even may not find a watercourse at the first place. If the terrain is shaped like a big bowl, then only the bottom point of the bowl is on the drainage network. This is undesirable. The fact that we are using a piecewise planar surface actually gives our definitions meanings, and makes them useful.  32  Chapter 5  D a t a Structure and Complexity of Drainage Processing In this chapter, we present the data structures used in and the computational complexities of the algorithms in terrain drainage information processing. The complexity is composed of several parts: complexity of building the T I N data structure; complexity of identifying local features; complexity of delineating global features, etc. If we are to handle queries, there are also the complexity of preprocessing and the complexity of queries themselves. We will address queries and their complexities in next chapter. Let n be the number of vertices in a TIN'. By Euler's relation, the number of edges is at most 3n, thus of 0(n).- The number of triangles is at most 2n, thus 0(n). each vertex, there are on average less than 6 adjacent vertices.  33  For  5.1  T I N Data Structure  Traditionally, a T I N is stored as a list of triangles, or as a graph, as indicated by Poiker in Chapter 39 of N C G I A Core Curriculum. Neither of the two has an explicit representation for the spatial ordering of adjacent edges of a vertex. The winged edge and quad edge[12] data structures, however, do explicitly represent the spatial ordering information. The operations on them are simple. There are only two primitive topology operations on quad edge data structure: "makeedge", which initializes a single edge, and "splice", which modifies the topology. A l l the complex operations can be implemented by using the two primitives. Quad edge data structure also has the additional advantage of representing the dual graph along with the original graph. For its simplicity and expressiveness, we choose to use the quad edge data structure to represent TINs. Given the input in the form of a list of 3-D surface points, we use a Delaunay Triangulation algorithm presented by Guibas et al. [12] to triangulate the x-y locations of the points. The algorithm uses quad edge data structure. The complexity of building a T I N is the same as that of the Delaunay Triangulation algorithm, which is  5.2  0(nlogn).  Local Terrain Features  To extract local terrain features (pits, peaks, passes), each vertex along with its adjacent vertices is examined. Let degree(v{) be the degree (in graph theoretic sense) of the ith vertex u,-, then the numbers of adjacent vertices, edges, and faces the extraction operation has to examine are all degree(vi). Furthermore, the complexity of the operation on a vertex is linear in the number of vertices, edges, and faces it has  34  to examine. Therefore the complexity of the operation on  is 0(degree(vi)).  The  total complexity of the operation on all the vertices is thus linear in the total degree of the graph: 0 ( ^ ™  = 1  degree(vi)).  By Euler's relation, we know the total degree of  the graph is at most 6n, in other words, YHi=i degree(v{) = 0(n).  Therefore the  complexity of the local feature extraction operation is 0(n). The lower bound of the computational complexity on local feature extraction is also O(n), therefore it is tight. Since each type of features is a distinct subset of the original vertices, the storage complexity for them is at most linear,  5.3 5.3.1  0(n).  Global Terrain Features Processing C o m p l e x i t y Is the Same as O u t p u t C o m p l e x i t y  A l l global features (drainage network, dividing lines, crests and watercourse traceups) starting points can be locally identified at T I N vertices in 0(n) time. The only exception are dividing lines which start at merge points of the drainage network. Since the drainage network is a forest, the number of merge points are linear in the number of watercourse origins, which is linear in the number of vertices n. Therefore, the number of starting points of dividing lines is 0(n).  As we already know, finding  starting points of dividing lines is just a by-product of generating drainage network. Watercourses and general separating lines each consist of connected line segments. Each line segment lies within a triangle. Its two end points lie on two edges of that triangle. To "extract" or "traverse" a watercourse or. a general separating line, we have to .go through all .of its. consisting line segments. A t each line segment, it takes constant 0(1) amount of processing. Also, it takes constant amount of space  35  to store a line segment. Thus, processing and storing a watercourse or a general separating line are linear in the number of its consisting segments. Let Si be a watercourse or a general separating line, |5,-| be the number of line segments in Si, Lemma 2 Extracting a watercourse or a general separating line Si takes © ( | 5 j | ) . Storing it also takes 0 ( | 5 | ) . t  We thus establish the following: after identifying all the local features, the time and space complexity of delineating drainage network and general separating lines is linear in the number of line segments comprising them.  In other words, the  processing complexity is the same as output complexity. We now set out to find the worst case and average case complexities of drainage networks and general separating lines, where "complexity" means the number of line segments.  5.3.2  U p p e r B o u n d on Drainage Network  Complexity  We start with drainage network complexity. The same analysis applies to general separating lines as well. We need to determine how many triangles (or, equivalently, triangulation edges) the drainage network can cross, noticing that watercourses are trickle paths of watercourse origins. We do this by showing that, first, there are at most 0(n)  water-  course origins; secondly, a drop of water trickling down along the surface can cross a particular triangulation edge at most 0(n) times; thirdly, there are at most 0(n) triangulation edges, therefore the total complexity of all watercourses is 0(n ). 3  first and the third facts are already established. We now show the second one.  36  The  Lemma 3 In a terrain represented by a TIN with n vertices, a trickle path meets any edge of the TIN at most 0{n) times. Proof: Break the trickle path into n subpaths by cutting at the elevation of each vertex of the T I N ; let p denote one of these subpaths. We will show that p hits each edge of the T I N at most twice. Subpath p visits a sequence of triangles and edges. Within each triangle, p's direction is determined by the steepest descent. Thus if p visits a triangle a second time, then it repeats the sequence that it visited before until it ends at the elevation of a vertex, as illustrated in Figure 5.1.  Figure 5.1: Repeating a Triangle Repeats the Subsequent Sequence Consider now the edge in the repeating sequence that has the shallowest slope, where the slope of a segment in 3-D is the ratio of the positive difference in its endpoints' z coordinates to the length of its projection in the x-y plane. If p visits this edge first-at- point p and lower at q, then we can derive a contradiction: The segment pq is shorter then the portion of p from p to q, denoted as p A q. p -4 q, however is steeper than pq, because p follows the steepest descent, thus the direction p follows on each triangle must be steeper than the preceding edge (so that p follows the face 37  of the triangle instead of the edge), which, in turn, is steeper than pq (pq being on the edge of the shallowest slope). Thus each segment on p A q is steeper than pq. A longer path (in projected length) from p having a steeper slope than pq, however, cannot end at q. Within a subpath, since the shallowest edge cannot be visited twice, and the sequence of visited triangles repeat itself, no other edge can be visited more the twice. Because there are n vertices, a single edge in T I N can be cut into at most n — 1 subpaths, thus a single edge in T I N can be visited at most 2n — 2, which is 0(n), times by a single trickle path. Q E D .  5.3.3  A Worst-case E x a m p l e  The above stated complexity 0(n ) 3  of the drainage network can actually be achieved,  as shown in Figure 5.2.  Figure 5.2: Worst Case Example for Drainage Network Complexity: the Pyramid Along one edge of the pyramid, form a scree ("rocky slope) of n triangles. On the top, carve n local channels that will catch water and start n rivers near the top of the scree.  38  Form n gently-sloping helical ducts that meet the rivers after the rivers meet the scree. Each duct begins with a nearly-flat top whose steepest descent angles slightly away from the pyramid face ( preventing itself from forming a local channel with the pyramid face), so that rivers that join the duct separately remain separate. A t a turn, the two top triangles of the duct share a vertex at the outside corner and a small, near-vertical triangle interpolates between them, forming separate waterfalls for separate rivers. The resulting terrain has 0(n)  vertices and n initial rivers each cross the n scree  triangles n times. Thus, Theorem 1 In a TIN with n vertices, the worst-case complexity of a watercourse is 6(n ), and of all watercourses is 9(n ) 2  3  This upper bound is tight. The worse-case complexity of storing global features is also  0(N ). 3  Theorem 1 applies to arbitrary triangulations. We suspect Delaunay Triangulation with its empty circle property and maximization of smallest angles will give us a less expensive upper-bound.  5.3.4  Average-case C o m p l e x i t y  Analysis of average-case complexity based on randomized surface elevations gives unrealistically high numbers of pits, peaks, and watercourses. For example, great number of pits and peaks are introduced by randomization, while in nature, the *numbersof pits and peaks are small' relative to the number of sample points. Table 5.1-shows .experimental results we have obtained from sample data on St. Mary's Lake, B . C . We can see that the network has slightly fewer line segments than there are vertices. 39  Table 5.1: Observed Drainage Network Complexity No. Vertices Network Size  5.3.5  450 353  1800 1427  4956 4134  11125 10579  Terrain Planar G r a p h  To obtain the planar graph we described in the previous chapter, we need to perform some extra processing.  The information we need can all be collected while we  generate the global features, without substantially increasing the running time of the algorithm. From there, we can generate the graph itself. The complexity is a function of the size of the graph: how many pits there are, how many crests there are, etc. Storage complexity is linear in the size of the graph. To indicate which crest belongs to an edge in the graph, we simply need a pointer to the crest.  40  Chapter 6  Drainage Queries Identifying local and global drainage features is useful in its own right, furthermore, it helps answer many interesting terrain drainage queries. In this chapter, we will examine what types of queries can be answered with the help of our data structure; what preprocessing steps we can take to reduce the complexity of the queries; and what the complexities of preprocessing and the queries themselves are. Some of the queries can be answered with little extra processing beyond generating the global features. We will discuss them first.  6.1 6.1.1  Simple Queries F i n d i n g Basin Boundaries  Theorem 2 A basin's boundary is formed by a closed chain of crests. With a planar graph of the terrain, finding basin boundaries amounts to finding all the incident edges of a face in the graph. We can do this by following the "righthand" rule: walk along the edges while keeping the right hand touching the wall.  41  Using quad edge data structure to store the graph makes following the "right-hand" rule even more trivial.  .  With known basin boundaries, we can compute the area and volume of basins. The area is defined as the projected area of the basin onto the x-y plane. The volume is defined as the the space bounded from below by the terrain surface and from above by a plane at the level of the basin's outlet (lowest point on the basin's boundary).  6.2  Finer Partition  In order to reduce the cost of watershed related queries, we need to generate more trace-up lines in addition to our global features. We add' a new type of trace-ups called "interior node trace-ups" starting at drainage forest's interior nodes. For each interior node, we add two trace-ups; each divides an incoming watercourse and an out going watercourse of the node, see Figure 6.1. W i t h these additions, we have obtained watershed boundaries of all drainage forest nodes. These additions along with drainage network and general separating lines form a finer partition of the terrain. We call each area bounded by these lines a partition unit.  6.3  Querying the Watershed of a Point  To find out the boundary of the watershed of a query point p, trace-ups have to be performed. The query operation takes several steps: first locate the triangle that contains p. This takes 0(y/n)  expected time if no other assisting data structure is  used along with the quad edge (we need to walk across the triangles). If other data structures are used, this may be improved to O(logn). The second step is to perform trace-ups, see Figure 6.2. Finally, the watershed  42  new trace-up line  watercourse  general separating line  Figure 6.1: New Trace-up Lines boundary of p can be delineated by using these trace-ups and the finer terrain partition boundaries.  p is not on watercourse  p is on a watercourse and a face  p is on a watercourse'and a local channel  Figure 6.2: Query Watershed of a Point  43  6.4  Querying the Size of the Watershed of a Point  If the query point is not on the drainage network, its watershed is zero in size. Otherwise, we need more preprocessing to achieve constant time query (excluding point location) using an interpolation scheme. Storing watershed sizes of drainage merge points does not provide enough information for performing accurate interpolation, since the interpolation scheme is not well defined between merge points. What we need is to store watershed sizes of as few selected points along the drainage network as possible, yet achieving constant time interpolation. Now we describe the criteria for selecting these points. We start by making a few more definitions.  6.4.1  Definitions  • A watercourse segment is a maximal portion of the drainage network that doesn't contain watercourse merge points, watercourse origins, or pits. Drainage network consisting of watercourse segments are called refined drainage network. • The drain of a point p, denoted drain(p), is the point where p's trickle path first meets the drainage network. If p is on the network, then drain(p) = p. • A catchment area of a watercourse segment s consists of all points p for which drain(p) is on s. It is the same as a partition unit of the finer partition in preceding discussions. •• - • A strip is formed as follows:-for every vertex v of-the T I N , we introduce a node at drain(y).  These nodes cut the watercourse segments into yet smaller  44  segments that we call the base segments. It is easily verified that all these nodes lie on local channels, see Figure 6.3. • For a base segment s, define the right strip to be the set of points p, such that drain(p) is on s, and p's trickle path comes from the right, where the front is in line with the direction of the flow. Define the left strip similarly.  Figure 6.3: Strips  6.4.2  L e m m a and Corollaries  Lemma 4 A strip is a region bounded by a base segment from a local channel, a segment of a local ridge, and two steepest descent paths. No vertices of the TIN lie inside a strip. We have two corollaries following directly from the lemma.  /  Corollary 1 Basins, source areas, and catchment areas are bounded by local ridges and steepest descent paths. Corollary 2. Watershed area,.projected area, or steady-state flow rate under uniform rainfall assumptions can be summarized by piecewise quadratic functions on the base segments of the refined drainage network. 45  Proof: For a query point q on the drainage network, the left and right strips of the base segment s containing q are partially inside the watershed of q. Every other strip is either entirely inside or outside of the watershed of q, and contributes a constant amount of area, projected area, or steady-state flow. The contribution from the left and right strips of base segment s, as shown before, can be summarized as a quadratic function at + bt + c, where t indicates 2  the position of q on s. We now give an argument for the quadratic behavior of the function: Imagine we move along the drainage network, the trace-ups starting at the current position would thus sweep the terrain. Within a triangle, the sweeping trace-ups are all parallel, since the steepest ascent direction is unique on a plane. The area bounded by the two parallel lines and the triangle's boundary is quadratic in the distance between the two parallel lines, provided none of the triangle's vertices lie between the two parallel lines (see Figure 6.4). The distance between the two parallel lines is linear in the distance between the intersections of an edge with the two parallel lines. The size of the swept area is quadratic in the distance swept. The interpolation parameters for the area change only when the trace-ups sweep across a vertex.  Figure 6.4: Sweeping Trace-up Line  46  This interpolation scheme can be extended to the whole trace-ups across multiple triangles: we can find the watershed sizes of drainage network points whose traceups sweep a vertex. We then pick a point on the drainage network that lies between each pair of consecutive points above described, and find its watershed area. These in-between points are necessary since we need three points to establish a quadratic interpolation. See Figure 6.5.  Figure 6.5: Interpolation Points The following algorithm for preprocessing proves to be efficient: Trace down from each vertex until hitting the drainage network. This takes  Ylt=i  where 5,- is a trace-up, and \Si\ is the size (number of line segments) of Si. For each hit point and a point between each pair of consecutive hit points, we perform a trace-up to calculate the point's watershed size. The work needed for tracing up and calculating the area of the strip between trace-ups is 0(|S,-|). Therefore, the total complexity is J2i=i  Storing the interpolation data takes 0(n)  storage.  Lemma 5 Preprocessing for constant time (excluding, point location operation) watershed size queries takes YLl=i  time, and O(n) space in the above algorithm.  47  Chapter 7  Conclusions We have delved into the particulars of our models and methods in the preceding chapters. We conclude with an overview of the whole problem.  7.1  Level of Abstraction  Like anything else, to solving drainage problems, we need to keep our goal and resources in mind. We claim that richer and important information is obtained through extra processing in our method.  7.1.1  Advantages  We recount the benefits of our approach as follows: • The model we have provides a common ground for discussions on terrain drainage characteristics.  The definitions are general rather than particular  to a specific representation such as T I N o r . D E M ! They are mathematical, thus precise. It is a step closer to reality compared with previous researches.  48  • We can generate all the standard terrain features such as pits, peaks, watercourses, ridge lines. We give more consistent answers to these problems, because we have a well defined model and generally applicable definitions, rather than define these features as a product of a particular algorithm. • We can easily extract the planar graph structure of terrain, which characterizes terrain drainage at the highest level.  This graph structure is an inherent  feature of the.terrain. It is know to researchers for decades. • We can accurately and efficiently answer watershed and catchment queries that are of interest to hydrology, hydraulic engineering and natural resource management. This is an important function that previous research has not fully addressed.  7.1.2  Disadvantages  The main disadvantages are: • Complex. Our methods are more complex than cell based methods. However, we note that the local-operator/line-connection approach of finding V shaped grid points and connecting them into drainage network is complex to specify completely, especially since the line-connection part is not well defined and may not work properly. Approaches like Mark's are simple in concept, but costly in computation (O(nlogn)). Our approach uses a more complete model than Mark's, but wins out in computation. , • Computational Cost. The computational cost is quite high in the theoretical worst-case. Because unlike D E M approaches, which never introduce new nodes or points besides the existing grid points, in our method we do need to add 49  nodes, intersection points between edges and trickle paths, etc.  Thus our  approach could potentially be very expensive, as indicated in the worst-case analysis, reaching 0 ( n ) . In tests of actual terrain data, the complexity of 3  global features in our method is linear 0(n) (excluding the T I N construction complexity).  7.2  Reflection of Reality  How well does our model reflect the physical process? It bears more resemblance to the processes of the physical world than previous methods. It makes water follow along the steepest descent of the surface.  Even  though it cannot account for dynamics, it does reflect the law of gravity to certain extent. Do the features we extract closely correspond to the features in the actual terrain? As mentioned in the introduction, the assumption that the surface, computed out of sample points possibly lying on the surface of water bodies, is the underlying dry surface water flows on is at variance with the nature of the sample data. But that's as good as we can do, since those are the only data we have. Furthermore, there are factors such as erosion, weather, sediment concentration, vegetation that affect the drainage process. We are not looking into any of those, although it is possible to incorporate them.  7.3  Future Research  Future research in terrain drainage will try to address some of the data structure and algorithm issues, make the model capable of handling special cases, and incorporate  50  different factors into our model.  7.3.1  M a k i n g the M o d e l M o r e Realistic  Currently, a pit is a true sink: water that goes in never comes out. Instead, we may want to make water fill to certain level, and have an absorption/evaporation model to keep the water level, or make water drain from the outlet.  7.3.2  D a t a Structures  Efficient data structures and algorithms need to be devised to process and store watercourses, general separating lines, planar terrain graph, and drainage forest.  7.3.3  H a n d l i n g Special Cases  In our model, we haven't explicitly addressed the problem of special cases. Some of the special cases are: • A t data level, there can be flat areas on the surface, where the direction of steepest descent is not well defined; there can be horizontal edges; there can be artifacts of sampling errors such as shallow pits, etc. We need to have a strategy to handle these cases. For example, consider flat areas as pits or lakes, and stop trickling down at such; use a filter to fill in small pits, make small pits overflow at the outlet, etc. • A t structure level, watercourses can go off the terrain boundary or come into -••I-  -the terrain boundary-from outside, and there can be partial basins. In such cases, we can never-.correctly trace all the drainage features because the necessary connections are outside of our data region. We need to devise methods to handle these cases. 51  7.3.4  A d d i n g Other Factors  We can add factors such as variable precipitation, water absorption of soil, accumulation of water at basins, sub-surface water flow.  7.3.5  M e r g i n g D a t a Sources  So far, we only used surface geometry data of the terrain. One possible additional data source is a river/stream map of the terrain. If we can get both the digital elevation data and a river/stream map, one way of merging the two pieces of information would be to overlay them, and adjust the elevations and/or x-y locations of the sample points with certain constraint, so that wherever there is a river in the river map, there will be one in the elevation model conforming to our definition of a river. Although there may be artifacts in adjusting the points, it is a first step in merging drainage information.  52  Bibliography [1] L . E . Band. Extraction of channel networks and topographic parameters from digital elevation data. In K . Beven and M . J . Kirkby, editors, Channel Network Hydrology, pages 13-42. John Wiley k Sons, New York, 1993. [2] J.R. Carter. Digital representation of topographic surfaces. Engineering and Remote Sensing, 54(11):1577-1580, 1988.  Photogrammetric  [3] A . Cayley. On contour and slope lines. Lond. Edin. Dublin Phil. Mag. and J. of Sci., 18(120):264-268, 1859. [4] Z. Chen and J . A . Guevara. Systematic selection of very important points(VIP) from digital terrain models for construction triangular irregular networ. In Proceedings of AutoCarto 8, pages 50-56, Falls Church, V A , 1987. [5] S. H . Collins. Terrain parameters directly from a digital terrain model. Canadian Surveyor, 29(5):507-518, 1975. [6] John C . Davis. Contouring algorithms. In Proc. Auto-Carto 2, pages 352-359. American Congress on Surveying and Mapping, 1975. [7] David H . Douglas. Experiments to locate ridges and channels to create a new type of digital elevation model. The Canadian Surveyer, 41(3):373-406, A u tumn 1986. [8] David H . Douglas and Thomas K . Peucker. Algorithms for the reduction of the number of points required to represent a line or its caricature. The Canadian Cartographer, 10(2):112-122, 1973. [9] L . De Floriani, B . Falcidieno, and C . Pienovi. Delauney-based representation of surfaces defined over arbitrarily shaped domains. Computer Vision, Graphics, and Image Processing, 32:127-140, 1985.  53  [10] A . U . Frank, B . Palmer, and V . B . Robinson. Formal methods for the accurate definition of some fundamental terms in physical geography. In Proc. 2nd Intl. Symp. Spatial Data Handling, pages 585-599, 1986. [11] John M . Gauch and Stephen M . Pizer. Multiresolution analysis of ridge and valleys in grey-scale images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 15(6):635-646, June 1993. [12] Leonidas Guibas and Jorge Stolfi. Primitives for the manipulation of general subdivisions and the computation of Voronoi Diagrams. ACM Transactions on Graphics, 4(2):74-123, April 1985. [13] S. K . Jenson and C . M . Trautwein. Methods and applications in surface depression analysis. In Proc. Auto-Carto 8, pages 137-144, 1987. [14] S.K. Jenson. Automated derivation of hydrologic basin characteristics from digital elevation model data. In Proc. Auto-Carto 7, pages 301-310, Falls Church, V A 22046, 1985. American Society of Photogrammetry and the American Congress on Surveying and Mapping. [15] S.K. Jenson and J . O . Domingue. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing, 54(11):1593-1600, Nov. 1988. [16] Norman L . Jones, Stephen G . Wright, and David R. Maidment. Watershed delineation with triangle-based terrain models. Journal of Hydraulic Engineering, 116(10):1232-1251, October 1990. [17] Jan J . Koenderink and Andrea J . van Doom. Local features of smooth shapes: Ridges and courses. In Geometric Methods in Computer Vision II, volume 2031, pages 2-13. SPIE, 1993. [18] Jan J . Koenderink and Andrea J . van Doom. Two-plus-one-dimensional differential geometry. Patton Recognition Letters, 15:439-443, May 1994. [19] M . P . Kumler. A n intensive comparison of triangulated irregular network (TINs) and digital elevation models (DEMs). Cartographica, 31(2):1, summer 1994. [20] I.S. Kweon and T . Kanade. Extracting topological terrain features from elevation maps. CVGIP: Image Understanding, 59(2):171-182, March 1994. [21] D . M . Mark. Topological properties of geographic surfaces: Applications in computer cartography. In G . Dutton, editor, First Int'l Adv. Study Symp. on Topological Data Structures for Geo. Info. Sys., volume 5. Harvard, 1978. 54  [22] D . M . Mark. Automated detection of drainage networks from digital elevation models. In Barry S. Wellar, editor, Automated cartography: international perspectives on achievements and challenges, pages 288-295. Proc. Auto-Carto 6, 1983. [23] D . Marks, J . Dozier, and J . Frew. Automated basin delineation from digital elevation data. Geo-Processing, 2:299-311, 1984. [24] J . C . Maxwell. On hills and dales. Lond. Edin. Dublin Phil. Mag. and J. of Sci., 40(269):421-425, 1870. [25] J . E . McCormack, M . N . Gahegan, S. A . Roberts, J . Hogg, and B . S. Hoyle. Feature-based derivation of drainage networks. IJGIS, 7(3):263-279, 1993. • [26] D . G . Morris and R. W . Flavin. A digital terrain model for hydrology. In Proc. 4th Intl. Symp. Spatial Data Handling, pages 250-262, 1990. [27] D . G . Morris and R . G . Heerdegen. Automatically derived catchment boundaries and channel networks and their hydrological applications. Geomorphology, 1:131-141, 1988. [28] David M . Mount. Storing the subdivision of a polyhedral surface. In Proceedings of the Second Annual Symposium on Computational Geometry, pages 150-158, Yorktown Heights, New York, June 1986. [29] John F . O'Callaghan and David M . Mark. The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing, 28:323-344,1984. [30] O . L . Palacios-Velez and B . Cuevas-Renaud. Automated river-course, ridge and basin delineation from digital elevation data. Journal of Hydrology, 86:299-314, 1986. [31] T . K . Peucker and D . H . Douglas. Detection of surface-specific points by local parallel processing of discrete terrain elevation data. CVGIP, 4:375-387, 1975. [32] T . K . Peucker, R. J . Fowler, J . J . Little, and D . M . Mark. The triangulated irregular network. In Amer. Soc. Photogrammetry Proc. Digital Terrain Models Symposium, pages 516-532, 1978. [33] T . K . Peucker and N . Chrisman. Cartographic data structures. American Cartographer, 2(l):55-69, 1975. [34] J . L . Pfaltz. Surface networks. Geographical Analysis, 8:77-93, 1976. 55  [35] Jianzhong Qian, Roger W . Ehrich, and James B . Campbell. D N E S Y S — A n expert system for automatic extraction of drainage networks from digital elevation data. IEEE Trans. Geosci. Remote Sens., 28(1):29-45, January 1990. [36] William W . Seemuller. The extraction of ordered vector drainage networks from elevation data. Computer Vision, Graphics, and Image Processing, 47:45-58, 1989. [37] Y . Shinagawa and T . L . Kunii. Constructing a Reeb Graph automatically from cross sections. IEEE Computer Graphics and Applications, 11(6):44^51, 1991. [38] A . T . Silfer, G . J . Kinn, and J . M . Hassett. A geographic information system utilizing the triangulated irregular network as a basis for hydrologic modeling. • In Proc. Auto-Carto 8, pages 129-136, 1987. [39] J . Snoeyink, M . van Kreveld, and S. Yu. The complexity of rivers in triangulated terrains. In Proceedings of Canadian Conference on Computational Geometry, Ottawa, Canada, 1996. [40] U.S. Geological Survey. U.S. Geological Survey Data User's Guide 5, chapter Digital Elevation Models, page 38. 1987.[41] S. Takahashi, T . Ikeda, Y . Shinagawa, T . L . Kunii, and M . Ueda. Algorithms for extracting correct critical points and constructing topological graphs from discrete geographical elevation data, 1995. [42] D . M . Theobald and M . F . Goodchild. Artifacts of tin-based surface flow modeling. In Proc. GIS/LIS'90,  pages 955-964, 1990.  [43] W . Warntz. The topology of a socio-economic terrain and spatial flow. Papers of the Regional Science Association, (17):47-61, 1966. [44] C . Werner. Several duality theorems for interlocking ridge and channel networks. Water Res. Res., 27(12):3237-3247, December 1991. [45] D . Wilcox and H . Moellering. Pass location to facilitate the direct extraction of warntz networks from grid digital elevation data. In Proc. Auto-Carto  12,  pages 22-31, 1995. [46] G . W . Wolf. Metric surface networks. In Proc. 4th Intl. Symp. Spatial Data Handling, pages 844-856, 1990.  56  [47] G . W . Wolf. Hydrologic applications of weighted surface networks. In Proc. 5th Intl. Symp. Spatial Data Handling, pages 567-579. I G U Commision on GIS, 1992. [48] S. Y u , M . van Kreveld, and J . Snoeyink. Drainage queries in TINs: From local to global and back again. In Proc. 7th Intl. Symp. Spatial Data Handling, Amsterdam, Netherlands, 1996.  57  Appendix A  Terrain Function The surface of a terrain can be represented as the graph in R  3  of a continuous real  function z = f(x,y) defined on a closed, compact, connected subset of the xy plane. For r £ R and r > 0, the forward directional derivative g(d) of f(x,y)  at point  (xo,yo) in the direction 6 £ [0,27r] (measured from the positive x direction) is: /0\ _ j .  5  (0)  =  5  fjxp + r cos (0), y + r sin (0)) - f(x , 0  0  y) 0  (2tt).  If the part of the surface (XQ, yo, /(xo, yo)) lies on is a plane, we derive:  g{6) =  kcos{6-6 ) 0  where k is the greatest slope value of the plane, and f? is the direction of the greatest 0  ascent. In a T I N , for a point (xo,yo), g(0) is defined on a partition of [0,27r], where the number of intervals (wrapped around intervals are considered as one interval) in the partition is the same as the number of triangle(s) incident upon the point. On each  58  interval, g(6) is a sinusoid with a possibly different coefficient k and phase angle 6Q. See Figure A . l .  There are four sinusoids indicating there are four planes incident upon the point  steepest descent Figure A . l : Sinusoidal Directional Derivative The analysis of drainage characteristics can be strictly based on the function g(0). A graph of g{6) will be very helpful in the analysis. For instance, the steepest descent direction is the direction 6 for which g{6) has the least negative value. A pit's g(6) is all positive. In T I N , g(6) of a point on a local ridge has two negative minima, and two maxima that lie on the division points of the intervals. A maximum ascent lies in the direction 6, where g{6) is a positive maximum. A maximum descent lies in the direction 6, where g(0) is a negative minimum..  59  Appendix B  Horizontal Slicing Plane Method Many of the relations between ascents and descents can be analyzed by a horizontal slicing plane method. In analyzing a point, imagine cutting the terrain surface with a horizontal plane which lies slightly above or below the point. The intersection of the terrain surface and the plane is polylines. We then project the edges that are below (or above) this plane onto the plane, we get a set of new polygons. If we project the steepest descents of faces onto this plane, we know that these projections are perpendicular to their corresponding faces' intersections with the plane. As an example, if we cut above a pit, we get a polygon, where the pit sees each vertex of the polygon. See Figure B . l . Because all the points on the polygon are of the same elevation, the further away the point is from the pit, the less g(6) is in the direction pointing from the pit to that point, ie. the smaller the slope is in that direction. Therefore the analysis of maximum and minimum of g{6) can be reduced to analyzing the distance between the intersection point and the point of interest.  60  By using these facts, we can easily show that a pit has at least and possibly only one local channel; between two incoming watercourses that are not separated by any outgoing watercourse, there can be only one maximum ascent.  Figure B . l : Horizontal Plane of a Pit  61  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items