{"id":135,"date":"2007-05-18T18:18:37","date_gmt":"2007-05-18T22:18:37","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2007\/05\/18\/upslope-area-part-3\/"},"modified":"2019-10-23T12:55:18","modified_gmt":"2019-10-23T16:55:18","slug":"upslope-area-part-3","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2007\/05\/18\/upslope-area-part-3\/","title":{"rendered":"Upslope area – Flow examples"},"content":{"rendered":"
\r\n

In Part 2<\/a> of the upslope area series, I showed an M-file, pixelFlow<\/tt>, that could calculate a pixel's flow direction and magnitude according to the D-infinity method of Tarboton.\r\n <\/p>\r\n

Let's try the function with a simple test case first, and then we'll try it with some DEM data.<\/p>

Z = peaks;\r\nmesh(peaks)<\/pre> 
pixelFlow<\/tt> works on a single pixel at a time, so we put it into a loop to calculate the flow direction, r<\/tt>, and the flow magnitude, s<\/tt>, for all the interior elements of Z<\/tt>.\r\n   <\/p>
[M,N] = size(Z);\r\nr = zeros(M-2, N-2);\r\ns = zeros(M-2, N-2);\r\nfor<\/span> i = 1:M-2\r\n    for<\/span> j = 1:N-2\r\n        [r(i,j), s(i,j)] = pixelFlow(Z,i+1,j+1);\r\n    end<\/span>\r\nend<\/span><\/pre>
Now display the peaks matrix as a grayscale image and superimpose a quiver plot.<\/p>
imshow(peaks,[],'initialmag'<\/span>,'fit'<\/span>)\r\nhold on<\/span>\r\n[x,y] = meshgrid(2:N-1,2:M-1);\r\nquiver(x, y, s.*cos(r), -s.*sin(r), 2, 'y'<\/span>)\r\naxis equal<\/span>\r\nhold off<\/span><\/pre> 
Zoom in on a low point:<\/p>
xlim([20 35]);\r\nylim([5 20])<\/pre> 
And a high point:<\/p>
xlim([30 40])\r\nylim([15 30])<\/pre> 
To get a real digital elevation model (DEM) to play with, I started at a page that the Mapping Toolbox development team fondly\r\n      calls Tech Support Note 2101: \"Accessing Geospatial Data on the Internet for the Mapping Toolbox.\"<\/a> That page directed me to GeoCommunity, one of the distribution points for free DEM data from the United States Geological Survey (USGS)<\/a>.  Downloading\r\n      DEM data set 1643091 and using the Mapping Toolbox function sdtsdemread<\/tt> gave me this data set:\r\n   <\/p>
s = load('upslope_sample_dem'<\/span>);\r\nimshow(s.Zm, [], 'initialmag'<\/span>, 'fit'<\/span>)<\/pre> 
pond = s.Zm(140:280, 160:230);\r\nimshow(pond, [], 'initialmag'<\/span>, 'fit'<\/span>)\r\ntitle('North Pond'<\/span>)<\/pre> 
Let's try pixelFlow<\/tt> on this data set.\r\n   <\/p>
[M,N] = size(pond);\r\nr = zeros(M-2, N-2);\r\ns = zeros(M-2, N-2);\r\nfor<\/span> i = 1:M-2\r\n    for<\/span> j = 1:N-2\r\n        [r(i,j), s(i,j)] = pixelFlow(pond,i+1,j+1);\r\n    end<\/span>\r\nend<\/span>\r\n\r\nimshow(pond,[],'initialmag'<\/span>,'fit'<\/span>)\r\nhold on<\/span>\r\n[x,y] = meshgrid(2:N-1,2:M-1);\r\nquiver(x, y, s.*cos(r), -s.*sin(r), 2, 'y'<\/span>)\r\naxis equal<\/span>\r\n\r\n% Zoom in on the southern end of the pond.<\/span>\r\nxlim([15 60])\r\nylim([75 105])\r\n\r\nhold off<\/span><\/pre> 
That seems to be working pretty well, but there's an unresolved problem to explore - what to do about plateaus.  Recall the\r\n      last few lines of pixelFlow<\/tt>:\r\n   <\/p>
dbtype pixelFlow<\/span> 47:51<\/span><\/pre>
\r\n47    else\r\n48        % The pixel is in a flat area or is a pit.  Flow direction angle is\r\n49        % unresolved.\r\n50        r = NaN;\r\n51    end\r\n\r\n<\/pre>
For example:<\/p>
Zp = repmat([1:3 3 3 3:6],5,1)<\/pre>
\r\nZp =\r\n\r\n     1     2     3     3     3     3     4     5     6\r\n     1     2     3     3     3     3     4     5     6\r\n     1     2     3     3     3     3     4     5     6\r\n     1     2     3     3     3     3     4     5     6\r\n     1     2     3     3     3     3     4     5     6\r\n\r\n<\/pre>
r = zeros(size(Zp)-2);\r\nfor<\/span> i = 1:size(r,1)\r\n    for<\/span> j = 1:size(r,2)\r\n        r(i,j) = pixelFlow(Zp, i+1, j+1);\r\n    end<\/span>\r\nend<\/span>\r\n\r\nr<\/pre>
\r\nr =\r\n\r\n    3.1416    3.1416       NaN       NaN       NaN    3.1416    3.1416\r\n    3.1416    3.1416       NaN       NaN       NaN    3.1416    3.1416\r\n    3.1416    3.1416       NaN       NaN       NaN    3.1416    3.1416\r\n\r\n<\/pre>
The pixel flow direction is to the left (r<\/tt> = pi) except for the middle pixels, for which pixelFlow<\/tt> could not compute a direction.\r\n   <\/p>\r\n   
I'll pursue the issue of plateaus next week.<\/p>
Near the center of this data set is North Pond<\/a>, close to where I live and about 3 miles southwest of the Boston Marathon starting line.  (The high spot west of the pond\r\n      is Peppercorn Hill, a great place for a Cub Scout hike.)\r\n   <\/p>