source: trunk/src/read_civxdata.m @ 879

Last change on this file since 879 was 809, checked in by g7moreau, 10 years ago
  • Add license
File size: 12.6 KB
Line 
1%'read_civxdata': reads civx data from netcdf files
2%------------------------------------------------------------------
3%
4% function [Field,VelTypeOut]=read_civxdata(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%                    .DijU; matrix of spatial derivatives (DijU(1,1,:)=DUDX,
29%                        DijU(1,2,:)=DUDY, Dij(2,1,:)=DVDX, DijU(2,2,:)=DVDY
30%
31% VelTypeOut: velocity type corresponding to the selected field: ='civ1','interp1','interp2','civ2'....
32%
33% INPUT:
34% filename: file name (string).
35% FieldNames =cell of field names to get, which can contain the strings:
36%             'ima_cor': image correlation, vec_c or vec2_C
37%             'vort','div','strain': requires velocity derivatives DUDX...
38%             'error': error estimate (vec_E or vec2_E)
39%             
40% VelType : character string indicating the types of velocity fields to read ('civ1','civ2'...)
41%            if vel_type=[] or'*', a  priority choice, given by vel_type_out{1,2}, is done depending
42%            if vel_type='filter'; a structured field is sought (filter2 in priority, then filter1)
43%
44% FUNCTIONS called:
45% 'varcivx_generator':, sets the names of vaiables to read in the netcdf file
46% 'nc2struct': reads a netcdf file
47
48%=======================================================================
49% Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France
50%   http://www.legi.grenoble-inp.fr
51%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
52%
53%     This file is part of the toolbox UVMAT.
54%
55%     UVMAT is free software; you can redistribute it and/or modify
56%     it under the terms of the GNU General Public License as published
57%     by the Free Software Foundation; either version 2 of the license,
58%     or (at your option) any later version.
59%
60%     UVMAT is distributed in the hope that it will be useful,
61%     but WITHOUT ANY WARRANTY; without even the implied warranty of
62%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
63%     GNU General Public License (see LICENSE.txt) for more details.
64%=======================================================================
65
66function [Field,VelTypeOut,errormsg]=read_civxdata(filename,FieldNames,VelType)
67errormsg='';
68Field=[];
69VelTypeOut=[];
70DataTest=nc2struct(filename,'ListGlobalAttribute','Conventions','CivStage');
71if isfield(DataTest,'Txt')
72    errormsg=['nc2struct / ' DataTest.Txt];
73    return
74elseif isequal(DataTest.Conventions,'uvmat/civdata')%test for new civ format
75     [Field,VelTypeOut,errormsg]=read_civdata(filename,FieldNames,VelType,DataTest.CivStage);
76      if ~isempty(errormsg),errormsg=['read_civdata / ' errormsg];end
77     return
78end
79   
80%% default input
81if ~exist('VelType','var')
82    VelType=[];
83end
84if isequal(VelType,'*')
85    VelType=[];
86end
87if ~exist('FieldNames','var')
88    FieldNames=[]; %default
89end
90
91%% reading data
92VelTypeOut=VelType;%default
93[var,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);%determine the names of constants and variables to read
94[Field,vardetect,ichoice]=nc2struct(filename,var);%read the variables in the netcdf file
95if isfield(Field,'Txt')
96    errormsg=['nc2struct / ' Field.Txt];
97    return
98end
99if vardetect(1)==0
100     errormsg=[ 'requested field not available in ' filename '/' VelType];
101     return
102end
103var_ind=find(vardetect);
104for ilist=1:length(FieldNames)
105        testinterp=isempty(regexp(FieldNames{ilist},'(^vec|^C)', 'once'));%test need for gridded data
106        if testinterp, break;end
107end   
108for ivar=1:length(var_ind)
109    Field.VarAttribute{ivar}.Role=role{var_ind(ivar)};
110    Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
111    if strcmp(role{var_ind(ivar)},'vector_x')
112        Field.VarAttribute{ivar}.FieldName=FieldNames;
113        if testinterp
114            Field.VarAttribute{ivar}.ProjModeRequest='interp_lin';
115        end
116    end
117end
118VelTypeOut=VelType;
119if ~isempty(ichoice)
120    VelTypeOut=vel_type_out_cell{ichoice};
121end
122
123%% adjust for Djui:
124if isfield(Field,'DjUi')
125    Field.ListVarName{end-3}='DjUi';
126    Field.VarDimName{end-3}=[Field.VarDimName{end-3} {'nb_coord'} {'nb_coord'}];
127    Field.ListVarName(end-2:end)=[];
128    Field.VarDimName(end-2:end)=[];
129    Field.VarAttribute(end-2:end)=[];
130end
131
132%% renaming for standard conventions
133Field.NbCoord=Field.nb_coord;
134Field.NbDim=Field.nb_dim;
135
136%% CivStage
137if isfield(Field,'patch2')&& isequal(Field.patch2,1)
138    Field.CivStage=6;
139elseif isfield(Field,'fix2')&& isequal(Field.fix2,1)
140    Field.CivStage=5;
141elseif isfield(Field,'civ2')&& isequal(Field.civ2,1)
142    Field.CivStage=4;
143elseif isfield(Field,'patch')&& isequal(Field.patch,1)
144    Field.CivStage=3;
145elseif isfield(Field,'fix')&& isequal(Field.fix,1)
146    Field.CivStage=2;
147else
148    Field.CivStage=1;
149end
150
151%determine the appropriate constant for time and dt for the PIV pair
152test_civ1=isequal(VelTypeOut,'civ1')||isequal(VelTypeOut,'interp1')||isequal(VelTypeOut,'filter1');
153test_civ2=isequal(VelTypeOut,'civ2')||isequal(VelTypeOut,'interp2')||isequal(VelTypeOut,'filter2');
154Field.Time=0; %default
155Field.TimeUnit='s';
156if test_civ1
157    if isfield(Field,'absolut_time_T0')
158        Field.Time=double(Field.absolut_time_T0);
159        Field.Dt=double(Field.dt);
160    else
161       errormsg='the input file is not civx';
162       Field.CivStage=0;
163       Field.Dt=0;
164    end
165elseif test_civ2
166    Field.Time=double(Field.absolut_time_T0_2);
167    Field.Dt=double(Field.dt2);
168else
169    errormsg='the input file is not civx';
170    Field.CivStage=0;
171    Field.Dt=0;
172end
173
174%% rescale fields to pixel coordinates
175if isfield(Field,'pixcmx')
176    Field.pixcmx=double(Field.pixcmx);
177    Field.pixcmy=double(Field.pixcmy);
178    Field.U=Field.U*Field.pixcmx;
179    Field.V=Field.V*Field.pixcmy;
180    Field.X=Field.X*Field.pixcmx;
181    Field.Y=Field.Y*Field.pixcmy;
182end
183if ~isequal(Field.Dt,0)
184    Field.U=Field.U*Field.Dt;%translate in px displacement
185    Field.V=Field.V*Field.Dt;
186    if isfield(Field,'DjUi')
187       Field.DjUi(:,1,1)=Field.Dt*Field.DjUi(:,1,1);
188       Field.DjUi(:,2,2)=Field.Dt*Field.DjUi(:,2,2);
189       Field.DjUi(:,1,2)=(Field.pixcmy/Field.pixcmx)*Field.Dt*Field.DjUi(:,1,2);
190       Field.DjUi(:,2,1)=(Field.pixcmx/Field.pixcmy)*Field.Dt*Field.DjUi(:,2,1);
191    end
192end
193
194%% update list of global attributes
195List=Field.ListGlobalAttribute;
196ind_remove=[];
197for ilist=1:length(List)
198    switch(List{ilist})
199        case {'patch2','fix2','civ2','patch','fix','dt','dt2','absolut_time_T0','absolut_time_T0_2','nb_coord','nb_dim','pixcmx','pixcmy'}
200            ind_remove=[ind_remove ilist];
201            Field=rmfield(Field,List{ilist});
202    end
203end
204List(ind_remove)=[];
205Field.ListGlobalAttribute=[{'NbCoord'},{'NbDim'} List {'Time','Dt','TimeUnit','CivStage','CoordUnit'}];
206Field.CoordUnit='pixel';
207
208
209%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
210% [var,role,units,vel_type_out]=varcivx_generator(FieldNames,vel_type)
211%INPUT:
212% FieldNames =cell of field names to get, which can contain the strings:
213%             'C': image correlation, vec_c or vec2_C
214%             'curl','div','strain': requires velocity derivatives DUDX...
215%             'error': error estimate (vec_E or vec2_E)
216%             
217% vel_type: character string indicating the types of velocity fields to read ('civ1','civ2'...)
218%            if vel_type=[] or'*', a  priority choice, given by vel_type_out{1,2}, is done depending
219%            if vel_type='filter'; a structured field is sought (filter2 in priority, then filter1)
220
221function [var,role,units,vel_type_out]=varcivx_generator(FieldNames,vel_type)
222
223%% default input values
224if ~exist('vel_type','var'),vel_type=[];end;
225if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
226if ~exist('FieldNames','var'),FieldNames={'C'};end;%default scalar
227if ischar(FieldNames), FieldNames={FieldNames}; end;
228
229%% select the priority order for automatic vel_type selection
230testder=0;
231for ilist=1:length(FieldNames)
232        testder=~isempty(regexp(FieldNames{ilist},'(^curl|^div|^strain)', 'once'));%test need for derivatives
233        if testder, break;end
234end     
235if isempty(vel_type) || isequal(vel_type,'*') %undefined velocity type (civ1,civ2...)
236    if testder
237         vel_type_out{1}='filter2'; %priority to filter2 for scalar reading, filter1 as second
238        vel_type_out{2}='filter1';
239    else
240        vel_type_out{1}='civ2'; %priority to civ2 for vector reading, civ1 as second priority     
241        vel_type_out{2}='civ1';
242    end
243elseif isequal(vel_type,'filter')
244        vel_type_out{1}='filter2'; %priority to filter2 for scalar reading, filter1 as second
245        vel_type_out{2}='filter1';
246        if ~testder
247            vel_type_out{3}='civ1';%civ1 as third priority if derivatives are not needed
248        end
249elseif testder
250    test_civ1=isequal(vel_type,'civ1')||isequal(vel_type,'interp1')||isequal(vel_type,'filter1');
251    if test_civ1
252        vel_type_out{1}='filter1'; %switch to filter for reading spatial derivatives
253    else
254        vel_type_out{1}='filter2';
255    end
256else   
257    vel_type_out{1}=vel_type;%imposed velocity field
258end
259vel_type_out=vel_type_out';
260
261%% determine names of netcdf variables to read
262var={'X','Y','Z','U','V','W','C','F','FF'};
263role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
264units={'pixel','pixel','pixel','pixel','pixel','pixel',[],[],[]};
265if testder
266    var=[var {'DjUi(:,1,1)','DjUi(:,1,2)','DjUi(:,2,1)','DjUi(:,2,2)'}];
267    role=[role {'tensor','tensor','tensor','tensor'}];
268    units=[units {'pixel','pixel','pixel','pixel'}];
269end
270for ilist=1:length(vel_type_out)
271    var=[var;varname1(vel_type_out{ilist},FieldNames)];
272end
273
274%------------------------------------------------------------------------ 
275%--- determine  var names to read
276function varin=varname1(vel_type,FieldNames)
277%------------------------------------------------------------------------
278testder=0;
279C1='';
280C2='';
281for ilist=1:length(FieldNames)
282    if ~isempty(FieldNames{ilist})
283    switch FieldNames{ilist}
284        case 'C' %image correlation corresponding to a vel vector
285            C1='vec_C';
286            C2='vec2_C';
287        case 'error'
288            C1='vec_E';
289            C2='vec2_E';
290        otherwise
291          testder=~isempty(regexp(FieldNames{ilist},'(^curl|^div|strain)', 'once'));%test need for derivatives
292    end
293    end
294end     
295switch vel_type
296    case 'civ1'
297        varin={'vec_X','vec_Y','vec_Z','vec_U','vec_V','vec_W',C1,'vec_F','vec_FixFlag'};
298    case 'interp1'
299        varin={'vec_patch_X','vec_patch_Y','','vec_patch0_U','vec_patch0_V','','','',''};
300    case 'filter1'
301        varin={'vec_patch_X','vec_patch_Y','','vec_patch_U','vec_patch_V','','','',''};
302    case 'civ2'
303        varin={'vec2_X','vec2_Y','vec2_Z','vec2_U','vec2_V','vec2_W',C2,'vec2_F','vec2_FixFlag'};
304    case 'interp2'
305        varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch0_U','vec2_patch0_V','vec2_patch0_W','','',''};
306    case 'filter2'
307        varin={'vec2_patch_X','vec2_patch_Y','vec2_patch_Z','vec2_patch_U','vec2_patch_V','vec2_patch0_W','','',''};
308end
309if testder
310     switch vel_type
311        case 'filter1'
312            varin=[varin {'vec_patch_DUDX','vec_patch_DVDX','vec_patch_DUDY','vec_patch_DVDY'}];
313        case 'filter2'
314            varin=[varin {'vec2_patch_DUDX','vec2_patch_DVDX','vec2_patch_DUDY','vec2_patch_DVDY'}];
315    end   
316end
Note: See TracBrowser for help on using the repository browser.