Efficient Algorithms for Large-Scale GIS Applications

1. Efficient Algorithms for Large-Scale GIS Applications Laura Toma Duke University

2. Why GIS? • How it all started.. • Duke Environmental researchers: • computing flow accumulation for Appalachian Mountains took 14 days (with 512MB memory) • 800km x 800km at 100m resolution  ~64 million points • GIS (Geographic Information Systems) • System that handles spatial data • Visualization, processing, queries, analysis… • Rich area of problems for Computer Science • Graphics, graph theory, computational geometry, scientific computing…

3. Indispensable tool Monitoring:keep an eye on the state of earth systems using satellites and monitoring stations (water, pollution, ecosystems, urban development,…) Modeling and simulation: predict consequences of human actions and natural processes Analysis and risk assessment:find the problem areas and analyse the possible causes (soil erosion, groundwater pollution,..) Planning and decision support: provide information and tools for better management of resources Precipitation in Tropical South America Lots of rain Dry Nitrogen in Chesapeake Bay High nitrogen concentrations GIS and the Environment

4. GIS and the Environment Bald Head Island Renourishment Sediment flow

5. Reality: Elevation of terrain is a continuous function of two variables h(x,y) Estimate, predict, simulate Flooding, pollution Erosion, deposition Vegetation structure …. GIS: DEM (Digital Elevation Model) is a set of sample points and their heights { x, y, hxy} Computations on Terrains Model and compute indices

6. 3 3 3 3 2 2 2 2 4 4 4 4 7 7 7 7 5 5 5 5 8 8 8 8 7 7 7 7 1 1 1 1 9 9 9 9 DEM Representations TIN Grid Contour lines Sample points

7. Modeling Flow on Terrains • What happens when it rains? • Predict areas susceptible to floods. • Predict location of streams. • Compute watersheds. • Flow is modeled using two basic attributes • Flow Direction (FD) • The direction water flows at a point • Flow Accumulation (FA) • Total amount of water that flows through a point (if water is distributed according to the flow directions)

8. DEM and Flow Accumulation[Panama]

11. Uses Flow direction and flow accumulation are used for: • Computing other hydrological attributes • river network • moisture indices • watersheds and watershed divides • Analysis and prediction of sediment and pollutant movement in landscapes. • Decision support in land management, flood and pollution prevention and disaster management

12. Remote sensing technology Massive amounts of terrain data Higher resolutions (1km, 100m, 30m, 10m, 1m,…) NASA-SRTM Mission launched in 2001 Acquired data for 80% of earth at 30m resolution 5TB USGS Most of US at 10m resolution LIDAR 1m res Massive Terrain Data

13. Example: LIDAR Terrain Data • Massive (irregular) point sets (1-10m resolution) • Relatively cheap and easy to collect Example: Jockey’s ridge (NC coast)

14. It’s Growing! • Appalachian Mountains Area if approx. 800 km x 800 km Sampled at: • 100m resolution:  64 million points (128MB) • 30m resolution:  640 (1.2GB) • 10m resolution:  6400 = 6.4 billion (12GB) • 1m resolution:  600.4 billion (1.2TB)

15. Computing on Massive Data • GRASS (open source GIS) • Killed after running for 17 days on a 6700 x 4300 grid (approx 50 MB dataset) • TARDEM (research, U. Utah) • Killed after running for 20 days on a 12000 x 10000 grid (appox 240MB dataset) • CPU utilization5%, 3GB swap file • ArcInfo (ESRI, commercial GIS) • Can handle the 240MB dataset • Doesn’t work for datasets bigger than 2GB

16. Outline • Introduction • Flow direction and flow accumulation • Definitions, assumptions, algorithm outline • Scalability to large terrains • I/O-model • I/O-efficient flow accumulation • TerraFlow • I/O-efficient graph algorithms • Conclusion and future work

