Spatial Data - MATLAB Recipes for Earth Sciences - page 288

for i = 1:8
flow(:,:,i) = flow(:,:,i).*(flow(:,:,i)>0)./upssum;
end
h e 2D matrix inflowsum will store the intermediate inl ow totals for each
iteration. h ese intermediate totals are then summed to reach a i gure for the
total accumulated l ow flowac at the end of each iteration. h e initial values
for inflowsum and flowac are obtained through upssum .
inflowsum = upssum;
flowac = upssum;
Another 3D matrix inflow is now needed, in which to store the total
intermediate inl ow from all neighbors:
inflow = grads*0;
Flow accumulation is terminated when there is no inl ow or, translating
this into MATLAB code, we use a conditional while loop that terminates
if sum(inflowsum(:))==0 . h e number of non-zero entries in inflowsum will
decrease during each loop. h is is achieved by alternately updating inflow
and inflowsum . Here, inflowsum is updated with the intermediate inflow of
the neighboring cells weighted by flow , under the condition that all of the
neighboring cells are contributing cells, i.e., where grads are positive. Where
not all neighboring cells are contributing cells, the intermediate inflowsum is
reduced, as also is inflow . h e l ow accumulation flowac increases through
consecutive summations of the intermediate inflowsum .
while sum(inflowsum(:))>0
for i = 1:8
inflow(:,:,i) = circshift(inflowsum,[N(i,1) N(i,2)]);
end
inflowsum = sum(inflow.*flow.*grads>0,3);
flowac = flowac + inflowsum;
end
We display the result as a pseudocolor map with log-scaled values (Fig.
7.15e).
h = pcolor(log(1+flowac));
colormap(flipud(jet)), colorbar
set(h,'LineStyle','none')
axis equal
title('Flow accumulation')
[r c] = size(flowac);
axis([1 c 1 r])
set(gca,'TickDir','out');