Ignore:
Timestamp:
May 2, 2016, 9:49:22 AM (8 years ago)
Author:
sommeria
Message:

merge_proj_polar corrected

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/series/merge_proj_polar.m

    r941 r943  
    1 %'merge_proj': concatene several fields from series, can project them on a regular grid in phys coordinates
     1%'merge_proj': concatene several fields from series, project on a polar grid
    22%------------------------------------------------------------------------
    33% function ParamOut=merge_proj(Param)
     
    104104HeadData.azimuth=azimuth_arclength;   
    105105thresh2=16; % square of the interpolation range
     106
    106107%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
    107108ParamOut=[]; %default output
     
    160161
    161162%% prepare output file content
    162     TimeData.ListGlobalAttribute={'Conventions','Project','CoordUnit','TimeUnit','ZPos'};
    163     TimeData.Conventions='uvmat';
    164     TimeData.Project='2016_Circumpolar';
    165     TimeData.CoordUnit='cm';
    166     TimeData.TimeUnit='s';
    167     TimeData.ZPos=0;
    168     TimeData.ListVarName={'time','radius','azimuth','U','V','curl','div'};
    169     TimeData.VarDimName={'time','radius','azimuth',{'time','radius','azimuth'},{'time','radius','azimuth'}...
    170                         {'time','radius','azimuth'},{'time','radius','azimuth'}};
    171                     TimeData.VarAttribute{1}.Role='';
    172                     TimeData.VarAttribute{2}.Role='';
    173                     TimeData.VarAttribute{3}.Role='';
    174                     TimeData.VarAttribute{4}.Role='vector_x';
    175                     TimeData.VarAttribute{5}.Role='vector_y';
    176                     TimeData.VarAttribute{6}.Role='scalar';
    177                     TimeData.VarAttribute{7}.Role='scalar';
    178                     TimeData.time=nan(1,NbField);
    179                     TimeData.radius=radius_shifted;
    180                     TimeData.azimuth=azimuth_arclength;
    181                     nby=numel(radius);
    182                     nbx=numel(azimuth);
    183                     TimeData.U=nan(NbField,nby,nbx);
    184                     TimeData.V=nan(NbField,nby,nbx);
    185                     TimeData.curl=nan(NbField,nby,nbx);
    186                     TimeData.div=nan(NbField,nby,nbx);
     163%     TimeData.ListGlobalAttribute={'Conventions','Project','CoordUnit','TimeUnit','ZPos'};
     164%     TimeData.Conventions='uvmat';
     165%     TimeData.Project='2016_Circumpolar';
     166%     TimeData.CoordUnit='cm';
     167%     TimeData.TimeUnit='s';
     168%     TimeData.ZPos=0;
     169%     TimeData.ListVarName={'time','radius','azimuth','U','V','curl','div'};
     170%     TimeData.VarDimName={'time','radius','azimuth',{'time','radius','azimuth'},{'time','radius','azimuth'}...
     171%                         {'time','radius','azimuth'},{'time','radius','azimuth'}};
     172%                     TimeData.VarAttribute{1}.Role='';
     173%                     TimeData.VarAttribute{2}.Role='';
     174%                     TimeData.VarAttribute{3}.Role='';
     175%                     TimeData.VarAttribute{4}.Role='vector_x';
     176%                     TimeData.VarAttribute{5}.Role='vector_y';
     177%                     TimeData.VarAttribute{6}.Role='scalar';
     178%                     TimeData.VarAttribute{7}.Role='scalar';
     179%                     TimeData.time=nan(1,NbField);
     180%                     TimeData.radius=radius_shifted;
     181%                     TimeData.azimuth=azimuth_arclength;
     182%                     nby=numel(radius);
     183%                     nbx=numel(azimuth);
     184%                     TimeData.U=nan(NbField,nby,nbx);
     185%                     TimeData.V=nan(NbField,nby,nbx);
     186%                     TimeData.curl=nan(NbField,nby,nbx);
     187%                     TimeData.div=nan(NbField,nby,nbx);
    187188                   
    188189%     if ~isempty(timeread)
     
    254255for iview =2:numel(XmlData)
    255256    if ~(isfield(XmlData{iview},'GeometryCalib')&& isequal(XmlData{iview}.GeometryCalib.SliceCoord,ProjObjectCoord))...
    256         disp('error: geometric calibration missing')
     257        disp('error: geometric calibration missing or inconsistent plane positions')
    257258        return
    258259    end
     
    285286
    286287% %% output file type
    287 % FileExtOut='.nc'; %netcdf output
    288 % if isempty(j1_series{1})
    289 %     NomTypeOut='_1';
    290 % else
    291 %     NomTypeOut='_1_1';
    292 % end
    293 % RootFileOut=RootFile{1};
    294 % for iview=2:NbView
    295 %     if ~strcmp(RootFile{iview},RootFile{1})
    296 %         RootFileOut='mproj';
    297 %         break
    298 %     end
    299 % end
     288if isempty(j1_series{1})
     289    NomTypeOut='_1';
     290else
     291    NomTypeOut='_1_1';
     292end
     293RootFileOut=RootFile{1};
     294for iview=2:NbView
     295    if ~strcmp(RootFile{iview},RootFile{1})
     296        RootFileOut='mproj';
     297        break
     298    end
     299end
    300300
    301301
     
    313313for index=1:NbField
    314314    disp(['index=' num2str(index)])
    315     disp(['ellapsed time ' num2str(toc(tstart)/60,2) ' minutes'])
     315    disp(['ellapsed time ' num2str(toc(tstart)/60,4) ' minutes'])
    316316    update_waitbar(WaitbarHandle,index/NbField)
    317317    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
     
    319319        return
    320320    end
    321    
    322 %     %% generating the name of the merged field
    323 %     i1=i1_series{1}(index);
    324 %     if ~isempty(i2_series{end})
    325 %         i2=i2_series{end}(index);
    326 %     else
    327 %         i2=i1;
    328 %     end
    329 %     j1=1;
    330 %     j2=1;
    331 %     if ~isempty(j1_series{1})
    332 %         j1=j1_series{1}(index);
    333 %         if ~isempty(j2_series{end})
    334 %             j2=j2_series{end}(index);
    335 %         else
    336 %             j2=j1;
    337 %         end
    338 %     end
    339 %    % OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomTypeOut,i1,i2,j1,j2);
    340 %     if ~CheckOverwrite && exist(OutputFile,'file')
    341 %         disp(['existing output file ' OutputFile ' already exists, skip to next field'])
    342 %         continue% skip iteration if the mode overwrite is desactivated and the result file already exists
    343 %     end
     321    
     322        %% generating the name of the merged field
     323        i1=i1_series{1}(index);
     324        if ~isempty(i2_series{end})
     325            i2=i2_series{end}(index);
     326        else
     327            i2=i1;
     328        end
     329        j1=1;
     330        j2=1;
     331        if ~isempty(j1_series{1})
     332            j1=j1_series{1}(index);
     333            if ~isempty(j2_series{end})
     334                j2=j2_series{end}(index);
     335            else
     336                j2=j1;
     337            end
     338        end
     339       OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,'.nc',NomTypeOut,i1,i2,j1,j2);
     340        if ~CheckOverwrite && exist(OutputFile,'file')
     341            disp(['existing output file ' OutputFile ' already exists, skip to next field'])
     342            continue% skip iteration if the mode overwrite is desactivated and the result file already exists
     343        end
    344344   
    345345    %% z position
    346346    ZIndex=mod(i1_series{1}(index)-1,NbSlice_calib{1})+1;%Zindex for phys transform
    347     ZPos=ProjObjectCoord(ZIndex,3);
     347    ZPosNew=ProjObjectCoord(ZIndex,3);
    348348    if index==1
    349     TimeData.ZPos=ZPos;
     349        ZPos=ZPosNew;
    350350    else
    351         if ZPos~=TimeData.ZPos
     351        if ZPosNew~=ZPos
    352352            disp('inconsistent z positions in the series')
    353353            return
    354354        end
    355355    end
    356     % radius of the topography section at z pos
    357     TopoRadius=0;
     356    % radius of the topography section at z position
    358357    ind_mask=[];
    359358    if ZPos<20
     
    372371            return
    373372        end
    374         disp([filecell{iview,index} ' read'])
    375373        ListVar=Data{iview}.ListVarName;
    376374        for ilist=1:numel(ListVar)
     
    438436    end
    439437   
     438   
    440439    %% time of the merged field: take the average of the different views
    441440    if ~isempty(time)
     
    456455        MergeData.div(ind_mask)=NaN;
    457456    end
    458     TimeData.time(index)=timeread;
    459     TimeData.U(index,:,:)=Unew;
    460     TimeData.V(index,:,:)=Vnew;
    461     TimeData.curl(index,:,:)=MergeData.curl;
    462     TimeData.div(index,:,:)=MergeData.div;
    463      
    464     %% recording the merged field
    465    
    466 %     MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','NbCoord','NbDim','CoordUnit','ZPos'};
    467 %     MergeData.Conventions='uvmat';
    468 %     MergeData.NbCoord=2;
    469 %     MergeData.NbDim=2;
    470 %     MergeData.CoordUnit=CoordUnit;
    471 %     MergeData.ZPos=ZPos;
    472 %     % time interval of PIV
    473 %     Dt=[];
    474 %     if isfield(Data{1},'Dt')&& isnumeric(Data{1}.Dt)
    475 %         Dt=Data{1}.Dt;
    476 %     end
    477 %     for iview =2:numel(Data)
    478 %         if ~(isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Dt))
    479 %             Dt=[];%dt not the same for all fields
    480 %         end
    481 %     end
    482 %     if ~isempty(timeread)
    483 %         MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
    484 %         MergeData.Time=timeread;
    485 %     end
    486 %     
    487 %     % time unit
    488 %     if isfield(Data{1},'TimeUnit')
    489 %         TimeUnit=Data{1}.TimeUnit;
    490 %         for iview =2:numel(Data)
    491 %             if ~(isfield(Data{iview},'TimeUnit')&& isequal(Data{iview}.TimeUnit,TimeUnit))
    492 %                 TimeUnit=[];%TimeUnit not the same for all fields
    493 %             end
    494 %         end
    495 %         if ~isempty(TimeUnit)
    496 %             MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'TimeUnit'}];
    497 %             MergeData.TimeUnit=TimeUnit;
    498 %         end
    499 %     end
    500 
    501 end
    502     error=struct2nc(OutputFile,TimeData);%save result file
    503     if isempty(error)
    504         disp(['output file ' OutputFile ' written'])
    505     else
    506         disp(error)
    507     end
     457    [npy,npx]=size(Unew);
     458   
     459    %% create the output file for the first iteration of the loop
     460    if index==1
     461        TimeData.ListGlobalAttribute={'Conventions','Project','CoordUnit','TimeUnit','ZPos','Time'};
     462        TimeData.Conventions='uvmat';
     463        TimeData.Project='2016_Circumpolar';
     464        TimeData.CoordUnit='cm';
     465        TimeData.TimeUnit='s';
     466        TimeData.ZPos=ZPos;
     467        TimeData.ListVarName={'radius','azimuth','U','V','curl','div'};
     468        TimeData.VarDimName={'radius','azimuth',{'radius','azimuth'},{'radius','azimuth'}...
     469            {'radius','azimuth'},{'radius','azimuth'}};
     470        TimeData.VarAttribute{1}.Role='';
     471        TimeData.VarAttribute{2}.Role='';
     472        TimeData.VarAttribute{3}.Role='vector_x';
     473        TimeData.VarAttribute{4}.Role='vector_y';
     474        TimeData.VarAttribute{5}.Role='scalar';
     475        TimeData.VarAttribute{6}.Role='scalar';
     476       
     477        TimeData.radius=radius_shifted;
     478        TimeData.azimuth=azimuth_arclength;
     479    end
     480       
     481        %% append data to the netcdf file for next iterations
     482        TimeData.Time=timeread;
     483       TimeData.U=Unew;
     484       TimeData.V=Vnew;
     485       TimeData.curl=MergeData.curl;
     486       TimeData.div=MergeData.div;
     487
     488            [error,ncid]=struct2nc(OutputFile,TimeData);%save result file
     489        if isempty(error)
     490            disp(['output file ' OutputFile ' written'])
     491        else
     492            disp(error)
     493        end
     494            ellapsed_time=toc(tstart);
     495    disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes'])
     496end
     497
    508498ellapsed_time=toc(tstart);
    509499disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes'])
Note: See TracChangeset for help on using the changeset viewer.