[958]  1  %'read_civdata': reads new civ data from netcdf files


 2  %


 3  %


 4  % function [Field,VelTypeOut]=read_civdata(FileName,FieldNames,VelType)


 5  %


 6  % OUTPUT:


 7  % Field: structure representing the selected field, containing


 8  % .Txt: (char string) error message if any


 9  % .ListGlobalAttribute: list of global attributes containing:


 10  % .NbCoord: number of vector components


 11  % .NbDim: number of dimensions (=2 or 3)


 12  % .Dt: time interval for the corresponding image pair


 13  % .Time: absolute time (average of the initial image pair)


 14  % .CivStage: =0,


 15  % =1, civ1 has been performed only


 16  % =2, fix1 has been performed


 17  % =3, pacth1 has been performed


 18  % =4, civ2 has been performed


 19  % =5, fix2 has been performed


 20  % =6, pacth2 has been performed


 21  % .CoordUnit: 'pixel'


 22  % .ListVarName: {'X' 'Y' 'U' 'V' 'F' 'FF'}


 23  % .X, .Y, .Z: set of vector coordinates


 24  % .U,.V,.W: corresponding set of vector components


 25  % .F: warning flags


 26  % .FF: error flag, =0 for good vectors


 27  % .C: scalar associated with velocity (used for vector colors)


 28  %


 29  % VelTypeOut: velocity type corresponding to the selected field: ='civ1','interp1','interp2','civ2'....


 30  %


 31  % INPUT:


 32  % FileName: file name (string).


 33  % FieldNames =cell of field names to get, which can contain the strings:


 34  % 'ima_cor': image correlation, vec_c or vec2_C


 35  % 'vort','div','strain': requires velocity derivatives DUDX...


 36  % 'error': error estimate (vec_E or vec2_E)


 37  %


 38  % VelType : character string indicating the types of velocity fields to read ('civ1','civ2'...)


 39  % if vel_type=[] or'*', a priority choice, given by vel_type_out{1,2}, is done depending


 40  % if vel_type='filter'; a structured field is sought (filter2 in priority, then filter1)


 41  %


 42  % FUNCTIONS called:


 43  % 'varcivx_generator':, sets the names of vaiables to read in the netcdf file


 44  % 'nc2struct': reads a netcdf file


 45 


 46  %=======================================================================


[1027]  47  % Copyright 20082018, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


[958]  48  % http://www.legi.grenobleinp.fr


 49  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


 50  %


 51  % This file is part of the toolbox UVMAT.


 52  %


 53  % UVMAT is free software; you can redistribute it and/or modify


 54  % it under the terms of the GNU General Public License as published


 55  % by the Free Software Foundation; either version 2 of the license,


 56  % or (at your option) any later version.


 57  %


 58  % UVMAT is distributed in the hope that it will be useful,


 59  % but WITHOUT ANY WARRANTY; without even the implied warranty of


 60  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


 61  % GNU General Public License (see LICENSE.txt) for more details.


 62  %=======================================================================


 63 


 64  function [Field,VelTypeOut,errormsg]=read_pivdata_fluidimage(FileName,FieldNames,VelType)


 65 


 66  %% default input


 67  if ~exist('VelType','var')


 68  VelType='';


 69  end


 70  if isequal(VelType,'*')


 71  VelType='';


 72  end


 73  if isempty(VelType)


 74  VelType='';


 75  end


 76  if ~exist('FieldNames','var')


 77  FieldNames=[]; %default


 78  end


 79  Field=[];


 80  VelTypeOut=VelType;


 81  errormsg='';


 82  if ~exist(FileName,'file')


 83  errormsg=['input file ' FileName ' does not exist'];


 84  return


 85  end


 86  if ischar(FieldNames), FieldNames={FieldNames}; end;


 87  ProjModeRequest='';


 88  for ilist=1:length(FieldNames)


 89  if ~isempty(FieldNames{ilist})


 90  switch FieldNames{ilist}


[987]  91  case {'U','V','norm(U,V)'}


[958]  92  if ~strcmp(FieldNames{1},'vec(U,V)')% if the scalar is not used as color of vectors


 93  ProjModeRequest='interp_lin';


 94  end


 95  case {'curl(U,V)','div(U,V)','strain(U,V)'}


 96  ProjModeRequest='interp_tps';


 97  end


 98  end


 99  end


 100 


 101  %% reading data


 102 


 103 


 104  % Data=nc2struct(FileName,'ListGlobalAttribute','CivStage');


 105  % if isfield(Data,'Txt')


 106  % erromsg=['error in read_civdata: ' Data.Txt];


 107  % return


 108  % end


 109  % % set the list of variables to read and their role


 110  %[varlist,role,VelTypeOut]=varcivx_generator(ProjModeRequest,VelType,Datasets.CivStage);


 111  % if isempty(varlist)


 112  % erromsg=['error in read_civdata: unknow velocity type ' VelType];


 113  % return


 114  % else


 115  Datasets=hdf5load(FileName);%read the variables in the netcdf file


 116  %[varlist,role,VelTypeOut]=varcivx_generator(ProjModeRequest,VelType,Datasets.CivStage);


 117  role={'coord_x','coord_y','vector_x','vector_y','ancillary','warnflag','errorflag'};


 118  vardetect=ones(size(role));


 119  Field.ListGlobalAttribute= {'Dt','Time','CivStage'};


 120  Field.Dt=1;


 121  Field.Time=0;


 122  if isfield(Datasets, 'piv1')


 123  Field.CivStage=6;


 124  else


 125  Field.CivStage=3;


 126  end


 127  VelTypeOut=VelType;


 128  switch VelType


 129  case {'civ1','filter1'}


 130  Data=Datasets.piv0;


 131  VelTypeOut='filter1';


 132  case {'civ2','filter2'}


 133  Data=Datasets.piv1;


 134  VelTypeOut='filter2';


 135  case ''


 136  if isfield(Datasets, 'piv1')


 137  Data=Datasets.piv1;


 138  VelTypeOut='filter2';


 139  else


 140  Data=Datasets.piv0;


 141  VelTypeOut='filter1';


 142  end


 143  end


 144  switch VelType


 145  case {'civ1','civ2'}


 146  Field.X= double(Data.xs);


 147  Field.Y= double(Data.ys);


 148  Field.U= double(Data.deltaxs);


 149  Field.U(isnan(Field.U)) = 0;


 150  Field.V= double(Data.deltays);


 151  Field.V(isnan(Field.V)) = 0;


 152  case {'filter1','filter2',''}