17. Flow Direction (FD) on Grids • Water flows downhill • follows the gradient • On grids: Approximated using 3x3 neighborhood • SFD (Single-Flow Direction): • FD points to the steepest downslope neighbor • MFD (Multiple-Flow direction) : • FD points to all downslope neighbors

18. Flow accumulation with MFD

19. Flow accumulation with SFD

20. Computing FD • Goal: compute FD for every cell in the grid (FD grid) • Algorithm: • For each cell compute SFD/MFD by inspecting 8 neighbor cells • Analysis: O(N) time for a grid of N cells • Is this all? • NO! flat areas: Plateas and sinks

21. FD on Flat Areas • …no obvious flow direction • Plateaus • Assign flow directions such that each cell flows towards the nearest spill point of the plateau • Sinks • Either catch the water inside the sink • Or route the water outside the sink using uphill flow directions • model steady state of water and remove (fill) sinks by simulating flooding: uniformly pouring water on terrain until steady state is reached • Assign uphill flow directions on the original terrain by assigning downhill flow directions on the flooded terrain

22. Flow Accumulation (FA) on Grids FA models water flow through each cell with “uniform rain” • Initially one unit of water in each cell • Water distributed from each cell to neighbors pointed to by its FD • Flow conservation: If several FD, distribute proportionally to height difference • Flow accumulation of cell is total flow through it Goal: compute FA for every cell in the grid (FA grid)

23. Computing FA • FD graph: • node for each cell • (directed) edge from cell a to b if FD of a points to b • FD graph must be acyclic • ok on slopes, be careful on plateaus • FD graph depends on the FD method used • SFD graph: a tree (or a set of trees) • MFD graph: a DAG (or a set of DAGs)

24. Computing FA: Plane Sweeping • Input: flow direction grid FD • Output: flow accumulation grid FA (initialized to 1) • Process cells in topological order. For each cell: • Read its flow from FA grid and its direction from FD grid • Update flow for downslope neighbors (all neighbors pointed to by cell flow direction) • Correctness • One sweep enough • Analysis • O(sort) + O(N) time for a grid of N cells • Note: Topological order means decreasing height order (since water flows downhill).

25. Dataset Size (log) Scalability Problem • We can compute FD and FA using simple O(N)-time algorithms • ..but.. for large sets..??

26. read/write head read/write arm track magnetic surface Scalability Problem: Why? • Most (GIS) programs assume data fits in memory • minimize only CPU computation • But.. Massive data does not fit in main memory! • OS places data on disk and moves data between memory and disk as needed • Disk systems try to amortize large access time by transferring large contiguous blocks of data • When processing massive data disk I/O is the bottleneck, rather than CPU time!

27. Disks are Slow “The difference in speed between modern CPU and disk technologies is analogous to the difference in speed in sharpening a pencil using a sharpener on one’s desk or by taking an airplane to the other side of the world and using a sharpener on someone else’s desk.” (D. Comer)

28. Scalability to Large Data • Example: reading an array from disk • Array size N = 10 elements • Disk block size = 2 elements • Memory size = 4 elements (2 blocks) 1 5 2 6 3 8 9 4 7 10 1 2 10 9 5 6 3 4 8 7 Algorithm 1: Loads 10 blocks Algorithm 2: Loads 5 blocks N blocks >> N/B blocks • Block size is large (32KB, 64KB) N >> N/B N = 256 x 106, B = 8000 , 1ms disk access time  NI/Os take 256 x 103 sec = 4266 min = 71 hr  N/BI/Os take 256/8 sec = 32 sec

29. I/O-model I/O-operation Read/write one block of data from/to disk I/O-complexity number of I/O-operations (I/Os) performed by the algorithm External memory or I/O-efficient algorithms: Minimize I/O-complexity • RAM model • CPU-operation • CPU-complexity • Number of CPU-operations performed by the algorithm • Internal memory algorithms: • Minimize CPU-complexity

30. I/O-Model D Parameters N = # elements in problem instance B = # elements that fit in disk block M = # elements that fit in main memory Fundamental bounds: • Scanning: scan(N) = • Sorting: sort(N) = Block I/O M P In practiceblock and main memory sizes are big

