Changeset 139 for altifloat

Timestamp:

03/27/15 15:19:37 (9 years ago)

Author:

jbrlod

Message:

add functionalities to mask the model experiment

Location:

altifloat/matlab_toolbox

Files:

• aviso2yao.m (modified) (5 diffs)
• interp_bck.m (modified) (2 diffs)
• script_plot.m (added)
• script_prepare_exp_modele.m (added)

altifloat/matlab_toolbox/aviso2yao.m

@@ -1 @@
-function [] = aviso2yao( aviso_netcdf, meshmask,outfile,outmask)
-% [] = aviso2yao( aviso_netcdf_file, meshmask,outfile,outmask)
+function [Yao_uv,mask1] = aviso2yao( aviso_netcdf, meshmask,outfile,outmask,convlon)
+% [Yao_uv,mask1] = aviso2yao( aviso_netcdf, meshmask,outfile,outmask)
 % Read AVISO file aviso_netcdf_file and save to be read by read_bck (time step=1)
 % Interpolate on meshmask (meshmask : two columns lon,lat)
@@ -11 +11 @@
   
 
-%%%%%%%%%% Usually use:
-%%%%%%%%%% Nlat=58; %Phi
-%%%%%%%%%% Nlong=87; %Lambda
-  %% aviso_netcdf='../../data/AVISO/dt_med_allsat_madt_uv_20130828_20140704.nc'
-  %% meshmask='../obs_float/meshgrid3.dat'
-  %% outfile='../obs_float/uv_bck_20130828.dat';
-  %% outmask='../obs_float/mask_20130828.dat';
+% Usually use:
+% Nlat=58; %Phi
+% Nlong=87; %Lambda
+% aviso_netcdf='../../data/AVISO/dt_med_allsat_madt_uv_20130828_20140704.nc'
+
+% meshmask='../obs_float/meshgrid3.dat'
+% outfile='../obs_float/uv_bck_20130828.dat';
+% outmask='../obs_float/mask_20130828.dat';
 
 write_u=1;
+if ~exist('convlon')
+    convlon=true;
+end
 
  
@@ -58 +62 @@
 [LAT,LON]=meshgrid(lat,lon); %C'est à l'envers à cause de matlab
 
+if convlon
 %Conversion des longitudes
 LON=LON-360;
+end
 
 %Conversion des vitesses inutile (malgré le champ scale_factor, ce que je trouve bizarre,
@@ -73 +79 @@
 
 mask1=mask(:);
+
+if (exist('outmask') & ~isempty(outmask))
 dlmwrite(outmask,mask1,'precision',1);
+end
 
-%%Masked velocity are zero
+%Masked velocity are zero
 Uint(mask==0)=0;
 Vint(mask==0)=0;
@@ -136 +145 @@
 Yao_uv=[ ones(Nlong*Nlat,1)  cd_lon(:) cd_lat(:)  u_yao1 v_yao1];
 
-if (write_u==1)
+if (write_u==1 & exist('outfile') & ~isempty(outfile))
 
 dlmwrite(outfile,Yao_uv,'\t');

altifloat/matlab_toolbox/interp_bck.m

@@ -1 @@
-function []=interp_bck(uv_b1,uv_b2,it1,it2,dt,outfile)
+function [UV_bck,mask_back]=interp_bck(uv,it,dt,meshmask,outfile,outmask)
  %[]=interp_bck(uv_b1,uv_b2,it1,it2,dt)
  %Time unit : hour
+ % Cells containing file names in nc format
+ % it : time step of each file
    
-   %Check if same coordinate
-     
-     %uv_b1='../obs_float/uv_bck_20130828.dat'
-     %uv_b2='../obs_float/uv_bck_20130828.dat'
-     %it1=1;
-     %it2=47;
-     %dt=2;
-     %outfile='../obs_float/uv_bk_24pdt_const.dat'
-       
-     uv1=load(uv_b1);
-     uv2=load(uv_b2);
-     delta=(uv1(:,2:3)-uv2(:,2:3)).^2;
+%Check if same coordinate
+
+%uv_b1='../obs_float/uv_bck_20130828.dat'
+%uv_b2='../obs_float/uv_bck_20130828.dat'
+%it1=1;
+%it2=47;
+%dt=2;
+%outfile='../obs_float/uv_bk_24pdt_const.dat'
+
+%The grid in degrees 
+meshgrid2=load(meshmask);
+Nlong=length(unique(meshgrid2(:,1)));
+Nlat=length(unique(meshgrid2(:,2)));
+
+UV_bck=[];
+mask_back=ones(Nlat*Nlong,1);
+
+
+for idf=1:length(uv)-1
+    
+    uv_b1=uv{idf};
+    uv_b2=uv{idf+1};
+    it1=it(idf);
+    it2=it(idf+1);
+    
+    %     uv1=load(uv_b1);
+    %     uv2=load(uv_b2);
+    [uv1,mask1]=aviso2yao(uv_b1,meshmask,[],[],false);
+    [uv2,mask2]=aviso2yao(uv_b1,meshmask,[],[],false);
+    
+    mask_back=mask_back.*mask1.*mask2;
+    
+    delta=(uv1(:,2:3)-uv2(:,2:3)).^2;
      if (sum(delta(:))>1e-8)
        fprintf(2,'Coordinates does not match\n');
@@ -31 +54 @@
 
        end
-       
        tt=repmat(1:length(t),size(uv1,1),1);
        tt=tt(:);
        
-       UV_bck=[tt,repmat(uv1(:,2:3),length(t),1),ui(:),vi(:)];
+       UV_bck=[UV_bck;[tt,repmat(uv1(:,2:3),length(t),1),ui(:),vi(:)]];
+
+end %for idf
        
        dlmwrite(outfile,UV_bck,'\t');
-     
+       dlmwrite(outmask,mask_back,'precision',1)
      end