I need to create a chiprofile with m/n 0.45, 0.50 and 0.60 for the whole drainage network and highlight the prolonged stream profile.

[negipr]

Dear Wolfgang,
I changed the mn value and created slope-area plot , but whenever I am trying to plot the chi-profile it is not showing the right plot.

DEM_FIL = fillsinks(DEM);
FD = FLOWobj(DEM_FIL,'preprocess','carve');
A = flowacc(FD);

DEM_resolution = DEM.refmat(3,2);
area_threshold = 1e6;
minApix = area_threshold/(DEM_resolution*DEM_resolution);

minApix = ceil(minApix);
A_W = A > minApix;
imageschs(DEM, log10(A_W));
minarea = 1e6;
W = A.*(A.cellsize.^2)>minarea;
STR_W = STREAMobj(FD, A_W);
STR = STREAMobj(FD, W)

min_str_gradient = 1e-3;
mg = imposemin(FD,DEM,min_str_gradient);
rivers_slope = gradient(STR,mg,'unit','tangent');
rivers_area = A.Z(STR.IXgrid).*(A.cellsize).^2;

DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false);

figure()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
grid on
hold on
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)

figure()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
hold on
DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false,'theta', 0.45);
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)

c = chiplot(STR,DEM,A);
c = chiplot(STR,DEM,A,'mnplot',false,'fitto',....
'ts','trunkstream',trunk(STR),'a0',1e6);

figure ()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
hold on
DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false,'theta', 0.5);
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)
c = chiplot(STR,DEM,A);
c = chiplot(STR,DEM,A,'mnplot',false,'fitto',....
'ts','trunkstream',trunk(STR),'a0',1e6);

figure ()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
hold on
DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false,'theta', 0.6);
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)
c = chiplot(STR,DEM,A);
c = chiplot(STR,DEM,A,'mnplot',false,'fitto',....
'ts','trunkstream',trunk(STR),'a0',1e6);

[wschwanghart]

Hi,

I am copying your code and comment it below:

% Don't use fillsinks before running FLOWobj. FLOWobj handles sinks internally
DEM_FIL = fillsinks(DEM);
FD = FLOWobj(DEM_FIL,'preprocess','carve');
A = flowacc(FD);

DEM_resolution = DEM.refmat(3,2);
area_threshold = 1e6;
minApix = area_threshold/(DEM_resolution*DEM_resolution);

minApix = ceil(minApix);
A_W = A > minApix;
imageschs(DEM, log10(A_W));
minarea = 1e6;
W = A.*(A.cellsize.^2)>minarea;
STR_W = STREAMobj(FD, A_W);
STR = STREAMobj(FD, W)

min_str_gradient = 1e-3;
mg = imposemin(FD,DEM,min_str_gradient);
rivers_slope = gradient(STR,mg,'unit','tangent');
rivers_area = A.Z(STR.IXgrid).*(A.cellsize).^2;

% Median may be more robust but sometimes causes problems if there are many slope values with 
% similar values. This may be a problem, if you use imposemin as done above, because it might 
% create many slope values of 1e-3
DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false);

figure()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
grid on
hold on
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)

figure()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
hold on
DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false,'theta', 0.45);
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)

% Chiplot is ok, but actually it is better to use chitransform together with plotdz. 
% This gives you more flexibility.
c = chiplot(STR,DEM,A);
c = chiplot(STR,DEM,A,'mnplot',false,'fitto',....
'ts','trunkstream',trunk(STR),'a0',1e6);

% Here's the same code with chitransform
mn = 0.5; % Determine mn using mnoptim or mnoptimvar
c = chitransform(STR,DEM,A,'mn',mn);
plotdz(STR,DEM,'distance',c)

% Include above lines in your code below. this should work

figure ()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
hold on
DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false,'theta', 0.5);
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)
c = chiplot(STR,DEM,A);
c = chiplot(STR,DEM,A,'mnplot',false,'fitto',....
'ts','trunkstream',trunk(STR),'a0',1e6);

figure ()
loglog(rivers_area, rivers_slope,'.','color',[0.7 0.7 0.7])
hold on
DEM_slopearea = slopearea(STR,mg,A,'areabinlocs','median','gradaggfun',...
'median','streamgradient','robust','plot',false,'theta', 0.6);
loglog(DEM_slopearea.a, DEM_slopearea.g, 'ks', 'markersize' ,5)
aeval = logspace(log10(min(DEM_slopearea.a)),log10(max(DEM_slopearea.a)),10);
geval = DEM_slopearea.ks(1)*aeval.^DEM_slopearea.theta;
loglog(aeval, geval,'k-','linewidth', 2)
c = chiplot(STR,DEM,A);
c = chiplot(STR,DEM,A,'mnplot',false,'fitto',....
'ts','trunkstream',trunk(STR),'a0',1e6);

[negipr]

Thank you for your response. yes, chi-transform has been worked out.
How can we highlight the trunk stream?

STR = klargestconncomps(STR, 1);

St = trunk(STR);
c = chitransform(STR,DEM,A,'mnplot',false,'fitto','ts','trunkstream',St,'a0',1e6)

[wschwanghart]

You can access node-attribute lists (such as c) with the function nal2nal. Hence, no need to calculate chi again.

c = chitransform(STR,DEM,A,'mn',mn);
plotdz(STR,DEM,'distance',c)
hold on
St = trunk(STR)
plotdz(St,DEM,'distance',nal2nal(St,STR,c),'LineWidth',2)
hold off

[negipr]

Thank you so much for your response.