wschwanghart/topotoolbox

over sampling number of stream orders

mberry28 opened this issue · 3 comments

As mentioned in a previous issue, I am working on plotting frequency vs. stream order and mean length vs stream order. As you know they tend to have slopes similar to this:
image

However, in using "networksegment.m" and extracting the length and segment values, the produced plots are either flat (all segments are ~similar lengths) or negatively sloped (higher order streams are shorter) in the stream order vs mean length. I have tested this on several landscapes with similar results.

I notice that, for example, networksegment.m is counting 3 segments of stream order 6. however, when it is plotted in map view, there is only 1 true segment of stream order 6; but there is 2 breaks caused by addition stream junctions, which seems to be splitting the true segment length. This would certainly underpredict the length of the segment of higher orders, but oversample the number of segments. Perhaps this is a feature and not a bug. Is there a workaround to this?

I haven't worked with the function networksegment, but I guess it's splitting streams at junctions. So for the Strahler stream ordering scheme, you would be splitting streams of the same order. To avoid this, try out the following:

minarea = 1e5;

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','c');
FA = flowacc(FD).*FD.cellsize.*FD.cellsize;
ST = STREAMobj(FD,'minarea',minarea./FA.cellsize.^2);

% Bifurcation ratio of stream network
MS = STREAMobj2mapstruct(ST);
so = [MS.streamorder];
ct = 0;
clear s sn
for i = min(so):max(so)
ix = so==i;
if nnz(ix)>0
ct = ct+1;
s(ct) = i;
sn(ct) = nnz(ix);
end
end
cumsn = fliplr(cumsum(fliplr(sn)));

figure
semilogy(s,cumsn,'ro')
xlabel('Stream order')
ylabel('N(SO>SO'')')

bf = mean(sn(1:end-1)./sn(2:end));
bfsig = std(sn(1:end-1)./sn(2:end));
fprintf('Bifurcation ratio = %1.2f +/- %1.2f\n',bf,bfsig/sqrt(numel(sn)-1))

My initial test looks like this improved the results, I will test this with the length and report back.

Thank you so much!

Hi
Below is my approach to getting the mean stream distance from your suggested work through. I didn't see a length listed in the structure files. I did a few tests to verify it was working correctly

Feel free to add it to the matlab scripts. Or not. Either way ;)
Thanks

for ii = 1:7; 
    so = zeros(6,1);
    %location of the folder
    tiff = strcat(files_dir(ii).folder,'/',files_dir(ii).name);   
 DEM           = GRIDobj(tiff);
    FD            = FLOWobj(DEM,'preprocess','c');
    FA            = flowacc(FD).*FD.cellsize.*FD.cellsize;
    ST            = STREAMobj(FD,'minarea',10);
    MS            = STREAMobj2mapstruct(ST);
    SO            = [MS.streamorder];
    ct            = 0;
    clear s sn
    
    for i = 1:6;
        ix = SO==i;     %find index of a particular stream order
        if nnz(ix)>0    %number of non zero elements
%            ct = ct+1;
%            s(ct) = i;
            so(i) = nnz(ix);
        end
        %index location of the specific segments
        idx = find(SO == i);
        lenSO = nnz(ix);
        D  = zeros(lenSO,1);
        for j = 1:lenSO;
            k = idx(j);
            %obtain x and y values of the segment, remove the NAN
            X   = MS(k).X(~(isnan(MS(k).X)));
            Y   = MS(k).Y(~(isnan(MS(k).Y)));
            D(j)   = sum(hypot(diff(X),diff(Y))); %find the distance
        end
            mean_lengths(i)  = mean(D);
            if mean_lengths(i) == 0; %if the mean length of a segment is zero, then there was no segment anyway
                so(i) = 0;
            end
            
    end
    cumsn = fliplr(cumsum(fliplr(so)));
end