Changeset 145 for altifloat

Timestamp:

04/06/15 19:15:40 (9 years ago)

Author:

jbrlod

Message:

prepare model dan

Location:

altifloat/matlab_toolbox

Files:

• script_perso.m (modified) (5 diffs)
• script_prepare_exp_modele.m (modified) (2 diffs)

altifloat/matlab_toolbox/script_perso.m

@@ -6 +6 @@
 %ut=load([ indir 'uzero.dat']);
 %vt=load([ indir 'vzero.dat']);
-robs=load([ indir 'rfloat_total.dat']);
-rsim=load([ indir 'rfloat_total_jum0.dat']);
-uvb=load([ indir 'uv_bck_dan.dat']);
-uv=load([indir 'uv_total_jum0.dat']);
+robs=load([ indir 'obs_jum.dat']);
+rsim=load([ indir 'rfloat_total_jum1.dat']);
+uvb=load([ indir 'uv_bck_dan_s.dat']);
+uvt=load([ indir 'uv_bck_dan.dat']);
+uvr=load([indir 'uv_total_jum1.dat']);
 meshg=load([indir 'meshgrid_dan.dat']);
 drifterdir='../altifloat/';
@@ -17 +18 @@
 iech=3;
 jech=3;
-LIM=[34.9 35.3 33.25 33.6];
+LIM=[34.8 35.3 33.2 33.6];
 %Format : t/s, ilon, ilat, val
 
@@ -48 +49 @@
 
 %it, ilon,ilat,u,v
-ntime=length(unique(uv(:,1)));
-Ur=reshape(uv(:,end-1),nlat,nlon,ntime);
-Vr=reshape(uv(:,end),nlat,nlon,ntime);
+ntime=length(unique(uvr(:,1)));
+Ur=reshape(uvr(:,end-1),nlat,nlon,ntime);
+Vr=reshape(uvr(:,end),nlat,nlon,ntime);
 
-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),3),mean(Vb(1:jech:end,1:iech:end,1),3),'b');
+ntime=length(unique(uvt(:,1)));
+Ut=reshape(uvt(:,end-1),nlat,nlon,ntime);
+Vt=reshape(uvt(:,end),nlat,nlon,ntime);
+
+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), ...
+          mean(Vr(1:jech:end,1:iech:end,1:end-1),3),'r');
+hold on
+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');
 
 %Ht=quiver(Ut,Vt,'g');
-hold on
+
 
 for j=1:nfloat     
   data= load([drifterdir fname{j}]);
   
-  rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]);
-  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
-  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
-  Hobs=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-b');
+
   
   rfloat=rsim(rsim(:,2)==j-1&rsim(:,end)>0,[end-1 end]);
@@ -69 +74 @@
   rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
   Hsim=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-r');
-  
+    rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]);
+  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
+  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
+  Hobs=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-m');
 
   %   plot(rfloat(:,2)+1,rfloat(:,1)+1,'.-m');
 
-  Hreal=plot(data.lon(1:4),data.lat(1:4),'.-k')
+  % Hreal=plot(data.lon(1:4),data.lat(1:4),'.-k')
   
 end
 
-Hb=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ur(1:jech:end,1:iech:end,1),3),mean(Vr(1:jech:end,1:iech:end,1),3),'b');
+
+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');
 axis(LIM);
 
@@ -86 +95 @@
 
 %legend([Hb Hr],'background','retrieved','Location','SouthEast');
-legend([Hobs Hsim],'real traj.','simulated traj');
+legend([Hb Hr Ht],'background.','retrieved','truth');
 
-print -dpng ../exp_forw_dan/simu_drifter_dan.png
+%print -dpng ../exp_forw_dan/simu_drifter_dan.png
 
