Gradient of a multiple line 2d plot,

Hi, i have the following code, which generates some plots:

%% Punktauswertung nach Nutzung

%%%Flussbett
clear all
load('Flussbett.mat')
C = table2array(Flussbett);
CC =C;
CC(CC==0)=nan;
c=1:17

%% Wasserhöhen
 for i=1:17
Wasserhoehe (i,:)= CC([i],[4:537])- CC(i,3);
 end
Wasserhoehe(Wasserhoehe<=0)=nan;

D_BHQ_WH= Wasserhoehe([c],[1:84]);
D_HQextrem_WH= Wasserhoehe([c],[85:238]);
D_HQ5000_WH= Wasserhoehe([c],[239:392]);
D_HQ10000_WH= Wasserhoehe([c],[393:534]);

A= [0:141];

%Maximale Wasserhöhen 
for i=1:17
HQ10000_MWH (i,1)= max(D_HQ10000_WH(i,[1:142]));
HQ5000_MWH (i,1)= max(D_HQ5000_WH(i,[1:154]));
HQextrem_MWH (i,1)= max(D_HQextrem_WH(i,[1:154]));
BHQ_MWH (i,1)= max(D_BHQ_WH(i,[1:84]));
end

%% Anstieg pro Zeitschritt [m/h]

Wasserhoehe(isnan(Wasserhoehe))=0;
for i=2:534
    Wasserhoehe_Wasseranstieg(:,i)= Wasserhoehe([c],i)- Wasserhoehe([c],i-1);
end

D_BHQ_WA= Wasserhoehe_Wasseranstieg([c],[1:84])
D_HQextrem_WA= Wasserhoehe_Wasseranstieg([c],[85:238])
D_HQ5000_WA= Wasserhoehe_Wasseranstieg([c],[239:392])
D_HQ10000_WA= Wasserhoehe_Wasseranstieg([c],[393:534])

%Maximaler Wasseranstieg
for i=1:17
D_HQ10000_MWA (i,1)= max(D_HQ10000_WA(i,[1:142]));
D_HQ5000_MWA (i,1)= max(D_HQ5000_WA(i,[1:154]));
D_HQextrem_MWA (i,1)= max(D_HQextrem_WA(i,[1:154]));
D_BHQ_MWA (i,1)= max(D_BHQ_WA(i,[1:84]));
end

%Maximaler Wasserrückgang
for i=1:17
D_HQ10000_MWR (i,1)= min(D_HQ10000_WA(i,[1:142]));
D_HQ5000_MWR (i,1)= min(D_HQ5000_WA(i,[1:154]));
D_HQextrem_MWR (i,1)= min(D_HQextrem_WA(i,[1:154]));
D_BHQ_MWR (i,1)= min(D_BHQ_WA(i,[1:84]));
end

%% Überflutungsdauer

Wasserhoehe(isnan(Wasserhoehe))=0;
D_BHQ_x= Wasserhoehe([c],[1:84]);
D_HQextrem_x= Wasserhoehe([c],[85:238]);
D_HQ5000_x= Wasserhoehe([c],[239:392]);
D_HQ10000_x= Wasserhoehe([c],[393:534]);

for i=1:17
D_HQ10000_FD (i,1)= nnz(D_HQ10000_x(i,[1:142]));
D_HQ5000_FD (i,1)= nnz(D_HQ5000_x(i,[1:154]));
D_HQextrem_FD (i,1)= nnz(D_HQextrem_x(i,[1:154]));
%D_BHQ_FD (i,1)= nnz(D_BHQ_x(i,[1:84]));
end

