{"id":133,"date":"2007-04-26T11:11:20","date_gmt":"2007-04-26T15:11:20","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2007\/04\/26\/upslope-area-part-2\/"},"modified":"2019-10-23T12:53:30","modified_gmt":"2019-10-23T16:53:30","slug":"upslope-area-part-2","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2007\/04\/26\/upslope-area-part-2\/","title":{"rendered":"Upslope area – D-infinity flow"},"content":{"rendered":"
\r\n

A few weeks ago<\/a> I promised a series on developing a MATLAB algorithm to compute upslope area. Now I'm ready to begin that project in earnest.\r\n <\/p>\r\n

First, though, I'd like to comment on the motivation for this series. Specifically, why look at a method from the area of\r\n hydrology in a blog about image processing algorithms? There are two reasons:\r\n <\/p>\r\n

\r\n
    \r\n
  • MATLAB users are in all disciplines of engineering and science, and the ones who use image processing methods use them for\r\n a variety of gridded data sets, not just for processing pictures. Applying image processing techniques to digital elevation\r\n models (DEMs), for example, is common.\r\n <\/li>\r\n
  • Several techniques originally developed for the computational analysis of DEMs have turned out to be broadly useful for processing\r\n other kinds of image data. Perhaps the most well-known example is the watershed transform.\r\n <\/li>\r\n <\/ul>\r\n <\/div>\r\n

    The first thing we need to think about is computing the flow direction<\/i> for a given pixel. That is, based on the height of a pixel and its neighbors, in which direction will water run off from\r\n the center of that pixel? The Tarboton 1997 paper says the earliest method is to say the flow direction is toward whichever of the eight neigbhors has the steepest downward\r\n slope. Because the flow direction can take on one of 8 possible values, this is called the D8 method.\r\n <\/p>\r\n

    Tarboton proposes a new flow direction method, called D-infinity because flow direction assigned to a pixel can take on any\r\n value between 0 and 2*pi. The method examines eight triangular facets for each pixel. Each facet has a vertex at the pixel\r\n center and two vertices at the centers of neighboring pixels.\r\n <\/p>\r\n

    <\/p>\r\n

    The flow direction for a particular facet is the direction of the steepest downward slope on the facet. The function facetFlow<\/tt> below computes this value for the east-northeast facet. Variable names in the code are based on this diagram:\r\n <\/p>\r\n

    <\/p>

    type facetFlow<\/span><\/pre>
    \r\nfunction [r,s] = facetFlow(e0, e1, e2, d1, d2)\r\n% facetFlow Facet flow direction.\r\n% [r, s] = facetFlow(e0, e1, e2, d1, d2) computes the facet flow direction\r\n% for an east-northeast pixel facet.  e0 is the height of the center pixel.\r\n% e1 is the height of the east neighbor.  e2 is the height of the northeast\r\n% neighbor.  d1 is the horizontal pixel center spacing.  d2 is the vertical\r\n% pixel center spacing.  d1 and d2 are optional; if omitted, a value of 1.0\r\n% is assumed.\r\n%\r\n% r is the facet flow direction in radians, and s is the magnitude of the\r\n% corresponding slope.\r\n\r\n% Steve Eddins\r\n% Copyright 2007 The MathWorks, Inc.\r\n\r\n% Reference: Tarboton, \"A new method for the determination of flow\r\n% directions and upslope areas in grid digital elevation models,\" Water\r\n% Resources Research, vol. 33, no. 2, pages 309-319, February 1997.\r\n\r\nif nargin < 4\r\n    d1 = 1;\r\n    d2 = 1;\r\nend\r\n\r\ns1 = (e0 - e1)\/d1;             % eqn (1)\r\ns2 = (e1 - e2)\/d2;             % eqn (2)\r\n\r\nr = atan2(s2, s1);             % eqn (3)\r\ns = hypot(s1, s2);             % eqn (3)\r\n\r\ndiagonal_angle    = atan2(d2, d1);\r\ndiagonal_distance = hypot(d1, d2);\r\n\r\nif r < 0\r\n    % eqn (4)\r\n    r = 0;\r\n    s = s1;\r\n\r\nelseif r > diagonal_angle\r\n    % eqn (5)\r\n    r = diagonal_angle;\r\n    s = (e0 - e2)\/diagonal_distance;\r\nend\r\n\r\n\r\n<\/pre>
    This same computation can be applied to all the facets via rotations or flips. Tarboton provides a lookup table of facet transformations\r\n      that can be used directly in the implementation.  The function pixelFlow<\/tt> below uses facetFlow<\/tt> and the lookup table to find both the flow direction and flow magnitude for a given pixel inside a DEM matrix. The flow direction\r\n      for the pixel is the flow direction with the highest magnitude for each of the eight facets.\r\n   <\/p>
    type pixelFlow<\/span><\/pre>
    \r\nfunction [r, s] = pixelFlow(E, i, j, d1, d2)\r\n% pixelFlow Downslope flow direction from a pixel in a DEM.\r\n% [r, s] = pixelFlow(E, i, j, d1, d2) computes the flow direction for a\r\n% single DEM pixel E(i, j).  d1 is the horizontal pixel center spacing.  d2\r\n% is the vertical pixel center spacing.  d1 and d2 are optional; if\r\n% omitted, a value of 1.0 is assumed.\r\n%\r\n% r is the facet flow direction in radians, and s is the magnitude of the\r\n% corresponding slope.\r\n\r\n% Steve Eddins\r\n% Copyright 2007 The MathWorks, Inc.\r\n\r\n% Reference: Tarboton, \"A new method for the determination of flow\r\n% directions and upslope areas in grid digital elevation models,\" Water\r\n% Resources Research, vol. 33, no. 2, pages 309-319, February 1997.\r\n\r\nif nargin < 4\r\n    d1 = 1;\r\n    d2 = 1;\r\nend\r\n\r\n% Table 1, page 311\r\ne0 = E(i, j);\r\ne1 = [E(i,j+1) E(i-1,j) E(i-1,j) E(i,j-1) ...\r\n      E(i,j-1) E(i+1,j) E(i+1,j) E(i,j+1)];\r\n\r\ne2 = [E(i-1,j+1) E(i-1,j+1) E(i-1,j-1) E(i-1,j-1) ...\r\n      E(i+1,j-1) E(i+1,j-1) E(i+1,j+1) E(i+1,j+1)];\r\n\r\nac = [0  1  1  2  2  3  3  4];\r\naf = [1 -1  1 -1  1 -1  1 -1];\r\n\r\nrp = zeros(1, 8);\r\nsp = zeros(1, 8);\r\n\r\nfor k = 1:8\r\n    [rp(k), sp(k)] = facetFlow(e0, e1(k), e2(k), d1, d2);\r\nend\r\n\r\n% Find the facet with the maximum down-slope.\r\n[s, k_max] = max(sp);\r\n\r\nif s > 0\r\n    % Maximum down-slope is positive.\r\n    r = (af(k_max) * rp(k_max)) + (ac(k_max) * pi \/ 2);\r\nelse\r\n    % The pixel is in a flat area or is a pit.  Flow direction angle is\r\n    % unresolved.\r\n    r = NaN;\r\nend\r\n\r\n    \r\n<\/pre>
    Note: At this point in algorithm development, I am primarily concerned with getting a correct implementation of the method\r\n      described in the paper. It's way too early yet to be worrying about performance.  Also, note the cross-references in my code\r\n      to specific equation and page numbers.  This is a really<\/b> good habit to get into.  (I learned that lesson the hard way.)\r\n   <\/p>\r\n   
    Next up: Try this pixel flow computation on some simple test cases, and start to think about what to do for plateaus.<\/p>