source: trunk/src/calc_field.m @ 582

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

bugs corrected for combining fields

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