%% Tabelle erstellen 
% load('Pegelpunkte1DSpatialJoinTN.mat')
% varNames={'OID', 'FKM', 'Lage','Landnutzung',  'X', 'Y', 'Höhe', 'Maximale Wasserhöhe', 'Maximaler Wasseranstieg', 'Maximaler Wasserrückgang', 'Überflutungsdauer'}
% %BHQ= table(PP_D.OID_,PP_D.FKM,PP_D.Lage,PP_D.X,PP_D.Y,PP_D.z,D_BHQ_MWH,D_BHQ_MWA,D_BHQ_MWR,'VariableNames',varNames)
% HQ_extrem= table(PP_D.OID_,PP_D.FKM,PP_D.Lage,Pegelpunkte1DSpatialJoinTN.objart_1, PP_D.X,PP_D.Y,PP_D.z,D_HQextrem_MWH,D_HQextrem_MWA,D_HQextrem_MWR,D_HQextrem_FD, 'VariableNames',varNames)
% HQ_5000= table(PP_D.OID_,PP_D.FKM,PP_D.Lage,Pegelpunkte1DSpatialJoinTN.objart_1,PP_D.X,PP_D.Y,PP_D.z,D_HQ5000_MWH,D_HQ5000_MWA,D_HQ5000_MWR,D_HQ5000_FD, 'VariableNames',varNames)
% HQ_10000= table(PP_D.OID_,PP_D.FKM,PP_D.Lage,Pegelpunkte1DSpatialJoinTN.objart_1,PP_D.X,PP_D.Y,PP_D.z,D_HQ10000_MWH,D_HQ10000_MWA,D_HQ10000_MWR,D_HQ10000_FD, 'VariableNames',varNames)


%% FRI
D_HQ10000_WH(isnan(D_HQ10000_WH))=0;
D_HQ5000_WH(isnan(D_HQ5000_WH))=0;
D_HQextrem_WH(isnan(D_HQextrem_WH))=0;
D_BHQ_WH(isnan(D_BHQ_WH))=0;

%Wassertiefe Ih
href=6
for t=1:84
  IhBHQ([c],t)= 1- D_BHQ_WH([c], [t]) /href;  
end
for t=1:154
Ih5000([c],t)= 1-D_HQ5000_WH([c], [t]) /href;
Ihextrem([c],t)= 1- D_HQextrem_WH([c], [t]) /href;
end

for  t=1:142
   Ih10000([c],t)= 1- D_HQ10000_WH([c], [t]) /href; 
end

% Ih10000(Ih10000 >=1)=nan;
% Ih5000(Ih5000 >=1)=nan;
% Ihextrem(Ihextrem >=1)=nan;
% IhBHQ(IhBHQ >=1)=nan;
% 
% Ih10000(Ih10000 <=-1)=0;
% Ih5000(Ih5000 <=-1)=0;
% Ihextrem(Ihextrem <=-1)=0;
% IhBHQ(IhBHQ <=-1)=0;
% 
%Wasseranstiegsrate
WARref= 2
for t=1:84
  IwarBHQ([c],t)= 1- D_BHQ_WA([c], [t]) /WARref;
end

for t=1:154
Iwar5000([c],t)= 1- D_HQ5000_WA([c], [t]) /WARref;
Iwarextrem([c],t)= 1- D_HQextrem_WA([c], [t]) /WARref;
end
for  t=1:142
   Iwar10000([c],t)= 1- D_HQ10000_WA([c], [t]) /WARref;
end
% Iwar10000(Iwar10000 >=1)=nan;
% Iwar5000(Iwar5000 >=1)=nan;
% Iwarextrem(Iwarextrem >=1)=nan;
% IwarBHQ(IwarBHQ >=1)=nan;
% Iwar10000(Iwar10000 >=-1)=0;
% Iwar5000(Iwar5000 >=-1)=0;
% Iwarextrem(Iwarextrem >=1)=0;
% IwarBHQ(IwarBHQ >=-1)=0;
% 
%FRI

for t=1:84
FRI_BHQ([c],t)= (3*IhBHQ([c],t)+2*IwarBHQ([c],t))/(2+3);
end
for t=1:154 
FRI_5000([c],t)= (3*Ih5000([c],t)+2*Iwar5000([c],t))/(2+3);
FRI_extrem([c],t)= (3*Ihextrem([c],t)+2*Iwarextrem([c],t))/(2+3);
end

for  t=1:142
 FRI_10000([c],t)= (3*Ih10000([c],t)+2*Iwar10000([c],t))/(2+3);
end

