source: trunk/src/calc_field.m @ 389

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

several bugs corrected: object managing, tps filter...
File size: 15.2 KB
Line 
1%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
2%---------------------------------------------------------------------
3% [DataOut,errormsg]=calc_field(FieldName,DataIn)
4%
5% OUTPUT:
6% Scal: matlab vector representing the scalar values (length nbvec defined by var_read)
7%      if no input, Scal=list of programmed scalar names (to put in menus)
8%      if only the field name is put as input, vec_A=type of scalar, which can be:
9%                   'discrete': related to the individual velocity vectors, not interpolated by patch
10%                   'vel': scalar calculated solely from velocity components
11%                   'der': needs spatial derivatives
12%                   'var': the scalar name directly corresponds to a field name in the netcdf files
13% error: error flag
14%      error = 0; OK
15%      error = 1; the prescribed scalar cannot be read or calculated from available fields
16%
17% INPUT:
18% FieldName: string representing the name of the field
19% DataIn: structure representing the field, as defined in check_field_srtructure.m
20%
21% FUNCTION related
22% varname_generator.m: determines the field names to read in the netcdf
23% file, depending on the scalar
24
25function [DataOut,errormsg]=calc_field(FieldName,DataIn,VelType,XI,YI)
26
27%list of defined scalars to display in menus (in addition to 'ima_cor').
28% a type is associated to each scalar:
29%              'discrete': related to the individual velocity vectors, not interpolated by patch
30%              'vel': calculated from velocity components, continuous field (interpolated with velocity)
31%              'der': needs spatial derivatives
32%              'var': the scalar name corresponds to a field name in the netcdf files
33% a specific variable name for civ1 and civ2 fields are also associated, if
34% the scalar is calculated from other fields, as explicited below
35
36list_field={'velocity';...%image correlation corresponding to a vel vector
37    'ima_cor';...%image correlation corresponding to a vel vector
38    'norm_vel';...%norm of the velocity
39    'vort';...%vorticity
40    'div';...%divergence
41    'strain';...%rate of strain
42    'u';... %u velocity component
43    'v';... %v velocity component
44    'w';... %w velocity component
45    'w_normal';... %w velocity component normal to the plane
46    'error'}; %error associated to a vector (for stereo or patch)
47errormsg=[]; %default error message
48if ~exist('FieldName','var')
49    DataOut=list_field;% gives the list of possible fields in the absence of input
50else
51    if ~exist('DataIn','var')
52        DataIn=[];
53    end
54    if ischar(FieldName)
55        FieldName={FieldName};%convert a string input to a cell with one string element
56    end
57    if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
58        nbcoord=3;
59    else
60        nbcoord=2;
61    end
62    ListVarName={};
63    ValueList={};
64    RoleList={};
65    units_cell={};
66   
67    %% interpolation with new civ data
68    if isfield(DataIn,'SubRange') &&(strcmp(VelType,'filter1')||strcmp(VelType,'filter2'))
69        XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
70        YMax=max(max(DataIn.SubRange(2,:,:)));
71        XMin=min(min(DataIn.SubRange(1,:,:)));
72        YMin=min(min(DataIn.SubRange(2,:,:)));
73        if ~(exist('XI','var')&&exist('YI','var'))
74            Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.Coord_tps(:,1,:)));%2D
75            xI=XMin:Mesh:XMax;
76            yI=YMin:Mesh:YMax;
77            [XI,YI]=meshgrid(xI,yI);% interpolation grid
78            XI=reshape(XI,[],1);
79            YI=reshape(YI,[],1);
80        end
81        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
82        for ilist=1:numel(DataOut.ListGlobalAttribute)
83            DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
84        end
85        DataOut.ListVarName={'coord_y','coord_x','FF'};
86        DataOut.VarDimName{1}='coord_y';
87        DataOut.VarDimName{2}='coord_x';
88        DataOut.coord_y=[yI(1) yI(end)];
89        DataOut.coord_x=[xI(1) xI(end)];
90        switch FieldName{1}
91            case {'velocity','u','v'}
92                DataOut.U=zeros(size(XI));
93                DataOut.V=zeros(size(XI));
94            case 'vort'
95                DataOut.vort=zeros(size(XI));
96            case 'div'
97                DataOut.div=zeros(size(XI));
98            case 'strain'
99                DataOut.strain=zeros(size(XI));
100        end
101        nbval=zeros(size(XI));
102        NbSubDomain=size(DataIn.SubRange,3);
103        for isub=1:NbSubDomain
104            if isfield(DataIn,'NbSites')
105                nbvec_sub=DataIn.NbSites(isub);
106            else
107                nbvec_sub=numel(find(DataIn.Indices_tps(:,isub)));
108            end
109            ind_sel=find(XI>=DataIn.SubRange(1,1,isub) & XI<=DataIn.SubRange(1,2,isub) & YI>=DataIn.SubRange(2,1,isub) & YI<=DataIn.SubRange(2,2,isub));
110            %rho smoothing parameter
111            epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites
112            ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
113            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
114            switch FieldName{1}
115                case {'velocity','u','v'}
116                    EM = tps_eval(epoints,ctrs);%kernels for calculating the velocity from tps 'sources'
117                case{'vort','div','strain'}
118                    [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%kernels for calculating the spatial derivatives from tps 'sources'
119            end
120            switch FieldName{1}
121                case 'velocity'
122                    ListFields={'U', 'V'};
123                    VarAttributes{1}.Role='vector_x';
124                    VarAttributes{2}.Role='vector_y';
125                    DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
126                    DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
127                case 'u'
128                    ListFields={'U'};
129                    VarAttributes{1}.Role='scalar';
130                    DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
131                case 'v'
132                    ListFields={'V'};
133                    VarAttributes{1}.Role='scalar';
134                    DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
135                case 'vort'
136                    ListFields={'vort'};
137                    VarAttributes{1}.Role='scalar';
138                    DataOut.vort(ind_sel)=DataOut.vort(ind_sel)+EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)-EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
139                case 'div'
140                    ListFields={'div'};
141                    VarAttributes{1}.Role='scalar';
142                    DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
143                case 'strain'
144                    ListFields={'strain'};
145                    VarAttributes{1}.Role='scalar';
146                    DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
147            end
148        end
149        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang     
150        DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
151        nbval(nbval==0)=1;
152        switch FieldName{1}
153              case {'velocity','u','v'}
154        DataOut.U=reshape(DataOut.U./nbval,numel(yI),numel(xI));
155        DataOut.V=reshape(DataOut.V./nbval,numel(yI),numel(xI));
156            case 'vort'
157        DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
158            case 'div'
159        DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
160            case 'strain'
161           DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI));
162        end
163        DataOut.ListVarName=[DataOut.ListVarName ListFields];
164        for ilist=3:numel(DataOut.ListVarName)
165            DataOut.VarDimName{ilist}={'coord_y','coord_x'};
166        end
167        DataOut.VarAttribute={[],[]};
168        DataOut.VarAttribute{3}.Role='errorflag';
169        DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
170    else
171       
172        %% civx data
173        DataOut=DataIn;
174        for ilist=1:length(FieldName)
175            if ~isempty(FieldName{ilist})
176                [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
177                ListVarName=[ListVarName VarName];
178                ValueList=[ValueList Value];
179                RoleList=[RoleList Role];
180                units_cell=[units_cell units];
181            end
182        end
183        %erase previous data (except coordinates)
184        for ivar=nbcoord+1:length(DataOut.ListVarName)
185            VarName=DataOut.ListVarName{ivar};
186            DataOut=rmfield(DataOut,VarName);
187        end
188        DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
189        if isfield(DataOut,'VarDimName')
190            DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
191        else
192            errormsg='element .VarDimName missing in input data';
193            return
194        end
195        DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
196        %append new data
197        DataOut.ListVarName=[DataOut.ListVarName ListVarName];
198        for ivar=1:length(ListVarName)
199            DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
200            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
201            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
202            eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
203        end
204    end
205end
206
207
208%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
209function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
210VarName={};
211ValCell={};
212Role={};
213units_cell={};
214if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
215    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
216else
217    units='pixel';
218end
219if isfield(DataIn,'U')
220    VarName=[VarName {'U'}];
221    ValCell=[ValCell {DataIn.U}];
222    Role=[Role {'vector_x'}];
223    units_cell=[units_cell {units}];
224end
225if isfield(DataIn,'V')
226    VarName=[VarName {'V'}];
227    ValCell=[ValCell {DataIn.V}];
228    Role=[Role {'vector_y'}];
229    units_cell=[units_cell {units}];
230end
231if isfield(DataIn,'W')
232    VarName=[VarName {'W'}];
233    ValCell=[ValCell {DataIn.W}];
234    Role=[Role {'vector_z'}];
235    units_cell=[units_cell {units}];
236end
237if isfield(DataIn,'F')
238    VarName=[VarName {'F'}];
239    ValCell=[ValCell {DataIn.F}];
240    Role=[Role {'warnflag'}];
241    units_cell=[units_cell {[]}];
242end
243if isfield(DataIn,'FF')
244    VarName=[VarName,{'FF'}];
245    ValCell=[ValCell {DataIn.FF}];
246    Role=[Role {'errorflag'}];
247    units_cell=[units_cell {[]}];
248end
249
250%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
251function [VarName,ValCell,Role,units]=ima_cor(DataIn)
252VarName={};
253ValCell={};
254Role={};
255units={};
256if isfield(DataIn,'C')
257    VarName{1}='C';
258    ValCell{1}=DataIn.C;
259    Role={'ancillary'};
260    units={[]};
261end
262
263%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
264function [VarName,ValCell,Role,units]=norm_vel(DataIn)
265VarName={};
266ValCell={};
267Role={};
268units={};
269if isfield(DataIn,'U') && isfield(DataIn,'V')
270    VarName{1}='norm_vel';
271    ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
272    if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
273        ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
274    end
275    ValCell{1}=sqrt(ValCell{1});
276    Role{1}='scalar';
277    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
278        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
279    else
280        units={'pixel'};
281    end
282end
283
284
285
286%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
287function [VarName,ValCell,Role,units]=vort(DataIn)
288VarName={};
289ValCell={};
290Role={};
291units={};
292if isfield(DataIn,'DjUi')
293    VarName{1}='vort';
294    ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
295    siz=size(ValCell{1});
296    ValCell{1}=reshape(ValCell{1},siz(1),1);
297    Role{1}='scalar';
298    if isfield(DataIn,'TimeUnit')
299        units={[DataIn.TimeUnit '-1']};
300    else
301        units={[]};
302    end
303end
304
305%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
306function [VarName,ValCell,Role,units]=div(DataIn)
307VarName={};
308ValCell={};
309Role={};
310units={};
311if isfield(DataIn,'DjUi')
312    VarName{1}='div';
313    ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
314    siz=size(ValCell{1});
315    ValCell{1}=reshape(ValCell{1},siz(1),1);
316    Role{1}='scalar';
317    if isfield(DataIn,'TimeUnit')
318        units={[DataIn.TimeUnit '-1']};
319    else
320        units={[]};
321    end
322end
323
324%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
325function [VarName,ValCell,Role,units]=strain(DataIn)
326VarName={};
327ValCell={};
328Role={};
329units={};
330if isfield(DataIn,'DjUi')
331    VarName{1}='strain';
332    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
333    siz=size(ValCell{1});
334    ValCell{1}=reshape(ValCell{1},siz(1),1);
335    Role{1}='scalar';
336    if isfield(DataIn,'TimeUnit')
337        units={[DataIn.TimeUnit '-1']};
338    else
339        units={[]};
340    end
341end
342
343%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
344function [VarName,ValCell,Role,units]=u(DataIn)
345VarName={};
346ValCell={};
347Role={};
348units={};
349if isfield(DataIn,'U')
350    VarName{1}='U';
351    ValCell{1}=DataIn.U;
352    Role{1}='scalar';
353    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
354        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
355    else
356        units={'pixel'};
357    end
358end
359
360%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
361function [VarName,ValCell,Role,units]=v(DataIn)
362VarName={};
363ValCell={};
364Role={};
365units={};
366if isfield(DataIn,'V')
367    VarName{1}='V';
368    ValCell{1}=DataIn.V;
369    Role{1}='scalar';
370    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
371        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
372    else
373        units={'pixel'};
374    end
375end
376
377%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
378function [VarName,ValCell,Role,units]=w(DataIn)
379VarName={};
380ValCell={};
381Role={};
382units={};
383if isfield(DataIn,'W')
384    VarName{1}='W';
385    ValCell{1}=DataIn.W;
386    Role{1}='scalar';%will remain unchanged by projection
387    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
388        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
389    else
390        units={'pixel'};
391    end
392end
393
394%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
395function [VarName,ValCell,Role,units]=w_normal(DataIn)
396VarName={};
397ValCell={};
398Role={};
399units={};
400if isfield(DataIn,'W')
401    VarName{1}='W';
402    ValCell{1}=DataIn.W;
403    Role{1}='vector_z';%will behave like a vector component  by projection
404    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
405        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
406    else
407        units={'pixel'};
408    end
409end
410
411%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
412function [VarName,ValCell,Role,units]=error(DataIn)
413VarName={};
414ValCell={};
415Role={};
416units={};
417if isfield(DataIn,'E')
418    VarName{1}='E';
419    ValCell{1}=DataIn.E;
420    Role{1}='ancillary'; %TODO CHECK units in actual fields
421    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
422        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
423    else
424        units={'pixel'};
425    end
426end
427
Note: See TracBrowser for help on using the repository browser.