31. Sorting • Mergesort illustrates often used features: • Main memory sized chunks (for N/M runs) • Multi-way merge (repeatedly merge M/B of them)

32. I/O-Efficient Priority Queue • Insert • insert in B0; when B0 full, empty it recursively on Bi • extract_min • extract_min from PQ; if PQ empty, fill it from Bi

33. Computing FAI/O-Analysis Algorithm:O(N) time • Process (sweep) cells in topological order. For each cell: • Read flow from FA grid and direction from FD grid • Update flow in FA grid for downslope neighbors • Problem: Cells in topological order distributed over the terrain  scattered access to FA gridandFD gridO(N) blocks

34. [ATV00] I/O-Efficient Flow Accumulation • Eliminating scattered accesses to FD grid • Store FD grid in topological order • Eliminating scattered accesses to FA grid • Obs: flow to neighbor cell is only needed at the time when the neighbor is processed: • Time when cell is processed Topological rank priority • Push flow by inserting flow increment in priority queue with priority equal to neighbor’s time • Flow of cell obtained using DeleteMin operations • Note: Augment each cell with priorities of 8 neighbors • Obs: Space (~9N) traded for I/O • Use I/O-efficient priority queue [A95,BK97] • O(N) operations in I/Os

35. TerraFlow • TerraFlow is our suite of programs for flow routing and flow accumulation on massive grids[ATV`00,AC&al`02] • Flow routing and flow accumulation modeled as graph problems and solved in optimal I/O bounds • Efficient • 2-1000 times faster on very large grids than existing software • Scalable • 1 billion elements!! (>2GB data) • Flexible • Allows multiple methods flow modeling http://www.cs.duke.edu/geo*/terraflow

36. Implementation and Platform • C++ • Uses TPIE (Transparent Parallel I/O Environment) • Library of I/O-efficient modules developed at Duke • http://www.cs.duke.edu/~tpie/ • Platform • TerraFlow, ArcInfo • 500MHz Alpha, FreeBSD 4.0, 1GB RAM • GRASS/TARDEM • 500MHz Intel PIII, FreeBSD/Windows, 1GB RAM

37. TerraFlow • GRASS cannot handle Hawaii dataset (killed after 17 days) • TARDEM cannot handle Cumberlands dataset (killed after 20 days) • Significant speedup over ArcInfo for large datasets • East-Coast TerraFlow: 8.7 Hours ArcInfo: 78 Hours • Washington state TerraFlow: 63 Hours ArcInfo: %

38. I/O-Efficient Graph Algorithms Graph G=(V,E) on disk • Basic graph (searching) problems • BFS, DFS, SSSP, topological sorting • ..are big open problems in the I/O-model! • Standard internal memory algorithms: O(E) I/Os • No I/O-efficient algorithms are known for any of these problems on general graphs! • Lower bound Ω (sort(V)), best known Ω (V/sqrt(B)) • O(sort(E)) algorithms for special classes of graphs • Trees, grid graphs, bounded-treewidth graphs, outerplanar graphs, planar graphs • Exploit existence of small separators or geometric structure

39. Dijkstra’s Algorithm in External Memory • Dijkstra’s algorithm • Use a priority queue to store <u, d(s,u)> • Repeat: DeleteMin(u) and for each adjacent edge (u,v) • if d(s,v) > d(s,u) + wuv then DecreaseKey(v, d(s,u) + wuv) • Analysis: O(E) I/Os • O(V+E/B) I/Os to load adjacent edges for each vertex • Use an I/O-efficient priority queue • O(sort(E)) I/Os for O(E)Insert/DeleteMin • DecreaseKey: O(1) I/Os to read key d(s,v) of v • O(E) + O(sort(E)) I/Os • Improved SSSP algorithm: • O(V + E/B log V) I/Os [KS’96]

