Changeset 943 for trunk/src/series/merge_proj_polar.m

Timestamp:

May 2, 2016, 9:49:22 AM (8 years ago)

Author:

sommeria

Message:

merge_proj_polar corrected

File:

• trunk/src/series/merge_proj_polar.m (modified) (10 diffs)

trunk/src/series/merge_proj_polar.m

@@ -1 @@
-%'merge_proj': concatene several fields from series, can project them on a regular grid in phys coordinates
+%'merge_proj': concatene several fields from series, project on a polar grid
 %------------------------------------------------------------------------
 % function ParamOut=merge_proj(Param)
@@ -104 +104 @@
 HeadData.azimuth=azimuth_arclength;    
 thresh2=16; % square of the interpolation range
+
 %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
 ParamOut=[]; %default output
@@ -160 +161 @@
 
 %% prepare output file content
-    TimeData.ListGlobalAttribute={'Conventions','Project','CoordUnit','TimeUnit','ZPos'};
-    TimeData.Conventions='uvmat';
-    TimeData.Project='2016_Circumpolar';
-    TimeData.CoordUnit='cm';
-    TimeData.TimeUnit='s';
-    TimeData.ZPos=0;
-    TimeData.ListVarName={'time','radius','azimuth','U','V','curl','div'};
-    TimeData.VarDimName={'time','radius','azimuth',{'time','radius','azimuth'},{'time','radius','azimuth'}...
-                        {'time','radius','azimuth'},{'time','radius','azimuth'}};
-                    TimeData.VarAttribute{1}.Role='';
-                    TimeData.VarAttribute{2}.Role='';
-                    TimeData.VarAttribute{3}.Role='';
-                    TimeData.VarAttribute{4}.Role='vector_x';
-                    TimeData.VarAttribute{5}.Role='vector_y';
-                    TimeData.VarAttribute{6}.Role='scalar';
-                    TimeData.VarAttribute{7}.Role='scalar';
-                    TimeData.time=nan(1,NbField);
-                    TimeData.radius=radius_shifted;
-                    TimeData.azimuth=azimuth_arclength;
-                    nby=numel(radius);
-                    nbx=numel(azimuth);
-                    TimeData.U=nan(NbField,nby,nbx);
-                    TimeData.V=nan(NbField,nby,nbx);
-                    TimeData.curl=nan(NbField,nby,nbx);
-                    TimeData.div=nan(NbField,nby,nbx);
+%     TimeData.ListGlobalAttribute={'Conventions','Project','CoordUnit','TimeUnit','ZPos'};
+%     TimeData.Conventions='uvmat';
+%     TimeData.Project='2016_Circumpolar';
+%     TimeData.CoordUnit='cm';
+%     TimeData.TimeUnit='s';
+%     TimeData.ZPos=0;
+%     TimeData.ListVarName={'time','radius','azimuth','U','V','curl','div'};
+%     TimeData.VarDimName={'time','radius','azimuth',{'time','radius','azimuth'},{'time','radius','azimuth'}...
+%                         {'time','radius','azimuth'},{'time','radius','azimuth'}};
+%                     TimeData.VarAttribute{1}.Role='';
+%                     TimeData.VarAttribute{2}.Role='';
+%                     TimeData.VarAttribute{3}.Role='';
+%                     TimeData.VarAttribute{4}.Role='vector_x';
+%                     TimeData.VarAttribute{5}.Role='vector_y';
+%                     TimeData.VarAttribute{6}.Role='scalar';
+%                     TimeData.VarAttribute{7}.Role='scalar';
+%                     TimeData.time=nan(1,NbField);
+%                     TimeData.radius=radius_shifted;
+%                     TimeData.azimuth=azimuth_arclength;
+%                     nby=numel(radius);
+%                     nbx=numel(azimuth);
+%                     TimeData.U=nan(NbField,nby,nbx);
+%                     TimeData.V=nan(NbField,nby,nbx);
+%                     TimeData.curl=nan(NbField,nby,nbx);
+%                     TimeData.div=nan(NbField,nby,nbx);
                     
 %     if ~isempty(timeread)
