wschwanghart / topotoolbox

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

stream orders with segment length 0 #10

Closed mberry28 closed 5 years ago

mberry28 commented 5 years ago

While attempting to build stream order vs. number of stream and stream order vs. mean stream length plots.

using the following method:

DEM   = GRIDobj(tiff);
FD    = FLOWobj(DEM,'preprocess','carve');
A     = flowacc(FD);
S     = STREAMobj(FD,'minarea',10);
segment  = networksegment(S,FD,DEM,A,0.39);
for j = 1:6
   so_idx           = find(segment.strahler==j);
   mean_lengths(j)  = mean(segment.length(so_idx));
end

This assumes that the indexes of segment.strahler correspond to the indexes of segment.length.

I found that several indexes, particularly of higher order segments have lengths of 0, yet it still registers that a segment of that length exists.

dirkscherler commented 5 years ago

This seems to be related to how the function networksegment (contributed by P. Steer) works. If you execute the following code, you will see that in the case of the Big Tujunga catchment, five segments exist with a length of zero. When displaying where they are, it seems they are confined to tributaries of only one pixel in length. In the function networksegment, segments are computed between different stream points of interest (streampoi). A one-pixel long tributary has a channel head that coincides with the b-confluence, hence you get a length of zero. In other words, the segments appear to be always either a pixelsize or sqrt(2)*pixelsize too short. Not perfect, but probably a minor issue in most cases. If you want to correct it, you would have to make sure that the links (segments) are allways connecting to junctions and not b-confluences. For the moment, you could simply delete entries with a length of zero.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); FD = FLOWobj(DEM,'preprocess','c'); A = flowacc(FD); S = STREAMobj(FD,'minarea',1000); segment = networksegment(S,FD,DEM,A,0.39); plotsegmentgeometry(S,segment) hold on

i0 = find(segment.length==0); IX = segment.IX(ix); [x,y] = ind2coord(DEM,IX); plot(x,y,'ro')

Cheers Dirk

mberry28 commented 5 years ago

Thanks! just wanted to bring attention to the issue