wschwanghart / topotoolbox

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

internal drainage stream extraction #9

Closed karimnorouzi closed 6 years ago

karimnorouzi commented 6 years ago

Hi, I am working on a catchment that drains internally. I want to use carving for the streams. I get an error. error I get is this:

fd = FLOWobj(dem, 'preprocess', 'carve','internaldrainage',true); Reference to non-existent field 'MaxIntIX'.

Error in FLOWobj (line 367) ixm = [STATS(r).MaxIntIX];

A fix is highly appreciated.

wschwanghart commented 6 years ago

Hi, thanks for highlighting this error. The 'internaldrainage', true argument can only be used together with the argument 'sink' set to a GRIDobj with sinks. I see that this is currently not sufficiently documented.

Here is an example code: DEM = GRIDobj('Your DEM.tif'); SINKcoords = [x y]; % <- coordinates of locations located in the sinks IX = coord2ind(DEM,x,y); SINKS = GRIDobj(DEM,'logical'); % create a logical grid SINKS.Z(IX) = true; % set sink coordinates to true FD = FLOWobj(DEM,'sinks',SINKS,'internaldrainage',true);

Hope this helps.

karimnorouzi commented 6 years ago

Maybe I am doing something wrong, but I had the sinks too.

dem = GRIDobj('lake1000mCarved.tif'); sinks = dem < 1270; % sinks are below this elevation. fd = FLOWobj(dem, 'preprocess', 'carve', 'sinks', sinks ,'internaldrainage',true); S = STREAMobj(fd,'minarea',50e6,'unit','map'); S.plotdz(dem) % shis still shows sinks in the river streams.

streams

what am I doing wrong ?

wschwanghart commented 6 years ago

Ah, I see. The terminology of preprocess option in FLOWobj are a bit confusing. The 'carve'-option uses an approach that is better suited to carve the DEM at a later stage. However, it actually doesn't carve the DEM that resides in your workspace. This must be done afterwards using the function FLOWobj/imposemin.

What about, if you try this:

dem = GRIDobj('lake1000mCarved.tif'); sinks = dem < 1270; % sinks are below this elevation. fd = FLOWobj(dem, 'preprocess', 'carve', 'sinks', sinks); dem = imposemin(fd,dem); S = STREAMobj(fd,'minarea',50e6,'unit','map'); S.plotdz(dem) % shis still shows sinks in the river streams.

Does this look better?

wschwanghart commented 6 years ago

Super, glad to here this.

karimnorouzi commented 6 years ago

Many thanks, It works quite fine now.