|
Return to GeoComputation 99 Index
Honglei Zhu and Kristin Schneider
Clark Labs, Clark University, 950 Main Street, Worcester, MA 01610 U.S.A.
E-Mail: hzhu@clarku.edu
A geographic information system (GIS) can be used to extract geomorphic information from a regular grid digital elevation model (DEM) or a triangulated irregular network (TIN) model for hydrological modeling. Besides natural surface depressions, other flat areas in a DEM or a TIN prevent us from obtaining the information that hydrologic modeling requires. This paper proposes algorithms for removing flat features in a TIN model. Flat triangles, flat channels and flat ridges are processed.
Remove flat triangles in a TIN model. Flow directions cannot be calculated if a triangle is flat. A TIN model must be flat-triangle-free before it is used to extract features such as channel networks or watershed boundaries. Two cases are considered in this study, and algorithms are developed.
TIN models generated from a discrete point data set: Flat triangles are processed according to the status of flat triangles’ neighboring triangles.
TIN model generated from a contour data set: Edges known as bridge and tunnel edges form flat triangles in a TIN when a TIN is generated from a contour data set. The contour data set is used to assist the process of flat triangle removal. By inserting a point at the middle of a bridge or tunnel edge, flat triangles can be removed from a TIN model. The elevations of added points are interpolated, and the TIN is adjusted.
Flat channel and ridge processing. Flat channels and ridges may prevent us from locating depression points, peaks, which can be used for depression removal process, and watershed information extraction process. Algorithms for processing flat channels and ridges are proposed.
The methodology is validated via a case study.
Regular grid digital elevation models (DEM) and triangulated irregular networks (TIN), among others, are widely used in a geographic information systems (GIS) to represent a digital terrain model. A GIS can be used to extract geomorphic information from a DEM or a TIN model for hydrologic modeling. Flat areas in such models prevent us from obtaining correct flow direction information that hydrologic modeling requires. Such features have to be removed before further processing. In the past, partly because of the relative ease access of a DEM and manipulation of grid data structures, DEMs have drawn more attention. Researchers have developed feature extraction algorithms based on DEM (Jenson and Trautwein, 1987; Jenson and Domingue, 1988, etc.). With the adoption of TIN models in GIS as topographic surfaces (Grayman, et al 1979; Mark, 1998), several papers have contributed to feature extraction based on a TIN model (Guercio and Soccodato 1996, Tachikawa, et al., 1996).
Few papers were found to process flat features based on a TIN model. This paper proposes algorithms for removing flat features in a TIN model. Flat triangles, flat channels and ridges are processed.
The Delaunay method commonly is used in a triangulation process to produce a TIN model. When constraints are adopted into such a process, a Delaunay triangulation is extended into a constrained Delaunay triangulation. In such triangulation processes, a flat triangle is generated if three vertices of a triangle have the same attributes, or elevation values.
Topographic information, such as slope and aspect, are primary features used in a hydrologic model. Flow directions cannot be calculated for flat triangle areas. There are two flat triangle removal processing. One is for a TIN created from a point data set, the other is a TIN created from a contour data set.
Neighboring triangles of a flat triangle are used to assist the removal process. A triangle has a neighboring triangle if it shares a common edge with another triangle. A triangle can have up to three neighboring triangles. For a flat triangle removal process, based on the number of neighboring triangles, and on the attributes of its neighboring triangles, the following algorithms are proposed.
The linear method is adopted to interpolate the elevation of a centroid point based on each neighboring triangle. The final elevation of the centroid point is the average of the three values calculated based on the three neighboring triangles.
In Figure 1, parameters used for a neighboring triangle’s interpolation are labeled. A middle point of a flat edge is used as a temporary point. It takes the flat triangle’s elevation, and is used to interpolate the centroid point’s elevation. The distance between point P4 and the middle point of a flat edge Pmid, and the distance between point P4 and the centroid point Pcenter of a flat triangle are used as weights in the interpolation.
Linear method is used to interpolate the elevation of the middle point.
The linear method is used to interpolate the elevation points at the middle points of the two split edges.
An experimental equation is used for adjusting the elevation value of the vertex. If the minimum elevation difference between points P4-P0 and P3-P0 is: dMinH = Min[( Hp4-Hp0), ( Hp3-Hp0)], then the vertex’s elevation is adjusted as: Hp0 = Hp0 ± dMinH/3.0.
For a flat triangle, if one of its three edges is a boundary edge of a TIN, this flat triangle only has two neighboring triangles. Two cases are considered in the flat triangle removal process.
If two edges were used as a boundary edge of a TIN, only the third edge has a neighboring triangle. Two cases are considered in the flat triangle removal process.
A similar linear interpolation method as illustrated in Figure 1 is used.
A similar linear interpolation method is used as illustrated in Figure 5.
A flat triangle is assigned a priority rating according to the following criteria, and then the prioritized flat triangles are processed.
A flat triangle without any flat neighboring triangles has the highest priority for processing. A flat triangle can have up to three non-flat neighboring triangles if there is no edge used as a boundary edge.
A flat triangle with two non-flat neighboring triangles has the second priority. The third edge of the flat triangle is shared with another flat triangle.
A flat triangle with only one non-flat neighboring triangle has the lowest priority. The other two edges are shared with two flat neighboring triangles.
The flat triangle removal process is an iterative process. Flat triangles with highest priorities always are processed before those with lower priorities. The priority status of a flat triangle is updated dynamically if there is any change. A flat triangle with three flat neighboring triangles is not processed until at lease one of its three neighboring triangles has been processed, and its priority has a turn to be processed. The iterations continue until no more flat triangles are detected in the TIN.
If a TIN is generated based on a contour data set, the information provided by the contour data can be used in the process of flat triangle removal. From a different perspective, knowing the presence of flat triangles ( Zhu, Eastman and Schneider, 1999) derives an edge known as a bridge or tunnel edge. Bridge and tunnel (B/T) edges are illustrated in Figure 10.
A B/T edge can be shared by two triangles, or used as a boundary edge. They are processed respectively. As illustrated in Figure 11, a point at the middle of a B/T edge is added, the B/T edge is split, and the TIN is adjusted. Figure 13 presents a comparison before and after the B/T edge removal process for a simple TIN data set generated from a contour data set.
The interpolation of an elevation value at the middle point of a B/T edge is conducted by using a parabolic method using an eight-direction elevation search. An eight-direction search process tries to find two intersections with contours from each direction (see Figure 12). Up to eight elevation values can be obtained. Three points are required to calculate an elevation value of which two intersection points are from one direction and the third from the opposite direction. The final elevation is the averaged value of eight interpolation values. Equation 1 is used to calculate one direction’s elevation value based on three intersection points.
Where hi, i = 0, 1, 2 are elevation values of three intersection points; Si, Sj, i, j=0,1,2 represent the distances between the intersection point and the added point. The furthest point in that direction is used as the original (see Figure 12); therefore, the second point is nearer to the added point, and the third point is on the other side of the added point. In Figure 12, detailed information for north direction’s interpolation is labeled. According to the above definition, we know S0≡0, while S1=P1P0 , and S2=P2P0.
After the end of processing, all of the flat triangles have been removed. A certain flow direction can be calculated for any point in a triangle. For a process of watershed related information extraction, a collection of basic surface features such as depression points, peaks, channel networks and ridge networks, among other features, must be extracted. Flat channel and flat ridge edges are additional features that may prevent us from obtaining such information. Algorithms for removing flat channels and flat ridges also are proposed.
Once there is a surface flow, based on a TIN model, a non-flat triangle must have at least one, at most two flow-in edges, this means the surface water flows into the triangle across the edge(s). And also it must have at least one, and at most two flow-out edges. When an edge is shared by two neighboring triangles and it has a flow-out attribute for both triangles, this edge is named as a channel edge. When an edge shared by two neighboring triangles and it has a flow-in attribute for both triangles, this edge is called a ridge edge. If a channel edge is flat, the flow direction of its accumulation can be calculated. If a ridge is flat, a ridge tracing procedure cannot reach a peak to form a complete watershed boundary. These phenomena are illustrated in Figure 14. And the algorithms for flat channel and flat ridge processing are illustrated in Figure 15.
If a flat edge is detected as a channel edge, a searching procedure for two lower connected vertices are started from both ending vertices of the flat channel edge. The search procedure is continued until two vertices are found which have a higher or lower elevation value than that of the flat channel edge.
If a flat edge is detected as a ridge edge, a searching procedure for two higher connected vertices are started from both ending vertices of the flat ridge edge. The search procedure is continued until two vertices are found which have a higher or lower elevation value than that of the flat ridge edge.
If there is only one flat edge detected, and no more flat edges connect with this flat edge, two cases are considered.
a. If the two ending vertices found in the searching are all higher or lower than the flat edge, a point will be inserted at the middle of the flat edge. A linear method is used for interpolation from both sides, and the final elevation of the added point is the average value (see (a) and (b) in Figure 15).
b. If one ending vertex is higher, the other is lower, the two vertices of the flat edge are adjusted using a linear method based on the two ending vertices (see (c) in Figure 15).
If two or more connected flat edges are detected in a channel or ridge network, vertices in the flat edge are adjusted using a linear method (see (d), (e) and (f) in Figure 15).
Figures 16 and 17 present two case studies based on two different data sets. Flat triangles, flat channels, and flat ridges were processed after applying the algorithms proposed by this study. The preliminary case study results have validated the algorithms and procedures proposed in the paper. After removing depressions, based on the primary information in a TIN model, extractions of channel networks and watershed information, and other information for hydrologic modeling, such as a surface flow accumulation and flow length can be extracted. Testing continues in order to evaluate the adequacy of the processed TIN model for deriving hydrological characteristics and surface process modeling.
This study has been done as part of a research and development project at the Clark Labs. Further research and development is being undertaken as well as the continued testing and validation of currently developed algorithms in different geographic regions.
The authors would like to give thanks to J. Ronald Eastman for his directives and valuable discussions.
Grayman Walter M., B.A. Woodridge, E.B. Long, and A.C. Vidra, 1979. Joint State/regional Environmental Planning Using the PEMSO/ADAPT Geographic Information System, AutoCarto 4, Proceedings of the International Symposium on Cartographic and Computing: Applications in Health and Environment, Vol. 1, November 4-8, Reston, VA. pp. 546-553.
Guercio, R. and F.M. Soccodato, 1996. GIS Procedures for automatic extraction of Geomorphological Attributes from TIN-DTM, HydroGIS 96: Application of Geographic Information Systems in Hydrology and Watershed Resources Management ( Proceedings of the Vienna Conference, April 1996), IAHS Publ. No. 235, pp. 175-182.
Jenson, S.K. and J.O. Domingue, 1988. Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis, Photogrammetric Engineering and Remote Sensing, 11, pp. 1,593-1,600.
Jones, N.L., S.G. Write, and D.R. Maidment, 1990. Watershed Delineation with Triangle-Based Terrain Models, Journal of Hydraulics Engineering, ASCE 116(16), pp. 1,232-1,251.
Mark, D.M., 1998. The History of Geographic Information Systems: Invention and Re-invention of Triangulated Irregular Networks (TINs), Proceedings of ASPRS-RTI, 1998 Annual Conference, March 30-April 3, pp. 284-289.
Tachikawa Y., T. Takasao, and M. Shiiba, 1996. TIN-based topographic modeling and runoff prediction using a basin geomorphic information system. HydroGIS 96: Application of Geographic Information Systems in Hydrology and Watershed Resources Management ( Proceedings of the Vienna Conference, April 1996), IAHS Publ. No. 235, pp. 225-232.
Zhu, H., J.R. Eastman, and K. Schneider, 1999. Constrained Delaunay Triangulation and TIN Optimization Using Contour Data, The Thirteenth International Conference on Applied Geologic Remote Sensing, Vancouver, British Columbia, Canada.
Figure1. Three neighboring triangles of a flat triangle are all higher or all lower than the flat triangle, a centroid point is interpolated and added into the TIN. The flat triangle is split into three triangles, the TIN data is adjusted. A ‘+’ indicates a higher vertex, while a ‘-‘ indicates a lower vertex. The flat triangle is filled with gray. Interpolation from one neighboring triangle is labeled. A linear method is adopted in the interpolation. The elevation of the centroid point is calculated based on points P4 and Pmid, and distances DS1 and DS2.
Figure 2: Two neighboring triangles are higher or lower than the flat triangle, and the third neighboring triangle is lower or higher, the edge shared with the third triangle is split by adding a middle point, related triangles are split, and the TIN is adjusted. A linear method is adopted in the interpolation. The elevation of the centroid point is calculated based on points P4 and P2, and distances DS1 and DS2. A ‘+’ indicates a higher vertex, while a ‘-‘ indicates a lower vertex. The flat triangle is filled with gray.
Figure 3. A flat triangle has three neighboring triangles, one is higher, one is lower, and the third one also is a flat triangle. The first two edges of the flat triangle are split. Three original triangles are split, and the TIN is adjusted. A linear method is adopted in the interpolation. The elevation of two adding points are calculated based on points P4, Pmid1, Pmid2 and P3, and distances DS1, DS2 and DS3. A ‘+’ indicates a higher vertex, a ‘-‘ indicates a lower vertex, and a ‘=’ indicates a vertex has the save elevation as the flat triangle. The flat triangle is filled with gray.
Figure 4. A flat triangle has three neighboring triangles, two of them are all higher or all lower than the flat triangle, and the third one is also a flat triangle. The vertex shared by the two non-flat neighboring triangles is modified. A ‘+’ indicates a higher vertex, a ‘-‘ indicates a lower vertex, and a ‘=’ indicates a vertex has the save elevation as the flat triangle. The flat triangle is filled with gray.
Figure 5. A flat triangle has three neighboring triangles, two of them also are flat triangles, and the third one is higher or lower than the flat triangle, then the edge shared with the non-flat triangle is split, and the TIN is adjusted. A ‘+’ indicates a higher vertex, a ‘-‘ indicates a lower vertex, and a ‘=’ indicates a vertex has the same elevation as the flat triangle. The flat triangle is filled with gray.
Figure 6. A flat triangle has two neighboring triangles, one is higher, and the other is lower. Two edges of the flat triangle are split. Three original triangles are split into seven triangles. A ‘+’ indicates a higher vertex, while a ‘-‘ indicates a lower vertex. The flat triangle is filled with gray.
Figure 7. A flat triangle has three neighboring triangles, two of them are all higher or all lower than the flat triangle, and the third one also is a flat triangle. The elevation value of the vertex shared by the two non-flat neighboring triangles is modified. A ‘+’ indicates a higher vertex, while a ‘-‘ indicates a lower vertex. The flat triangle is filled with gray.
Figure 8: The flat triangle has only one neighboring triangle, the other two edges are boundary edges of a TIN. If the neighboring is lower than the flat triangle, a centroid point is added into the flat triangle, the flat triangle is split into three triangles. A ‘-‘ indicates a lower vertex. The flat triangle is filled with gray.
Figure 9: The flat triangle has only one neighboring triangle, the other two edges are boundary edges of a TIN. If the neighboring is higher than the flat triangle, a point is added at the middle of the shared edge, the elevation value is interpolated. A ‘+’ indicates a higher vertex. The flat triangle is filled with gray.
Figure 10. B/T edges in a contour based TIN. The left contour represents a valley while the right contour describes a hill. Flat triangles are formed due to the acceptance of B/T edges during a triangulation.
Figure 11. B/T edge processing. (a) A shared B/T edge is split into two edges by inserting a point at its middle. Four triangles are formed to replace the previous two triangles. (b) a non-shared B/T edge is split in the same way. Two triangles are formed to replace the previous flat triangle.
Figure 12. Eight-direction search for the interpolation of a middle point in a B/T edge. Parameters are labeled for north direction’s (direction 0) interpolation. The furthest point in that direction is used as the original; therefore, the second point is nearer to the added point, and the third point is on the other side of the added point. According to the above definition, we know S0≡0, while S1=P1P0 , and S2=P2P0.
Figure 13. TIN optimization: eliminating B/T edges. (a) a TIN with B/T edges before optimization. (b) a TIN after B/T edge removal process.
Figure 14. Process of flat channels and flat ridges. (a) a flat ridge. (b) a flat channel.
Figure 15. Profile views of flat channels and ridges, and algorithms of a removal process. (a) A single flat channel edge. A middle point is added, TIN is adjusted. (b) A single flat edge is only as a ridge edge. A middle point is added, and the TIN is adjusted. (c) A single flat edge is detected. One ending point is higher than the flat edge, and the other is lower. (d) Two or more flat edges are detected in a ridge network. Vertices in the flat edges are adjusted. The two ending points are both lower than the flat edges. (e) Two or more flat edges are detected in a channel network. Vertices in the flat edges are adjusted. The two ending points are both higher than the flat edges. (f) Two or more flat edges are detected in a channel or ridge network, there are more than two flat points. One ending point is higher than the flat edge, and the other is lower.
(a) (b)
Figure 16. A comparison before and after flat feature processing—a TIN generated based on a discrete point data set. (a) Before flat feature processing. (b) After flat feature processing.
(a) (b)
Figure 17. A comparison before and after flat feature processing—a TIN generated based on a contour data set. (a) Before flat feature processing. (b) After flat feature processing.