@@ -254 +255 @@
 for iview =2:numel(XmlData)
     if ~(isfield(XmlData{iview},'GeometryCalib')&& isequal(XmlData{iview}.GeometryCalib.SliceCoord,ProjObjectCoord))...
-        disp('error: geometric calibration missing')
+        disp('error: geometric calibration missing or inconsistent plane positions')
         return
     end
@@ -285 +286 @@
 
 % %% output file type
-% FileExtOut='.nc'; %netcdf output
-% if isempty(j1_series{1})
-%     NomTypeOut='_1';
-% else
-%     NomTypeOut='_1_1';
-% end
-% RootFileOut=RootFile{1};
-% for iview=2:NbView
-%     if ~strcmp(RootFile{iview},RootFile{1})
-%         RootFileOut='mproj';
-%         break
-%     end
-% end
+if isempty(j1_series{1})
+    NomTypeOut='_1';
+else
+    NomTypeOut='_1_1';
+end
+RootFileOut=RootFile{1};
+for iview=2:NbView
+    if ~strcmp(RootFile{iview},RootFile{1})
+        RootFileOut='mproj';
+        break
+    end
+end
 
 
@@ -313 +313 @@
 for index=1:NbField
     disp(['index=' num2str(index)])
-    disp(['ellapsed time ' num2str(toc(tstart)/60,2) ' minutes'])
+    disp(['ellapsed time ' num2str(toc(tstart)/60,4) ' minutes'])
     update_waitbar(WaitbarHandle,index/NbField)
     if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
@@ -319 +319 @@
         return
     end
-   
-%     %% generating the name of the merged field
-%     i1=i1_series{1}(index);
-%     if ~isempty(i2_series{end})
-%         i2=i2_series{end}(index);
-%     else
-%         i2=i1;
-%     end
-%     j1=1;
-%     j2=1;
-%     if ~isempty(j1_series{1})
-%         j1=j1_series{1}(index);
-%         if ~isempty(j2_series{end})
-%             j2=j2_series{end}(index);
-%         else
-%             j2=j1;
-%         end
-%     end
-%    % OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomTypeOut,i1,i2,j1,j2);
-%     if ~CheckOverwrite && exist(OutputFile,'file')
-%         disp(['existing output file ' OutputFile ' already exists, skip to next field'])
-%         continue% skip iteration if the mode overwrite is desactivated and the result file already exists
-%     end
+    
+        %% generating the name of the merged field
+        i1=i1_series{1}(index);
+        if ~isempty(i2_series{end})
+            i2=i2_series{end}(index);
+        else
+            i2=i1;
+        end
+        j1=1;
+        j2=1;
+        if ~isempty(j1_series{1})
+            j1=j1_series{1}(index);
+            if ~isempty(j2_series{end})
+                j2=j2_series{end}(index);
+            else
+                j2=j1;
+            end
+        end
+       OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,'.nc',NomTypeOut,i1,i2,j1,j2);
+        if ~CheckOverwrite && exist(OutputFile,'file')
+            disp(['existing output file ' OutputFile ' already exists, skip to next field'])
+            continue% skip iteration if the mode overwrite is desactivated and the result file already exists
+        end
     
     %% z position
     ZIndex=mod(i1_series{1}(index)-1,NbSlice_calib{1})+1;%Zindex for phys transform
-    ZPos=ProjObjectCoord(ZIndex,3);
+    ZPosNew=ProjObjectCoord(ZIndex,3);
     if index==1
-    TimeData.ZPos=ZPos;
+        ZPos=ZPosNew;
     else
-        if ZPos~=TimeData.ZPos
+        if ZPosNew~=ZPos
             disp('inconsistent z positions in the series')
             return
         end
     end
-    % radius of the topography section at z pos
-    TopoRadius=0;
+    % radius of the topography section at z position
     ind_mask=[];
     if ZPos<20
@@ -372 +371 @@
             return
         end
