source: trunk/src/calc_field.m @ 516

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

improvement of calc-field and combination of two fields

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