[1003]  153  Field.X= double(Data.ixvecs_final);


 154  Field.Y= double(Data.iyvecs_final);


 155  Field.U= double(Data.deltaxs_final);


[958]  156  Field.U(isnan(Field.U)) = 0;


[1003]  157  Field.V= double(Data.deltays_final);


[958]  158  Field.V(isnan(Field.V)) = 0;


 159  end


 160  Field.ListVarName={'X' 'Y' 'U' 'V' 'C' 'F' 'FF'};


 161  Field.VarDimName={'nb_vec' 'nb_vec' 'nb_vec' 'nb_vec' 'nb_vec' 'nb_vec' 'nb_vec'};


 162  Field.C = double(Data.correls_max);


 163  Field.F=zeros(size(Field.U)); Field.F(Data.errors.keys + 1)=1; % !!! convention matlab vs python


 164  Field.FF=zeros(size(Field.U)); Field.FF(Data.errors.keys + 1)=1;


 165  % if vardetect(1)==0


 166  % errormsg=[ 'requested field not available in ' FileName '/' VelType ': need to run patch'];


 167  % return


 168  % end


 169  switch VelTypeOut


[987]  170  case {'civ1','filter1'}


[958]  171  if isfield(Field,'Patch1_SubDomain')


 172  Field.SubDomain=Field.Patch1_SubDomain;


 173  Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'SubDomain'}];


 174  end


 175  if isfield(Field,'Civ1_Dt')


 176  Field.Dt=Field.Civ1_Dt;


 177  end


 178  if isfield(Field,'Civ1_Time')


 179  Field.Time=Field.Civ1_Time;


 180  end


[987]  181  case {'civ2','filter2'}


[958]  182  if isfield(Field,'Patch2_SubDomain')


 183  Field.SubDomain=Field.Patch2_SubDomain;


 184  Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'SubDomain'}];


 185  end


 186  if isfield(Field,'Civ2_Dt')


 187  Field.Dt=Field.Civ2_Dt;


 188  end


 189  if isfield(Field,'Civ2_Time')


 190  Field.Time=Field.Civ2_Time;


 191  end


 192  end


 193  %Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'Dt','Time'}];


 194  ivar_U_tps=[];


 195  ivar_V_tps=[];


 196  var_ind=find(vardetect);


 197  for ivar=1:numel(var_ind)


 198  Field.VarAttribute{ivar}.Role=role{var_ind(ivar)};


 199  %Field.VarAttribute{ivar}.Mesh=0.025;%typical mesh for histograms O.025 pixel (used in series)


 200  Field.VarAttribute{ivar}.ProjModeRequest=ProjModeRequest;


 201  if strcmp(role{var_ind(ivar)},'vector_x')


 202  Field.VarAttribute{ivar}.FieldName=FieldNames;


 203  ivar_U=ivar;


 204  end


 205  if strcmp(role{var_ind(ivar)},'vector_x_tps')


 206  Field.VarAttribute{ivar}.FieldName=FieldNames;


 207  ivar_U_tps=ivar;


 208  end


 209  if strcmp(role{var_ind(ivar)},'vector_y')


 210  Field.VarAttribute{ivar}.FieldName=FieldNames;


 211  ivar_V=ivar;


 212  end


 213  if strcmp(role{var_ind(ivar)},'vector_y_tps')


 214  Field.VarAttribute{ivar}.FieldName=FieldNames;


 215  ivar_V_tps=ivar;


 216  end


 217  end


 218  if ~isempty(ivar_U_tps)


 219  Field.VarAttribute{ivar_U}.VarIndex_tps=ivar_U_tps;


 220  end


 221  if ~isempty(ivar_V_tps)


 222  Field.VarAttribute{ivar_V}.VarIndex_tps=ivar_V_tps;


 223  end


 224 


 225  %% update list of global attributes


 226  Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'NbCoord','NbDim','TimeUnit','CoordUnit'}];


 227  Field.NbCoord=2;


 228  Field.NbDim=2;


 229  Field.TimeUnit='s';


 230  Field.CoordUnit='pixel';

