# Upslope area – Flow examples

In Part 2 of the upslope area series, I showed an M-file, `pixelFlow`, that could calculate a pixel's flow direction and magnitude according to the D-infinity method of Tarboton.

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

Z = peaks; mesh(peaks)

`pixelFlow` works on a single pixel at a time, so we put it into a loop to calculate the flow direction, `r`, and the flow magnitude, `s`, for all the interior elements of `Z`.

[M,N] = size(Z); r = zeros(M-2, N-2); s = zeros(M-2, N-2); for i = 1:M-2 for j = 1:N-2 [r(i,j), s(i,j)] = pixelFlow(Z,i+1,j+1); end end

Now display the peaks matrix as a grayscale image and superimpose a quiver plot.

imshow(peaks,[],'initialmag','fit') hold on [x,y] = meshgrid(2:N-1,2:M-1); quiver(x, y, s.*cos(r), -s.*sin(r), 2, 'y') axis equal hold off

Zoom in on a low point:

xlim([20 35]); ylim([5 20])

And a high point:

xlim([30 40]) ylim([15 30])

To get a real digital elevation model (DEM) to play with, I started at a page that the Mapping Toolbox development team fondly
calls Tech Support Note 2101: "Accessing Geospatial Data on the Internet for the Mapping Toolbox." That page directed me to GeoCommunity, one of the distribution points for free DEM data from the United States Geological Survey (USGS). Downloading
DEM data set 1643091 and using the Mapping Toolbox function `sdtsdemread` gave me this data set:

s = load('upslope_sample_dem'); imshow(s.Zm, [], 'initialmag', 'fit')

Near the center of this data set is North Pond, close to where I live and about 3 miles southwest of the Boston Marathon starting line. (The high spot west of the pond is Peppercorn Hill, a great place for a Cub Scout hike.)

pond = s.Zm(140:280, 160:230); imshow(pond, [], 'initialmag', 'fit') title('North Pond')

Let's try `pixelFlow` on this data set.

[M,N] = size(pond); r = zeros(M-2, N-2); s = zeros(M-2, N-2); for i = 1:M-2 for j = 1:N-2 [r(i,j), s(i,j)] = pixelFlow(pond,i+1,j+1); end end imshow(pond,[],'initialmag','fit') hold on [x,y] = meshgrid(2:N-1,2:M-1); quiver(x, y, s.*cos(r), -s.*sin(r), 2, 'y') axis equal % Zoom in on the southern end of the pond. xlim([15 60]) ylim([75 105]) hold off

That seems to be working pretty well, but there's an unresolved problem to explore - what to do about plateaus. Recall the
last few lines of `pixelFlow`:

dbtype pixelFlow 47:51

47 else 48 % The pixel is in a flat area or is a pit. Flow direction angle is 49 % unresolved. 50 r = NaN; 51 end

For example:

Zp = repmat([1:3 3 3 3:6],5,1)

Zp = 1 2 3 3 3 3 4 5 6 1 2 3 3 3 3 4 5 6 1 2 3 3 3 3 4 5 6 1 2 3 3 3 3 4 5 6 1 2 3 3 3 3 4 5 6

r = zeros(size(Zp)-2); for i = 1:size(r,1) for j = 1:size(r,2) r(i,j) = pixelFlow(Zp, i+1, j+1); end end r

r = 3.1416 3.1416 NaN NaN NaN 3.1416 3.1416 3.1416 3.1416 NaN NaN NaN 3.1416 3.1416 3.1416 3.1416 NaN NaN NaN 3.1416 3.1416

The pixel flow direction is to the left (`r` = pi) except for the middle pixels, for which `pixelFlow` could not compute a direction.

I'll pursue the issue of plateaus next week.

**범주:**- Upslope area

## 댓글

댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.