source: trunk/src/calc_field.m @ 397

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

civ_matlab and patch improved, changes in the management of interpolation (still in progress).
adapatation to movies (use of VideoReader?)
File size: 15.9 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,Coord_interp)
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') && isfield(DataIn,'Coord_tps')&& exist('Coord_interp','var')
69        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
70        for ilist=1:numel(DataOut.ListGlobalAttribute)
71            DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
72        end
73        DataOut.ListVarName={'coord_y','coord_x','FF'};
74        DataOut.VarDimName{1}='coord_y';
75        DataOut.VarDimName{2}='coord_x';
76        XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
77        YMax=max(max(DataIn.SubRange(2,:,:)));
78        XMin=min(min(DataIn.SubRange(1,:,:)));
79        YMin=min(min(DataIn.SubRange(2,:,:)));
80        %         if ~exist('Coord_interp','var')
81        %             Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.Coord_tps(:,1,:)));%2D
82        %             xI=XMin:Mesh:XMax;
83        %             yI=YMin:Mesh:YMax;
84        %             [XI,YI]=meshgrid(xI,yI);% interpolation grid
85        %             Coord_interp(:,1)=reshape(XI,[],1);
86        %             Coord_interp(:,2)=reshape(YI,[],1);
87        %             DataOut.coord_y=[yI(1) yI(end)];
88        %             DataOut.coord_x=[xI(1) xI(end)];
89        %      end
90        check_der=0;
91        check_val=0;
92        for ilist=1:length(FieldName)
93            switch FieldName{ilist}
94                case 'velocity'
95                    check_val=1;
96                    DataOut.U=zeros(size(Coord_interp,1));
97                    DataOut.V=zeros(size(Coord_interp,1));
98                case{'vort','div','strain'}% case of spatial derivatives
99                    check_der=1;
100                    DataOut.(FieldName{ilist})=zeros(size(Coord_interp,1));
101                otherwise % case of a scalar
102                    check_val=1;
103                    DataOut.(FieldName{ilist})=zeros(size(Coord_interp,1));
104            end
105        end
106        nbval=zeros(size(Coord_interp,1));
107        NbSubDomain=size(DataIn.SubRange,3);
108        for isub=1:NbSubDomain
109            if isfield(DataIn,'NbSites')
110                nbvec_sub=DataIn.NbSites(isub);
111            else
112                nbvec_sub=numel(find(DataIn.Indices_tps(:,isub)));
113            end
114            ind_sel=find(Coord_interp>=DataIn.SubRange(:,1,isub) & Coord_interp<=DataIn.SubRange(:,2,isub));
115            %rho smoothing parameter
116            %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
117            %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
118            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
119            if check_val
120                EM = tps_eval(Coord_interp(ind_sel),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
121            end
122            if check_der
123                [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
124            end
125            for ilist=1:length(FieldName)
126                switch FieldName{ilist}
127                    case 'velocity'
128                        ListFields={'U', 'V'};
129                        VarAttributes{1}.Role='vector_x';
130                        VarAttributes{2}.Role='vector_y';
131                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
132                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
133                    case 'u'
134                        ListFields={'U'};
135                        VarAttributes{1}.Role='scalar';
136                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
137                    case 'v'
138                        ListFields={'V'};
139                        VarAttributes{1}.Role='scalar';
140                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
141                    case 'vort'
142                        ListFields={'vort'};
143                        VarAttributes{1}.Role='scalar';
144                        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);
145                    case 'div'
146                        ListFields={'div'};
147                        VarAttributes{1}.Role='scalar';
148                        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);
149                    case 'strain'
150                        ListFields={'strain'};
151                        VarAttributes{1}.Role='scalar';
152                        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);
153                end
154            end
155            DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
156            DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
157            nbval(nbval==0)=1;
158            switch FieldName{1}
159                case {'velocity','u','v'}
160                    DataOut.U=reshape(DataOut.U./nbval,numel(yI),numel(xI));
161                    DataOut.V=reshape(DataOut.V./nbval,numel(yI),numel(xI));
162                case 'vort'
163                    DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
164                case 'div'
165                    DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
166                case 'strain'
167                    DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI));
168            end
169            DataOut.ListVarName=[DataOut.ListVarName ListFields];
170            for ilist=3:numel(DataOut.ListVarName)
171                DataOut.VarDimName{ilist}={'coord_y','coord_x'};
172            end
173            DataOut.VarAttribute={[],[]};
174            DataOut.VarAttribute{3}.Role='errorflag';
175            DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
176        end
177    else
178       
179        %% civx data
180        DataOut=DataIn;
181        for ilist=1:length(FieldName)
182            if ~isempty(FieldName{ilist})
183                [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
184                ListVarName=[ListVarName VarName];
185                ValueList=[ValueList Value];
186                RoleList=[RoleList Role];
187                units_cell=[units_cell units];
188            end
189        end
190        %erase previous data (except coordinates)
191        for ivar=nbcoord+1:length(DataOut.ListVarName)
192            VarName=DataOut.ListVarName{ivar};
193            DataOut=rmfield(DataOut,VarName);
194        end
195        DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
196        if isfield(DataOut,'VarDimName')
197            DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
198        else
199            errormsg='element .VarDimName missing in input data';
200            return
201        end
202        DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
203        %append new data
204        DataOut.ListVarName=[DataOut.ListVarName ListVarName];
205        for ivar=1:length(ListVarName)
206            DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
207            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
208            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
209            DataOut.(ListVarName{ivar})=ValueList{ivar};
210        end
211    end
212end
213
214
215%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
216function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
217VarName={};
218ValCell={};
219Role={};
220units_cell={};
221if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
222    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
223else
224    units='pixel';
225end
226if isfield(DataIn,'U')
227    VarName=[VarName {'U'}];
228    ValCell=[ValCell {DataIn.U}];
229    Role=[Role {'vector_x'}];
230    units_cell=[units_cell {units}];
231end
232if isfield(DataIn,'V')
233    VarName=[VarName {'V'}];
234    ValCell=[ValCell {DataIn.V}];
235    Role=[Role {'vector_y'}];
236    units_cell=[units_cell {units}];
237end
238if isfield(DataIn,'W')
239    VarName=[VarName {'W'}];
240    ValCell=[ValCell {DataIn.W}];
241    Role=[Role {'vector_z'}];
242    units_cell=[units_cell {units}];
243end
244if isfield(DataIn,'F')
245    VarName=[VarName {'F'}];
246    ValCell=[ValCell {DataIn.F}];
247    Role=[Role {'warnflag'}];
248    units_cell=[units_cell {[]}];
249end
250if isfield(DataIn,'FF')
251    VarName=[VarName,{'FF'}];
252    ValCell=[ValCell {DataIn.FF}];
253    Role=[Role {'errorflag'}];
254    units_cell=[units_cell {[]}];
255end
256
257%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
258function [VarName,ValCell,Role,units]=ima_cor(DataIn)
259VarName={};
260ValCell={};
261Role={};
262units={};
263if isfield(DataIn,'C')
264    VarName{1}='C';
265    ValCell{1}=DataIn.C;
266    Role={'ancillary'};
267    units={[]};
268end
269
270%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
271function [VarName,ValCell,Role,units]=norm_vel(DataIn)
272VarName={};
273ValCell={};
274Role={};
275units={};
276if isfield(DataIn,'U') && isfield(DataIn,'V')
277    VarName{1}='norm_vel';
278    ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
279    if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
280        ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
281    end
282    ValCell{1}=sqrt(ValCell{1});
283    Role{1}='scalar';
284    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
285        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
286    else
287        units={'pixel'};
288    end
289end
290
291
292
293%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
294function [VarName,ValCell,Role,units]=vort(DataIn)
295VarName={};
296ValCell={};
297Role={};
298units={};
299if isfield(DataIn,'DjUi')
300    VarName{1}='vort';
301    ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
302    siz=size(ValCell{1});
303    ValCell{1}=reshape(ValCell{1},siz(1),1);
304    Role{1}='scalar';
305    if isfield(DataIn,'TimeUnit')
306        units={[DataIn.TimeUnit '-1']};
307    else
308        units={[]};
309    end
310end
311
312%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
313function [VarName,ValCell,Role,units]=div(DataIn)
314VarName={};
315ValCell={};
316Role={};
317units={};
318if isfield(DataIn,'DjUi')
319    VarName{1}='div';
320    ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
321    siz=size(ValCell{1});
322    ValCell{1}=reshape(ValCell{1},siz(1),1);
323    Role{1}='scalar';
324    if isfield(DataIn,'TimeUnit')
325        units={[DataIn.TimeUnit '-1']};
326    else
327        units={[]};
328    end
329end
330
331%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
332function [VarName,ValCell,Role,units]=strain(DataIn)
333VarName={};
334ValCell={};
335Role={};
336units={};
337if isfield(DataIn,'DjUi')
338    VarName{1}='strain';
339    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
340    siz=size(ValCell{1});
341    ValCell{1}=reshape(ValCell{1},siz(1),1);
342    Role{1}='scalar';
343    if isfield(DataIn,'TimeUnit')
344        units={[DataIn.TimeUnit '-1']};
345    else
346        units={[]};
347    end
348end
349
350%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
351function [VarName,ValCell,Role,units]=u(DataIn)
352VarName={};
353ValCell={};
354Role={};
355units={};
356if isfield(DataIn,'U')
357    VarName{1}='U';
358    ValCell{1}=DataIn.U;
359    Role{1}='scalar';
360    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
361        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
362    else
363        units={'pixel'};
364    end
365end
366
367%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
368function [VarName,ValCell,Role,units]=v(DataIn)
369VarName={};
370ValCell={};
371Role={};
372units={};
373if isfield(DataIn,'V')
374    VarName{1}='V';
375    ValCell{1}=DataIn.V;
376    Role{1}='scalar';
377    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
378        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
379    else
380        units={'pixel'};
381    end
382end
383
384%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
385function [VarName,ValCell,Role,units]=w(DataIn)
386VarName={};
387ValCell={};
388Role={};
389units={};
390if isfield(DataIn,'W')
391    VarName{1}='W';
392    ValCell{1}=DataIn.W;
393    Role{1}='scalar';%will remain unchanged by projection
394    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
395        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
396    else
397        units={'pixel'};
398    end
399end
400
401%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
402function [VarName,ValCell,Role,units]=w_normal(DataIn)
403VarName={};
404ValCell={};
405Role={};
406units={};
407if isfield(DataIn,'W')
408    VarName{1}='W';
409    ValCell{1}=DataIn.W;
410    Role{1}='vector_z';%will behave like a vector component  by projection
411    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
412        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
413    else
414        units={'pixel'};
415    end
416end
417
418%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
419function [VarName,ValCell,Role,units]=error(DataIn)
420VarName={};
421ValCell={};
422Role={};
423units={};
424if isfield(DataIn,'E')
425    VarName{1}='E';
426    ValCell{1}=DataIn.E;
427    Role{1}='ancillary'; %TODO CHECK units in actual fields
428    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
429        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
430    else
431        units={'pixel'};
432    end
433end
434
Note: See TracBrowser for help on using the repository browser.