Changeset 940


Ignore:
Timestamp:
Apr 12, 2016, 4:01:13 PM (9 years ago)
Author:
sommeria
Message:

merge_proj_polar made

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series.m

    r937 r940  
    14241424    return
    14251425end
    1426 if isfield(Param,'InputFields')&& isequal(Param.InputFields.FieldName,'get_field...')
     1426if isfield(Param,'InputFields')&& isfield(Param.InputFields,'FieldName')&& isequal(Param.InputFields.FieldName,'get_field...')
    14271427    errormsg='input field name(s) not defined, select get_field...';
    14281428    return
    … …  
    22942294
    22952295%% Detect the types of input files and set menus and default options in 'FieldName'
    2296 if FieldNameRequest && numel(iview_netcdf)>=1
    2297     set(handles.InputFields,'Visible','on')
    2298     %if CheckList==0        % not civ input made
    2299         if isfield(SeriesData.FileInfo{iview_netcdf(1)},'ListVarName')
     2296if (FieldNameRequest || VelTypeRequest) && numel(iview_netcdf)>=1
     2297    set(handles.InputFields,'Visible','on')   
     2298    if FieldNameRequest && isfield(SeriesData.FileInfo{iview_netcdf(1)},'ListVarName')
     2299        set(handles.FieldName,'Visible','on')
    23002300        ListVarName=SeriesData.FileInfo{iview_netcdf(1)}.ListVarName;
    23012301        ind_var=get(handles.FieldName,'Value');%indices of previously selected variables
    … …  
    23542354            set(handles.FieldName_1,'Visible','off')
    23552355        end
    2356         end
    2357 %     else
    2358 %         set(handles_coord,'Visible','off')% no coord display for civ data
    2359 %     end
     2356    else
     2357        set(handles.FieldName,'Visible','off')
     2358    end
    23602359    set(handles.FieldName,'String',FieldList)
    23612360else

  • trunk/src/series/merge_proj_polar.m

    r939 r940  
    6363    ParamOut.AllowInputSort='on';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    6464    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
    65     ParamOut.NbSlice='off'; %nbre of slices ('off' by default)
     65    ParamOut.NbSlice='on'; %nbre of slices ('off' by default)
    6666    ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    67     ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
     67    ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
    6868    ParamOut.FieldTransform = 'on';%can use a transform function
    6969    ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only)
    … …  
    9090% calculate the positions on which to interpolate
    9191radius_ref=450;% radius of the mountain top
    92 radius_shifted=-120:2:120;% radius shifted by the radius of the origin at the topography summit
     92radius_shifted=-130:2:130;% radius shifted by the radius of the origin at the topography summit
    9393radius=radius_ref+radius_shifted;%radius from centre of the tank
    94 azimuth_arclength=(-120:2:500);%azimuth in arc length at origin position
     94azimuth_arclength=(-150:2:400);%azimuth in arc length at origin position
    9595azimuth=pi/2-azimuth_arclength/radius_ref;%azimuth in radian
    9696[Radius,Azimuth]=meshgrid(radius,azimuth);
    … …  
    103103HeadData.radius=radius_shifted;
    104104HeadData.azimuth=azimuth_arclength;   
    105 
     105thresh2=16; % square of the interpolation range
    106106%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
    107107ParamOut=[]; %default output
    … …  
    117117    RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
    118118    WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
    119 end
    120 
    121 %% define the directory for result file (with path=RootPath{1})
    122 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
    123 
    124 if ~isfield(Param,'InputFields')
    125     Param.InputFields.FieldName='';
    126119end
    127120
    … …  
    150143NbField=NbField_j*NbField_i; %total number of fields
    151144
     145%% define the name for result file (with path=RootPath{1})
     146OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
     147OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1',i1_series{1}(1));
     148CheckOverwrite=1;%default
     149if isfield(Param,'CheckOverwrite')
     150    CheckOverwrite=Param.CheckOverwrite;
     151end
     152if ~CheckOverwrite && exist(OutputFile,'file')
     153    disp(['existing output file ' OutputFile ' already exists, skip to next field'])
     154    return% skip iteration if the mode overwrite is desactivated and the result file already exists
     155end
     156
     157if ~isfield(Param,'InputFields')
     158    Param.InputFields.FieldName='';
     159end
     160
     161%% 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);
     187                   
     188%     if ~isempty(timeread)
     189%         MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
     190%         MergeData.Time=timeread;
     191%     end
     192   
     193    % time unit
     194%     if isfield(Data{1},'TimeUnit')
     195%         TimeUnit=Data{1}.TimeUnit;
     196%         for iview =2:numel(Data)
     197%             if ~(isfield(Data{iview},'TimeUnit')&& isequal(Data{iview}.TimeUnit,TimeUnit))
     198%                 TimeUnit=[];%TimeUnit not the same for all fields
     199%             end
     200%         end
     201%         if ~isempty(TimeUnit)
     202%             MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'TimeUnit'}];
     203%             MergeData.TimeUnit=TimeUnit;
     204%         end
     205%     end
     206
     207
    152208%% determine the file type on each line from the first input file
    153209ImageTypeOptions={'image','multimage','mmreader','video'};
    … …  
    228284end
    229285
    230 %% output file type
    231 FileExtOut='.nc'; %netcdf output
    232 if isempty(j1_series{1})
    233     NomTypeOut='_1';
    234 else
    235     NomTypeOut='_1_1';
    236 end
    237 RootFileOut=RootFile{1};
    238 for iview=2:NbView
    239     if ~strcmp(RootFile{iview},RootFile{1})
    240         RootFileOut='mproj';
    241         break
    242     end
    243 end
     286% %% 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
    244300
    245301
    … …  
    253309%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    254310tstart=tic; %used to record the computing time
    255 CheckOverwrite=1;%default
    256 if isfield(Param,'CheckOverwrite')
    257     CheckOverwrite=Param.CheckOverwrite;
    258 end
     311
     312
    259313for index=1:NbField
    260314    update_waitbar(WaitbarHandle,index/NbField)
    … …  
    263317        return
    264318    end
    265    
    266     %% generating the name of the merged field
    267     i1=i1_series{1}(index);
    268     if ~isempty(i2_series{end})
    269         i2=i2_series{end}(index);
     319   
     320%     %% generating the name of the merged field
     321%     i1=i1_series{1}(index);
     322%     if ~isempty(i2_series{end})
     323%         i2=i2_series{end}(index);
     324%     else
     325%         i2=i1;
     326%     end
     327%     j1=1;
     328%     j2=1;
     329%     if ~isempty(j1_series{1})
     330%         j1=j1_series{1}(index);
     331%         if ~isempty(j2_series{end})
     332%             j2=j2_series{end}(index);
     333%         else
     334%             j2=j1;
     335%         end
     336%     end
     337%    % OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomTypeOut,i1,i2,j1,j2);
     338%     if ~CheckOverwrite && exist(OutputFile,'file')
     339%         disp(['existing output file ' OutputFile ' already exists, skip to next field'])
     340%         continue% skip iteration if the mode overwrite is desactivated and the result file already exists
     341%     end
     342   
     343    %% z position
     344    ZIndex=mod(i1_series{1}(index)-1,NbSlice_calib{1})+1;%Zindex for phys transform
     345    ZPos=ProjObjectCoord(ZIndex,3);
     346    if index==1
     347    TimeData.ZPos=ZPos;
    270348    else
    271         i2=i1;
    272     end
    273     j1=1;
    274     j2=1;
    275     if ~isempty(j1_series{1})
    276         j1=j1_series{1}(index);
    277         if ~isempty(j2_series{end})
    278             j2=j2_series{end}(index);
    279         else
    280             j2=j1;
    281         end
    282     end
    283     OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFileOut,FileExtOut,NomTypeOut,i1,i2,j1,j2);
    284     if ~CheckOverwrite && exist(OutputFile,'file')
    285         disp(['existing output file ' OutputFile ' already exists, skip to next field'])
    286         continue% skip iteration if the mode overwrite is desactivated and the result file already exists
    287     end
    288    
    289     %% z position
    290        ZIndex=mod(i1_series{1}(index)-1,NbSlice_calib{1})+1;%Zindex for phys transform
    291        ZPos=ProjObjectCoord(ZIndex,3);
    292        % radius of the topography section at z pos
    293        TopoRadius=0;
    294        ind_mask=[];
    295        if ZPos<20
    296     TopoRadius=40*sin(acos((20+ZPos)/40));
    297        ind_mask=(XI'.*XI'+YI'.*YI')<TopoRadius*TopoRadius;% indidces of data to mask
    298        end
    299        
     349        if ZPos~=TimeData.ZPos
     350            disp('inconsistent z positions in the series')
     351            return
     352        end
     353    end
     354    % radius of the topography section at z pos
     355    TopoRadius=0;
     356    ind_mask=[];
     357    if ZPos<20
     358        TopoRadius=40*sin(acos((20+ZPos)/40));
     359        ind_mask=(XI'.*XI'+YI'.*YI')<TopoRadius*TopoRadius;% indidces of data to mask
     360    end
     361   
    300362    %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    301363    Data=cell(1,NbView);%initiate the set Data
    … …  
    308370            return
    309371        end
    310        
     372        disp([filecell{iview,index} 'read'])
    311373        ListVar=Data{iview}.ListVarName;
    312374        for ilist=1:numel(ListVar)
    … …  
    336398        [DataOut,VarAttribute,errormsg]=calc_field_tps(Data{iview}.Coord_tps,Data{iview}.NbCentre,Data{iview}.SubRange,...
    337399            cat(3,Data{iview}.U_tps,Data{iview}.V_tps),FieldNames,cat(3,XI,YI));
     400        % set to NaN interpolation points which are too far from any initial data (more than 2 CoordMesh)
     401        Coord=permute(Data{iview}.Coord_tps,[1 3 2]);
     402        Coord=reshape(Coord,size(Coord,1)*size(Coord,2),2);
     403        if exist('scatteredInterpolant','file')%recent Matlab versions
     404            F=scatteredInterpolant(Coord,Coord(:,1),'nearest');
     405            G=scatteredInterpolant(Coord,Coord(:,2),'nearest');
     406        else
     407            F=TriScatteredInterp(Coord,Coord(:,1),'nearest');
     408            G=TriScatteredInterp(Coord,Coord(:,2),'nearest');
     409        end
     410        Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point
     411        Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point
     412        Dist=Distx.*Distx+Disty.*Disty;
    338413        ListVarName=(fieldnames(DataOut))';
     414        VarDimName=cell(size(ListVarName));
    339415        ProjData{iview}=HeadData;
    340416        ProjData{iview}.ListVarName= [ProjData{iview}.ListVarName ListVarName];
    … …  
    343419        for ivar=1:numel(ListVarName)
    344420            ProjData{iview}.VarDimName{ivar+2}={'radius','azimuth'};
    345             ProjData{iview}.(ListVarName{ivar})=(DataOut.(ListVarName{ivar}))';
    346         end
    347         %% mask
    348         %         if Param.CheckMask && ~isempty(MaskData{iview})
    349         %              [Data{iview},errormsg]=mask_proj(Data{iview},MaskData{iview});
    350         %         end
     421            VarName=ListVarName{ivar};
     422            if ~isempty(thresh2)
     423                DataOut.(VarName)(Dist>thresh2)=NaN;% put to NaN interpolated positions further than RangeInterp from initial data
     424            end
     425            ProjData{iview}.(VarName)=(DataOut.(VarName))';
     426        end
     427       
    351428    end
    352429    %%%%%%%%%%%%%%%% END LOOP ON VIEWS %%%%%%%%%%%%%%%%
    … …  
    370447    %% rotating the velocity vectors to the local axis of the polatr coordinates
    371448    Unew=MergeData.U.*sin(Azimuth')-MergeData.V.*cos(Azimuth');
    372     MergeData.V=MergeData.U.*cos(Azimuth')+MergeData.V.*sin(Azimuth');
    373     MergeData.U=Unew;
     449    Vnew=MergeData.U.*cos(Azimuth')+MergeData.V.*sin(Azimuth');
    374450    if ~isempty(ind_mask)
    375         MergeData.U(ind_mask)=NaN;
    376         MergeData.V(ind_mask)=NaN;
     451        Unew(ind_mask)=NaN;
     452        Vnew(ind_mask)=NaN;
    377453        MergeData.curl(ind_mask)=NaN;
    378454        MergeData.div(ind_mask)=NaN;
    379     end
    380    
     455    end
     456    TimeData.time(index)=timeread;
     457    TimeData.U(index,:,:)=Unew;
     458    TimeData.V(index,:,:)=Vnew;
     459    TimeData.curl(index,:,:)=MergeData.curl;
     460    TimeData.div(index,:,:)=MergeData.div;
     461     
    381462    %% recording the merged field
    382463   
    383     MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','NbCoord','NbDim','CoordUnit','ZPos'};
    384     MergeData.Conventions='uvmat';
    385     MergeData.NbCoord=2;
    386     MergeData.NbDim=2;
    387     MergeData.CoordUnit=CoordUnit;
    388     MergeData.ZPos=ZPos;
    389     % time interval of PIV
    390     Dt=[];
    391     if isfield(Data{1},'Dt')&& isnumeric(Data{1}.Dt)
    392         Dt=Data{1}.Dt;
    393     end
    394     for iview =2:numel(Data)
    395         if ~(isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Dt))
    396             Dt=[];%dt not the same for all fields
    397         end
    398     end
    399     if ~isempty(timeread)
    400         MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
    401         MergeData.Time=timeread;
    402     end
    403 
    404     % time unit
    405     if isfield(Data{1},'TimeUnit')
    406         TimeUnit=Data{1}.TimeUnit;
    407         for iview =2:numel(Data)
    408             if ~(isfield(Data{iview},'TimeUnit')&& isequal(Data{iview}.TimeUnit,TimeUnit))
    409                 TimeUnit=[];%TimeUnit not the same for all fields
    410             end
    411         end
    412         if ~isempty(TimeUnit)
    413             MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'TimeUnit'}];
    414             MergeData.TimeUnit=TimeUnit;
    415         end
    416     end
    417     error=struct2nc(OutputFile,MergeData);%save result file
     464%     MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','NbCoord','NbDim','CoordUnit','ZPos'};
     465%     MergeData.Conventions='uvmat';
     466%     MergeData.NbCoord=2;
     467%     MergeData.NbDim=2;
     468%     MergeData.CoordUnit=CoordUnit;
     469%     MergeData.ZPos=ZPos;
     470%     % time interval of PIV
     471%     Dt=[];
     472%     if isfield(Data{1},'Dt')&& isnumeric(Data{1}.Dt)
     473%         Dt=Data{1}.Dt;
     474%     end
     475%     for iview =2:numel(Data)
     476%         if ~(isfield(Data{iview},'Dt')&& isequal(Data{iview}.Dt,Dt))
     477%             Dt=[];%dt not the same for all fields
     478%         end
     479%     end
     480%     if ~isempty(timeread)
     481%         MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'Time'}];
     482%         MergeData.Time=timeread;
     483%     end
     484%     
     485%     % time unit
     486%     if isfield(Data{1},'TimeUnit')
     487%         TimeUnit=Data{1}.TimeUnit;
     488%         for iview =2:numel(Data)
     489%             if ~(isfield(Data{iview},'TimeUnit')&& isequal(Data{iview}.TimeUnit,TimeUnit))
     490%                 TimeUnit=[];%TimeUnit not the same for all fields
     491%             end
     492%         end
     493%         if ~isempty(TimeUnit)
     494%             MergeData.ListGlobalAttribute=[MergeData.ListGlobalAttribute {'TimeUnit'}];
     495%             MergeData.TimeUnit=TimeUnit;
     496%         end
     497%     end
     498
     499end
     500    error=struct2nc(OutputFile,TimeData);%save result file
    418501    if isempty(error)
    419502        disp(['output file ' OutputFile ' written'])
    … …  
    421504        disp(error)
    422505    end
    423 end
    424 
    425506ellapsed_time=toc(tstart);
    426507disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes'])
    … …  
    439520NbView=length(Data);
    440521if NbView==1% if there is only one field, just reproduce it in MergeData
    441     return 
     522    return
    442523end
    443524
    … …  
    484565                            errormsg='sizes of the input matrices do not agree, need to interpolate on a common grid using a projection object';
    485566                            return
    486                         else                             
     567                        else
    487568                            MergeData.(VarName)=MergeData.(VarName) +double(Data{iview}.(VarName));%add data
    488569                            NbAver=NbAver + ~check_bad;% add 1 for good data, 0 else
    … …  
    493574                end
    494575        end
    495 %         if isempty(FFName)
    496 %             FFName='FF';
    497 %         end
    498 %         MergeData.(FFName)(NbAver~=0)=0;% flag to 1 undefined summed data
    499 %         MergeData.(FFName)(NbAver==0)=1;% flag to 1 undefined summed data
    500     end
    501 end
    502 
    503 
    504    
     576        %         if isempty(FFName)
     577        %             FFName='FF';
     578        %         end
     579        %         MergeData.(FFName)(NbAver~=0)=0;% flag to 1 undefined summed data
     580        %         MergeData.(FFName)(NbAver==0)=1;% flag to 1 undefined summed data
     581    end
     582end
     583
     584
     585   
Note: See TracChangeset for help on using the changeset viewer.