source: trunk/src/calc_field.m @ 405

Last change on this file since 405 was 405, checked in by sommeria, 14 years ago

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