D8 Algorithm for Hydrological Models

Step 1: Find flow direction (find min. potential of all adjacent cells)
Step 2: Follow flow direction to sum the cumulative # cells flowing into
a given "minimum" cell.
Step 3: Find the maximum likely channel location in each E-W direction
for a given number of channels.

NOTE: The user should make sure that the upstream portion of the
potential matrix Psi corresponds to first row of the matrix.

USAGES
[channel,I] = d8alg(Psi)
[channel,I] = d8alg(Psi,numChans)
[channel,I] = d8alg(Psi,Accum, numChans)
[channel,I] = d8alg(Psi,Accum,numChans,'plot')

INPUT:
Psi: An MxN array giving the hydraulic potential over a grid.
Accum: The accumulation matrix. This weights each cell with a fluid
volume. Must be same size as Psi.
numChans: The number of channels to be chosen per column.
plotopt: Gives an option to plot images of the potential and the channel
locations, is a string 'plot'

OUTPUT:
channel: The channel locations are specified by a matrix in which all the
elements without channels are 0 and all the elements with a
channel are equal to the accumulation.
I: the index matrix giving the channel locations.

-------------------------------------------------------------------------
EXAMPLE: See example on d-8 tutorial, provided by the following tutorial:
http://www.niu.edu/landform/nonlinearRules.htm

%make potential function
Psi =[14,13,13,14,11;
11,12,14,13,10;
11,10,11,11, 9;
8, 9, 11,10, 8;
6, 8, 10, 7, 9];

Acc = ones(size(Psi));

%make accumulation function. Compare to website accumulation function.
%Note the last row of the website prediction contains an error. One arrow
%in the predicted flow direction is pointed the wrong way.

[channel,I] = d8alg(Psi,Acc,2,'plot');

%The predicted channel locations match with what you would expect from
%the accumulation plot?