forked from topicster/iMHEA_scripts
-
Notifications
You must be signed in to change notification settings - Fork 1
/
iMHEA_Pair.m
192 lines (174 loc) · 4.94 KB
/
iMHEA_Pair.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
function [Indices,Climate,PM,QM,FDC1,FDC2,IDC1,IDC2] = iMHEA_Pair(Date1,P1,Q1,A1,Date2,P2,Q2,A2)
%iMHEA Pair-catchments comparisons.
% [Indices,Climate] = iMHEA_Pair(Date1,P1,Q1,A1,Date2,P2,Q2,A2) processes
% data from paired catchments and produce indices ands plots.
%
% Input:
% Date = dd/mm/yyyy hh:mm:ss [date format].
% P = Precipitation [mm].
% Q = Discharge [l/s].
% A = Catchment area [km2].
%
% Output:
% Indices = Matrix of hydrological indices from streamflow.
% Climate = Matrix of climate indices from precipitation.
% PM = Monthly precipitation (mm) per month number [Jan=1, Dec=12].
% QM = Monthly Mean Daily flow (l/s) per month number [Jan=1, Dec=12].
% IDCi = Maximum Intensity - Duration Curve [mm/h v time] for catchment i.
% FDCi = Flow Duration Curve [l/s v %] for catchment i.
% Comparative plots.
%
% Boris Ochoa Tocachi
% Imperial College London
% Created in June, 2014
% Last edited in November, 2017
%% PROCESS
% Initialise variables
Climate = zeros(13,2);
Indices = zeros(57,2);
PM = zeros(12,2);
QM = zeros(12,2);
% Calculate indices for Discharge and Precipitation.
[Climate(:,1),Indices(:,1)] = iMHEA_IndicesTotal(Date1,P1,Q1,A1);
[Climate(:,2),Indices(:,2)] = iMHEA_IndicesTotal(Date2,P2,Q2,A2);
[~,PM(:,1),IDC1,~,QM(:,1),FDC1,~,~,~,~,CumP1,CumQ1,DP1,DQ1] = iMHEA_Indices(Date1,P1,Q1,A1);
[~,PM(:,2),IDC2,~,QM(:,2),FDC2,~,~,~,~,CumP2,CumQ2,DP2,DQ2] = iMHEA_Indices(Date2,P2,Q2,A2);
%% PLOT PAIR COMPARISON
figure
subplot(3,1,3)
plot(Date1,Q1/A1,'k',Date2,Q2/A2,'r')
xlabel('Date')
ylabel('Discharge (l/s/km2)')
title('Input Discharge')
Xlim1 = get(gca,'XLim');
legend('Streamflow 1','Streamflow 2',...
'Location','NorthWest')
box on
subplot(3,1,1)
hold on
bar(Date1,P1,'k','EdgeColor','k')
xlabel('Date')
ylabel('Precipitation (mm)')
title('Input Precipitation')
Xlim2 = get(gca,'XLim');
legend('Rainfall 1',...
'Location','SouthWest')
box on
subplot(3,1,2)
hold on
bar(Date2,P2,'r','EdgeColor','r')
xlabel('Date')
ylabel('Precipitation (mm)')
title('Input Precipitation')
Xlim3 = get(gca,'XLim');
legend('Rainfall 2',...
'Location','SouthWest')
box on
Xlim = [min([Xlim1(1),Xlim2(1),Xlim3(1)]),max([Xlim1(2),Xlim2(2),Xlim3(2)])];
subplot(3,1,1)
set(gca,'YDir','reverse','XLim',Xlim)
subplot(3,1,2)
set(gca,'YDir','reverse','XLim',Xlim)
subplot(3,1,3)
set(gca,'XLim',Xlim)
figure
subplot(2,1,1)
hold on
plot(CumP1(:,1),CumP1(:,2),'-k',CumP2(:,1),CumP2(:,2),'-r')
plot(CumQ1(:,1),CumQ1(:,2),'--k',CumQ2(:,1),CumQ2(:,2),'--r','LineWidth',1.5)
title('Cumulative comparison')
legend('Rainfall 1','Rainfall 2','Streamflow 1','Streamflow 2',...
'Location','NorthWest')
xlabel('Date')
ylabel('Cumulative variables')
box on
subplot(2,1,2)
hold on
plot(CumP1(:,2),CumQ1(:,2),'-k')
plot(CumP2(:,2),CumQ2(:,2),'--r')
title('Double Mass Plot')
xlabel('Precipitation')
ylabel('Discharge')
legend('Catchment 1','Catchment 2',...
'Location','NorthWest')
box on
figure
semilogx(IDC1(:,1),IDC1(:,2),'-k')
hold on
semilogx(IDC2(:,1),IDC2(:,2),'--r')
xlabel('Duration (min)')
ylabel('Maximum precipitation intensity (mm/h)')
title('Maximum Intensity-Duration Curve')
legend('Catchment 1','Catchment 2')
box on
figure
semilogy(FDC1(:,1),FDC1(:,2),'-k')
hold on
semilogy(FDC2(:,1),FDC2(:,2),'--r')
xlabel('Exceedance probability')
ylabel('Discharge (l/s/km2)')
title('Flow Duration Curve')
legend('Catchment 1','Catchment 2')
box on
figure
subplot(2,1,1)
plot((0:12)',[PM(12,1);PM(:,1)],'-k',(0:12)',[QM(12,1);QM(:,1)],'--k')
xlabel('Date')
ylabel('Variable (mm)')
legend('Precipitation (mm)','Discharge (mm)')
title('Monthly Data')
title('Catchment 1')
box on
subplot(2,1,2)
plot((0:12)',[PM(12,2);PM(:,2)],'-r',(0:12)',[QM(12,2);QM(:,2)],'--r')
xlabel('Date')
ylabel('Variable (mm)')
legend('Precipitation (mm)','Discharge (mm)')
title('Monthly Data')
title('Catchment 2')
box on
dateP1 = datetime(DP1(:,1),'ConvertFrom','datenum');
dateP2 = datetime(DP2(:,1),'ConvertFrom','datenum');
dateQ1 = datetime(DQ1(:,1),'ConvertFrom','datenum');
dateQ2 = datetime(DQ2(:,1),'ConvertFrom','datenum');
figure
subplot(4,1,1)
bar(dateP1,DP1(:,2),'k','EdgeColor','k')
xlabel('Date')
ylabel('Precipitation (mm)')
set(gca,'YDir','reverse','XLim',Xlim);
title('Daily Precipitation in Catchment 1')
legend('Rainfall',...
'Location','South')
box on
subplot(4,1,2)
hold on
plot(dateQ1,DQ1(:,2),dateQ1,DQ1(:,3),dateQ1,DQ1(:,4))
xlabel('Date')
ylabel('Discharge (l/s/km2)')
set(gca,'XLim',Xlim)
legend('Total flow','Baseflow','Stormflow',...
'Location','North')
set(gca,'XLim',Xlim);
title('Daily Discharge in Catchment 1')
box on
subplot(4,1,3)
bar(dateP2,DP2(:,2),'r','EdgeColor','r')
xlabel('Date')
ylabel('Precipitation (mm)')
set(gca,'YDir','reverse','XLim',Xlim);
title('Daily Precipitation in Catchment 2')
legend('Rainfall',...
'Location','South')
box on
subplot(4,1,4)
hold on
plot(dateQ2,DQ2(:,2),dateQ2,DQ2(:,3),dateQ2,DQ2(:,4))
xlabel('Date')
ylabel('Discharge (l/s/km2)')
set(gca,'XLim',Xlim)
legend('Total flow','Baseflow','Stormflow',...
'Location','North')
set(gca,'XLim',Xlim);
title('Daily Discharge in Catchment 2')
box on