source: trunk/src/calc_field.m @ 575

Last change on this file since 575 was 575, checked in by sommeria, 11 years ago

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