wschwanghart / topotoolbox

A MATLAB software for the analysis of digital elevation models -
https://topotoolbox.wordpress.com/
153 stars 89 forks source link

over sampling number of stream orders #11

Closed mberry28 closed 5 years ago

mberry28 commented 5 years ago

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?

dirkscherler commented 5 years ago

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))

mberry28 commented 5 years ago

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

Thank you so much!

mberry28 commented 5 years ago

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