source: trunk/src/calc_field.m @ 404

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

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