{"id":156,"date":"2007-07-27T14:00:46","date_gmt":"2007-07-27T18:00:46","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2007\/07\/27\/upslope-area-part-7\/"},"modified":"2019-10-23T13:44:05","modified_gmt":"2019-10-23T17:44:05","slug":"upslope-area-part-7","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2007\/07\/27\/upslope-area-part-7\/","title":{"rendered":"Upslope area – Recursive calculation procedure"},"content":{"rendered":"
\r\n

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

  Procedure DPAREA(i, j)\r\n    if AREA(i, j) is known\r\n      then\r\n        no action\r\n      else\r\n        AREA(i, j) = 1 (the area of a single pixel)\r\n        for each neighbor (location in, jn)\r\n          p = proportion of neighbor (in, jn) that drains to\r\n              pixel (i, j) based on angle\r\n          if (p > 0) then\r\n            call DPAREA(in, jn) (this is the recursive call to\r\n                 calculate the area for the neighbor)\r\n            AREA(i, j) = AREA(i, j) + p x AREA(in, jn)\r\n  Return<\/pre>
I had to read the surrounding text a couple of times to figure out exactly how the author intends to calculate p<\/i>.  Consider the flow direction for a particular neighbor, (in, jn).  If the flow direction points directly at pixel (i, j),\r\n      then the corresponding weight p<\/i> 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\r\n      the direction toward (i, j), then the weight p<\/i> 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<\/i> varies proportionally between 1 and 0.\r\n   <\/p>\r\n   
For example, consider the computation of upslope area for pixel number 5 in this set of 9 pixels:<\/p>
   1 4 7\r\n   2 5 8\r\n   3 6 9<\/pre>
If the flow direction for pixel 2 is zero radians, its flow is pointing directly at pixel 5, so the corresponding weight is\r\n      1.0.  If the flow direction for pixel 8 is pi\/4 radians, then its flow is pointing directly at pixel 4. The corresponding\r\n      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\r\n      for the pixel 5 computation is 0.5.\r\n   <\/p>\r\n   
It also took me a few minutes to figure out how to compute the radial difference between two angles, properly taking into\r\n      account the 2*pi periodicity of angles.  (For example, the radial difference between pi\/8 and 2*pi - pi\/8 should be pi\/4.)\r\n       Here's what I came up with:\r\n   <\/p>
   radial_difference = abs(mod(theta1 - theta2 + pi), 2*pi) - pi)<\/pre>
(Does anyone else have a better way to compute this quantity?)  To calculate the weight for a given neighbor of pixel (i,j),\r\n      first determine the angle that points directly from that neighbor to the pixel (i, j).  Call that angle theta_c.  Then find\r\n      the flow direction for that neighbor; call it theta_f.  Now compute the radial difference:\r\n   <\/p>
   radial_difference = abs(mod(theta_c - theta_f + pi), 2*pi) - pi)<\/pre>
And finally compute the weight:<\/p>
   p = max(1 - (radial_difference * 4 \/ pi), 0)<\/pre>
So, are we ready to start writing a recursive MATLAB function based on DPAREA above? Nope! I don't actually want to use a\r\n      recursive formulation at all.\r\n   <\/p>\r\n   
Here's a different way to think about it:<\/p>
  AREA(i,j) = 1 + p1*AREA(i1,j1) + p2*AREA(i2,j2) + ... + p8*AREA(i8,j8)<\/pre>
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.\r\n       Although the system is very large (N-by-N, where N is the number of pixels), it is also very sparse, because each equation\r\n      involves no more than nine of the unknowns.\r\n   <\/p>\r\n   
So we are getting very close to calculating upslope area. We \"just\" have to set up an extremely large sparse system of equations\r\n      and then solve it.\r\n   <\/p>\r\n   
Next time I'll tackle that.<\/p>