Gradient of a multiple line 2d plot,

5 views (last 30 days)
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

Accepted Answer

Mathieu NOE
Mathieu NOE on 14 Jun 2022
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'));
  4 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!