Terrain Analysis by yPsYn3x


									Advanced Terrain Analysis
            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


  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

                                         5                                  1

                                         6                                 8

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

                          Contributing Area using

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,

“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,
                                       • 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
                                  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 contributing neighbors

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
        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
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)
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
                             Reverse Accumulation
        1.35    1.8          0

         0.9    1.8          0

         0      1.8         0

         0     1.2          0.6
         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
            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
• TauDEM – ArcMap toolbar using Visual
  Basic/C++ (implementation of
  functionality not available in ArcGIS)
• Visual Basic Programming of simple grid
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
       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

                                          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

Private Function runareadinf(Optional toadd As Boolean = False) As
  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
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         30   30
            30   30    30   30
            30   30    30   30
   Summing the distances down from
            each grid cell


30+42.4               72.4

   Summing the distances down from
            each grid cell

30+42.4                  60   72.4

30+30+42.4               90   100.4 114.8

                        120   132.4 144.8 174.8 204.8

 Number of additions N  P
           Recursive Approach


30+42.4                 72.4

           Recursive Approach


30+42.4                60    72.4

30+72.4                90    100.4 114.8

                       120   132.4 144.8 174.8 204.8

           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
       Call recursive procedure on the neighbor,
          DistCalc(in, jn)
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
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
• 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

ESRI, (1999), ArcObjects Developers Guide:
  ArcInfo 8, ESRI Press, Redlands,
Zeiler, M., (2001), Exploring ArcObjects.
  Vol 1. Applications and Cartography. Vol
  2. Geographic Data Management, ESRI,
  Redlands, CA.
Are there any questions ?

                    AREA 2


                  AREA 1


To top