source: trunk/src/read_civdata.m @ 609

Last change on this file since 609 was 594, checked in by sommeria, 12 years ago

bugs corrected. One step further for using compiled fcts and cluster with series

File size: 9.5 KB
Line 
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
43% FUNCTIONS called:
44% 'varcivx_generator':, sets the names of vaiables to read in the netcdf file
45% 'nc2struct': reads a netcdf file
46
47function [Field,VelTypeOut,errormsg]=read_civdata(filename,FieldNames,VelType)
48
49%% default input
50if ~exist('VelType','var')
51    VelType='';
52end
53if isequal(VelType,'*')
54    VelType='';
55end
56if isempty(VelType)
57    VelType='';
58end
59if ~exist('FieldNames','var')
60    FieldNames=[]; %default
61end
62errormsg='';
63if ischar(FieldNames), FieldNames={FieldNames}; end;
64ProjModeRequest='';
65for ilist=1:length(FieldNames)
66    if ~isempty(FieldNames{ilist})
67        switch FieldNames{ilist}
68            case{'U','V','norm(U,V)'}
69                ProjModeRequest='interp_lin';
70            case {'curl(U,V)','div(U,V)','strain(U,V)'}
71                ProjModeRequest='interp_tps';
72        end
73    end
74end
75
76%% reading data
77Data=nc2struct(filename,'ListGlobalAttribute','CivStage');
78[varlist,role,VelTypeOut]=varcivx_generator(ProjModeRequest,VelType,Data.CivStage);
79if isempty(varlist)
80    erromsg=['error in read_civdata: unknow velocity type ' VelType];
81    return
82else
83    [Field,vardetect]=nc2struct(filename,varlist);%read the variables in the netcdf file
84end
85if isfield(Field,'Txt')
86    errormsg=Field.Txt;
87    return
88end
89if vardetect(1)==0
90     errormsg=[ 'requested field not available in ' filename '/' VelType ': need to run patch'];
91     return
92end
93switch VelTypeOut
94    case{'civ1','filter1'}
95        if isfield(Field,'Patch1_SubDomain')
96            Field.SubDomain=Field.Patch1_SubDomain;
97            Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'SubDomain'}];
98        end
99        if isfield(Field,'Civ1_Dt')
100        Field.Dt=Field.Civ1_Dt;
101        end
102        if isfield(Field,'Civ1_Time')
103        Field.Time=Field.Civ1_Time;
104        end
105    case{'civ2','filter2'}
106        if isfield(Field,'Patch2_SubDomain')
107            Field.SubDomain=Field.Patch2_SubDomain;
108            Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'SubDomain'}];
109        end
110        if isfield(Field,'Civ2_Dt')
111        Field.Dt=Field.Civ2_Dt;
112        end
113        if isfield(Field,'Civ2_Time')
114        Field.Time=Field.Civ2_Time;
115        end
116end
117Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'Dt','Time'}];
118ivar_U_tps=[];
119ivar_V_tps=[];
120var_ind=find(vardetect);
121for ivar=1:numel(var_ind)
122    Field.VarAttribute{ivar}.Role=role{var_ind(ivar)};
123    if strcmp(role{var_ind(ivar)},'vector_x')
124        Field.VarAttribute{ivar}.FieldName=FieldNames;
125        Field.VarAttribute{ivar}.ProjModeRequest=ProjModeRequest;
126        ivar_U=ivar;
127    end
128    if strcmp(role{var_ind(ivar)},'vector_x_tps')
129        Field.VarAttribute{ivar}.FieldName=FieldNames;
130        Field.VarAttribute{ivar}.ProjModeRequest=ProjModeRequest;
131        ivar_U_tps=ivar;
132    end
133    if strcmp(role{var_ind(ivar)},'vector_y')
134        Field.VarAttribute{ivar}.FieldName=FieldNames;
135        Field.VarAttribute{ivar}.ProjModeRequest=ProjModeRequest;
136        ivar_V=ivar;
137    end
138    if strcmp(role{var_ind(ivar)},'vector_y_tps')
139        Field.VarAttribute{ivar}.FieldName=FieldNames;
140        Field.VarAttribute{ivar}.ProjModeRequest=ProjModeRequest;
141        ivar_V_tps=ivar;
142    end
143    Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
144end
145if ~isempty(ivar_U_tps)
146    Field.VarAttribute{ivar_U}.VarIndex_tps=ivar_U_tps;
147end
148if ~isempty(ivar_V_tps)
149    Field.VarAttribute{ivar_V}.VarIndex_tps=ivar_V_tps;
150end
151
152%% update list of global attributes
153Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'NbCoord','NbDim','TimeUnit','CoordUnit'}];
154Field.NbCoord=2;
155Field.NbDim=2;
156Field.TimeUnit='s';
157Field.CoordUnit='pixel';
158
159%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
160% [var,role,units,vel_type_out]=varcivx_generator(FieldNames,vel_type)
161%INPUT:
162% FieldNames =cell of field names to get, which can contain the strings:
163%             'ima_cor': image correlation, vec_c or vec2_C
164%             'vort','div','strain': requires velocity derivatives DUDX...
165%             'error': error estimate (vec_E or vec2_E)
166%             
167% vel_type: character string indicating the types of velocity fields to read ('civ1','civ2'...)
168%            if vel_type=[] or'*', a  priority choice is done, civ2 considered better than civ1 )
169
170function [var,role,vel_type_out,errormsg]=varcivx_generator(ProjModeRequest,vel_type,CivStage)
171
172%% default input values
173if ~exist('vel_type','var'),vel_type='';end;
174if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
175errormsg='';
176
177%% select the priority order for automatic vel_type selection
178if strcmp(vel_type,'civ2') && strcmp(ProjModeRequest,'derivatives')
179    vel_type='filter2';
180elseif strcmp(vel_type,'civ1') && strcmp(ProjModeRequest,'derivatives')
181    vel_type='filter1';
182end
183if isempty(vel_type)||strcmp(vel_type,'*')
184    switch CivStage
185        case {6} %filter2 available
186            vel_type='filter2';
187        case {4,5}% civ2 available but not filter2
188            if strcmp(ProjModeRequest,'derivatives')% derivatives needed
189                vel_type='filter1';
190            else
191                vel_type='civ2';
192            end
193        case 3
194            vel_type='filter1';
195        case {1,2}% civ1 available but not filter1
196            vel_type='civ1';
197    end
198end
199
200var={};
201switch vel_type
202    case 'civ1'
203        var={'X','Y','Z','U','V','W','C','F','FF';...
204            'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
205        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
206    %    units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
207    case 'filter1'
208        var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbCentres','NbCentres';...
209            'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U_smooth','Civ1_V_smooth','Civ1_W','Civ1_C','Civ1_F','Civ1_FF',...
210            'Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps','Civ1_W_tps','Civ1_SubRange','Civ1_NbCentres','Civ1_NbSites'};
211        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',...
212            'vector_y_tps','vector_z_tps','ancillary','ancillary','ancillary'};
213     %  rmq: NbSites obsolete replaced by NbCentres, kept for consistency with previous data
214    case 'civ2'
215        var={'X','Y','Z','U','V','W','C','F','FF';...
216            'Civ2_X','Civ2_Y','Civ2_Z','Civ2_U','Civ2_V','Civ2_W','Civ2_C','Civ2_F','Civ2_FF'};
217        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
218    case 'filter2'
219        var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbCentres','NbCentres';...
220            'Civ2_X','Civ2_Y','Civ2_Z','Civ2_U_smooth','Civ2_V_smooth','Civ2_W','Civ2_C','Civ2_F','Civ2_FF',...
221            'Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps','','Civ2_SubRange','Civ2_NbCentres','Civ2_NbSites'};
222        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',...
223            'vector_y_tps','vector_z_tps','ancillary','ancillary','ancillary'};
224end
225if ~strcmp(ProjModeRequest,'interp_tps')
226    var=var(:,1:9);%suppress tps if not needed
227end
228vel_type_out=vel_type;
229
230
231
232
233
Note: See TracBrowser for help on using the repository browser.