Gradient of a multiple line 2d plot,

1 次查看(过去 30 天)
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

采纳的回答

Mathieu NOE
Mathieu NOE 2022-6-14
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 个评论

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Annotations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by