Get most upstream nodes that fulfill condition

[Richard-Ott]

I want to locate the most upstream points in a stream network where several conditions from rasters are fulfilled. Below is some dummy code to illustrate the problem:

locrel  = localtopography(DEM, 1000);
elev_min = 300; % minimum elevation of outlets;
elev_max = 700;
lorel_threshold = 100;
stream_elevation = getnal(S,DEM);
stream_relief = getnal(S,locrel);
IX_correct_elevation = and(stream_elevation < elev_max, stream_elevation > elev_min)  ; 
IX_correct_relief = stream_relief < locrel_threshold;
IX_candidates = and(IX_correct_elevation, IX_correct_relief);

Now how would I find the most upstream IX_candidate along every flowpath efficiently? Seems easy but I couldn't come up with a simple way that didn't involve looping over flowpaths. And I'd like to avoid looping over flowpaths because this needs to be computed for a large area.

Ideas appreciated!

[wschwanghart]

Hi Richard,
I think there's no function right now to do this (but we could add one, if needed). However, here's a solution:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM);
S   = STREAMobj(FD,'minarea',300);

locrel  = localtopography(DEM, 1000);
elev_min = 300; % minimum elevation of outlets;
elev_max = 700;
locrel_threshold = 100;
stream_elevation = getnal(S,DEM);
stream_relief = getnal(S,locrel);
IX_correct_elevation = and(stream_elevation < elev_max, stream_elevation > elev_min)  ; 
IX_correct_relief = stream_relief < locrel_threshold;
IX_candidates = and(IX_correct_elevation, IX_correct_relief);

% Now find the most upstream candidates. 
hasUpstream = false(size(S.x));
for r = 1:numel(S.ix)
    if hasUpstream(S.ix(r)) || IX_candidates(S.ix(r)) 
        hasUpstream(S.ixc(r)) = true;
    end
end
IX_candidates_up = IX_candidates;
IX_candidates_up(hasUpstream) = false;

% Let's check the results
plotdz(S,DEM)
hold on
z = getnal(S,DEM);
plot(S.distance(IX_candidates_up),z(IX_candidates_up),'sk',...
    'MarkerFaceColor','r')

The approach requires one scan in downstream direction and is reasonably fast. Hope this solves the issue.
Cheers, Wolfgang

[Richard-Ott]

Danke Wolfgang!!!

[wschwanghart]

Gerne ;-)