Exercise 4: GIS/DEM Programming

Hydro GIS Short Course 
University of Padua 
Spring 2000

Prepared by David Tarboton, Utah State University.
The purpose of this exercise is to learn how to write computer programs to work with GRID DEMs.

To get started you will need the following:


Develop a FORTRAN or C program to compute the distance to the outlet from each grid cell.  The inputs should be:
  • Name of dem flow direction grid file to work with (e.g. reydemp.asc), recognizing that it is only the flow directions that are needed to compute distances.  This will need to be exported from ArcView grid reydemp.
  • X and Y coordinates of outlet point.
  • Suggested algorithm

  • Read the DEM flow direction grid.  This is a set of integer values 1 to 8 indicating flow direction
  • Read the outlet coordinates. Convert these to row and column references.
  • Initialize a distance to outlet grid as -1.
  • Start from the outlet point
  • Set the distance to 0.
  • Call a recursive procedure at the outlet point, DISTANCE(rowout,colout)
  • end

  • Recursive Procedure DISTANCE(i,j)
    do for each neighbor (location in, jn)

    end do

    A recursive procedure is one that calls itself.  This was not possible in early versions of FORTRAN, but it is now.  It is possible in C.

    To understand the logic behind this consider the following set of flow directions:

    The outlet is the red cell (3,1) at the bottom left. The program starts at this cell and sets its distance to 0.  The program then examines each neighbor of this cell.  The first neighbor examined is in direction k=1, cell (3,2).  This cell has direction dir(3,2)=7, and dir(3,2)-k = 5.  This is not 4 so the neighbor does not drain to the outlet.  The program then examines the neighbor in direction k=2, cell (2,2).  This cell has direction dir(2,2)=6, and dir(2,2)-k=4, indicating that this neighbor does drain to the outlet.  The flow distance from this neighbor is then evaluated as dist(2,2)=dist(3,3)+distance from (3,3) to (2,2) = 0 + 1.4 = 1.4.  The procedure is now called at cell (2,2).  This temporarily suspends operation of the procedure at cell (3,3).  All necessary variables are placed on a stack for later use.  At cell (2,2) the procedure  finds grid cell (2,3) that drains to it, so evaluates the distance from it as dist(2,3) = dist(2,2) + distance from (2,3) to (2,2) = 1.4 + 1 = 2.4.  The procedure is now called at cell (2,3), temporarily suspending operation for the cell(2,2), placing all necessary variables on a stack.  Cell (2,3) has no grid cells draining to it, so the procedure is not called again and ends, returning control to the suspended operation at cell (2,2). The necessary variables saved on a stack are retrieved from the stack.  The terminology for this is that the recursion is unravelling.  No more cells draining to cell (2,2) are found, so the procedure ends, returning control to the suspended operation at cell (3,1).  The procedure at cell (3,1) next examines the neighbor in direction k=3 and so on.  It ends up tracing up through cells (2,1) then (1,2) and (1,3) before unravelling and ending.  Note that only the cells in blue, that actually drain to cell (3,1) are visited.  Note that this recursive algorithm appears to be in danger of getting into infinite loops, because it calls itself repeatedly.  For it to terminate the structure has to be such that the function can end, and the nesting of calls unravel.  This is the case for the algorithm above as long as there are NO LOOPS.  If cell (3,1) drained to cell (3,2), and this drained to (2,2) which drained back to (3,1) an infinite loop would result.  Therefore it is very important that the directions be determined properly.  The programs flood and d8 have put a lot of effort in to this, but they are still not foolproof when existing channels that may be inconsistent with DEM elevations are used.

    Some FORTRAN tips
    Declare the flow direction grid a short integer to save space

    integer*2 dir(irows,icols)
    Then read it using the following
    call sfgridread(filename,dir,nx,ny,dx,dy,bndbox,csize,xndv,irows)

    Use data arrays d1 and d2 to define neighbors

    integer*2 d1(8), d2(8)
    data d1/0,-1,-1,-1,0,1,1,1/
    data d2/1,1,0,-1,-1,-1,0,1/
    Then a loop over each neighbor can be done as
    do k = 1,8
    The test as to whether a neighbor in direction k drains towards the cell i,j can be done as
    if(abs(k-dir(i+d1(k),j+d2(k))) .eq. 4)then
    The logic of the above is that since there are eight possible directions a differerence of 4 in absolute value means the oposite direction.
    Distances between a cell i,j and the neighbor in the k direction can be calculated as
    distk=sqrt((dy*d1(k))**2+(dx*d2(k))**2)    !  Pythagoras

    Summary of material to turn in

    1. a layout showing the distance to the outlet computed using your program
    2. a listing of your source code
    Ok, you're done!

    These materials may be used for study, research, and education, but please credit the authors and the Utah Water Research Laboratory, Utah State University. All commercial rights reserved. Copyright 2000 Utah State University.