source: trunk/src/calc_field.m @ 399

Last change on this file since 399 was 399, checked in by sommeria, 13 years ago

implementation of thin plate interpolation (proj on planes with mode 'filter'), rationalisation of variable formats in civ_matlab
File size: 15.7 KB
Line 
1%'calc_field': defines fields (velocity, vort, div...) from civx data 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
25
26function [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
37list_field={'velocity';...%image correlation corresponding to a vel vector
38    'ima_cor';...%image correlation corresponding to a vel vector
39    'norm_vel';...%norm of the velocity
40    'vort';...%vorticity
41    'div';...%divergence
42    'strain';...%rate of strain
43    'u';... %u velocity component
44    'v';... %v velocity component
45    'w';... %w velocity component
46    'w_normal';... %w velocity component normal to the plane
47    'error'}; %error associated to a vector (for stereo or patch)
48errormsg=[]; %default error message
49if ~exist('FieldList','var')
50    DataOut=list_field;% gives the list of possible fields in the absence of input
51else
52    if ~exist('DataIn','var')
53        DataIn=[];
54    end
55    if ischar(FieldList)
56        FieldList={FieldList};%convert a string input to a cell with one string element
57    end
58    if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
59        nbcoord=3;
60    else
61        nbcoord=2;
62    end
63    ListVarName={};
64    ValueList={};
65    RoleList={};
66    units_cell={};
67   
68    %% interpolation with new civ data
69    if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps')&& exist('Coord_interp','var')
70        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
71        for ilist=1:numel(DataOut.ListGlobalAttribute)
72            DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
73        end
74        DataOut.ListVarName={'coord_y','coord_x','FF'};
75        DataOut.VarDimName{1}='coord_y';
76        DataOut.VarDimName{2}='coord_x';
77        XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
78        YMax=max(max(DataIn.SubRange(2,:,:)));
79        XMin=min(min(DataIn.SubRange(1,:,:)));
80        YMin=min(min(DataIn.SubRange(2,:,:)));
81        check_der=0;
82        check_val=0;
83        nb_sites=size(Coord_interp,1);
84        nb_coord=size(Coord_interp,2);
85        for ilist=1:length(FieldList)
86            switch FieldList{ilist}
87                case 'velocity'
88                    check_val=1;
89                    DataOut.U=zeros(nb_sites,1);
90                    DataOut.V=zeros(nb_sites,1);
91                case{'vort','div','strain'}% case of spatial derivatives
92                    check_der=1;
93                    DataOut.(FieldList{ilist})=zeros(nb_sites,1);
94%                 otherwise % case of a scalar
95%                     check_val=1;
96%                     DataOut.(FieldList{ilist})=zeros(size(Coord_interp,1));
97            end
98        end
99        nbval=zeros(nb_sites,1);
100        NbSubDomain=size(DataIn.SubRange,3);
101        %DataIn.Coord_tps=DataIn.Coord_tps(1:end-3,:,:);% suppress the 3 zeros used to fit with the dimensions of variables
102        for isub=1:NbSubDomain
103            nbvec_sub=DataIn.NbSites(isub);
104            check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');
105            ind_sel=find(sum(check_range,2)==nb_coord);
106            %rho smoothing parameter
107            %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
108            %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
109            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
110            if check_val
111                EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
112            end
113            if check_der
114                [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'
115            end
116            for ilist=1:length(FieldList)
117                switch FieldList{ilist}
118                    case 'velocity'
119                        ListFields={'U', 'V'};
120                        VarAttributes{1}.Role='vector_x';
121                        VarAttributes{2}.Role='vector_y';
122                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
123                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
124                    case 'u'
125                        ListFields={'U'};
126                        VarAttributes{1}.Role='scalar';
127                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
128                    case 'v'
129                        ListFields={'V'};
130                        VarAttributes{1}.Role='scalar';
131                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
132                    case 'vort'
133                        ListFields={'vort'};
134                        VarAttributes{1}.Role='scalar';
135                        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);
136                    case 'div'
137                        ListFields={'div'};
138                        VarAttributes{1}.Role='scalar';
139                        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);
140                    case 'strain'
141                        ListFields={'strain'};
142                        VarAttributes{1}.Role='scalar';
143                        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);
144                end
145            end
146            DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
147%            DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
148            nbval(nbval==0)=1;
149%             switch FieldList{1}
150%                 case {'velocity','u','v'}
151%                     DataOut.U=reshape(DataOut.U./nbval,numel(yI),numel(xI));
152%                     DataOut.V=reshape(DataOut.V./nbval,numel(yI),numel(xI));
153%                 case 'vort'
154%                     DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
155%                 case 'div'
156%                     DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
157%                 case 'strain'
158%                     DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI));
159%             end
160            DataOut.ListVarName=[DataOut.ListVarName ListFields];
161            for ilist=3:numel(DataOut.ListVarName)
162                DataOut.VarDimName{ilist}={'coord_y','coord_x'};
163            end
164            DataOut.VarAttribute={[],[]};
165            DataOut.VarAttribute{3}.Role='errorflag';
166            DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
167        end
168    else
169       
170        %% civx data
171        DataOut=DataIn;
172        for ilist=1:length(FieldList)
173            if ~isempty(FieldList{ilist})
174                [VarName,Value,Role,units]=feval(FieldList{ilist},DataIn);%calculate field with appropriate function named FieldList{ilist}
175                ListVarName=[ListVarName VarName];
176                ValueList=[ValueList Value];
177                RoleList=[RoleList Role];
178                units_cell=[units_cell units];
179            end
180        end
181        %erase previous data (except coordinates)
182        for ivar=nbcoord+1:length(DataOut.ListVarName)
183            VarName=DataOut.ListVarName{ivar};
184            DataOut=rmfield(DataOut,VarName);
185        end
186        DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
187        if isfield(DataOut,'VarDimName')
188            DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
189        else
190            errormsg='element .VarDimName missing in input data';
191            return
192        end
193        DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
194        %append new data
195        DataOut.ListVarName=[DataOut.ListVarName ListVarName];
196        for ivar=1:length(ListVarName)
197            DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
198            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
199            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
200            DataOut.(ListVarName{ivar})=ValueList{ivar};
201        end
202    end
203end
204
205
206%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
207function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
208VarName={};
209ValCell={};
210Role={};
211units_cell={};
212if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
213    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
214else
215    units='pixel';
216end
217if isfield(DataIn,'U')
218    VarName=[VarName {'U'}];
219    ValCell=[ValCell {DataIn.U}];
220    Role=[Role {'vector_x'}];
221    units_cell=[units_cell {units}];
222end
223if isfield(DataIn,'V')
224    VarName=[VarName {'V'}];
225    ValCell=[ValCell {DataIn.V}];
226    Role=[Role {'vector_y'}];
227    units_cell=[units_cell {units}];
228end
229if isfield(DataIn,'W')
230    VarName=[VarName {'W'}];
231    ValCell=[ValCell {DataIn.W}];
232    Role=[Role {'vector_z'}];
233    units_cell=[units_cell {units}];
234end
235if isfield(DataIn,'F')
236    VarName=[VarName {'F'}];
237    ValCell=[ValCell {DataIn.F}];
238    Role=[Role {'warnflag'}];
239    units_cell=[units_cell {[]}];
240end
241if isfield(DataIn,'FF')
242    VarName=[VarName,{'FF'}];
243    ValCell=[ValCell {DataIn.FF}];
244    Role=[Role {'errorflag'}];
245    units_cell=[units_cell {[]}];
246end
247
248%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
249function [VarName,ValCell,Role,units]=ima_cor(DataIn)
250VarName={};
251ValCell={};
252Role={};
253units={};
254if isfield(DataIn,'C')
255    VarName{1}='C';
256    ValCell{1}=DataIn.C;
257    Role={'ancillary'};
258    units={[]};
259end
260
261%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
262function [VarName,ValCell,Role,units]=norm_vel(DataIn)
263VarName={};
264ValCell={};
265Role={};
266units={};
267if isfield(DataIn,'U') && isfield(DataIn,'V')
268    VarName{1}='norm_vel';
269    ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
270    if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
271        ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
272    end
273    ValCell{1}=sqrt(ValCell{1});
274    Role{1}='scalar';
275    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
276        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
277    else
278        units={'pixel'};
279    end
280end
281
282
283
284%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
285function [VarName,ValCell,Role,units]=vort(DataIn)
286VarName={};
287ValCell={};
288Role={};
289units={};
290if isfield(DataIn,'DjUi')
291    VarName{1}='vort';
292    ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
293    siz=size(ValCell{1});
294    ValCell{1}=reshape(ValCell{1},siz(1),1);
295    Role{1}='scalar';
296    if isfield(DataIn,'TimeUnit')
297        units={[DataIn.TimeUnit '-1']};
298    else
299        units={[]};
300    end
301end
302
303%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
304function [VarName,ValCell,Role,units]=div(DataIn)
305VarName={};
306ValCell={};
307Role={};
308units={};
309if isfield(DataIn,'DjUi')
310    VarName{1}='div';
311    ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
312    siz=size(ValCell{1});
313    ValCell{1}=reshape(ValCell{1},siz(1),1);
314    Role{1}='scalar';
315    if isfield(DataIn,'TimeUnit')
316        units={[DataIn.TimeUnit '-1']};
317    else
318        units={[]};
319    end
320end
321
322%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
323function [VarName,ValCell,Role,units]=strain(DataIn)
324VarName={};
325ValCell={};
326Role={};
327units={};
328if isfield(DataIn,'DjUi')
329    VarName{1}='strain';
330    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
331    siz=size(ValCell{1});
332    ValCell{1}=reshape(ValCell{1},siz(1),1);
333    Role{1}='scalar';
334    if isfield(DataIn,'TimeUnit')
335        units={[DataIn.TimeUnit '-1']};
336    else
337        units={[]};
338    end
339end
340
341%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
342function [VarName,ValCell,Role,units]=u(DataIn)
343VarName={};
344ValCell={};
345Role={};
346units={};
347if isfield(DataIn,'U')
348    VarName{1}='U';
349    ValCell{1}=DataIn.U;
350    Role{1}='scalar';
351    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
352        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
353    else
354        units={'pixel'};
355    end
356end
357
358%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
359function [VarName,ValCell,Role,units]=v(DataIn)
360VarName={};
361ValCell={};
362Role={};
363units={};
364if isfield(DataIn,'V')
365    VarName{1}='V';
366    ValCell{1}=DataIn.V;
367    Role{1}='scalar';
368    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
369        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
370    else
371        units={'pixel'};
372    end
373end
374
375%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
376function [VarName,ValCell,Role,units]=w(DataIn)
377VarName={};
378ValCell={};
379Role={};
380units={};
381if isfield(DataIn,'W')
382    VarName{1}='W';
383    ValCell{1}=DataIn.W;
384    Role{1}='scalar';%will remain unchanged by projection
385    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
386        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
387    else
388        units={'pixel'};
389    end
390end
391
392%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
393function [VarName,ValCell,Role,units]=w_normal(DataIn)
394VarName={};
395ValCell={};
396Role={};
397units={};
398if isfield(DataIn,'W')
399    VarName{1}='W';
400    ValCell{1}=DataIn.W;
401    Role{1}='vector_z';%will behave like a vector component  by projection
402    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
403        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
404    else
405        units={'pixel'};
406    end
407end
408
409%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
410function [VarName,ValCell,Role,units]=error(DataIn)
411VarName={};
412ValCell={};
413Role={};
414units={};
415if isfield(DataIn,'E')
416    VarName{1}='E';
417    ValCell{1}=DataIn.E;
418    Role{1}='ancillary'; %TODO CHECK units in actual fields
419    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
420        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
421    else
422        units={'pixel'};
423    end
424end
425
Note: See TracBrowser for help on using the repository browser.