FRI_BHQ(FRI_BHQ<=0)=0;
FRI_extrem(FRI_extrem<=0)=0;
FRI_5000(FRI_5000<=0)=0;
FRI_10000(FRI_10000<=0)=0;
FRI_BHQ(FRI_BHQ>=1)=1;
FRI_extrem(FRI_extrem>=1)=1;
FRI_5000(FRI_5000>=1)=1;
FRI_10000(FRI_10000>=1)=1;

A= [0:28];
B=[0:141];
X=[0:153]
Y=[0:83]

for i=1:17
figure (1)
plot(B ,FRI_10000([i],[1:142]));
grid on 
hold on 
title('FRI HQ10000');
xlabel('Time [h]');
ylabel('FRI []');
set(gcf, 'Position', get(0, 'Screensize'));
end

for i=1:17
figure (2)
plot(X ,FRI_5000([i],[1:154]));
hold on 
grid on 
title('FRI HQ5000')
xlabel('Time [h]');
ylabel('FRI []');
set(gcf, 'Position', get(0, 'Screensize'));
end
for i=1:17
figure (3)
plot(X ,FRI_extrem([i],[1:154]));
hold on 
grid on 
title('FRI HQextrem')
xlabel('Time [h]');
ylabel('FRI []');
set(gcf, 'Position', get(0, 'Screensize'));
end

for i=1:17
plot(Y ,FRI_BHQ([i],[1:84]))
hold on 
grid on 
title('FRI BHQ')
xlabel('Time [h]');
ylabel('FRI []');
set(gcf, 'Position', get(0, 'Screensize'));
end

I want to calculate the gradient at the end of the lines and I dont want the little jumps to count.

Attached a little image of what i want for a value.

Thanks for your help

Answer 1

[Mathieu NOE]

hello

I see you have a passion for for loops !

maybe you can do yours plots without a for loop this time .

Added some code for the end slope computation and display :

I don't put the beginning of your code , this is one example for the last plot :

% plot FRI_BHQ
figure (4) 
    % end slope calculation on last 10 samples
    dx_mean = mean(diff(Y)); % 
    FRI_BHQ_end = FRI_BHQ(:,end-9:end); % last 10 samples
    dy_mean = mean(diff(FRI_BHQ_end,[],2),2); % diff on 2nd dimension , mean on 2nd dimension
    slope = dy_mean./dx_mean; % by definition
    % draw it 
    Y_draw = Y(:,end-9:end); % last 10 samples
    FRI_BHQ_draw = FRI_BHQ(:,end-9)+slope*(Y_draw-Y_draw(1)); % 10 samples before end point + slope = line

% plot
plot(Y ,FRI_BHQ,Y_draw,FRI_BHQ_draw,'k','linewidth',2.5)
grid on 
title('FRI BHQ')
xlabel('Time [h]');
ylabel('FRI []');
set(gcf, 'Position', get(0, 'Screensize'));

[Frederik Reese]

Thanks. Great job. How can I extend the line until it touches the 1 on the y axis ?

[Mathieu NOE]

here you are

% plot FRI_BHQ
figure (4) 
    % end slope calculation on last 10 samples
    dx_mean = mean(diff(Y)); % 
    FRI_BHQ_end = FRI_BHQ(:,end-9:end); % last 10 samples
    dy_mean = mean(diff(FRI_BHQ_end,[],2),2); % diff on 2nd dimension , mean on 2nd dimension
    slope = dy_mean./dx_mean; % by definition
    % draw it 
%     Y_draw = Y(:,end-9:end); % last 10 samples
    Y_draw = Y(:,end-9):200 ; % from en-9 sample to extension max until it reaches y = 1 
    FRI_BHQ_draw = FRI_BHQ(:,end-9)+slope*(Y_draw-Y_draw(1)); % 10 samples before end point + slope = line

% plot
plot(Y ,FRI_BHQ,Y_draw,FRI_BHQ_draw,'k','linewidth',2.5)
grid on 
title('FRI BHQ')
xlabel('Time [h]');
ylabel('FRI []');
set(gcf, 'Position', get(0, 'Screensize'));

[Frederik Reese]

thanks

[Mathieu NOE]

my pleasure !