Ignore:
Timestamp:
Apr 3, 2013, 10:21:53 AM (11 years ago)
Author:
sommeria
Message:

various bug corrections. Steps further for civ_series (still development needed)

File:
1 edited

Legend:

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

    r596 r599  
    4545if isstruct(Param) && isequal(Param.Action.RUN,0)
    4646    ParamOut.AllowInputSort='off';...% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    47     ParamOut.WholeIndexRange='off';...% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
    48     ParamOut.NbSlice='on'; ...%nbre of slices ('off' by default)
    49     ParamOut.VelType='two';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    50     ParamOut.FieldName='two';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
    51     ParamOut.FieldTransform = 'on';...%can use a transform function
    52     ParamOut.ProjObject='on';...%can use projection object(option 'off'/'on',
    53     ParamOut.Mask='off';...%can use mask option   (option 'off'/'on', 'off' by default)
    54     ParamOut.OutputDirExt='.tseries';%set the output dir extension
    55 return
     47        ParamOut.WholeIndexRange='off';...% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
     48        ParamOut.NbSlice='on'; ...%nbre of slices ('off' by default)
     49        ParamOut.VelType='two';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
     50        ParamOut.FieldName='two';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
     51        ParamOut.FieldTransform = 'on';...%can use a transform function
     52        ParamOut.ProjObject='on';...%can use projection object(option 'off'/'on',
     53        ParamOut.Mask='off';...%can use mask option   (option 'off'/'on', 'off' by default)
     54        ParamOut.OutputDirExt='.tseries';%set the output dir extension
     55    return
    5656end
    5757
     
    7878% filecell{iview,fileindex}:
    7979%        iview: line in the table corresponding to a given file series
    80 %        fileindex: file index within  the file series, 
    81 % i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j 
     80%        fileindex: file index within  the file series,
     81% i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j
    8282% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
    8383%%%%%%%%%%%%
     
    9090nbfield_i=size(i1_series{1},2); %nb of fields for the i index
    9191nbfield=nbfield_j*nbfield_i; %total number of fields
    92 nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices) 
     92nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices)
    9393nbfield=nbfield_i*NbSlice; %total number of fields after adjustement
    9494
    95 %determine the file type on each line from the first input file 
     95%determine the file type on each line from the first input file
    9696ImageTypeOptions={'image','multimage','mmreader','video'};
    9797NcTypeOptions={'netcdf','civx','civdata'};
     
    117117    if diff_time>0
    118118        displ_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)],checkrun)
    119     end   
     119    end
    120120    time=time(1,:);% choose the time data from the first sequence
    121121end
     
    130130
    131131%%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
    132  % EDIT FROM HERE
     132% EDIT FROM HERE
    133133
    134134%% check the validity of  ctinput file types
     
    137137elseif CheckNc{1}
    138138    FileExtOut='.nc';% write result as .nc files for netcdf inputs
    139 else 
     139else
    140140    displ_uvmat('ERROR',['invalid file type input ' FileType{1}],checkrun)
    141141    return
    142142end
    143143if nbview==2 && ~isequal(CheckImage{1},CheckImage{2})
    144         displ_uvmat('ERROR','input must be two image series or two netcdf file series',checkrun)
     144    displ_uvmat('ERROR','input must be two image series or two netcdf file series',checkrun)
    145145    return
    146146end
     
    196196end
    197197
    198 %% LOOP ON SLICES
    199198nbmissing=0; %number of undetected files
    200 for i_slice=1:NbSlice
    201     index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
    202     nbfile=0;
    203     nbmissing=0;
    204    
    205     %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    206     for index=index_slice 
    207         if checkrun
    208             stopstate=get(Param.RUNHandle,'BusyAction');
    209             update_waitbar(Param.WaitbarHandle,index/nbfield)
    210         else
    211             stopstate='queue';
    212         end
    213         if isequal(stopstate,'queue')% enable STOP command
    214             Data=cell(1,nbview);%initiate the set Data;
    215             nbtime=0;
    216             dt=[];
    217             %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    218             for iview=1:nbview
    219                 % reading input file(s)
    220                 [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
    221                 if ~isempty(errormsg)
    222                     errormsg=['time_series / read_field / ' errormsg];
    223                     display(errormsg)
    224                     break
     199% for i_slice=1:NbSlice
     200%index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
     201nbfile=0;
     202nbmissing=0;
     203
     204%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
     205for index=1:nbfield
     206    if checkrun
     207        stopstate=get(Param.RUNHandle,'BusyAction');
     208        update_waitbar(Param.WaitbarHandle,index/nbfield)
     209    else
     210        stopstate='queue';
     211    end
     212    if ~isequal(stopstate,'queue')% enable STOP command
     213        break
     214    end
     215    Data=cell(1,nbview);%initiate the set Data;
     216    nbtime=0;
     217    dt=[];
     218    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
     219    for iview=1:nbview
     220        % reading input file(s)
     221        [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
     222        if ~isempty(errormsg)
     223            errormsg=['time_series / read_field / ' errormsg];
     224            display(errormsg)
     225            break
     226        end
     227        if ~isempty(NbSlice_calib)
     228            Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
     229        end
     230    end
     231    if isempty(errormsg)
     232        Field=Data{1}; % default input field structure
     233        % coordinate transform (or other user defined transform)
     234        if ~isempty(transform_fct)
     235            switch nargin(transform_fct)
     236                case 4
     237                    if length(Data)==2
     238                        Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
     239                    else
     240                        Field=transform_fct(Data{1},XmlData{1});
     241                    end
     242                case 3
     243                    if length(Data)==2
     244                        Field=transform_fct(Data{1},XmlData{1},Data{2});
     245                    else
     246                        Field=transform_fct(Data{1},XmlData{1});
     247                    end
     248                case 2
     249                    Field=transform_fct(Data{1},XmlData{1});
     250                case 1
     251                    Field=transform_fct(Data{1});
     252            end
     253        end
     254       
     255        % calculate tps coefficients if needed
     256        if isfield(Param.ProjObject,'ProjMode')&& strcmp(Param.ProjObject.ProjMode,'interp_tps')
     257            Field=tps_coeff_field(Field,check_proj_tps);
     258        end
     259       
     260        %field projection on an object
     261        if Param.CheckObject
     262            [Field,errormsg]=proj_field(Field,Param.ProjObject);
     263            if ~isempty(errormsg)
     264                msgbox_uvmat('ERROR',['time_series / proj_field / ' errormsg])
     265                return
     266            end
     267        end
     268        nbfile=nbfile+1;
     269       
     270        % initiate the time series at the first iteration
     271        if nbfile==1
     272            % stop program if the first field reading is in error
     273            if ~isempty(errormsg)
     274                displ_uvmat('ERROR',['time_series / sub_field / ' errormsg],checkrun)
     275                return
     276            end
     277            DataOut=Field;%default
     278            DataOut.NbDim=Field.NbDim+1; %add the time dimension for plots
     279            nbvar=length(Field.ListVarName);
     280            if nbvar==0
     281                displ_uvmat('ERROR','no input variable selected',checkrun)
     282                return
     283            end
     284            testsum=2*ones(1,nbvar);%initiate flag for action on each variable
     285            if isfield(Field,'VarAttribute') % look for coordinate and flag variables
     286                for ivar=1:nbvar
     287                    if length(Field.VarAttribute)>=ivar && isfield(Field.VarAttribute{ivar},'Role')
     288                        var_role=Field.VarAttribute{ivar}.Role;%'role' of the variable
     289                        if isequal(var_role,'errorflag')
     290                            displ_uvmat('ERROR','do not handle error flags in time series',checkrun)
     291                            return
     292                        end
     293                        if isequal(var_role,'warnflag')
     294                            testsum(ivar)=0;  % not recorded variable
     295                            eval(['DataOut=rmfield(DataOut,''' Field.ListVarName{ivar} ''');']);%remove variable
     296                        end
     297                        if isequal(var_role,'coord_x')| isequal(var_role,'coord_y')|...
     298                                isequal(var_role,'coord_z')|isequal(var_role,'coord')
     299                            testsum(ivar)=1; %constant coordinates, record without time evolution
     300                        end
     301                    end
     302                    % check whether the variable ivar is a dimension variable
     303                    DimCell=Field.VarDimName{ivar};
     304                    if ischar(DimCell)
     305                        DimCell={DimCell};
     306                    end
     307                    if numel(DimCell)==1 && isequal(Field.ListVarName{ivar},DimCell{1})%detect dimension variables
     308                        testsum(ivar)=1;
     309                    end
    225310                end
    226                 if ~isempty(NbSlice_calib)
    227                     Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
     311            end
     312            for ivar=1:nbvar
     313                if testsum(ivar)==2
     314                    eval(['DataOut.' Field.ListVarName{ivar} '=[];'])
    228315                end
    229316            end
    230             if isempty(errormsg)
    231             Field=Data{1}; % default input field structure
    232             % coordinate transform (or other user defined transform)
    233             if ~isempty(transform_fct)
    234                 switch nargin(transform_fct)
    235                     case 4
    236                         if length(Data)==2
    237                             Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
    238                         else
    239                             Field=transform_fct(Data{1},XmlData{1});
     317            DataOut.ListVarName=[{'Time'} DataOut.ListVarName];
     318        end
     319       
     320        % add data to the current field
     321        for ivar=1:length(Field.ListVarName)
     322            VarName=Field.ListVarName{ivar};
     323            VarVal=Field.(VarName);
     324            if testsum(ivar)==2% test for recorded variable
     325                if isempty(errormsg)
     326                    if isequal(Param.ProjObject.ProjMode,'inside')% take the average in the domain for 'inside' mode
     327                        if isempty(VarVal)
     328                            displ_uvmat('ERROR',['empty result at frame index ' num2str(i1_series{iview}(index))],checkrun)
     329                            return
    240330                        end
    241                     case 3
    242                         if length(Data)==2
    243                             Field=transform_fct(Data{1},XmlData{1},Data{2});
    244                         else
    245                             Field=transform_fct(Data{1},XmlData{1});
    246                         end
    247                     case 2
    248                         Field=transform_fct(Data{1},XmlData{1});
    249                     case 1
    250                         Field=transform_fct(Data{1});
     331                        VarVal=mean(VarVal,1);
     332                    end
     333                    VarVal=shiftdim(VarVal,-1); %shift dimension
     334                    DataOut.(VarName)=cat(1,DataOut.(VarName),VarVal);%concanete the current field to the time series
     335                else
     336                    DataOut.(VarName)=cat(1,DataOut.(VarName),0);% put each variable to 0 in case of input reading error
    251337                end
    252             end
    253            
    254             % calculate tps coefficients if needed
    255             if isfield(Param.ProjObject,'ProjMode')&& strcmp(Param.ProjObject.ProjMode,'interp_tps')
    256                 Field=tps_coeff_field(Field,check_proj_tps);
    257             end
    258            
    259             %field projection on an object
    260             if Param.CheckObject
    261                 [Field,errormsg]=proj_field(Field,Param.ProjObject);
    262                 if ~isempty(errormsg)
    263                     msgbox_uvmat('ERROR',['time_series / proj_field / ' errormsg])
     338            elseif testsum(ivar)==1% variable representing fixed coordinates
     339                VarInit=DataOut.(VarName);
     340                if isempty(errormsg) && ~isequal(VarVal,VarInit)
     341                    displ_uvmat('ERROR',['time series requires constant coordinates ' VarName],checkrun)
    264342                    return
    265343                end
    266344            end
    267             nbfile=nbfile+1;
    268            
    269             % initiate the time series at the first iteration
    270             if nbfile==1
    271                 % stop program if the first field reading is in error
    272                 if ~isempty(errormsg)
    273                     displ_uvmat('ERROR',['time_series / sub_field / ' errormsg],checkrun)
    274                     return
    275                 end
    276                 DataOut=Field;%default
    277                 DataOut.NbDim=Field.NbDim+1; %add the time dimension for plots
    278                 nbvar=length(Field.ListVarName);
    279                 if nbvar==0
    280                     displ_uvmat('ERROR','no input variable selected',checkrun)
    281                     return
    282                 end
    283                 testsum=2*ones(1,nbvar);%initiate flag for action on each variable
    284                 if isfield(Field,'VarAttribute') % look for coordinate and flag variables
    285                     for ivar=1:nbvar
    286                         if length(Field.VarAttribute)>=ivar && isfield(Field.VarAttribute{ivar},'Role')
    287                             var_role=Field.VarAttribute{ivar}.Role;%'role' of the variable
    288                             if isequal(var_role,'errorflag')
    289                                 displ_uvmat('ERROR','do not handle error flags in time series',checkrun)
    290                                 return
    291                             end
    292                             if isequal(var_role,'warnflag')
    293                                 testsum(ivar)=0;  % not recorded variable
    294                                 eval(['DataOut=rmfield(DataOut,''' Field.ListVarName{ivar} ''');']);%remove variable
    295                             end
    296                             if isequal(var_role,'coord_x')| isequal(var_role,'coord_y')|...
    297                                     isequal(var_role,'coord_z')|isequal(var_role,'coord')
    298                                 testsum(ivar)=1; %constant coordinates, record without time evolution
    299                             end
    300                         end
    301                         % check whether the variable ivar is a dimension variable
    302                         DimCell=Field.VarDimName{ivar};
    303                         if ischar(DimCell)
    304                             DimCell={DimCell};
    305                         end
    306                         if numel(DimCell)==1 && isequal(Field.ListVarName{ivar},DimCell{1})%detect dimension variables
    307                             testsum(ivar)=1;
    308                         end
    309                     end
    310                 end
    311                 for ivar=1:nbvar
    312                     if testsum(ivar)==2
    313                         eval(['DataOut.' Field.ListVarName{ivar} '=[];'])
    314                     end
    315                 end
    316                 DataOut.ListVarName=[{'Time'} DataOut.ListVarName];
    317             end
    318            
    319             % add data to the current field
    320             for ivar=1:length(Field.ListVarName)
    321                 VarName=Field.ListVarName{ivar};
    322                 VarVal=Field.(VarName);
    323                 if testsum(ivar)==2% test for recorded variable
    324                     if isempty(errormsg)
    325                         if isequal(Param.ProjObject.ProjMode,'inside')% take the average in the domain for 'inside' mode
    326                             if isempty(VarVal)
    327                                 displ_uvmat('ERROR',['empty result at frame index ' num2str(i1_series{iview}(index))],checkrun)
    328                                 return
    329                             end
    330                             VarVal=mean(VarVal,1);
    331                         end
    332                         VarVal=shiftdim(VarVal,-1); %shift dimension
    333                         DataOut.(VarName)=cat(1,DataOut.(VarName),VarVal);%concanete the current field to the time series
    334                     else
    335                         DataOut.(VarName)=cat(1,DataOut.(VarName),0);% put each variable to 0 in case of input reading error
    336                     end
    337                 elseif testsum(ivar)==1% variable representing fixed coordinates
    338                     VarInit=DataOut.(VarName);
    339                     if isempty(errormsg) && ~isequal(VarVal,VarInit)
    340                         displ_uvmat('ERROR',['time series requires constant coordinates ' VarName],checkrun)
    341                         return
    342                     end
    343                 end
    344             end
    345            
    346             % record the time:
    347             if isempty(time)% time not set by xml filer(s)
    348                 if isfield(Data{1},'Time')
    349                     DataOut.Time(nbfile,1)=Field.Time;
    350                 else
    351                     DataOut.Time(nbfile,1)=index;%default
    352                 end
    353             else % time from ImaDoc prevails  TODO: correct
    354                 DataOut.Time(nbfile,1)=time(index);%
    355             end
    356            
    357             % record the number of missing input fields
    358             if ~isempty(errormsg)
    359                 nbmissing=nbmissing+1;
    360                 display(['index=' num2str(index) ':' errormsg])
    361             end
    362             end
    363         end
    364     end
    365     %%%%%%% END OF LOOP WITHIN A SLICE
     345        end
     346       
     347        % record the time:
     348        if isempty(time)% time not set by xml filer(s)
     349            if isfield(Data{1},'Time')
     350                DataOut.Time(nbfile,1)=Field.Time;
     351            else
     352                DataOut.Time(nbfile,1)=index;%default
     353            end
     354        else % time from ImaDoc prevails  TODO: correct
     355            DataOut.Time(nbfile,1)=time(index);%
     356        end
     357       
     358        % record the number of missing input fields
     359        if ~isempty(errormsg)
     360            nbmissing=nbmissing+1;
     361            display(['index=' num2str(index) ':' errormsg])
     362        end
     363    end
    366364   
    367     %remove time for global attributes if exists
    368     Time_index=find(strcmp('Time',DataOut.ListGlobalAttribute));
    369     if ~isempty(Time_index)
    370         DataOut.ListGlobalAttribute(Time_index)=[];
    371     end
    372     DataOut.Conventions='uvmat';
    373     for ivar=1:numel(DataOut.ListVarName)
    374         VarName=DataOut.ListVarName{ivar};
    375         eval(['DataOut.' VarName '=squeeze(DataOut.' VarName ');']) %remove singletons
    376     end
    377    
    378     % add time dimension
    379     for ivar=1:length(Field.ListVarName)
    380         DimCell=Field.VarDimName(ivar);
    381         if testsum(ivar)==2%variable used as time series
    382             DataOut.VarDimName{ivar}=[{'Time'} DimCell];
    383         elseif testsum(ivar)==1
    384             DataOut.VarDimName{ivar}=DimCell;
    385         end
    386     end
    387     indexremove=find(~testsum);
    388     if ~isempty(indexremove)
    389         DataOut.ListVarName(1+indexremove)=[];
    390         DataOut.VarDimName(indexremove)=[];
    391         if isfield(DataOut,'Role') && ~isempty(DataOut.Role{1})%generaliser aus autres attributs
    392             DataOut.Role(1+indexremove)=[];
    393         end
    394     end
    395    
    396     %shift variable attributes
    397     if isfield(DataOut,'VarAttribute')
    398         DataOut.VarAttribute=[{[]} DataOut.VarAttribute];
    399     end
    400     DataOut.VarDimName=[{'Time'} DataOut.VarDimName];
    401     DataOut.Action=Param.Action;%name of the processing programme
    402     test_time=diff(DataOut.Time)>0;% test that the readed time is increasing (not constant)
    403     if ~test_time
    404         DataOut.Time=1:filecounter;
    405     end
    406    
    407     % display nbmissing
    408     if ~isequal(nbmissing,0)
    409         displ_uvmat('WARNING',[num2str(nbmissing) ' files skipped: missing files or bad input, see command window display'],checkrun)
    410     end
    411    
    412     %name of result file
    413     OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(1),i1_series{1}(end),i_slice,[]);
    414     errormsg=struct2nc(OutputFile,DataOut); %save result file
    415     if isempty(errormsg)
    416         display([OutputFile ' written'])
    417     else
    418         displ_uvmat('ERROR',['error in Series/struct2nc: ' errormsg],checkrun)
    419     end
    420 end
     365end
     366%%%%%%% END OF LOOP WITHIN A SLICE
     367
     368%remove time for global attributes if exists
     369Time_index=find(strcmp('Time',DataOut.ListGlobalAttribute));
     370if ~isempty(Time_index)
     371    DataOut.ListGlobalAttribute(Time_index)=[];
     372end
     373DataOut.Conventions='uvmat';
     374for ivar=1:numel(DataOut.ListVarName)
     375    VarName=DataOut.ListVarName{ivar};
     376    eval(['DataOut.' VarName '=squeeze(DataOut.' VarName ');']) %remove singletons
     377end
     378
     379% add time dimension
     380for ivar=1:length(Field.ListVarName)
     381    DimCell=Field.VarDimName(ivar);
     382    if testsum(ivar)==2%variable used as time series
     383        DataOut.VarDimName{ivar}=[{'Time'} DimCell];
     384    elseif testsum(ivar)==1
     385        DataOut.VarDimName{ivar}=DimCell;
     386    end
     387end
     388indexremove=find(~testsum);
     389if ~isempty(indexremove)
     390    DataOut.ListVarName(1+indexremove)=[];
     391    DataOut.VarDimName(indexremove)=[];
     392    if isfield(DataOut,'Role') && ~isempty(DataOut.Role{1})%generaliser aus autres attributs
     393        DataOut.Role(1+indexremove)=[];
     394    end
     395end
     396
     397%shift variable attributes
     398if isfield(DataOut,'VarAttribute')
     399    DataOut.VarAttribute=[{[]} DataOut.VarAttribute];
     400end
     401DataOut.VarDimName=[{'Time'} DataOut.VarDimName];
     402DataOut.Action=Param.Action;%name of the processing programme
     403test_time=diff(DataOut.Time)>0;% test that the readed time is increasing (not constant)
     404if ~test_time
     405    DataOut.Time=1:filecounter;
     406end
     407
     408% display nbmissing
     409if ~isequal(nbmissing,0)
     410    displ_uvmat('WARNING',[num2str(nbmissing) ' files skipped: missing files or bad input, see command window display'],checkrun)
     411end
     412
     413%name of result file
     414OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(1),i1_series{1}(end),i_slice,[]);
     415errormsg=struct2nc(OutputFile,DataOut); %save result file
     416if isempty(errormsg)
     417    display([OutputFile ' written'])
     418else
     419    displ_uvmat('ERROR',['error in Series/struct2nc: ' errormsg],checkrun)
     420end
     421
    421422
    422423%% plot the time series (the last one in case of multislices)
     
    425426    haxes=axes;
    426427    plot_field(DataOut,haxes)
    427        
     428   
    428429    %% display the result file using the GUI get_field
    429430    hget_field=findobj(allchild(0),'name','get_field');
Note: See TracChangeset for help on using the changeset viewer.