So far in this series I've talked about calculating pixel flow directions and handling plateaus, but I haven't yet discussed the actual upslope area calculation. The Tarboton paper presents a recursive formulation of the upslope area calculation (from page 313):

  Procedure DPAREA(i, j)
    if AREA(i, j) is known
      then
        no action
      else
        AREA(i, j) = 1 (the area of a single pixel)
        for each neighbor (location in, jn)
          p = proportion of neighbor (in, jn) that drains to
              pixel (i, j) based on angle
          if (p > 0) then
            call DPAREA(in, jn) (this is the recursive call to
                 calculate the area for the neighbor)
            AREA(i, j) = AREA(i, j) + p x AREA(in, jn)
  Return

I had to read the surrounding text a couple of times to figure out exactly how the author intends to calculate p. Consider the flow direction for a particular neighbor, (in, jn). If the flow direction points directly at pixel (i, j), then the corresponding weight p is 1. In other words, pixel (i, j) gets all of the flow from (in, jn). If the flow direction is pi/4 or greater away from the direction toward (i, j), then the weight p is 0. Pixel (i, j) gets none of the flow from (in, jn). In between an angular difference of 0 and pi/4, the weight p varies proportionally between 1 and 0.

For example, consider the computation of upslope area for pixel number 5 in this set of 9 pixels:

   1 4 7
   2 5 8
   3 6 9

If the flow direction for pixel 2 is zero radians, its flow is pointing directly at pixel 5, so the corresponding weight is 1.0. If the flow direction for pixel 8 is pi/4 radians, then its flow is pointing directly at pixel 4. The corresponding weight for pixel 5 is 0.0. If the flow direction is pi/8 radians, its flow is pointing between pixels 5 and 4, and the weight for the pixel 5 computation is 0.5.

It also took me a few minutes to figure out how to compute the radial difference between two angles, properly taking into account the 2*pi periodicity of angles. (For example, the radial difference between pi/8 and 2*pi - pi/8 should be pi/4.) Here's what I came up with:

   radial_difference = abs(mod(theta1 - theta2 + pi), 2*pi) - pi)

(Does anyone else have a better way to compute this quantity?) To calculate the weight for a given neighbor of pixel (i,j), first determine the angle that points directly from that neighbor to the pixel (i, j). Call that angle theta_c. Then find the flow direction for that neighbor; call it theta_f. Now compute the radial difference:

   radial_difference = abs(mod(theta_c - theta_f + pi), 2*pi) - pi)

And finally compute the weight:

   p = max(1 - (radial_difference * 4 / pi), 0)

So, are we ready to start writing a recursive MATLAB function based on DPAREA above? Nope! I don't actually want to use a recursive formulation at all.

Here's a different way to think about it:

  AREA(i,j) = 1 + p1*AREA(i1,j1) + p2*AREA(i2,j2) + ... + p8*AREA(i8,j8)

If we write down this equation for each pixel in the image, we'll end up with a system of N linear equations in N unknowns. Although the system is very large (N-by-N, where N is the number of pixels), it is also very sparse, because each equation involves no more than nine of the unknowns.

So we are getting very close to calculating upslope area. We "just" have to set up an extremely large sparse system of equations and then solve it.

Next time I'll tackle that.

Get the MATLAB code

Published with MATLAB® 7.4