source: trunk/src/calc_field.m @ 380

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

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