source: trunk/src/calc_field.m @ 491

Last change on this file since 491 was 491, checked in by sommeria, 9 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
RevLine 
[404]1
2
[124]3%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
[8]4%---------------------------------------------------------------------
[399]5% [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)
[8]6%
[124]7% OUTPUT:
8% Scal: matlab vector representing the scalar values (length nbvec defined by var_read)
[123]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
[8]12%                   'vel': scalar calculated solely from velocity components
[124]13%                   'der': needs spatial derivatives
[8]14%                   'var': the scalar name directly corresponds to a field name in the netcdf files
15% error: error flag
[89]16%      error = 0; OK
17%      error = 1; the prescribed scalar cannot be read or calculated from available fields
18%
[8]19% INPUT:
[399]20% FieldList: cell array of strings representing the name(s) of the field(s) to calculate
[8]21% DataIn: structure representing the field, as defined in check_field_srtructure.m
[399]22% Coord_interp(:,nb_coord) optional set of coordinates to interpolate the field (use with thin plate shell)
[89]23%
[8]24% FUNCTION related
[124]25% varname_generator.m: determines the field names to read in the netcdf
26% file, depending on the scalar
[8]27
[399]28function [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)
[8]29
30%list of defined scalars to display in menus (in addition to 'ima_cor').
[124]31% a type is associated to each scalar:
[8]32%              'discrete': related to the individual velocity vectors, not interpolated by patch
[124]33%              'vel': calculated from velocity components, continuous field (interpolated with velocity)
34%              'der': needs spatial derivatives
[8]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
[123]37% the scalar is calculated from other fields, as explicited below
38
[404]39%% list of field options implemented
40FieldOptions={'velocity';...%image correlation corresponding to a vel vector
[124]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)
[8]51errormsg=[]; %default error message
[399]52if ~exist('FieldList','var')
[404]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}
[491]69        case {'u','v','velocity','norm_vel','ima_cor'}
[404]70            check_grid=1;% needs a regular grid
71        case{'vort','div','strain'}% needs spatial derivatives spatial derivatives
72            check_der=1;
[491]73%         case {'velocity','norm_vel','ima_cor'};
[404]74        otherwise
75            check_calc(ilist)=0;
76    end
77end
78FieldList=FieldList(check_calc==1);
[405]79% if isempty(FieldList)
80%     DataOut=DataIn;
81%     return
82% end
[404]83if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
84    nbcoord=3;
[8]85else
[404]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)
[406]95    %reproduce global attributes
96    DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;
[404]97    for ilist=1:numel(DataOut.ListGlobalAttribute)
98        DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
[8]99    end
[406]100
[405]101    %create a default grid if needed
102    if  ~exist('Coord_interp','var')
[406]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,:,:)));
[405]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;
[491]121%         npx=length(coord_x);
122%         npy=length(coord_y);
[406]123        DataOut.coord_x=[XMin XMax];
124        DataOut.coord_y=[YMin YMax];
[405]125        [XI,YI]=meshgrid(coord_x,coord_y);
[491]126%         XI=reshape(XI,[],1);
127%         YI=reshape(YI,[],1);
128        Coord_interp=cat(3,XI,YI);%[XI YI];
[405]129    end
[491]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);
[406]134    %initialise output
[404]135    nb_sites=size(Coord_interp,1);
136    nb_coord=size(Coord_interp,2);
[406]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}=[];
[404]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);
[389]150        end
[404]151    end
[406]152   
153    % interpolate data in each subdomain
[404]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)
[405]162        if check_grid
[404]163            EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
[246]164        end
[404]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
[412]168        ListVar={};
[399]169        for ilist=1:length(FieldList)
[412]170            var_count=numel(ListVar);
[404]171            switch FieldList{ilist}
172                case 'velocity'
[412]173                    ListVar=[ListVar {'U', 'V'}];
[411]174                    VarAttribute{var_count+1}.Role='vector_x';
175                    VarAttribute{var_count+2}.Role='vector_y';
[404]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'
[412]179                    ListVar=[ListVar {'u'}];
[411]180                    VarAttribute{var_count+1}.Role='scalar';
[406]181                    DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
[404]182                case 'v'
[412]183                    ListVar=[ListVar {'v'}];
[411]184                    VarAttribute{var_count+1}.Role='scalar';
[406]185                    DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
[405]186                case 'norm_vel'
[412]187                    ListVar=[ListVar {'norm_vel'}];
[411]188                    VarAttribute{var_count+1}.Role='scalar';
189                    U=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
[404]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'
[412]193                    ListVar=[ListVar {'vort'}];
[411]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);
[404]196                case 'div'
[412]197                    ListVar=[ListVar {'div'}];
[411]198                    VarAttribute{var_count+1}.Role='scalar';
[404]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'
[412]201                    ListVar=[ListVar {'strain'}];
[411]202                    VarAttribute{var_count+1}.Role='scalar';
[404]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);
[246]204            end
205        end
[404]206    end
[411]207    DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
[412]208    nbval(nbval==0)=1;% to avoid division by zero for averaging
[411]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
[412]214    DataOut.ListVarName=[DataOut.ListVarName ListVar];
215    for ifield=1:numel(ListVar)
[411]216        VarDimName{ifield}={'coord_y','coord_x'};
[412]217        DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
218        DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
[411]219    end
220    DataOut.FF=reshape(DataOut.FF,npy,npx);
221    DataOut.VarDimName=[DataOut.VarDimName VarDimName];     
222    DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
[404]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];
[246]234        end
[124]235    end
[404]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
[8]257end
258
[246]259
[404]260
[8]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}];
[124]281    Role=[Role {'vector_y'}];
[8]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';
[124]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')
[8]331        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
[124]332    else
[8]333        units={'pixel'};
[124]334    end
335end
[8]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
[124]356end
[8]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
[124]375end
[8]376
377%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
378function [VarName,ValCell,Role,units]=strain(DataIn)
379VarName={};
380ValCell={};
381Role={};
382units={};
383if isfield(DataIn,'DjUi')
[124]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);
[156]388    Role{1}='scalar';
[124]389    if isfield(DataIn,'TimeUnit')
[8]390        units={[DataIn.TimeUnit '-1']};
[124]391    else
[8]392        units={[]};
[124]393    end
394end
[8]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;
[124]456    Role{1}='vector_z';%will behave like a vector component  by projection
[8]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.