Changeset 905 for trunk/src/proj_field.m


Ignore:
Timestamp:
May 30, 2015, 8:44:20 PM (8 years ago)
Author:
sommeria
Message:

projection with tps corrected + minor changes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r896 r905  
    522522function  [ProjData,errormsg] = proj_line(FieldData, ObjectData)
    523523%-----------------------------------------------------------------
     524
     525%% prepare heading for the projected field
    524526[ProjData,errormsg]=proj_heading(FieldData,ObjectData);%transfer global attributes
    525527if ~isempty(errormsg)
     
    544546ListIndex={};
    545547
     548%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
     549[CellInfo,NbDim,errormsg]=find_field_cells(FieldData);
     550if ~isempty(errormsg)
     551    errormsg=['error in proj_field/proj_line:' errormsg];
     552    return
     553end
     554CellInfo=CellInfo(NbDim==2); %keep only the 2D cells
     555%%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation
     556cell_select=true(size(CellInfo));
     557
     558for icell=1:length(CellInfo)
     559    if isfield(CellInfo{icell},'ProjModeRequest')
     560        if ~strcmp(CellInfo{icell}.ProjModeRequest, ProjMode)
     561            cell_select(icell)=0;
     562        end
     563        if strcmp(ProjMode,'interp_tps')&& ~strcmp(CellInfo{icell}.CoordType,'tps')
     564            cell_select(icell)=0;
     565        end
     566    end
     567end
     568if isempty(find(cell_select))
     569    errormsg=[' invalid projection mode ''' ProjMode ''': use ''interp_tps'' to interpolate spatial derivatives'];
     570    return
     571end
     572CellInfo=CellInfo(cell_select);
    546573
    547574%% projection line: object types selected from  proj_field='line','polyline','polygon','rectangle','ellipse':
     
    633660    end
    634661end
    635 
    636 
    637 %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
    638 [CellInfo,NbDim,errormsg]=find_field_cells(FieldData);
    639 if ~isempty(errormsg)
    640     errormsg=['error in proj_field/proj_line:' errormsg];
    641     return
    642 end
    643 CellInfo=CellInfo(NbDim==2); %keep only the 2D cells
    644 %%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation
    645662
    646663%% loop on variable cells with the same space dimension 2
     
    12871304                    Dist=Distx.*Distx+Disty.*Disty;
    12881305                    for ivar=1:numel(VarVal)
    1289                         VarVal{ivar}(Dist>4*ProjData.CoordMesh)=NaN;% put to NaN interpolated positions too far from initial data
     1306                        VarVal{ivar}(Dist>16*ProjData.CoordMesh)=NaN;% % put to NaN interpolated positions further than 4 meshes from initial data
    12901307                    end 
    12911308                   
     
    13251342               
    13261343                % set to NaN interpolation points which are too far from any initial data (more than 2 CoordMesh)
    1327                     if exist('scatteredInterpolant','file')%recent Matlab versions
    1328                         F=scatteredInterpolant(coord_X, coord_Y,coord_X,'nearest');
    1329                         G=scatteredInterpolant(coord_X, coord_Y,coord_Y,'nearest');
    1330                     else
    1331                         F=TriScatteredInterp([coord_X coord_Y],coord_X,'nearest');
    1332                         G=TriScatteredInterp([coord_X coord_Y],coord_Y,'nearest');
    1333                     end
    1334                     Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point
    1335                     Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point
    1336                     Dist=Distx.*Distx+Disty.*Disty;
    1337                     for ivar=1:numel(VarVal)
    1338                         VarVal{ivar}(Dist>4*ProjData.CoordMesh)=NaN;% put to NaN interpolated positions too far from initial data
    1339                     end 
    1340                    
    1341                
     1344                Coord=permute(Coord,[1 3 2]);
     1345                Coord=reshape(Coord,size(Coord,1)*size(Coord,2),2);
     1346                if exist('scatteredInterpolant','file')%recent Matlab versions
     1347                    F=scatteredInterpolant(Coord,Coord(:,1),'nearest');
     1348                    G=scatteredInterpolant(Coord,Coord(:,2),'nearest');
     1349                else
     1350                    F=TriScatteredInterp(Coord,Coord(:,1),'nearest');
     1351                    G=TriScatteredInterp(Coord,Coord(:,2),'nearest');
     1352                end
     1353                Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point
     1354                Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point
     1355                Dist=Distx.*Distx+Disty.*Disty;
    13421356                ListVarName=(fieldnames(DataOut))';
    13431357                VarDimName=cell(size(ListVarName));
     
    13451359                    VarName=ListVarName{ilist};
    13461360                    ProjData.(VarName)=DataOut.(VarName);
     1361                    ProjData.(VarName)(Dist>16*ProjData.CoordMesh)=NaN;% put to NaN interpolated positions further than 4 meshes from initial data
    13471362                    VarDimName{ilist}={'coord_y','coord_x'};
    13481363                end
Note: See TracChangeset for help on using the changeset viewer.