VIEWS: 30 PAGES: 42 POSTED ON: 11/29/2011
Advanced Terrain Analysis Concepts Specific Catchment Area The D Surface Flow Model Topmodel Terrain based calculations of saturated areas and soil moisture deficit Generalized flow accumulation functions Flow algebra Extending ArcGIS Scripts and Programming Drainage area can be concentrated or dispersed (specific catchment area) representing concentrated or dispersed flow. Area defining concentrated contributing area at P Contour width b P Specific catchment Flow path originating area is A/b at divide with dispersed contributing area A Specific catchment area a is the upslope area per unit contour length [m2/m m] Stream line Specific Catchment Area a = A/b Contour line Unit contour length b Contributing area A The D Algorithm Proportion Steepest direction flowing to downslope neighboring Proportion flowing to grid cell 4 is neighboring grid cell 3 1/(1+2) is 2/(1+2) 3 2 4 2 1 Flow direction. 5 1 6 8 7 Tarboton, D. G., (1997), "A New Method for the Determination of Flow Directions and Contributing Areas in Grid Digital Elevation Models," Water Resources Research, 33(2): 309-319.) (http://www.engineering.usu.edu/cee/faculty/dtarb/dinf.pdf) Contributing Area using D Contributing Area using D8 TOPMODEL Beven, K., R. Lamb, P. Quinn, R. Romanowicz and J. Freer, (1995), "TOPMODEL," Chapter 18 in Computer Models of Watershed Hydrology, Edited by V. P. Singh, Water Resources Publications, Highlands Ranch, Colorado, p.627-668. “TOPMODEL is not a hydrological modeling package. It is rather a set of conceptual tools that can be used to reproduce the hydrological behaviour of catchments in a distributed or semi- distributed way, in particular the dynamics of surface or subsurface contributing areas.” Hydrological processes within a catchment are complex, involving: • Macropores • Heterogeneity • Fingering flow • Local pockets of saturation The general tendency of water to flow downhill is however subject to macroscale conceptualization Map of saturated areas showing expansion during a single rainstorm. The solid black shows the saturated area at the beginning of the rain; the lightly shaded area is saturated Seasonal variation in pre-storm by the end of the storm and is the area over saturated area [from Dunne and which the water table had risen to the ground Leopold, 1978] surface. [from Dunne and Leopold, 1978] TOPMODEL Key Ideas • Surface saturation and soil moisture deficits based on topography – Slope – Specific Catchment Area – Topographic Convergence • Partial contributing area concept • Saturation from below (Dunne) runoff generation mechanism Topmodel - Assumptions Specific catchment area a [m2/m m] • The soil profile at each point has a finite (per unit contour length) capacity to transport water laterally downslope based on slope and transmissivity, qcap=T tan. • The actual lateral discharge is proportional to specific catchment area, q=Ra • Relative wetness at a point and depth to water table is determined by comparing q and qcap • The dynamics of the saturated zone can be approximated by successive steady tan state representations involving the depth to water table adjusting to accommodate q that results in expansion and contraction of the variable source area contributing to saturation excess runoff TOPMODEL Slope and GIS Specific Catchment Area Wetness Index ln(a/tan) from Raster Calculator. Average, l = 6.9 Numerical Example Given Compute • Ko=10 m/hr z 0.46 m • R=0.0002 m/h • f=5 m-1 • Qb = 0.8 m3/s • l=6.9 1 a tan l z z ln • a (from GIS) • T=2 m2/hr f • ne = 0.2 Raster calculator -( [ln(sca/S)] - 6.9)/5+0.46 Depth to saturation z Flat (0.5%) -3 - 0 (7.8%) 0 - 0.1 (2.5%) 0.1 - 0.2 (4.0%) 0.2 - 0.5 (29%) 0.5 - 1 (56%) 1 - 1.5 (0.2%) Generalized Flow Accumulation Continuum Definition. A[ r ( x )] r (x )dx Contributing Area A[.] is a functional operator that takes as input a spatial field r(x), and the topographic flow direction field (not denoted) and produces a field A(x) representing the accumulation of r(x) up to each point x. Numerical evaluation A[r(x)] = A(i,j) = r( i, j)2+ p A(i , jk ) k k contributing neighbors k pk is the proportion of flow from neighbor k contributing to the grid cell (i,j). p k =1 is required to ensure „conservation‟. Flow directions must not have loops. Downslope Influence The contributing area only of points in a target set y. I(x|y) says what the contribution from points y are at each mapped point x. I(x|y)=A[i(x|y)] Grid cell y 1 0 0 0.5 0.5 0 0 1 0 0 0.6 0.4 0 0.6 0.4 Influence function of grid cell y Useful for example to track where sediment or contaminant moves Upslope Dependence. Quantifies the amount a point x contributes to the point or zone y. The inverse of the influence function D(x|y) = I(y|x) 0.6 0 0 0.3 0.3 0 0 0.6 0 0 1 0 0 1 0 Grid cells y Dependence function of grid cells y Useful for example to track where a contaminant may come from Decaying Accumulation A decayed accumulation operator DA[.] takes as input a mass loading field m(x) expressed at each grid location as m(i, j) that is assumed to move with the flow field but is subject to first order decay in moving from cell to cell. The output is the accumulated mass at each location DA(x). The accumulation of m at each grid cell can be numerically evaluated DA[m(x)] = DA(i, j) = m(i, j)2 + pkd(ik , jk )DA(ik , jk ) k contributing neighbors Here d(x) = d(i ,j) is a decay multiplier giving the fractional (first order) reduction in mass in moving from grid cell x to the next downslope cell. If travel (or residence) times t(x) associated with flow Useful for a tracking between cells are available d(x) may be evaluated as exp(lt( x )) where l is a contaminant or compound first order decay parameter. subject to decay or attenuation Concentration Limited Accumulation An unlimited supply of a substance that is loaded into flow at a concentration or solubility threshold Csol. The set of points y, delineating the area of the substance supply are mapped using the (0,1) indicator field i(x;y). 1. Flow (weighted accumulation) Q(x)=A[w(x)] 2. Supply area concentration at threshold If i(x; y) = 1 C(x) = Csol L(x) = Csol Q(x) 3. Remaining locations by load accumulation and dilution L(x) = L(i, j) = pk L(ik , jk ) k contributing neighbors C(x) = L(x)/Q(x) Useful for a tracking a contaminant released or partitioned to flow at a fixed threshold concentration Transport limited accumulation Erodability Transport Capacity Transport Flux, Deposition, e.g. E = a0.7 S0.6 e.g. Tcap = a2 S2 T D Tout min( E Tin , Tcap ) D E Tin Tout Useful for modeling erosion and sediment delivery, the spatial dependence of sediment delivery ratio and contaminant that adheres to sediment Reverse Accumulation 1.35 1.8 0 0.9 1.8 0 0 1.8 0 0 1.2 0.6 0.4 0 0.8 0.6 0.8 0.6 Reverse accumulation of field weights indicated in red Useful for destabilization sensitivity in landslide hazard assessment with Bob Pack Extending ArcGIS via programming Why Programming Automation of repetitive tasks Implementation of functionality not available Programming functionality Scripts (AML, VB, Python, Avenue) Interfaces for application programmers Model Builder ArcObjects COM Integration Three Views of GIS Geodatabase view: Structured data sets that represent geographic information in terms of a generic GIS data model. Geovisualization view: A GIS is a set of intelligent maps and other views that shows features and feature relationships on the earth's surface. "Windows into the database" to support queries, analysis, and editing of the information. Geoprocessing view: Information transformation tools that derives new geographic data sets from existing data sets. adapted from www.esri.com Examples • TauDEM – ArcMap toolbar using Visual Basic/C++ (implementation of functionality not available in ArcGIS) • Visual Basic Programming of simple grid calculations TauDEM Software Functionality Pit removal (standard flooding approach) Flow directions and slope D8 (standard) D (Tarboton, 1997, WRR 33(2):309) Flat routing (Garbrecht and Martz, 1997, JOH 193:204) Drainage area (D8 and D) Network and watershed delineation Support area threshold/channel maintenance coefficient (Standard) Combined area-slope threshold (Montgomery and Dietrich, 1992, Science, 255:826) Local curvature based (using Peuker and Douglas, 1975, Comput. Graphics Image Proc. 4:375) Threshold/drainage density selection by stream drop analysis (Tarboton et al., 1991, Hyd. Proc. 5(1):81) Wetness index and distance to streams Water Quality Functions TauDEM in ArcGIS Visual Basic Visual Basic Standalone ESRI ArcGIS 8.x+9.x MapWindow plugin command line Toolbar applications C++ COM DLL interface Available from Fortran (legacy) http://www.engineering.usu.edu/dtarb/ TauDEM C++ library components Data Access libraries (grid, shape, image, dbf) ESRI gridio API (Spatial analyst) Data Vector shape ASCII text Binary direct ESRI binary formats files grid access grid grid Implementation Details Spatial Analyst includes a C programming API (Application Programming Interface) that allows you to read and write ESRI grid data sets directly. Excerpt from gioapi.h / * GetWindowCell - Get a cell within the window for a layer, * Client must interpret the type of the output 32 Bit Ptr * to be the type of the layer being read from. * * PutWindowCell - Put a cell within the window for a layer. * Client must ensure that the type of the input 32 Bit Ptr * is the type of the layer being read from. * */ int GetWindowCell(int channel, int rescol, int resrow, CELLTYPE *cell); int PutWindowCell(int channel, int col, int row, CELLTYPE cell); C++ COM Methods used to implement functionality using Microsoft Visual C++ STDMETHODIMP CtkTauDEM::Areadinf(BSTR angfile, BSTR scafile, long x, long y, int doall, BSTR wfile, int usew, int contcheck, long *result) { USES_CONVERSION; //needed to convert from BSTR to Char* or String *result = area( OLE2A(angfile), OLE2A(scafile), x,y,doall, OLE2A(wfile), usew, contcheck); return S_OK; } Visual Basic for the GUI and ArcGIS linkage Private TarDEM As New tkTauDEM … Private Function runareadinf(Optional toadd As Boolean = False) As Boolean Dim i As Long runareadinf = False i = TarDEM.Areadinf(tdfiles.ang, tdfiles.sca, 0, 0, 1, "", 0, 1) If TDerror(i) Then Exit Function If toadd Then AddMap tdfiles.sca, 8 End If runareadinf = True End Function Using TauDEM - Exercise Calculating the distance to a drain point (e.g. watershed outlet or gage) • Introduce VB scripting • Introduce Recursion • Isolate the Watershed draining to a point • Could be done with the Flow Path function (maybe), but a programmed solution is instructive and can be generalized to other applications Distance from each grid cell to outlet along flow path Write program to do this 27.5 27.5 27.4 27.5 31.2 38 44 52 29.9 28.6 28.7 31.3 35.5 39.8 43.7 47.3 36.6 32.4 32.5 33.4 36.9 42.8 48.3 51.1 46.5 41.5 38.8 36.5 39.1 46.7 52.7 59.3 53.7 50.3 47.6 42.9 40.3 44.7 52.7 64.5 60.1 58.4 56.9 51.7 45.4 43.6 50.7 63.1 30 Distances 30 30 30 30 30 30 30 30 30 30 30 30 30 Summing the distances down from each grid cell 30 30 30+42.4 72.4 102.4 30+30+42.4 132.4 30+30+42.4+30 Summing the distances down from each grid cell 0 30 30 30+42.4 60 72.4 30+30+42.4 90 100.4 114.8 120 132.4 144.8 174.8 204.8 30+30+42.4+30 Number of additions N P Recursive Approach 30 30 30+42.4 72.4 102.4 30+72.4 132.4 30+102.4 Recursive Approach 0 30 30 30+42.4 60 72.4 30+72.4 90 100.4 114.8 120 132.4 144.8 174.8 204.8 30+102.4 Number of additions N Recursive Procedure DistCalc (i,j) For each neighbor (location in, jn) If neighbor (in, jn) drains to cell (i,j) Distance to outlet = Distance from me to outlet + Distance from neighbor to me D(in,jn)=D(i,j)+(in,jn) Call recursive procedure on the neighbor, DistCalc(in, jn) endif end for loop This requires (assumes) that the distance for cell i,j is initialized or has been calculated to start the process Programming 4 3 2 the calculation 5 1 of distance to 6 7 8 the outlet Direction encoding 1 2 3 1 72.4 102.4 7 6 5 2 30 42.4 72.4 7 6 5 3 0 6 7 7 Distances to outlet Visual Basic Implementation 'RECURSIVE DISTANCE CALCULATION FUNCTION Sub DistCalc(i, j) Dim k As Integer, inb As Long, jnb As Long For k = 1 To 8 ' for each neighbor inb = i + di(k) jnb = j + dj(k) If (inb >= 0 And inb < nrow And jnb >= 0 And jnb < ncol) Then ' guard ' against out of domain If pPixels(jnb, inb) > 0 Then ' guard against no data If (pPixels(jnb, inb) - 4 = k Or pPixels(jnb, inb) + 4 = k) Then ' Here we have a grid cell that drains back to the grid cell we are at dPixels(jnb, inb) = dPixels(j, i) + dd(k) DistCalc inb, jnb ' Call the function for the neighbor pixel End If End If End If Next k End Sub Steps for distance to outlet program • Read the outlet coordinates • Read the DEM flow direction grid. This is a set of integer values 1 to 8 indicating flow direction • Initialize a distance to outlet grid with a no data value • Convert outlet to row and column references • Start from the outlet point. Set the distance to 0. • Call the DistCalc function at the outlet. Distance Calculation VBA Exercise Visual Basic Programming in ArcMAP References ESRI, (1999), ArcObjects Developers Guide: ArcInfo 8, ESRI Press, Redlands, California. Zeiler, M., (2001), Exploring ArcObjects. Vol 1. Applications and Cartography. Vol 2. Geographic Data Management, ESRI, Redlands, CA. Are there any questions ? AREA 2 3 AREA 1 12