Changeset 871 for trunk/src/series/aver_stat.m
 Timestamp:
 Feb 16, 2015, 12:15:23 AM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/aver_stat.m
r869 r871 71 71 ParamOut.OutputDirExt='.stat';%set the output dir extension 72 72 ParamOut.OutputFileMode='NbSlice';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice 73 % check the existence of the first file in the series 73 % check for selection of a projection object 74 hseries=findobj(allchild(0),'Tag','series');% handles of the GUI series 75 if ~isfield(Param,'ProjObject') 76 answer=msgbox_uvmat('INPUT_YN','use a projection object?'); 77 if strcmp(answer,'Yes') 78 hhseries=guidata(hseries); 79 set(hhseries.CheckObject,'Visible','on') 80 set(hhseries.CheckObject,'Value',1) 81 Param.CheckObject=1; 82 series('CheckObject_Callback',hseries,[],hhseries); %file input with xml reading in uvmat, show the image in phys coordinates 83 end 84 end 85 % introduce bin size for histograms 86 if Param.CheckObject 87 SeriesData=get(hseries,'UserData'); 88 if ismember(SeriesData.ProjObject.ProjMode,{'inside','outside'}) 89 answer=msgbox_uvmat('INPUT_TXT','set bin size for histograms (or keep ''auto'' by default)?','auto'); 90 ParamOut.ActionInput.VarMesh=str2double(answer); 91 end 92 end 93 % check the existence of the first and last file in the series 74 94 first_j=[]; 75 95 if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end … … 87 107 LastFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},... 88 108 Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2); 89 if ~exist( FirstFileName,'file')109 if ~exist(LastFileName,'file') 90 110 msgbox_uvmat('WARNING',['the last input file ' LastFileName ' does not exist']) 91 111 end … … 170 190 % EDIT FROM HERE 171 191 172 %% check the validity of input file types 192 %% check the validity of input file types and set the output file type 173 193 if CheckImage{1} 174 194 FileExtOut='.png'; % write result as .png images for image inputs … … 182 202 disp_uvmat('ERROR','input must be two image series or two netcdf file series',checkrun) 183 203 return 204 end 205 if isfield(Param,'ProjObject') && ~strcmp(Param.ProjObject.Type,'plane') 206 FileExtOut='.nc';% write result as .nc files (even for image input) 184 207 end 185 208 … … 212 235 end 213 236 end 214 % for i_slice=1:NbSlice 215 % index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice 216 nbfiles=0; 217 % nbmissing=0;218 237 nbfiles=0;%counter of the successfully read files (bad files are skipped) 238 VarMesh=NaN; 239 if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})&& isfield(Param.ActionInput,'VarMesh')%case of histograms 240 VarMesh=Param.ActionInput.VarMesh; 241 end 219 242 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% 220 243 for index=1:NbField … … 243 266 if isempty(errormsg) 244 267 Field=Data{1}; % default input field structure 268 nbfiles=nbfiles+1; %increment the file counter 245 269 %% coordinate transform (or other user defined transform) 246 270 if ~isempty(transform_fct) … … 265 289 end 266 290 267 %% calculate tps coefficients if needed 268 if isfield(Param,'ProjObject')&&isfield(Param.ProjObject,'ProjMode')&& strcmp(Param.ProjObject.ProjMode,'interp_tps') 269 Field=tps_coeff_field(Field,check_proj_tps); 270 end 271 272 %field projection on an object 291 %% field projection on an object 273 292 if Param.CheckObject 274 [Field,errormsg]=proj_field(Field,Param.ProjObject); 293 if strcmp(Param.ProjObject.ProjMode,'interp_tps') 294 Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed 295 end 296 [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh); 275 297 if ~isempty(errormsg) 276 298 disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun) … … 278 300 end 279 301 end 280 nbfiles=nbfiles+1; 281 302 282 303 %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% 283 304 if nbfiles==1 %first field … … 286 307 time_1=Field.Time(1); 287 308 end 288 DataOut=Field;%default 289 DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files 290 errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable 291 for ivar=1:numel(Field.ListVarName) 292 VarName=Field.ListVarName{ivar}; 293 DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero 294 NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero 295 for iivar=1:length(Field.VarAttribute) 296 if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')... 297 && strcmp(Field.VarAttribute{iivar}.Role,'errorflag') 298 errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar 309 DataOut=Field;%outcome reproduces the first (projected) field by default 310 DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files 311 if ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms 312 for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram) 313 VarName=Field.ListVarName{ivar}; 314 if isfield(Data{1},VarName) 315 DataOut.(VarName)=Field.(VarName); 316 DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName))); 317 VarMesh=DataOut.(VarName)(2)DataOut.(VarName)(1); 299 318 end 300 319 end 301 end 302 DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result 303 DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result 304 DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result 320 disp(['mesh for histogram = ' num2str(VarMesh)]) 321 else 322 errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable 323 for ivar=1:numel(Field.ListVarName) 324 VarName=Field.ListVarName{ivar}; 325 DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero 326 NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero 327 for iivar=1:length(Field.VarAttribute) 328 if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')... 329 && strcmp(Field.VarAttribute{iivar}.Role,'errorflag') 330 errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar 331 end 332 end 333 end 334 DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result 335 DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result 336 DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result 337 end 305 338 end %current field 306 339 for ivar=1:length(DataOut.ListVarName) 307 VarName= Field.ListVarName{ivar};340 VarName=DataOut.ListVarName{ivar}; 308 341 sizmean=size(DataOut.(VarName)); 309 342 siz=size(Field.(VarName)); 310 if ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean) 343 if ismember(Param.ProjObject.ProjMode,{'inside','outside'}) 344 if isfield(Data{1},VarName) 345 MaxValue=max(DataOut.(VarName));% current max of histogram absissa 346 MinValue=min(DataOut.(VarName));% current min of histogram absissa 347 % VarMesh=Field.VarAttribute{ivar}.Mesh; 348 MaxIndex=round(MaxValue/VarMesh); 349 MinIndex=round(MinValue/VarMesh); 350 MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field 351 MinIndex_new=round(min(Field.(VarName)/VarMesh)); 352 if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one 353 DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values 354 DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_newMaxIndex)]; % append the new histo values 355 end 356 if MinIndex_new <= MinIndex1 357 DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex1) DataOut.(VarName)];% insert the new variable values 358 DataOut.([VarName 'Histo'])=[zeros(1,MinIndexMinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values 359 ind_start=1; 360 else 361 ind_start=MinIndex_newMinIndex+1; 362 end 363 DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_newMinIndex_new)=... 364 DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_newMinIndex_new)+Field.([VarName 'Histo']); 365 end 366 elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean) 311 367 disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project on a grid'],checkrun) 312 368 return … … 328 384 end 329 385 %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%% 330 331 for ivar=1:length(Field.ListVarName) 332 VarName=Field.ListVarName{ivar}; 333 DataOut.(VarName)=DataOut.(VarName)./NbData.(VarName); % normalize the mean 386 if ~ismember(Param.ProjObject.ProjMode,{'inside','outside'}) 387 for ivar=1:length(Field.ListVarName) 388 VarName=Field.ListVarName{ivar}; 389 DataOut.(VarName)=DataOut.(VarName)./NbData.(VarName); % normalize the mean 390 end 334 391 end 335 392 nbmissing=NbFieldnbfiles; 336 393 if nbmissing~=0 337 disp_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skip ted'],checkrun)394 disp_uvmat('WARNING',[num2str(nbmissing) ' input files are missing or skipped'],checkrun) 338 395 end 339 396 if isempty(time) % time is read from files … … 352 409 %% writing the result file 353 410 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j); 354 if CheckImage{1}%case of images411 if strcmp(FileExtOut,'.png') %case of images 355 412 if isequal(FileInfo{1}.BitDepth,16)(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16)) 356 413 DataOut.A=uint16(DataOut.A); … … 361 418 end 362 419 disp([OutputFile ' written']); 363 else %case of netcdf inputfile , determine global attributes420 else %case of netcdf file , determine global attributes 364 421 errormsg=struct2nc(OutputFile,DataOut); %save result file 365 422 if isempty(errormsg)
Note: See TracChangeset
for help on using the changeset viewer.