1 | close all |
---|
2 | clear all |
---|
3 | indir='../exp_nemed/'; |
---|
4 | %indir='../exp_twin_aviso/'; |
---|
5 | %indir='../../code_leila/altifloat/obs_float/'; |
---|
6 | %indir='../obs_float/'; |
---|
7 | %ut=load([ indir 'uzero.dat']); |
---|
8 | %vt=load([ indir 'vzero.dat']); |
---|
9 | robs=load([ indir 'obs.dat']); |
---|
10 | rsim=load([ indir 'rfloat_nodiv.dat']); |
---|
11 | uvb=load([ indir 'uv_back_aviso.dat']); |
---|
12 | uvt=load([ indir 'uv_back_aviso.dat']); |
---|
13 | uvr=load([indir 'uv_tot_nodiv.dat']); |
---|
14 | uvgeo=load([indir 'uv_geos_nodiv.dat']); |
---|
15 | meshg=load([indir 'meshgrid_aviso.dat']); |
---|
16 | Div_ana_file=[indir 'Div_ana.dat']; |
---|
17 | Div_fg_file=[indir 'Div_fg.dat']; |
---|
18 | |
---|
19 | %drifterdir='../altifloat/'; |
---|
20 | %fname={'a300234011043340.mat','a300234011120770.mat','a300234011249440.mat'}; |
---|
21 | %X1=ut(:,end-3); |
---|
22 | %Y1=ut(:,end-2); |
---|
23 | iech=1; |
---|
24 | jech=1; |
---|
25 | %LIM=[34.5 36 32.5 35]; |
---|
26 | LON = [32.5 35.6]; |
---|
27 | LAT = [32.5 35.5]; |
---|
28 | LIM= [LON LAT]; |
---|
29 | %Format : t/s, ilon, ilat, val |
---|
30 | |
---|
31 | %nlon=87; |
---|
32 | %nlat=58; |
---|
33 | lon=unique(meshg(:,1)); |
---|
34 | lat=unique(meshg(:,2)); |
---|
35 | nlon=length(lon); |
---|
36 | nlat=length(lat); |
---|
37 | |
---|
38 | [Lon,Lat]=meshgrid(lon',lat); |
---|
39 | |
---|
40 | %X1=X1+2; |
---|
41 | %Y1=Y1+2 |
---|
42 | |
---|
43 | %Ut=reshape(ut(:,end),nlat,nlon); |
---|
44 | %Vt=reshape(vt(:,end),nlat,nlon); |
---|
45 | |
---|
46 | |
---|
47 | |
---|
48 | nfloat=length(unique(robs(:,2))); |
---|
49 | |
---|
50 | |
---|
51 | |
---|
52 | %it, ilon,ilat,u,v |
---|
53 | ntime=length(unique(uvb(:,1))); |
---|
54 | Ub=reshape(uvb(:,end-1),nlat,nlon,ntime); |
---|
55 | Vb=reshape(uvb(:,end),nlat,nlon,ntime); |
---|
56 | |
---|
57 | ntime=length(unique(uvt(:,1))); |
---|
58 | Ut=reshape(uvt(:,end-1),nlat,nlon,ntime); |
---|
59 | Vt=reshape(uvt(:,end),nlat,nlon,ntime); |
---|
60 | |
---|
61 | |
---|
62 | %it, ilon,ilat,u,v |
---|
63 | ntime=length(unique(uvr(:,1))); |
---|
64 | Ur=reshape(uvr(:,end-1),nlat,nlon,ntime); |
---|
65 | Vr=reshape(uvr(:,end),nlat,nlon,ntime); |
---|
66 | |
---|
67 | if exist('uvgeo') |
---|
68 | ntime=length(unique(uvgeo(:,1))); |
---|
69 | Ugeo=reshape(uvgeo(:,end-1),nlat,nlon,ntime); |
---|
70 | Vgeo=reshape(uvgeo(:,end),nlat,nlon,ntime); |
---|
71 | end |
---|
72 | |
---|
73 | |
---|
74 | Hr=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ur(1:jech:end,1:iech:end,1:end-1),3), ... |
---|
75 | mean(Vr(1:jech:end,1:iech:end,1:end-1),3),'r'); |
---|
76 | hold on |
---|
77 | Hb=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ub(1:jech:end,1:iech:end,1:end-1),3),mean(Vb(1:jech:end,1:iech:end,1:end-1),3),'b'); |
---|
78 | |
---|
79 | %Ht=quiver(Ut,Vt,'g'); |
---|
80 | |
---|
81 | |
---|
82 | for j=1:nfloat |
---|
83 | |
---|
84 | |
---|
85 | |
---|
86 | rfloat=rsim(rsim(:,2)==j-1&rsim(:,end)>0,[end-1 end]); |
---|
87 | rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); |
---|
88 | rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); |
---|
89 | Hsim=plot(rfloat_lon,rfloat_lat,'.-r'); |
---|
90 | rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]); |
---|
91 | rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); |
---|
92 | rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); |
---|
93 | Hobs=plot(rfloat_lon,rfloat_lat,'.-m'); |
---|
94 | |
---|
95 | % plot(rfloat(:,2)+1,rfloat(:,1)+1,'.-m'); |
---|
96 | |
---|
97 | % Hreal=plot(data.lon(1:4),data.lat(1:4),'.-k') |
---|
98 | |
---|
99 | end |
---|
100 | |
---|
101 | |
---|
102 | Ht=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ut(1:jech:end,1:iech:end,1:end-1),3),mean(Vt(1:jech:end,1:iech:end,1:end-1),3),'m'); |
---|
103 | axis(LIM); |
---|
104 | |
---|
105 | |
---|
106 | |
---|
107 | %Ajout des trajectoires réelles |
---|
108 | |
---|
109 | |
---|
110 | |
---|
111 | %legend([Hb Hr],'background','retrieved','Location','SouthEast'); |
---|
112 | legend([Hb Hr Ht],'background.','retrieved','truth','Location','SouthEast'); |
---|
113 | |
---|
114 | %print -dpng ../exp_forw_dan/simu_drifter_dan.png |
---|
115 | |
---|
116 | |
---|
117 | if (false) |
---|
118 | delta_b=sqrt(mean((Ut(:,:,1:end-1)-Ub(:,:,1:end-1)),3).^2+mean((Vt(:,:,1:end-1)-Vb(:,:,1:end-1)),3).^2); |
---|
119 | delta_r=sqrt(mean((Ut(:,:,1:end-1)-Ur(:,:,1:end-1)),3).^2+mean((Vt(:,:,1:end-1)-Vr(:,:,1:end-1)),3).^2); |
---|
120 | delta_inc=sqrt(mean((Ub(:,:,1:end-1)-Ur(:,:,1:end-1)),3).^2+mean((Vb(:,:,1:end-1)-Vr(:,:,1:end-1)),3).^2); |
---|
121 | |
---|
122 | cscale=[0 5e-4]; |
---|
123 | |
---|
124 | figure(2) |
---|
125 | subplot(2,2,1) |
---|
126 | imagesc(lon,lat,delta_b,cscale); |
---|
127 | axis xy;axis(LIM); |
---|
128 | hold on |
---|
129 | for j=1:3 |
---|
130 | data= load([drifterdir fname{j}]); |
---|
131 | rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]); |
---|
132 | rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); |
---|
133 | rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); |
---|
134 | plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m'); |
---|
135 | end |
---|
136 | title('Backround error'); |
---|
137 | |
---|
138 | |
---|
139 | subplot(2,2,2) |
---|
140 | imagesc(lon,lat,delta_r,cscale); |
---|
141 | axis xy;axis(LIM); |
---|
142 | hold on |
---|
143 | for j=1:3 |
---|
144 | data= load([drifterdir fname{j}]); |
---|
145 | rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]); |
---|
146 | rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); |
---|
147 | rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); |
---|
148 | plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m'); |
---|
149 | end |
---|
150 | |
---|
151 | title('Analysis error'); |
---|
152 | |
---|
153 | |
---|
154 | |
---|
155 | subplot(2,2,4) |
---|
156 | imagesc(lon,lat,delta_inc); |
---|
157 | axis xy;axis(LIM); |
---|
158 | colorbar |
---|
159 | hold on |
---|
160 | for j=1:3 |
---|
161 | data= load([drifterdir fname{j}]); |
---|
162 | rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]); |
---|
163 | rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); |
---|
164 | rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); |
---|
165 | plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m'); |
---|
166 | end |
---|
167 | title('Increment'); |
---|
168 | |
---|
169 | |
---|
170 | %print -f1 -dpng ../exp_forw_dan/twin_exp_uv_field.png |
---|
171 | %print -f2 -dpng ../exp_forw_dan/twin_exp_errors.png |
---|
172 | end |
---|
173 | |
---|
174 | |
---|
175 | %% Erreur le long du trajet |
---|
176 | if(false) |
---|
177 | Id=0; |
---|
178 | for id=Id |
---|
179 | figure(5+id) |
---|
180 | I=find(robs(:,2)==id-1); |
---|
181 | Xq=robs(I,end)+1; |
---|
182 | Yq=robs(I,end-1)+1; |
---|
183 | Tq=robs(I,1)+1; |
---|
184 | |
---|
185 | err_b=interp3(1:nlon,1:nlat,1:ntime,delta_b,Xq,Yq,Tq); |
---|
186 | err_r=interp3(1:nlon,1:nlat,1:ntime,delta_r,Xq,Yq,Tq); |
---|
187 | incr=interp3(1:nlon,1:nlat,1:ntime,delta_inc,Xq,Yq,Tq); |
---|
188 | end %id |
---|
189 | end %if false |
---|
190 | |
---|
191 | %% Divergence |
---|
192 | Div_ana=load(Div_ana_file); |
---|
193 | Div_fg=load(Div_fg_file); |
---|
194 | |
---|
195 | Div2_ana=reshape(Div_ana(:,end),nlat,nlon,ntime-1); |
---|
196 | Div2_fg=reshape(Div_fg(:,end),nlat,nlon,ntime-1); |
---|
197 | |
---|
198 | figure(3) |
---|
199 | subplot(3,1,1) |
---|
200 | imagesc(lon,lat,mean(Div2_fg,3),[-2 2]*1e-5); |
---|
201 | axis xy; |
---|
202 | axis(LIM); |
---|
203 | title('Mean divergence Background'); |
---|
204 | colorbar |
---|
205 | |
---|
206 | |
---|
207 | figure(3) |
---|
208 | subplot(3,1,3) |
---|
209 | imagesc(lon,lat,mean(Div2_ana,3),[-2 2]*1e-5); |
---|
210 | axis xy; |
---|
211 | axis(LIM); |
---|
212 | title('Mean divergence Analyse'); |
---|
213 | colorbar |
---|
214 | |
---|
215 | subplot(3,1,2); |
---|
216 | imagesc(lon,lat,mean(Div2_ana,3)-mean(Div2_fg,3),[-2 2]*1e-5); |
---|
217 | axis xy; |
---|
218 | axis(LIM); |
---|
219 | title('Mean difference'); |
---|
220 | colorbar |
---|