Changeset 940 for trunk/src/series/merge_proj_polar.m
 Timestamp:
 Apr 12, 2016, 4:01:13 PM (8 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/merge_proj_polar.m
r939 r940 63 63 ParamOut.AllowInputSort='on';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) 64 64 ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) 65 ParamOut.NbSlice='o ff'; %nbre of slices ('off' by default)65 ParamOut.NbSlice='on'; %nbre of slices ('off' by default) 66 66 ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default) 67 ParamOut.FieldName='o ne';% 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) 68 68 ParamOut.FieldTransform = 'on';%can use a transform function 69 69 ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only) … … 90 90 % calculate the positions on which to interpolate 91 91 radius_ref=450;% radius of the mountain top 92 radius_shifted=1 20:2:120;% radius shifted by the radius of the origin at the topography summit92 radius_shifted=130:2:130;% radius shifted by the radius of the origin at the topography summit 93 93 radius=radius_ref+radius_shifted;%radius from centre of the tank 94 azimuth_arclength=(1 20:2:500);%azimuth in arc length at origin position94 azimuth_arclength=(150:2:400);%azimuth in arc length at origin position 95 95 azimuth=pi/2azimuth_arclength/radius_ref;%azimuth in radian 96 96 [Radius,Azimuth]=meshgrid(radius,azimuth); … … 103 103 HeadData.radius=radius_shifted; 104 104 HeadData.azimuth=azimuth_arclength; 105 105 thresh2=16; % square of the interpolation range 106 106 %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%% 107 107 ParamOut=[]; %default output … … 117 117 RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series 118 118 WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series 119 end120 121 %% define the directory for result file (with path=RootPath{1})122 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files123 124 if ~isfield(Param,'InputFields')125 Param.InputFields.FieldName='';126 119 end 127 120 … … 150 143 NbField=NbField_j*NbField_i; %total number of fields 151 144 145 %% define the name for result file (with path=RootPath{1}) 146 OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files 147 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},'.nc','_1',i1_series{1}(1)); 148 CheckOverwrite=1;%default 149 if isfield(Param,'CheckOverwrite') 150 CheckOverwrite=Param.CheckOverwrite; 151 end 152 if ~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 155 end 156 157 if ~isfield(Param,'InputFields') 158 Param.InputFields.FieldName=''; 159 end 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 152 208 %% determine the file type on each line from the first input file 153 209 ImageTypeOptions={'image','multimage','mmreader','video'}; … … 228 284 end 229 285 230 % % output file type231 FileExtOut='.nc'; %netcdf output232 if isempty(j1_series{1})233 NomTypeOut='_1';234 else235 NomTypeOut='_1_1';236 end237 RootFileOut=RootFile{1};238 for iview=2:NbView239 if ~strcmp(RootFile{iview},RootFile{1})240 RootFileOut='mproj';241 break242 end243 end286 % %% 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 244 300 245 301 … … 253 309 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% 254 310 tstart=tic; %used to record the computing time 255 CheckOverwrite=1;%default 256 if isfield(Param,'CheckOverwrite') 257 CheckOverwrite=Param.CheckOverwrite; 258 end 311 312 259 313 for index=1:NbField 260 314 update_waitbar(WaitbarHandle,index/NbField) … … 263 317 return 264 318 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; 270 348 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 300 362 %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% 301 363 Data=cell(1,NbView);%initiate the set Data … … 308 370 return 309 371 end 310 372 disp([filecell{iview,index} 'read']) 311 373 ListVar=Data{iview}.ListVarName; 312 374 for ilist=1:numel(ListVar) … … 336 398 [DataOut,VarAttribute,errormsg]=calc_field_tps(Data{iview}.Coord_tps,Data{iview}.NbCentre,Data{iview}.SubRange,... 337 399 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; 338 413 ListVarName=(fieldnames(DataOut))'; 414 VarDimName=cell(size(ListVarName)); 339 415 ProjData{iview}=HeadData; 340 416 ProjData{iview}.ListVarName= [ProjData{iview}.ListVarName ListVarName]; … … 343 419 for ivar=1:numel(ListVarName) 344 420 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 351 428 end 352 429 %%%%%%%%%%%%%%%% END LOOP ON VIEWS %%%%%%%%%%%%%%%% … … 370 447 %% rotating the velocity vectors to the local axis of the polatr coordinates 371 448 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'); 374 450 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; 377 453 MergeData.curl(ind_mask)=NaN; 378 454 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 381 462 %% recording the merged field 382 463 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 499 end 500 error=struct2nc(OutputFile,TimeData);%save result file 418 501 if isempty(error) 419 502 disp(['output file ' OutputFile ' written']) … … 421 504 disp(error) 422 505 end 423 end424 425 506 ellapsed_time=toc(tstart); 426 507 disp(['total ellapsed time ' num2str(ellapsed_time/60,2) ' minutes']) … … 439 520 NbView=length(Data); 440 521 if NbView==1% if there is only one field, just reproduce it in MergeData 441 return 522 return 442 523 end 443 524 … … 484 565 errormsg='sizes of the input matrices do not agree, need to interpolate on a common grid using a projection object'; 485 566 return 486 else 567 else 487 568 MergeData.(VarName)=MergeData.(VarName) +double(Data{iview}.(VarName));%add data 488 569 NbAver=NbAver + ~check_bad;% add 1 for good data, 0 else … … 493 574 end 494 575 end 495 % if isempty(FFName)496 % FFName='FF';497 % end498 % MergeData.(FFName)(NbAver~=0)=0;% flag to 1 undefined summed data499 % MergeData.(FFName)(NbAver==0)=1;% flag to 1 undefined summed data500 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 582 end 583 584 585
Note: See TracChangeset
for help on using the changeset viewer.