forked from mikedurand/SWOTAprimeCalcs
-
Notifications
You must be signed in to change notification settings - Fork 0
/
TestAprimeCalcsPepsi.m
175 lines (140 loc) · 5.15 KB
/
TestAprimeCalcsPepsi.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
%% Computing cross-sectional area changes from height and width observations
% Mike Durand, [email protected]
%
% May 11, 2018
clear all; close all
addpath('./PepsiSrc');
%% 3. Setup: Read data and create synthetic observations
% So far, I've looked at example cases on the Sacramento Upstream (cf=13)
% Garonne Downstream (cf=4), and Seine (cf=14) from the Pepsi 1.
%%
cf=10;
pathtoncfiles='./Pepsi1/';
Files=dir([pathtoncfiles '*.nc']);
numbfiles=size(Files,1);
ReadRiver;
FilterRiver;
R=Rivers(cf).Reaches;
disp(['Working with ' Rivers(cf).Name ' data.'])
tObs=ceil(Rivers(cf).orbit./86400)';
nObsCycle=length(tObs);
clear Rivers
[nR,nt]=size(R.W);
t=1:nt;
% SWOT observation
nCycle=floor(nt/21);
nReg=3; %number of regions
% for i=2:nCycle
% tObs=[tObs tObs(1:nObsCycle)+(i-1)*21];
% end
tObs=t;
nObs=length(tObs);
ReDoFits=true;
if ReDoFits
stdH=0.5;
stdW=10;
Hobs=R.H(:,tObs)+stdH.*randn(nR,nObs);
Wobs=R.W(:,tObs)+stdW.*randn(nR,nObs);
save('HWData.mat','Hobs','Wobs','stdH','stdW')
else
load HWData
end
%% 4. Visualize height-width relationship & summary statistics
Hbar=median(Hobs,2);
Wbar=median(Wobs,2);
m_zz=nan(2,2,nR); %m_zz is the W-H observation covariance matrix for each reach
for i=1:nR
m_zz(:,:,i)=cov(Wobs(i,:),Hobs(i,:));
end
r=4;
disp(['Working with reach ' num2str(r)])
figure(1)
plot(R.H(r,tObs),R.W(r,tObs),'o','LineWidth',2);
set(gca,'FontSize',14)
xlabel('Water surface elevation, m')
ylabel('River width, m')
title(['Fig. 1: "True" Height-width data for reach ' num2str(r)])
grid on
%% 5. Fitting funtions to the width & height data
[p,xbreak,iout,hval,wval,UseReg] = FitData(ReDoFits,nR,Wobs,stdW,Hobs,stdH,nReg);
figure(2)
errorbar(Hobs(r,:),Wobs(r,:),stdW*ones(nObs,1),stdW*ones(nObs,1),stdH*ones(nObs,1),stdH*ones(nObs,1),'LineStyle','none','LineWidth',2,'Marker','o'); hold on
plot(hval{r},wval{r},'LineWidth',2)
a=axis;
for i=1:length(xbreak{r})
plot(xbreak{r}(i)*ones(2,1),a(3:4),'g-')
end
plot(Hobs(r,iout{r}),Wobs(r,iout{r}),'kx','MarkerSize',12,'LineWidth',2)
hold off;
set(gca,'FontSize',14)
xlabel('Water surface elevation, m')
ylabel('River width, m')
title(['Fig. 2: Height-width fits & obs. for reach ' num2str(r)])
grid on
%% 6. Performing the Aprime calculations
dAHbar=CalculatedAEIV(Hbar(r),Wbar(r),xbreak{r},p{r},nReg,0,stdW^2,stdH^2,m_zz(:,:,r),nObs); %sample calculation
rMax=find(UseReg{r}, 1, 'last' );
Htest=linspace(xbreak{r}(1)+.001,xbreak{r}(rMax+1)-.001,100);
for i=1:length(Htest)
[dAtest(i)] = CalculatedAEIV(Htest(i),[],xbreak{r},p{r},nReg,dAHbar,stdW^2,stdH^2,m_zz(:,:,r),nObs); %sample calculation
end
[~,iSort]=sort(R.H(r,:));
[Hs,iUnique]=unique(R.H(r,:),'sorted');
As=R.A(r,iUnique);
AHbarhat_true=interp1(Hs,As,Hbar(r));
epsilonAbar=AHbarhat_true-median(R.A(r,:));
disp(['epsilon_Abar=' num2str(epsilonAbar) ' m^2'])
disp(['Abar=' num2str(median(R.A(r,:))) ' m^2'])
dAtrue=R.A(r,:)-median(R.A(r,:)); %true dA for comparison
figure(3)
plot(R.H(r,tObs),dAtrue(tObs),'o',Htest,dAtest,Hbar(r),0,'ko','LineWidth',2)
set(gca,'FontSize',14)
xlabel('Height, m')
ylabel('A'', m^2')
legend('True','SWOT calibrated','Hbar,0','Location','Best')
grid on;
title(['Fig. 3: Illustration of A''(H) for reach ' num2str(r)])
%% 7. Timeseries estimation of $A\prime$ and preliminary uncertainty calculations
for i=1:nObs
[dAhat(i),What(i),Hhat(i),dAunc(i)] = CalculatedAEIV(Hobs(r,i),Wobs(r,i),xbreak{r},p{r},nReg,dAHbar,stdW^2,stdH^2,m_zz(:,:,r),nObs);
end
figure(4)
plot(tObs,dAtrue(tObs),'s-','LineWidth',2); hold on;
errorbar(tObs,dAhat,dAunc,'s-','LineStyle','none','LineWidth',2,'MarkerSize',10); hold off;
set(gca,'FontSize',14)
xlabel('Obs Time')
ylabel('A'', m^2')
title(['Fig. 4: A'' timeseries for reach ' num2str(r)])
legend('True','SWOT')
Err.AprimeErrStd=std(dAhat-dAtrue(tObs));
Err.AprimeAvgUnc=mean(dAunc);
%% 8. Taking a quick look at the constrained widths and heights
figure(5)
plot([min(Wobs(r,~iout{r})) max(Wobs(r,~iout{r}))],[min(Wobs(r,~iout{r})) max(Wobs(r,~iout{r}))],'k-');
hold on;
plot(R.W(r,tObs(~iout{r})),Wobs(r,~iout{r}),'ro','LineWidth',2,'MarkerSize',10);
plot(R.W(r,tObs(~iout{r})),What(~iout{r}),'bx','LineWidth',2,'MarkerSize',10);
hold off
set(gca,'FontSize',14)
xlabel('True Width, m')
ylabel('SWOT Width, m')
legend('1:1','Observations','Constrained Obs','Location','Best')
title(['Fig. 5: Observed and constrained width for reach ' num2str(r)])
grid on;
Err.WidthObsStdDev=std(R.W(r,tObs(~iout{r}))-Wobs(r,~iout{r})); %close to stdW
Err.WidthConstrainedStdDev=std(R.W(r,tObs(~iout{r}))-What(~iout{r}));
Err.HeightObsStdDev=std(R.H(r,tObs(~iout{r}))-Hobs(r,~iout{r})); %close to stdW
Err.HeightConstrainedStdDev=std(R.H(r,tObs(~iout{r}))-Hhat(~iout{r}));
disp(Err)
figure(6)
plot([min(Hobs(r,~iout{r})) max(Hobs(r,~iout{r}))],[min(Hobs(r,~iout{r})) max(Hobs(r,~iout{r}))],'k-');
hold on;
plot(R.H(r,tObs(~iout{r})),Hobs(r,~iout{r}),'ro','LineWidth',2,'MarkerSize',10);
plot(R.H(r,tObs(~iout{r})),Hhat(~iout{r}),'bx','LineWidth',2,'MarkerSize',10);
hold off
set(gca,'FontSize',14)
xlabel('True Height, m')
ylabel('SWOT Height, m')
legend('1:1','Observations','Constrained Obs','Location','Best')
title(['Fig. 6: Observed and constrained width for reach ' num2str(r)])
grid on;