Changeset 943


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

merge_proj_polar corrected

Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/geometry_calib.m

    r926 r943  
    331331        GeometryCalib.SliceCoord=Z_plane'*[0 0 1];
    332332    end
     333   
    333334else
    334335    GeometryCalib=[];
    335336    index=1;
    336337end
     338
    337339
    338340%------------------------------------------------------------------------
  • trunk/src/nc2struct.m

    r924 r943  
    276276                    end
    277277                end
    278                 Data.(VarName)=double(netcdf.getVar(nc,var_index(ivar)-1,ind_vec,ind_size)); %read the variable data
     278                Data.(VarName)=netcdf.getVar(nc,var_index(ivar)-1,ind_vec,ind_size); %read the variable data
    279279                Data.(VarName)=squeeze(Data.(VarName));%remove singeton dimension
    280280            else
  • trunk/src/plot_field.m

    r932 r943  
    958958            colormap(map);
    959959        else
    960             colormap('default'); % default matlab colormap ('jet')
     960            colormap('jet'); % default matlab colormap ('jet')
    961961        end
    962962       
  • 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'])
  • trunk/src/struct2nc.m

    r938 r943  
    11% 'struct2nc': create a netcdf file from a Matlab structure
    22%---------------------------------------------------------------------
    3 % errormsg=struct2nc(flname,Data)
     3% errormsg=struct2nc(flname,Data,action)
    44%
    55% OUTPUT:
     
    1616%       (optional) .VarAttribute: cell array of structures of the form .VarAttribute{ivar}.key=value, defining an attribute key name and value for the variable #ivar
    1717%      (requested) .Var1, .Var2....: variables (Matlab arrays) with names listed in .ListVarName
     18% action: if ='keep_open', don't close the file after creation, keep it open for further data input
    1819
    1920%=======================================================================
     
    3536%=======================================================================
    3637
    37 function [errormsg,nc]=struct2nc(flname,Data,action)
     38function [errormsg,nc]=struct2nc(flname,Data,action,ListDimName,DimValue,VarDimIndex)
    3839nc=[];
    39 % if ~ischar(flname)
    40 %     errormsg='invalid input for the netcf file name';
    41 %     return
    42 % end
     40if ~ischar(flname)
     41    errormsg='invalid input for the netcf file name';
     42    return
     43end
    4344if ~exist('Data','var')
    4445     errormsg='no data  input for the netcdf file';
    4546    return
    4647end
    47 if ~exist('action','var')
    48     action='one_input'; %fill the file with data and close it
    49 end
    5048
    5149
    5250%% check the validity of the input field structure
     51if ~ (exist('action','var') && strcmp(action,'keep_open'))
    5352[errormsg,ListDimName,DimValue,VarDimIndex]=check_field_structure(Data);
    5453if ~isempty(errormsg)
     
    5655    return
    5756end
     57end
    5858ListVarName=Data.ListVarName;
    5959
    6060%% create the netcdf file with name flname in format NETCDF4
    61 if ischar(flname)
     61% if ischar(flname)
    6262    FilePath=fileparts(flname);
    6363    if ~strcmp(FilePath,'') && ~exist(FilePath,'dir')
     
    6969    cmode = bitor(cmode, netcdf.getConstant('CLOBBER'));
    7070    nc = netcdf.create(flname, cmode);
    71 else
    72     nc=flname;
    73 end
     71% else
     72%     nc=flname;
     73% end
    7474
    7575%% write global constants
     
    113113for ivar=1:length(ListVarName)
    114114    if isfield(Data,ListVarName{ivar})
    115         VarClass=class(Data.(ListVarName{ivar}));
     115        VarClass{ivar}=class(Data.(ListVarName{ivar}));
    116116        VarType='';
    117         switch VarClass
     117        switch VarClass{ivar}
    118118            case {'single','double'}
    119119                VarType='nc_float'; % store all floating reals as single
     
    146146
    147147%% fill the variables with input data
    148 for ivar=1:length(ListVarName)
    149     if ~isnan(varid(ivar))
    150         VarVal=Data.(ListVarName{ivar});
    151         %varval=values of the current variable
    152         VarDimName=Data.VarDimName{ivar};
    153         if ischar(VarDimName)
    154             VarDimName={VarDimName};
    155         end
    156         siz=size(VarVal);
    157         testrange=(numel(VarDimName)==1 && strcmp(VarDimName{1},ListVarName{ivar}) && numel(VarVal)==2);% case of a coordinate defined on a regular mesh by the first and last values.
    158         testline=isequal(length(siz),2) && isequal(siz(1),1)&& isequal(siz(2), DimValue(VarDimIndex{ivar}));%matlab vector
    159         %testcolumn=isequal(length(siz),2) && isequal(siz(1), DimValue(VarDimIndex{ivar}))&& isequal(siz(2),1);%matlab column vector
    160         if testline || testrange
    161             if testrange
    162                 VarVal=linspace(VarVal(1),VarVal(2),DimValue(VarDimIndex{ivar}));% restitute the whole array of coordinate values
    163             end
    164             netcdf.putVar(nc,varid(ivar), double(VarVal'));
    165         else
    166             netcdf.putVar(nc,varid(ivar), double(VarVal));
    167         end     
    168     end
    169 end
    170 if strcmp(action,'one_input')
    171 netcdf.close(nc)
     148if ~(exist('action','var') && strcmp(action,'keep_open'))
     149    for ivar=1:length(ListVarName)
     150        if ~isnan(varid(ivar))
     151            VarVal=Data.(ListVarName{ivar});
     152            if strcmp(VarClass{ivar},'single')
     153                VarVal=single(VarVal);
     154            end
     155            %varval=values of the current variable
     156            VarDimName=Data.VarDimName{ivar};
     157            if ischar(VarDimName)
     158                VarDimName={VarDimName};
     159            end
     160            siz=size(VarVal);
     161            testrange=(numel(VarDimName)==1 && strcmp(VarDimName{1},ListVarName{ivar}) && numel(VarVal)==2);% case of a coordinate defined on a regular mesh by the first and last values.
     162            testline=isequal(length(siz),2) && isequal(siz(1),1)&& isequal(siz(2), DimValue(VarDimIndex{ivar}));%matlab vector
     163            %testcolumn=isequal(length(siz),2) && isequal(siz(1), DimValue(VarDimIndex{ivar}))&& isequal(siz(2),1);%matlab column vector
     164            if testline || testrange
     165                if testrange
     166                    VarVal=linspace(VarVal(1),VarVal(2),DimValue(VarDimIndex{ivar}));% restitute the whole array of coordinate values from the first and last values
     167                end
     168                netcdf.putVar(nc,varid(ivar), VarVal');
     169            else
     170                netcdf.putVar(nc,varid(ivar), VarVal);
     171            end
     172        end
     173    end
     174end
     175 netcdf.close(nc)
    172176[success,errormsg] = fileattrib(flname ,'+w');% allow writing access for the group of users
    173 end
    174 
    175177%'check_field_structure': check the validity of the field struture representation consistant with the netcdf format
    176178%------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.