40. SSSP on Grid Graphs [ATV’00] Grid graph O(N) vertices, O(N) edges • Previous bound: O(N) I/Os • Lemma The portion of δ(s,t)between intersection points with boundaries of any subgrid is the shortest path within the subgrid

41. SSSP on Grid Graphs [ATV’00] Idea: Compute shortest paths locally in each subgrid then compute the shortest way to combine them together • Divide grid into subgrids of size M (assumeM > B2) • Replace each subgrid with complete graph on its boundary vertices • Edge weight: shortest path between the two boundary vertices within the subgrid • Reduced graph GR O(N/B) vertices, O(N) edges

42. SSSP on Grid Graphs [ATV’00] Algorithm • Compute SSSP on GRfrom s to all boundary vertices • Find SSSP from s to all interior vertices: for any subgrid σ, for any t in σ δ(s,t) = min v in Bnd(σ){δ(s,v) + δ σ(v,t)} • Correctness: • easy to show using Lemma • Analysis: O(sort(N)) I/Os • Dijkstra’s algorithm using I/O efficient priority queue and graph blocking

44. DFS O(sort(N)) I/Os [AMTZ’01] O(sort(N)) I/Os [ABT’00] O(sort(N)) I/Os [ABT’00] ε-separators BFS SSSP Results on Planar graphs Planar graph G with N vertices • Separators can be computed in O(sort(N)) I/Os • I/O-efficient reductions [ABT’00, AMTZ’01]  BFS, DFS, SSSP in O(sort(N)) I/Os

45. SSSP on Planar Graphs • Similar with grid graphs. AssumeM > B2, bounded degree • Assume graph is separated • O(N/B2) subgraphs, O(B2) vertices each, S=O(N/B) separators • each subgraph adjacent to O(B) separators

46. SSSP on Planar Graphs • Reduced graph GR • S = O(N/B) vertices • O(N/B2) x O(B2) = O(N) edges • Compute SSSP on GR • Dijkstra’s algorithm and I/O-efficient priority queue • Keep list L S = {d(s,v) | v in S} • For each vertex in S read from L Sits O(B) adjacent boundary vertices O(N/B x B) = O(N) I/Os • ..butO(N/B2)boundary sets • Store L Ssuch that vertices in same boundary set are consecutive • Each boundary set is accessed once by its O(B) adjacent vertices  O(N/B2x B) = O(N/B)I/Os

47. On I/O-Efficient DFS • General graphs • Internal memory algorithm O(V+E) I/Os • Improved upper bound: O(V + E/B log V) I/Os [KS’96] • Lower bound Ω (sort(V)) • DFS is a big open problem!! • Note: PRAM DFS is P-complete • DFS on planar graphs • O(sort(N) log N) I/Os algorithm [AMTZ’01] • O(sort(N)) I/Os using DFS to BFS reduction [AMTZ’01]

48. A Divide-and-Conquer DFS Algorithm on Planar Graphs • Based on PRAM algorithm [Smith86] • Algorithm • Compute a cycle separator C • Compute DFS recursively in the connected components Gi of G\C • Attach the DFS trees of Gionto the cycle • I/O-analysis • O(log N) recursive steps • O(sort(N)) I/Os per step O(sort(N) log N) I/Os in total

49. DFS to BFS Reduction on Planar Graphs Idea: Partition the faces of G into levels around a source face containing s and grow DFS level-by-level • Levels can be obtained from BFS in dual graph • Denote • Gi = union of the boundaries of faces at level <= i • Ti= DFS tree of Gi • Hi = Gi \ G i-1 • Algorithm: Compute a spanning forest of Hi and attach it onto T i-1 • Structure of levels is simple • The bicomps of the Hi are the boundary cycles of Gi • Glueing onto T i-1 is simple • A spanning tree is a DFS tree if and only if it has no cross edges

50. DFS to BFS Reduction on Planar Graphs Idea: Partition the faces of G into levels around a source face containing s and grow DFS level-by-level