-        disp([filecell{iview,index} ' read'])
         ListVar=Data{iview}.ListVarName;
         for ilist=1:numel(ListVar)
@@ -438 +436 @@
     end
     
+    
     %% time of the merged field: take the average of the different views
     if ~isempty(time)
@@ -456 +455 @@
         MergeData.div(ind_mask)=NaN;
     end
-    TimeData.time(index)=timeread;
-    TimeData.U(index,:,:)=Unew;
-    TimeData.V(index,:,:)=Vnew;
-    TimeData.curl(index,:,:)=MergeData.curl;
-    TimeData.div(index,:,:)=MergeData.div;
-     
-    %% recording the merged field
-    
-%     MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','NbCoord','NbDim','CoordUnit','ZPos'};
-%     MergeData.Conventions='uvmat';
-%     MergeData.NbCoord=2;
-%     MergeData.NbDim=2;
-%     MergeData.CoordUnit=CoordUnit;
-%     MergeData.ZPos=ZPos;
-%     % time interval of PIV
-%     Dt=[];
-%     if isfield(Data{1},'Dt')&& isnumeric(Data{1}.Dt)
-%         Dt=Data{1}.Dt;
-%     end
-%     for iview =2:numel(Data)
-%         if ~(isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Dt))
-%             Dt=[];%dt not the same for all fields
-%         end
-%     end
-%     if ~isempty(timeread)
-%         MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
-%         MergeData.Time=timeread;
-%     end
-%     
-%     % time unit
-%     if isfield(Data{1},'TimeUnit')
-%         TimeUnit=Data{1}.TimeUnit;
-%         for iview =2:numel(Data)
-%             if ~(isfield(Data{iview},'TimeUnit')&& isequal(Data{iview}.TimeUnit,TimeUnit))
-%                 TimeUnit=[];%TimeUnit not the same for all fields
-%             end
-%         end
-%         if ~isempty(TimeUnit)
-%             MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'TimeUnit'}];
-%             MergeData.TimeUnit=TimeUnit;
-%         end
-%     end
-
-end
-    error=struct2nc(OutputFile,TimeData);%save result file
-    if isempty(error)
-        disp(['output file ' OutputFile ' written'])
-    else
-        disp(error)
-    end
+    [npy,npx]=size(Unew);
+    
+    %% create the output file for the first iteration of the loop
+    if index==1
+        TimeData.ListGlobalAttribute={'Conventions','Project','CoordUnit','TimeUnit','ZPos','Time'};
+        TimeData.Conventions='uvmat';
+        TimeData.Project='2016_Circumpolar';
+        TimeData.CoordUnit='cm';
+        TimeData.TimeUnit='s';
+        TimeData.ZPos=ZPos;
+        TimeData.ListVarName={'radius','azimuth','U','V','curl','div'};
+        TimeData.VarDimName={'radius','azimuth',{'radius','azimuth'},{'radius','azimuth'}...
+            {'radius','azimuth'},{'radius','azimuth'}};
+        TimeData.VarAttribute{1}.Role='';
+        TimeData.VarAttribute{2}.Role='';
+        TimeData.VarAttribute{3}.Role='vector_x';
+        TimeData.VarAttribute{4}.Role='vector_y';
+        TimeData.VarAttribute{5}.Role='scalar';
+        TimeData.VarAttribute{6}.Role='scalar';
+        
+        TimeData.radius=radius_shifted;
+        TimeData.azimuth=azimuth_arclength;
+    end
+        
+        %% append data to the netcdf file for next iterations
+        TimeData.Time=timeread;
+       TimeData.U=Unew;
+       TimeData.V=Vnew;
+       TimeData.curl=MergeData.curl;
+       TimeData.div=MergeData.div;
+
+            [error,ncid]=struct2nc(OutputFile,TimeData);%save result file
+        if isempty(error)
+            disp(['output file ' OutputFile ' written'])
+        else
+            disp(error)
+        end
+            ellapsed_time=toc(tstart);
+    disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes'])
+end
+
 ellapsed_time=toc(tstart);
 disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes'])