Ignore:
Timestamp:
Aug 27, 2012, 4:38:41 PM (12 years ago)
Author:
sommeria
Message:

new conventions for find_field_cells .

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/find_field_cells.m

    r527 r530  
    22%    group the variables  into 'fields' with common dimensions
    33%------------------------------------------------------------------------
    4 % function  [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data)
     4% function  [CellInfo,NbDim,errormsg]=find_field_cells(Data)
    55%
    66% OUTPUT:
    7 % CellVaxIndex: cell whose elements are arrays of indices in the list data.ListVarName 
    8 %              CellvarIndex{i} represents a set of variables with the same dimensions
    9 % NbDim: array with the length of CellVarIndex, giving its  space dimension
    10 % CellVarType: cell array of structures with fields
    11 %      .coord_x, y, z: indices (in .ListVarname) of variables representing  unstructured coordinates x, y, z
    12 %      .vector_x,_y,_z: indices of variables giving the vector components x, y, z
    13 %      .warnflag: index of warnflag
    14 %      .errorflag: index of error flag
    15 %      .ancillary: indices of ancillary variables
    16 %      .image   : B/W image, (behaves like scalar)
    17 %      .color : color image, the last index, which is not a coordinate variable, represent the 3 color components rgb
    18 %      .discrete: like scalar, but set of data points without continuity, represented as dots in a usual plot, instead of continuous lines otherwise
    19 %      .scalar: scalar field (default)
    20 %      .coord: vector of indices of coordinate variables corresponding to matrix dimensions
    21 %
     7% CellInfo: cell of structures describing field cells
     8%     .CoordType:  type of coordinates for each field cell = 'scattered','grid','tps';
     9%     .CoordIndex: array of the indices of the variables representing the coordinates (in the order z,y,x)
     10%     .CoordSize: array of the nbre of values for each  coordinate in a grid, nbre of points in the unstructured case
     11%     .NbSite_tps:
     12%     .SubRange_tps
     13%     .VarIndex: arrays of the variable indices in the field cell
     14%     .VarIndex_ancillary: indices of ancillary variables
     15%              _color : color image, the last index, which is not a coordinate variable, represent the 3 color components rgb
     16%              _discrete: like scalar, but set of data points without continuity, represented as dots in a usual plot, instead of continuous lines otherwise
     17%              _errorflag: index of error flag
     18%              _image   : B/W image, (behaves like scalar)
     19%              _vector_x,_y,_z: indices of variables giving the vector components x, y, z
     20%              _warnflag: index of warnflag   
    2221%      .FieldRequest= 'interp_lin', 'interp_tps' indicate whether lin interpolation  or derivatives (tps) is needed to calculate the requested field
    23 %      .FieldNames = cell of fields to calculate from the fied cell
    24 %
     22%      .Operation = operation to be performed to finalise the field cell after projection
     23%      .SubCheck=0 /1 indicate that the field must be substracted (second  entry in uvmat)
     24% NbDim: array with the length of CellVarIndex, giving the space dimension of each field cell
    2525% errormsg: error message
    2626%   
    … …  
    5454%    group the variables  into 'fields' with common dimensions
    5555
    56 function [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data)
    57 CellVarIndex={};
    58 
    59 CellVarType=[];
     56function [CellInfo,NbDim,errormsg]=find_field_cells(Data)
     57% CellVarIndex={};
     58% CellInfo={}
     59NbDim=0;
     60% CellVarType=[];
    6061errormsg=[];
     62if ~isfield(Data,'ListVarName'), errormsg='the list of variables .ListVarName is missing';return;end
     63if ~isfield(Data,'VarDimName'), errormsg='the list of dimensions .VarDimName is missing';return;end
    6164nbvar=numel(Data.ListVarName);%number of field variables
     65if ~isequal(numel(Data.VarDimName),nbvar), errormsg='.ListVarName and .VarDimName have unequal length';return;end
     66if isfield(Data,'ListDimName')&& isfield(Data,'DimValue')&&isequal(numel(Data.ListDimName),numel(Data.DimValue))
     67    check_dim=1;% dimensions of data defined, data not needed for this function
     68else
     69    check_dim=0;
     70    for ilist=1:numel(Data.ListVarName)
     71        if ~isfield(Data,Data.ListVarName{ilist})
     72            errormsg=['missing variable ' Data.ListVarName{ilist}];
     73            return
     74        end
     75    end
     76end
    6277icell=0;
    6378
    64 NbDim=[];
    65 VarDimIndex=[];
    66 VarDimName={};
    67 if ~isfield(Data,'VarDimName')
    68     errormsg='missing .VarDimName';
    69     return
    70 end
     79% NbDim=[];
     80% VarDimIndex=[];
     81% VarDimName={};
    7182
    7283%% role of variables and list of requested operations
     84%ListRole={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','vector_x_tps','vector_y_tps','warnflag','errorflag',...
     85%   'ancillary','image','color','discrete','scalar','coord_tps'};% rmq vector_x_tps and vector_y_tps to be replaced by vector_x and vector_y
    7386Role=num2cell(blanks(nbvar));%initialize a cell array of nbvar blanks
    7487FieldRequest=regexprep(Role,' ',''); % fieldRequest set to '' by default
    … …  
    93106end
    94107
    95 %% loop on the list of variables, group them by common dimensions
    96 CellVarType=cell(1,length(CellVarIndex));
    97 for ivar=1:nbvar
    98     if ischar(Data.VarDimName{ivar})
    99         Data.VarDimName{ivar}=Data.VarDimName(ivar);%transform char chain into cell
    100     end
    101     DimCell=Data.VarDimName{ivar}; %dimensions associated with the variable #ivar
    102     testnewcell=1;
    103     for icell_prev=1:numel(CellVarIndex)%detect whether the dimensions of ivar fit with an existing cell
    104         PrevVarIndex=CellVarIndex{icell_prev};%list of variable indices in cell # icell_prev
    105         PrevDimCell=Data.VarDimName{PrevVarIndex(1)};%list of corresponding variable names
    106         if isequal(PrevDimCell,DimCell)
    107             CellVarIndex{icell_prev}=[CellVarIndex{icell_prev} ivar];% add variable index #ivar to the cell #icell_prev
    108             testnewcell=0; %existing cell detected
    109             break
    110         end
    111     end
    112     if testnewcell
    113         icell=icell+1;
    114         CellVarIndex{icell}=ivar;%put the current variable index in the new cell
    115         NbDim(icell)=numel(DimCell);%default   
    116         CellVarType{icell}=[];
    117     end
    118     if ~isempty(FieldRequest{ivar})
    119        CellVarType{icell}.FieldRequest=FieldRequest{ivar};
    120     end
    121     if ~isempty(Operation{ivar})
    122        CellVarType{icell}.Operation=Operation{ivar};
    123     end
    124     if CheckSub(ivar)
    125     CellVarType{icell}.CheckSub=1;
    126     end
    127 end
    128 
    129 %% find dimension variables
    130 checksinglecell=cellfun(@numel,CellVarIndex)==1 & NbDim==1;% find isolated cells with a single dimension
    131 ind_dim_var_cell=find(checksinglecell);
    132 %CoordType(ind_dim_var_cell)='dim_var';% to be used in output
    133 %VarDimIndex=cell(size(ind_dim_var_cell));
    134 VarDimName=cell(size(ind_dim_var_cell));
    135 for icoord=1:numel(ind_dim_var_cell)
    136     VarDimIndex(icoord)=CellVarIndex{ind_dim_var_cell(icoord)};
    137     VarDimName{icoord}=Data.VarDimName{VarDimIndex(icoord)}{1};
    138 end
    139 
    140 %% find the spatial dimensions and vector components
    141 ListRole={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','vector_x_tps','vector_y_tps','warnflag','errorflag',...
    142     'ancillary','image','color','discrete','scalar','coord_tps'};% rmq vector_x_tps and vector_y_tps to be replaced by vector_x and vector_y
    143 
    144 for ilist=1:numel(ListRole)
    145     VarType.(ListRole{ilist})=find(strcmp(ListRole{ilist},Role));
    146 end
    147 %look for tps coordinates
    148 if ~isempty(VarType.coord_tps)
    149     VarType.subrange_tps=[];
    150     VarType.nbsites_tps=[];
    151     select=zeros(1,numel(VarType.coord_tps));
    152     for ifield=1:numel(VarType.coord_tps)
    153         DimCell=Data.VarDimName{VarType.coord_tps(ifield)};
    154         if numel(DimCell)==3
    155             for ivardim=1:numel(Data.VarDimName)
    156                 if strcmp(Data.VarDimName{ivardim},DimCell{3})
    157                     VarType.nbsites_tps=[VarType.nbsites_tps ivardim];
    158                     select(ifield)=select(ifield)+1;
    159                 elseif strcmp(Data.VarDimName{ivardim}{1},DimCell{2}) && strcmp(Data.VarDimName{ivardim}{3},DimCell{3})
    160                     VarType.subrange_tps=[VarType.subrange_tps ivardim];
    161                     select(ifield)=select(ifield)+1;
     108%% find scattered (unstructured) coordinates
     109ivar_coord_x=find(strcmp('coord_x',Role));
     110% VarDimCell=cell(numel(ivar_coord_x));
     111check_select=zeros(1,nbvar);
     112CellInfo=cell(1,numel(ivar_coord_x));
     113NbDim=zeros(1,numel(ivar_coord_x));
     114for icell=1:numel(ivar_coord_x)
     115    DimCell=Data.VarDimName{ivar_coord_x(icell)};
     116    if ischar(DimCell),DimCell={DimCell};end
     117    check_cell=zeros(numel(DimCell),nbvar);
     118    for idim=1:numel(DimCell)
     119        for ivar=1:nbvar
     120            check_cell(idim,ivar)=max(strcmp(DimCell{idim},Data.VarDimName{ivar}));
     121        end
     122    end
     123    check_cell=sum(check_cell,1)==numel(DimCell);%logical array=1 for variables belonging to the current cell
     124    VarIndex=find(check_cell);
     125    if ~(numel(VarIndex)==1 && numel(DimCell)==1)% exclude case of isolated coord_x variable (treated later)
     126        if ~(numel(VarIndex)==1 && numel(DimCell)>1)% a variable is associated to coordinate
     127            CellInfo{icell}.CoordIndex=ivar_coord_x(icell);
     128            % size of coordinate var
     129            if check_dim
     130                for idim=1:numel(DimCell)
     131                 check_index= strcmp(DimCell{idim},Data.ListDimName);
     132                 CellInfo{icell}.CoordSize(idim)=Data.DimValue(check_index);
     133                end
     134                CellInfo{icell}.CoordSize=prod(CellInfo{icell}.CoordSize);
     135            else
     136                CellInfo{icell}.CoordSize=numel(Data.(Data.ListVarName{ivar_coord_x(icell)}));
     137            end
     138            ind_y=find(strcmp('coord_y',Role(VarIndex)));
     139            if numel(VarIndex)==2||isempty(ind_y)% no variable, except possibly y
     140                NbDim(icell)=1;
     141            else
     142                CellInfo{icell}.CoordType='scattered';
     143                ind_z=find(strcmp('coord_z',Role(VarIndex)));
     144                if numel(VarIndex)==3||isempty(ind_z)% no z variable, except possibly as a fct z(x,y)
     145                    CellInfo{icell}.CoordIndex=[VarIndex(ind_y) CellInfo{icell}.CoordIndex];
     146                    NbDim(icell)=2;
     147                else
     148                    CellInfo{icell}.CoordIndex=[VarIndex(ind_z) CellInfo{icell}.CoordIndex];
     149                    NbDim(icell)=3;
    162150                end
    163151            end
    164152        end
    165     end
    166     VarType.coord_tps=VarType.coord_tps(select==2);
    167     VarType.subrange_tps=VarType.subrange_tps(select==2);
    168     VarType.nbsites_tps=VarType.nbsites_tps(select==2);
    169 end
    170 
    171 index_remove=[];
    172 for icell=1:length(CellVarIndex)
    173     if checksinglecell(icell)
    174         continue
    175     end
    176     VarIndex=CellVarIndex{icell};%set of variable indices with the same dim
    177     check_remove=0;
    178     for ifield=1:numel(VarType.coord_tps)
    179         if isequal(VarIndex,VarType.coord_tps(ifield))||isequal(VarIndex,VarType.subrange_tps(ifield))||isequal(VarIndex,VarType.nbsites_tps(ifield))
    180             index_remove=[index_remove icell];% removes Coord_tps as field cell
    181             check_remove=1;
    182         end
    183     end
    184    
    185     if ~check_remove
    186         for ilist=1:numel(ListRole)
    187             CellVarType{icell}.(ListRole{ilist})=VarIndex(find(strcmp(ListRole{ilist},Role(VarIndex))));
    188         end
    189         DimCell=Data.VarDimName{VarIndex(1)};% list of dimensions for each variable in the cell #icell
    190         if numel(CellVarType{icell}.coord_x)>1 || numel(CellVarType{icell}.coord_y)>1 || numel(CellVarType{icell}.coord_z)>1
    191             errormsg='multiply defined coordinates  in the same cell';
    192             return
    193         end
    194         % case of x cordinate marked as a dimension variable (var name=dimension name)
    195         if isempty(CellVarType{icell}.coord_x)
    196             var_dim_index=find(strcmp(DimCell{1},Data.ListVarName(VarIndex)));
    197             if ~isempty(var_dim_index)
    198                 CellVarType{icell}.coord_x=VarIndex(var_dim_index);
     153        CellInfo{icell}.VarIndex=VarIndex;
     154        check_select=check_select|check_cell;
     155    end
     156end
     157
     158%% look for tps coordinates
     159ivar_remain=find(~check_select);
     160check_coord_tps= strcmp('coord_tps',Role(~check_select));
     161ivar_tps=ivar_remain(check_coord_tps);
     162for icell_tps=1:numel(ivar_tps)
     163    DimCell=Data.VarDimName{ivar_tps(icell_tps)};
     164    icell=numel(CellInfo)+icell_tps;
     165    CellInfo{icell}.CoordIndex=ivar_tps(icell);
     166    CellInfo{icell}.VarIndex_subrange_tps=[];
     167    CellInfo{icell}.VarIndex_nbsites_tps=[];
     168    if numel(DimCell)==3
     169        VarDimName=Data.VarDimName(~check_select);
     170        for ivardim=1:numel(VarDimName)
     171            if strcmp(VarDimName{ivardim},DimCell{3})
     172                CellInfo{icell}.NbSite_tps= ivar_remain(ivardim);
     173                check_cell(ivar_remain(ivardim))=1;% nbre of sites for each tps subdomain
     174            elseif strcmp(VarDimName{ivardim}{1},DimCell{2}) && strcmp(VarDimName{ivardim}{3},DimCell{3})
     175                CellInfo{icell}.SubRange_tps=ivar_remain(ivardim);
     176                check_cell(ivar_remain(ivardim))=1;% subrange definiton for tps
     177            elseif strcmp(VarDimName{ivardim}{1},DimCell{1}) && strcmp(VarDimName{ivardim}{2},DimCell{3})
     178                check_cell(ivar_remain(ivardim))=1;% variable
    199179            end
    200         end         
    201         if numel(CellVarType{icell}.errorflag)>1
    202             errormsg='multiply defined error flag in the same cell';
    203             return
    204         end
    205         if numel(CellVarType{icell}.warnflag)>1
    206             errormsg='multiply defined warning flag in the same cell';
    207             return
    208         end
    209         test_coord=0;
    210         % look for unstructured coordinates
    211         if numel(VarIndex)>1
    212             if ~isempty(CellVarType{icell}.coord_z)
    213                 NbDim(icell)=3;
    214                 test_coord=1;
    215             elseif ~isempty(CellVarType{icell}.coord_y)
    216                 NbDim(icell)=2;
    217                 test_coord=1;
    218             elseif ~isempty(CellVarType{icell}.coord_x)
    219                 NbDim(icell)=1;
    220                 test_coord=1;
    221             elseif numel(DimCell)==1
    222                 NbDim(icell)=0;% set of data without coordinates
    223             end
    224         end
    225         % look for coordinates variables
    226         coord=zeros(1,numel(DimCell));%default
    227         if  ~test_coord && ~isempty(VarDimName)
    228             for idim=1:numel(DimCell)   %loop on the dimensions of the variables in cell #icell
    229                 ind_coord=find(strcmp(DimCell{idim},VarDimName));
    230                 if ~isempty(ind_coord)
    231                     coord(idim)=VarDimIndex(ind_coord);
    232                 end
    233             end
    234             NbDim(icell)=numel(find(coord));
    235         end
    236         CellVarType{icell}.coord=coord;
    237         %look for tps data
    238         if ~isempty(VarType.coord_tps)
    239             for ilist=1:numel(VarType.coord_tps)
    240             tps_dimnames=Data.VarDimName{VarType.coord_tps(ilist)};
    241             if length(tps_dimnames)==3 && strcmp(tps_dimnames{1},DimCell{1}) && strcmp(tps_dimnames{3},DimCell{2})
    242                 CellVarIndex{icell}=[CellVarIndex{icell} VarType.coord_tps(ilist) VarType.nbsites_tps(ilist) VarType.subrange_tps(ilist)];
    243                 CellVarType{icell}.coord_tps=VarType.coord_tps(ilist);
    244                 CellVarType{icell}.nbsites_tps=VarType.nbsites_tps(ilist);
    245                 CellVarType{icell}.subrange_tps=VarType.subrange_tps(ilist);
    246                 if isfield(Data,'ListDimName')
    247                     dim_index=find(strcmp(tps_dimnames{2},Data.ListDimName));
    248                     NbDim(icell)=Data.DimValue(dim_index);
    249                 else
    250                 NbDim(icell)=size(Data.(Data.ListVarName{VarType.coord_tps(ilist)}),2);
    251                 end
    252             end
    253             end
    254         end
    255     end
    256 end
    257 if ~isempty(index_remove)
    258     CellVarIndex(index_remove)=[];
    259     CellVarType(index_remove)=[];
    260     NbDim(index_remove)=[];
    261 end
     180        end
     181    end
     182    CellInfo{icell}.CoordType='tps';
     183    CellInfo{icell}.VarIndex=find(check_cell);
     184    if check_dim
     185        check_index_1= strcmp(DimCell{1},Data.ListDimName);
     186        check_index_2= strcmp(DimCell{2},Data.ListDimName);
     187        check_index_3= strcmp(DimCell{3},Data.ListDimName);
     188        NbDim(icell)=Data.DimValue(check_index_2);
     189        CellInfo{icell}.CoordSize=Data.DimValue(check_index_1)*Data.DimValue(check_index_3);
     190    else
     191        NbDim(icell)=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),2);
     192        CellInfo{icell}.CoordSize=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),1)*size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),3);
     193    end
     194    check_select=check_select|check_cell;
     195end
     196
     197%% look for dimension variables and corresponding gridded data
     198ivar_remain=find(~check_select);
     199VarDimName=Data.VarDimName(~check_select);%dimensions of remaining variables
     200check_coord= cellfun(@numel,VarDimName)==1|cellfun(@ischar,VarDimName)==1;% find variables with a single dimension
     201ListCoordIndex=ivar_remain(check_coord);
     202ListCoordName=Data.ListVarName(ListCoordIndex);
     203ListDimName=Data.VarDimName(ListCoordIndex);
     204%remove redondant values
     205check_keep=logical(ones(size(ListDimName)));
     206for idim=1:numel(ListDimName)
     207    prev_ind=strcmp(ListDimName{idim},ListDimName(1:idim-1));
     208    if ~isempty(prev_ind)
     209        if strcmp(ListCoordName{idim},ListDimName{idim}) %coordinate variable
     210            check_keep(prev_ind)=0;
     211        else
     212           check_keep(idim)=0;
     213        end
     214    end
     215end
     216ListCoordIndex=ListCoordIndex(check_keep);
     217ListCoordName=ListCoordName(check_keep);
     218ListDimName=ListDimName(check_keep);
     219
     220CoordSize=[];
     221for ilist=1:numel(ListCoordIndex)
     222    if iscell(ListDimName{ilist})
     223        ListDimName(ilist)=ListDimName{ilist};%transform cell to string
     224    end
     225    if check_dim% if the list of dimensions is directly defined
     226        check_index= strcmp(ListDimName{ilist},Data.ListDimName);
     227        DimValue=Data.DimValue(check_index);
     228    else
     229        DimValue=numel(Data.(ListCoordName{ilist}));
     230    end
     231    if DimValue==2% case of uniform grid coordinate defined by lower and upper bounds only
     232        ListDimName{ilist}=ListCoordName{ilist};% look for dimensions with name equal to coordinate for
     233        if check_dim
     234            check_index= strcmp(ListCoordName{ilist},Data.ListDimName);
     235            CoordSize(ilist)=Data.DimValue(check_index);
     236        else
     237            CoordSize(ilist)=numel(Data.(ListCoordName{ilist}));
     238        end
     239    else
     240        CoordSize(ilist)=DimValue;
     241    end
     242end
     243NewCellInfo={};
     244NewCellDimIndex={};
     245NewNbDim=[];
     246for ivardim=1:numel(VarDimName) % loop on the list of remaining variables
     247    DimCell=VarDimName{ivardim};% dimension names of the current variable
     248    if ischar(DimCell), DimCell={DimCell}; end %transform char to cell if needed
     249    DimIndices=[];
     250    for idim=1:numel(DimCell)
     251        ind_dim=find(strcmp(DimCell{idim},ListDimName));%find the dim index in the list of coord dimensions
     252        if ~isempty(ind_dim)
     253            DimIndices=[DimIndices ind_dim]; %update the list of coord dimensions included in DimCell
     254        end
     255    end
     256    check_previous=0;
     257    for iprev=1:numel(NewCellInfo)
     258        if isequal(DimIndices,NewCellDimIndex{iprev})
     259            check_previous=1;
     260            NewCellInfo{iprev}.VarIndex=[NewCellInfo{iprev}.VarIndex ivar_remain(ivardim)];%append the current variable index to the found field cell
     261            break
     262        end
     263    end
     264    if ~check_previous
     265        nbcell=numel(NewCellInfo)+1;
     266        NewCellDimIndex{nbcell}=DimIndices;
     267        NewCellInfo{nbcell}.VarIndex=ivar_remain(ivardim);% create a new field cell with the current variable index
     268        NewNbDim(nbcell)=numel(DimIndices);
     269        NewCellInfo{nbcell}.CoordType='grid';
     270        NewCellInfo{nbcell}.CoordSize=CoordSize;
     271        NewCellInfo{nbcell}.CoordIndex=ListCoordIndex(DimIndices);
     272    end
     273end
     274NbDim=[NbDim NewNbDim];
     275CellInfo=[CellInfo NewCellInfo];
     276
     277%% suppress empty cells
     278check_empty=cellfun(@isempty,CellInfo);
     279%check_empty=cellfun(@isempty,CellVarIndex);
     280CellInfo(check_empty)=[];
     281
     282% CellVarIndex(check_empty)=[];
     283NbDim(check_empty)=[];
     284% CoordType(check_empty)=[];
     285% VarRole(check_empty)=[];
     286
     287%% document roles of non-coordinate variables
     288% ListRole={'vector_x','vector_y','vector_z','vector_x_tps','vector_y_tps','warnflag','errorflag',...
     289%    'ancillary','image','color','discrete','scalar'};% except coord,coord_x,_y,_z,Coord_tps already taken, into account
     290for icell=1:numel(CellInfo)
     291    VarIndex=CellInfo{icell}.VarIndex;
     292    for ivar=VarIndex
     293        if isfield(CellInfo{icell},['VarIndex_' Role{ivar}])
     294            CellInfo{icell}.(['VarIndex_' Role{ivar}])=[CellInfo{icell}.(['VarIndex_' Role{ivar}]) ivar];
     295        else
     296            CellInfo{icell}.(['VarIndex_' Role{ivar}])= ivar;
     297        end
     298        if ~isempty(FieldRequest{ivar})
     299            CellInfo{icell}.FieldRequest=FieldRequest{ivar};
     300        end
     301        if ~isempty(Operation{ivar})
     302            CellInfo{icell}.Operation=Operation{ivar};
     303        end
     304        if CheckSub(ivar)==1
     305            CellInfo{icell}.CheckSub=1;
     306        end
     307    end
     308end
Note: See TracChangeset for help on using the changeset viewer.