+
+
+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);
+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);
+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);
+
+cscale=[0 5e-4];
+
+figure(2)
+subplot(2,2,1)
+imagesc(lon,lat,delta_b,cscale);
+axis xy;axis(LIM);
+hold on
+for j=1:3
+   data= load([drifterdir fname{j}]);
+  rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]);
+  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
+  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
+plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m');
+end
+title('Backround error');
+
+
+subplot(2,2,2)
+imagesc(lon,lat,delta_r,cscale);
+axis xy;axis(LIM);
+hold on
+for j=1:3
+   data= load([drifterdir fname{j}]);
+  rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]);
+  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
+  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
+plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m');
+end
+
+title('Analysis error');
+
+
+
+subplot(2,2,4)
+imagesc(lon,lat,delta_inc);
+axis xy;axis(LIM);
+colorbar
+hold on
+for j=1:3
+   data= load([drifterdir fname{j}]);
+  rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]);
+  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
+  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
+plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m');
+end
+title('Increment');
+
+
+%print -f1 -dpng ../exp_forw_dan/twin_exp_uv_field.png
+%print -f2 -dpng ../exp_forw_dan/twin_exp_errors.png
+
+%% Erreur le long du trajet
+if(false)
+Id=0;
+for id=Id
+  figure(5+id)
+  I=find(robs(:,2)==id-1);
+  Xq=robs(I,end)+1;
+  Yq=robs(I,end-1)+1;
+  Tq=robs(I,1)+1;
+  
+  err_b=interp3(1:nlon,1:nlat,1:ntime,delta_b,Xq,Yq,Tq);
+  err_r=interp3(1:nlon,1:nlat,1:ntime,delta_r,Xq,Yq,Tq);
+  incr=interp3(1:nlon,1:nlat,1:ntime,delta_inc,Xq,Yq,Tq);
+end %id
+end %if false

altifloat/matlab_toolbox/script_prepare_exp_modele.m

@@ -21 +21 @@
 interp_bck(aviso_netcdf,it,dt,meshfile,outfile,outmask);
 
+%% Lissage
+outfile_s='../exp_forw_dan/uv_bck_dan_s.dat';
+
+
+uvb=load(outfile);
+meshg=load(meshfile);
+lon=unique(meshg(:,1));
+lat=unique(meshg(:,2));
+nlon=length(lon);
+nlat=length(lat);
+ntime=length(unique(uvb(:,1)));
+Ub=reshape(uvb(:,end-1),nlat,nlon,ntime);
+Vb=reshape(uvb(:,end),nlat,nlon,ntime);
+
+Ub_s=nan*ones(size(Ub));
+Vb_s=nan*ones(size(Vb));
+
+K=ones(95,95);
+K=K./sum(K(:));
+
+for it=1:ntime
+Ub_s(:,:,it)=conv2(Ub(:,:,it),K,'same');
+Vb_s(:,:,it)=conv2(Vb(:,:,it),K,'same');
+end
+
+uvb_s=uvb;
+uvb_s(:,end-1) = Ub_s(:);
+uvb_s(:,end)   = Vb_s(:);
+dlmwrite(outfile_s,uvb_s,'\t');
+
 %% STEP THREE : Prepare init
 drifterdir='../altifloat/';
@@ -40 +70 @@
 dlmwrite(['../exp_forw_dan/' initfile],MM,'\t');
 
-%STEP OBS : modify obs file
+%% STEP OBS : modify obs file
+init=load(['../exp_forw_dan/' initfile]);
+obs=load(['../exp_forw_dan/rfloat_total.dat']);
 
+ilon=find(init(:,1)==1);
+ilat=find(init(:,1)==2);
+
+init=[init(ilon,2)-1 init(ilon,3)-1 init(ilat,4) init(ilon,4)];
+
+istep=5;
+I=find(obs(:,end)>1e-8);
+obs=obs(I,:);
+
+obs=[init;obs];
+obs=sortrows(obs,[2 1]);
+
+t=unique(obs(:,1));
+tech=t(1:istep:end);
+
+I=find(ismember(obs(:,1),tech));
+obs=obs(I,:);
+dlmwrite(['../exp_forw_dan/obs_jum.dat'],obs,'\t');
+
+%% Create bck
+aviso_netcdf={[dir_uvfiles 'uveta_cyco1_2013092812.nc']...
+              [dir_uvfiles 'uveta_cyco1_2013092912.nc'];};
+
+it=1:24:1+24*(length(aviso_netcdf)-1);
+
+%meshmask='../exp_forw_dan/meshgrid_dan.dat';
+outfile='../exp_forw_dan/uv_bck2_dan.dat';
+outmask='../exp_forw_dan/mask_dan2.dat';
+
+dt=1;
+
+interp_bck(aviso_netcdf,it,dt,meshfile,outfile,outmask);
+