Terrain Analysis by yPsYn3x

VIEWS: 30 PAGES: 42

• pg 1
```									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
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[.]
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)
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
 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.

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
Private TarDEM As New tkTauDEM
…

Boolean
Dim i As Long
i = TarDEM.Areadinf(tdfiles.ang, tdfiles.sca, 0, 0, 1, "", 0, 1)
If TDerror(i) Then Exit Function
End If
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

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 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

```
To top