source: trunk/src/calc_field.m @ 491

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

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