To get started you will need the following:

• Computer with FORTRAN or C Compiler that you know how to use
• Example Fortran code or C code for working with DEM's.  This reads a digital elevation model grid, does some checking for gaps then writes the result.
• FORTRAN ASCII grid input and output functions, qgio.f.
• C gridio.c, gridio_null.c, gridio.h and gioapi.h files.

Exercise

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)

If neighbor (in, jn) drains to cell (i,j)
Distance from (in, jn) is distance from (i,j) + distance between cells (accounting for possible diagonals)
Call recursive procedure on the neighbor, DISTANCE(in, jn)
endif

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
e.g.

parameter(irows=1000,icols=1000)
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

in=i+d1(k)
jn=j+d2(k)
...

enddo

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

...

endif

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