Changeset 491 for trunk/src/proj_field.m
- Timestamp:
- Jul 11, 2012, 2:56:40 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/proj_field.m
r445 r491 970 970 ProjData.ListVarName={}; 971 971 ProjData.VarDimName={}; 972 ProjData.VarAttribute={}; 972 973 if ~isequal(DX,0)&& ~isequal(DY,0) 973 974 ProjData.Mesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots … … 1004 1005 VarIndex=CellVarIndex{icell};% indices of the selected variables in the list FieldData.ListVarName 1005 1006 VarType=VarTypeCell{icell}; 1006 ivar_X=VarType.coord_x;1007 ivar_Y=VarType.coord_y;1008 ivar_Z=VarType.coord_z;1007 % ivar_X=VarType.coord_x; 1008 % ivar_Y=VarType.coord_y; 1009 % ivar_Z=VarType.coord_z; 1009 1010 ivar_U=VarType.vector_x; 1010 1011 ivar_V=VarType.vector_y; 1012 if ~isempty(VarType.vector_x_tps)&&~isempty(VarType.vector_y_tps) 1013 ivar_U=VarType.vector_x_tps; 1014 ivar_V=VarType.vector_y_tps; 1015 end 1011 1016 ivar_W=VarType.vector_z; 1012 % ivar_C=VarType.scalar ; 1013 ivar_Anc=VarType.ancillary; 1014 test_anc=zeros(size(VarIndex)); 1015 test_anc(ivar_Anc)=ones(size(ivar_Anc)); 1016 ivar_F=VarType.warnflag; 1017 ivar_FF=VarType.errorflag; 1018 check_unstructured_coord=(~isempty(ivar_X) && ~isempty(ivar_Y))||~isempty(VarType.coord_tps); 1017 % ivar_Anc=VarType.ancillary; 1018 % test_anc=zeros(size(VarIndex)); 1019 % test_anc(ivar_Anc)=ones(size(ivar_Anc)); 1020 % ivar_F=VarType.warnflag; 1021 % ivar_FF=VarType.errorflag; 1022 1023 %type of coordinates 1024 if ~isempty(VarType.coord_x) && ~isempty(VarType.coord_y) 1025 CoordType='unstructured'; 1026 elseif ~isempty(VarType.coord_tps) 1027 CoordType='tps'; 1028 else 1029 CoordType='structured'; 1030 end 1031 1032 %dimensions 1019 1033 DimCell=FieldData.VarDimName{VarIndex(1)}; 1020 1034 if ischar(DimCell) 1021 1035 DimCell={DimCell};%name of dimensions 1022 1036 end 1037 coord_z=0;%default 1023 1038 1024 %% case of input fields with unstructured coordinates 1025 coord_z=0;%default 1026 if check_unstructured_coord 1027 if ~isempty(ivar_X)% case of tps excluded. TODO similar procedure 1028 XName=FieldData.ListVarName{ivar_X}; 1029 YName=FieldData.ListVarName{ivar_Y}; 1030 coord_x=FieldData.(XName); 1031 coord_y=FieldData.(YName); 1032 if length(ivar_Z)==1 1033 ZName=FieldData.ListVarName{ivar_Z}; 1034 coord_z=FieldData.(ZName); 1039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1040 switch CoordType 1041 1042 %% case of input fields with unstructured coordinates 1043 case 'unstructured' 1044 if strcmp(ObjectData.ProjMode,'filter') 1045 continue %skip for filter (needs tps) 1046 end 1047 coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x}); 1048 coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y}); 1049 if ~isempty(VarType.coord_z) 1050 coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z}); 1035 1051 end 1036 1052 … … 1038 1054 coord_x=coord_x-ObjectData.Coord(1,1); 1039 1055 coord_y=coord_y-ObjectData.Coord(1,2); 1040 if ~isempty( ivar_Z)1056 if ~isempty(VarType.coord_z) 1041 1057 coord_z=coord_z-ObjectData.Coord(1,3); 1042 1058 end 1043 1059 1044 1060 % selection of the vectors in the projection range (3D case) 1045 if length(ivar_Z)==1&& width > 01061 if ~isempty(VarType.coord_z) && width > 0 1046 1062 %components of the unitiy vector normal to the projection plane 1047 1063 fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane … … 1105 1121 coord_X=coord_X(indcut); 1106 1122 coord_Y=coord_Y(indcut); 1107 if length(ivar_Z)==11123 if ~isempty(VarType.coord_z) 1108 1124 coord_Z=coord_Z(indcut); 1109 1125 end 1110 1126 end 1111 end 1112 % different cases of projection 1113 switch ObjectData.ProjMode 1114 case 'projection' 1115 for ivar=VarIndex %transfer variables to the projection plane 1116 VarName=FieldData.ListVarName{ivar}; 1117 if ivar==ivar_X %x coordinate 1118 eval(['ProjData.' VarName '=coord_X;']) 1119 elseif ivar==ivar_Y % y coordinate 1120 eval(['ProjData.' VarName '=coord_Y;']) 1121 elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced) 1122 eval(['ProjData.' VarName '=FieldData.' VarName ';']) 1123 end 1124 if isempty(ivar_Z) || ivar~=ivar_Z 1125 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1126 ProjData.VarDimName=[ProjData.VarDimName DimCell]; 1127 nbvar=nbvar+1; 1128 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar 1129 ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; 1127 1128 % different cases of projection 1129 switch ObjectData.ProjMode 1130 case 'projection' 1131 for ivar=VarIndex %transfer variables to the projection plane 1132 VarName=FieldData.ListVarName{ivar}; 1133 if ivar==VarType.coord_x %x coordinate 1134 ProjData.(VarName)=coord_X; 1135 elseif ivar==VarType.coord_y % y coordinate 1136 ProjData.(VarName)=coord_Y; 1137 elseif isempty(VarType.coord_z) || ivar~=VarType.coord_z % other variables (except Z coordinate wyhich is not reproduced) 1138 ProjData.(VarName)=FieldData.(VarName); 1130 1139 end 1131 end 1132 end 1133 case 'interp'%interpolate data on a regular grid 1140 if isempty(VarType.coord_z) || ivar~=VarType.coord_z 1141 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1142 ProjData.VarDimName=[ProjData.VarDimName DimCell]; 1143 nbvar=nbvar+1; 1144 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar 1145 ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; 1146 end 1147 end 1148 end 1149 case 'interp'%interpolate data on a regular grid 1150 coord_x_proj=XMin:DX:XMax; 1151 coord_y_proj=YMin:DY:YMax; 1152 DimCell={'coord_y','coord_x'}; 1153 ProjData.ListVarName={'coord_y','coord_x'}; 1154 ProjData.VarDimName={'coord_y','coord_x'}; 1155 nbcoord=2; 1156 ProjData.coord_y=[YMin YMax]; 1157 ProjData.coord_x=[XMin XMax]; 1158 % if isempty(ivar_X), ivar_X=0; end; 1159 % if isempty(ivar_Y), ivar_Y=0; end; 1160 % if isempty(ivar_Z), ivar_Z=0; end; 1161 % if isempty(ivar_U), ivar_U=0; end; 1162 % if isempty(ivar_V), ivar_V=0; end; 1163 % if isempty(ivar_W), ivar_W=0; end; 1164 % if isempty(ivar_F), ivar_F=0; end; 1165 % if isempty(ivar_FF), ivar_FF=0; end; 1166 % if ~isequal(ivar_FF,0) 1167 if ~isempty(VarType.errorflag) 1168 VarName_FF=FieldData.ListVarName{VarType.errorflag}; 1169 indsel=find(FieldData.(VarName_FF)==0); 1170 coord_X=coord_X(indsel); 1171 coord_Y=coord_Y(indsel); 1172 end 1173 1174 FF=zeros(1,length(coord_y_proj)*length(coord_x_proj)); 1175 testFF=0; 1176 %for ivar=VarIndex % loop on field variables to project 1177 for ivar=[VarType.scalar VarType.vector_x VarType.vector_y] 1178 VarName=FieldData.ListVarName{ivar}; 1179 % if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF ||... 1180 % ( numel(test_anc)>=ivar && test_anc(ivar)==1)) 1181 ivar_new=ivar_new+1; 1182 ProjData.ListVarName=[ProjData.ListVarName {VarName}]; 1183 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1184 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar 1185 ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar}; 1186 end 1187 %if ~isequal(ivar_FF,0) 1188 if ~isempty(VarType.errorflag) 1189 FieldData.(VarName)=FieldData.(VarName)(indsel); 1190 end 1191 ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj'); 1192 varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj)); 1193 FFlag= isnan(varline); %detect undefined values NaN 1194 indnan=find(FFlag); 1195 if~isempty(indnan) 1196 varline(indnan)=zeros(size(indnan)); 1197 ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj)); 1198 FF(indnan)=ones(size(indnan)); 1199 testFF=1; 1200 end 1201 if ivar==ivar_U 1202 ivar_U=ivar_new; 1203 end 1204 if ivar==ivar_V 1205 ivar_V=ivar_new; 1206 end 1207 if ivar==ivar_W 1208 ivar_W=ivar_new; 1209 end 1210 % end 1211 end 1212 if testFF 1213 ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); 1214 ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; 1215 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1216 ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; 1217 end 1218 end 1219 1220 %% case of tps interpolation (applies only in filter mode) 1221 case 'tps' 1222 if strcmp(ObjectData.ProjMode,'filter') 1134 1223 coord_x_proj=XMin:DX:XMax; 1135 1224 coord_y_proj=YMin:DY:YMax; 1136 DimCell={'coord_y','coord_x'}; 1137 ProjData.ListVarName={'coord_y','coord_x'}; 1138 ProjData.VarDimName={'coord_y','coord_x'}; 1139 nbcoord=2; 1225 np_x=numel(coord_x_proj); 1226 np_y=numel(coord_y_proj); 1227 [XI,YI]=meshgrid(coord_x_proj,coord_y_proj'); 1228 XI=XI+ObjectData.Coord(1,1); 1229 YI=YI+ObjectData.Coord(1,2); 1230 DataOut=calc_field(FieldData.FieldList,FieldData,cat(3,XI,YI)); 1231 ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName]; 1232 ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName]; 1233 ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute]; 1234 for ilist=1:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise 1235 VarName=DataOut.ListVarName{ilist}; 1236 % ProjData.(VarName)=DataOut.(VarName); 1237 if ilist>=3 1238 ProjData.(VarName)=reshape(DataOut.(VarName),np_y,np_x); 1239 end 1240 end 1241 ProjData.coord_x=[XMin XMax]; 1140 1242 ProjData.coord_y=[YMin YMax]; 1141 ProjData.coord_x=[XMin XMax]; 1142 if isempty(ivar_X), ivar_X=0; end; 1143 if isempty(ivar_Y), ivar_Y=0; end; 1144 if isempty(ivar_Z), ivar_Z=0; end; 1145 if isempty(ivar_U), ivar_U=0; end; 1146 if isempty(ivar_V), ivar_V=0; end; 1147 if isempty(ivar_W), ivar_W=0; end; 1148 if isempty(ivar_F), ivar_F=0; end; 1149 if isempty(ivar_FF), ivar_FF=0; end; 1150 if ~isequal(ivar_FF,0) 1151 VarName_FF=FieldData.ListVarName{ivar_FF}; 1152 indsel=find(FieldData.(VarName_FF)==0); 1153 coord_X=coord_X(indsel); 1154 coord_Y=coord_Y(indsel); 1155 end 1156 1157 FF=zeros(1,length(coord_y_proj)*length(coord_x_proj)); 1158 testFF=0; 1159 for ivar=VarIndex 1160 VarName=FieldData.ListVarName{ivar}; 1161 if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF ||... 1162 ( numel(test_anc)>=ivar && test_anc(ivar)==1)) 1163 ivar_new=ivar_new+1; 1164 ProjData.ListVarName=[ProjData.ListVarName {VarName}]; 1165 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1166 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar 1167 ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar}; 1243 end 1244 1245 %% case of input fields defined on a structured grid 1246 case 'structured' 1247 1248 VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions 1249 eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions 1250 DimValue(DimValue==1)=[];%remove singleton dimensions 1251 NbDim=numel(DimValue);%update number of space dimensions 1252 nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate 1253 if NbDim>=3 1254 if NbDim>3 1255 errormsg='matrices with more than 3 dimensions not handled'; 1256 return 1257 else 1258 if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate 1259 nbcolor=DimValue(3); 1260 DimValue(3)=[]; %number of 'color' components updated 1261 NbDim=2;% space dimension set to 2 1262 end 1263 end 1264 end 1265 AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection) 1266 AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection) 1267 if testangle% TODO modify name also in case of origin shift in x or y 1268 AYProjName='Y'; 1269 AXProjName='X'; 1270 count=0; 1271 %modify coordinate names if they are already used 1272 while ~(isempty(find(strcmp('AXName',ProjData.ListVarName),1)) && isempty(find(strcmp('AYName',ProjData.ListVarName),1))) 1273 count=count+1; 1274 AYProjName=[AYProjName '_' num2str(count)]; 1275 AXProjName=[AXProjName '_' num2str(count)]; 1276 end 1277 else 1278 AYProjName=AYName;% (name preserved on projection) 1279 AXProjName=AXName;%name of input y coordinate (name preserved on projection) 1280 end 1281 ListDimName=FieldData.VarDimName{VarIndex(1)}; 1282 ProjData.ListVarName=[ProjData.ListVarName {AYProjName} {AXProjName}]; %TODO: check if it already exists in Projdata (several cells) 1283 ProjData.VarDimName=[ProjData.VarDimName {AYProjName} {AXProjName}]; 1284 Coord_z=[]; 1285 Coord_y=[]; 1286 Coord_x=[]; 1287 1288 for idim=1:NbDim %loop on space dimensions 1289 test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default 1290 ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension 1291 if ~isequal(ivar,0)% a variable corresponds to the dimension #idim 1292 eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field 1293 if numel(Coord{idim})==2 %input array defined on a regular grid 1294 DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim); 1295 else 1296 DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field 1297 DCoord_min(idim)=min(DCoord); 1298 DCoord_max=max(DCoord); 1299 % test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise 1300 if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000) 1301 msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']) 1302 return 1168 1303 end 1169 if ~isequal(ivar_FF,0) 1170 FieldData.(VarName)=FieldData.(VarName)(indsel); 1304 test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity 1305 end 1306 test_direct(idim)=(DCoord_min(idim)>0); 1307 else % no variable associated with the dimension #idim, the coordinate value is set equal to the matrix index by default 1308 Coord_i_str=['Coord_' num2str(idim)]; 1309 DCoord_min(idim)=1;%default 1310 Coord{idim}=[0.5 DimValue(idim)-0.5]; 1311 test_direct(idim)=1; 1312 end 1313 end 1314 if DY==0 1315 DY=abs(DCoord_min(NbDim-1)); 1316 end 1317 npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol 1318 if DX==0 1319 DX=abs(DCoord_min(NbDim)); 1320 end 1321 npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol 1322 for idim=1:NbDim 1323 if test_interp(idim) 1324 DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri 1325 end 1326 end 1327 Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY); 1328 test_direct_y=test_direct(NbDim-1); 1329 Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX); 1330 test_direct_x=test_direct(NbDim); 1331 DAX=DCoord_min(NbDim); 1332 DAY=DCoord_min(NbDim-1); 1333 minAX=min(Coord_x); 1334 maxAX=max(Coord_x); 1335 minAY=min(Coord_y); 1336 maxAY=max(Coord_y); 1337 xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1); 1338 ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2); 1339 xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame 1340 ycor_new=-xcorner*sin_om+ycorner*cos_om; 1341 if ~testXMax 1342 XMax=max(xcor_new); 1343 end 1344 if ~testXMin 1345 XMin=min(xcor_new); 1346 end 1347 if ~testYMax 1348 YMax=max(ycor_new); 1349 end 1350 if ~testYMin 1351 YMin=min(ycor_new); 1352 end 1353 DXinit=(maxAX-minAX)/(DimValue(NbDim)-1); 1354 DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1); 1355 if DX==0 1356 DX=DXinit; 1357 end 1358 if DY==0 1359 DY=DYinit; 1360 end 1361 if NbDim==3 1362 DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1); 1363 if ~test_direct(1) 1364 DZ=-DZ; 1365 end 1366 Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1)); 1367 test_direct_z=test_direct(1); 1368 end 1369 npX=floor((XMax-XMin)/DX+1); 1370 npY=floor((YMax-YMin)/DY+1); 1371 if test_direct_y 1372 coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line 1373 else 1374 coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line 1375 end 1376 if test_direct_x 1377 coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line 1378 else 1379 coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line 1380 end 1381 % case with no interpolation 1382 if isequal(ProjMode,'projection') && (~testangle || test90y || test90x) 1383 if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax 1384 ProjData=FieldData;% no change by projection 1385 else 1386 indY=NbDim-1; 1387 if test_direct(indY) 1388 min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1; 1389 max_indy=floor((YMax-Coord{indY}(1))/DYinit)+1; 1390 Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1); 1391 Ybound(2)=Coord{indY}(1)+DYinit*(max_indy-1); 1392 else 1393 min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1; 1394 max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1; 1395 Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1); 1396 Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1); 1397 end 1398 if test_direct(NbDim)==1 1399 min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1; 1400 max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1; 1401 Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); 1402 Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); 1403 else 1404 min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1; 1405 max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1; 1406 Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); 1407 Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); 1408 end 1409 min_indy=max(min_indy,1);% deals with margin (bound lower than the first index) 1410 min_indx=max(min_indx,1); 1411 1412 if test90y 1413 ind_new=[3 2 1]; 1414 DimCell={AYProjName,AXProjName}; 1415 % DimValue=DimValue(ind_new); 1416 iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1; 1417 for ivar=VarIndex 1418 VarName=FieldData.ListVarName{ivar}; 1419 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1420 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1421 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1422 eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new);'])% permute x and z indices for 90 degree rotation 1423 eval(['ProjData.' VarName '=squeeze(ProjData.' VarName '(iz,:,:));'])% select the z index iz 1171 1424 end 1172 ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj'); 1173 varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj)); 1174 FFlag= isnan(varline); %detect undefined values NaN 1175 indnan=find(FFlag); 1176 if~isempty(indnan) 1177 varline(indnan)=zeros(size(indnan)); 1178 ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj)); 1179 FF(indnan)=ones(size(indnan)); 1180 testFF=1; 1425 eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1426 eval(['ProjData.' AXProjName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates 1427 else 1428 if NbDim==3 1429 DimCell(1)=[]; %suppress z variable 1430 DimValue(1)=[]; 1431 if test_direct(1) 1432 iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1; 1433 else 1434 iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1; 1435 end 1181 1436 end 1182 if ivar==ivar_U 1183 ivar_U=ivar_new; 1437 max_indy=min(max_indy,DimValue(1));%introduce bounds in y and x indices 1438 max_indx=min(max_indx,DimValue(2)); 1439 for ivar=VarIndex% loop on non coordinate variables 1440 VarName=FieldData.ListVarName{ivar}; 1441 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1442 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1443 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar 1444 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; 1445 end 1446 if NbDim==3 1447 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']); 1448 else 1449 eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']); 1450 end 1184 1451 end 1185 if ivar==ivar_V 1186 ivar_V=ivar_new; 1187 end 1188 if ivar==ivar_W 1189 ivar_W=ivar_new; 1190 end 1191 end 1192 end 1193 if testFF 1194 ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); 1195 ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; 1196 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1197 ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; 1198 end 1199 case 'filter' 1200 if ~isempty(VarType.coord_tps) 1201 %Coord_tps=FieldData.ListVarName{VarType.coord_tps}; 1202 %TODO: possibly translate and rotate coordiantes translate initial coordinates 1203 coord_x_proj=XMin:DX:XMax; 1204 coord_y_proj=YMin:DY:YMax; 1205 np_x=numel(coord_x_proj); 1206 np_y=numel(coord_y_proj); 1207 [XI,YI]=meshgrid(coord_x_proj,coord_y_proj'); 1208 XI=reshape(XI,[],1)+ObjectData.Coord(1,1); 1209 YI=reshape(YI,[],1)+ObjectData.Coord(1,2); 1210 ProjData=calc_field(FieldData.FieldList,FieldData,[XI YI]); 1211 for ilist=3:length(ProjData.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise 1212 VarName=ProjData.ListVarName{ilist}; 1213 ProjData.(VarName)=reshape(ProjData.(VarName),np_y,np_x); 1214 end 1215 ProjData.coord_x=[XMin XMax]; 1216 ProjData.coord_y=[YMin YMax]; 1217 end 1218 end 1219 %% case of input fields defined on a structured grid 1220 else 1221 VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions 1222 eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions 1223 DimValue(DimValue==1)=[];%remove singleton dimensions 1224 NbDim=numel(DimValue);%update number of space dimensions 1225 nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate 1226 if NbDim>=3 1227 if NbDim>3 1228 errormsg='matrices with more than 3 dimensions not handled'; 1229 return 1230 else 1231 if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate 1232 nbcolor=DimValue(3); 1233 DimValue(3)=[]; %number of 'color' components updated 1234 NbDim=2;% space dimension set to 2 1235 end 1236 end 1237 end 1238 AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection) 1239 AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection) 1240 if testangle% TODO modify name also in case of origin shift in x or y 1241 AYProjName='Y'; 1242 AXProjName='X'; 1243 count=0; 1244 %modify coordinate names if they are already used 1245 while ~(isempty(find(strcmp('AXName',ProjData.ListVarName),1)) && isempty(find(strcmp('AYName',ProjData.ListVarName),1))) 1246 count=count+1; 1247 AYProjName=[AYProjName '_' num2str(count)]; 1248 AXProjName=[AXProjName '_' num2str(count)]; 1249 end 1250 else 1251 AYProjName=AYName;% (name preserved on projection) 1252 AXProjName=AXName;%name of input y coordinate (name preserved on projection) 1253 end 1254 ListDimName=FieldData.VarDimName{VarIndex(1)}; 1255 ProjData.ListVarName=[ProjData.ListVarName {AYProjName} {AXProjName}]; %TODO: check if it already exists in Projdata (several cells) 1256 ProjData.VarDimName=[ProjData.VarDimName {AYProjName} {AXProjName}]; 1257 Coord_z=[]; 1258 Coord_y=[]; 1259 Coord_x=[]; 1260 1261 for idim=1:NbDim %loop on space dimensions 1262 test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default 1263 ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension 1264 if ~isequal(ivar,0)% a variable corresponds to the dimension #idim 1265 eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field 1266 if numel(Coord{idim})==2 %input array defined on a regular grid 1267 DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim); 1268 else 1269 DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field 1270 DCoord_min(idim)=min(DCoord); 1271 DCoord_max=max(DCoord); 1272 % test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise 1273 if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000) 1274 msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']) 1275 return 1276 end 1277 test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity 1278 end 1279 test_direct(idim)=(DCoord_min(idim)>0); 1280 else % no variable associated with the dimension #idim, the coordinate value is set equal to the matrix index by default 1281 Coord_i_str=['Coord_' num2str(idim)]; 1282 DCoord_min(idim)=1;%default 1283 Coord{idim}=[0.5 DimValue(idim)-0.5]; 1284 test_direct(idim)=1; 1285 end 1286 end 1287 if DY==0 1288 DY=abs(DCoord_min(NbDim-1)); 1289 end 1290 npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol 1291 if DX==0 1292 DX=abs(DCoord_min(NbDim)); 1293 end 1294 npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol 1295 for idim=1:NbDim 1296 if test_interp(idim) 1297 DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri 1298 end 1299 end 1300 Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY); 1301 test_direct_y=test_direct(NbDim-1); 1302 Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX); 1303 test_direct_x=test_direct(NbDim); 1304 DAX=DCoord_min(NbDim); 1305 DAY=DCoord_min(NbDim-1); 1306 minAX=min(Coord_x); 1307 maxAX=max(Coord_x); 1308 minAY=min(Coord_y); 1309 maxAY=max(Coord_y); 1310 xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1); 1311 ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2); 1312 xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame 1313 ycor_new=-xcorner*sin_om+ycorner*cos_om; 1314 if ~testXMax 1315 XMax=max(xcor_new); 1316 end 1317 if ~testXMin 1318 XMin=min(xcor_new); 1319 end 1320 if ~testYMax 1321 YMax=max(ycor_new); 1322 end 1323 if ~testYMin 1324 YMin=min(ycor_new); 1325 end 1326 DXinit=(maxAX-minAX)/(DimValue(NbDim)-1); 1327 DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1); 1328 if DX==0 1329 DX=DXinit; 1330 end 1331 if DY==0 1332 DY=DYinit; 1333 end 1334 if NbDim==3 1335 DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1); 1336 if ~test_direct(1) 1337 DZ=-DZ; 1338 end 1339 Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1)); 1340 test_direct_z=test_direct(1); 1341 end 1342 npX=floor((XMax-XMin)/DX+1); 1343 npY=floor((YMax-YMin)/DY+1); 1344 if test_direct_y 1345 coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line 1346 else 1347 coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line 1348 end 1349 if test_direct_x 1350 coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line 1351 else 1352 coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line 1353 end 1354 % case with no interpolation 1355 if isequal(ProjMode,'projection') && (~testangle || test90y || test90x) 1356 if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax 1357 ProjData=FieldData;% no change by projection 1358 else 1359 indY=NbDim-1; 1360 if test_direct(indY) 1361 min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1; 1362 max_indy=floor((YMax-Coord{indY}(1))/DYinit)+1; 1363 Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1); 1364 Ybound(2)=Coord{indY}(1)+DYinit*(max_indy-1); 1365 else 1366 min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1; 1367 max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1; 1368 Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1); 1369 Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1); 1370 end 1371 if test_direct(NbDim)==1 1372 min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1; 1373 max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1; 1374 Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); 1375 Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); 1376 else 1377 min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1; 1378 max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1; 1379 Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); 1380 Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); 1381 end 1382 min_indy=max(min_indy,1);% deals with margin (bound lower than the first index) 1383 min_indx=max(min_indx,1); 1384 1385 if test90y 1386 ind_new=[3 2 1]; 1387 DimCell={AYProjName,AXProjName}; 1388 % DimValue=DimValue(ind_new); 1389 iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1; 1452 eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1453 eval(['ProjData.' AXProjName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates 1454 end 1455 end 1456 else % case with rotation and/or interpolation 1457 if NbDim==2 %2D case 1458 [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates 1459 XIMA=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-Y*sin(PlaneAngle(3));%corresponding coordinates in the original image 1460 YIMA=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+Y*cos(PlaneAngle(3)); 1461 XIMA=(XIMA-minAX)/DXinit+1;% image index along x 1462 YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y 1463 XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' 1464 YIMA=reshape(round(YIMA),1,npX*npY); 1465 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 1466 if isequal(ObjectData.ProjMode,'filter') 1467 npx_filter=ceil(abs(DX/DAX)); 1468 npy_filter=ceil(abs(DY/DAY)); 1469 Mfilter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter); 1470 test_filter=1; 1471 else 1472 test_filter=0; 1473 end 1474 eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates 1475 eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates 1390 1476 for ivar=VarIndex 1391 1477 VarName=FieldData.ListVarName{ivar}; 1478 if test_interp(1) || test_interp(2)%interpolate on a regular grid 1479 eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST 1480 end 1481 %filter the field (image) if option 'filter' is used 1482 if test_filter 1483 Aclass=class(FieldData.A); 1484 eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');']) 1485 if ~isequal(Aclass,'double') 1486 eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values 1487 end 1488 end 1489 eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line 1490 %ind_in=find(flagin); 1491 ind_out=find(~flagin); 1492 ICOMB=(XIMA-1)*DimValue(1)+YIMA; 1493 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 1494 vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); 1495 for icolor=1:nbcolor 1496 vec_B(ind_out,icolor)=zeros(size(ind_out)); 1497 end 1392 1498 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1393 1499 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1394 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1395 eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new);'])% permute x and z indices for 90 degree rotation 1396 eval(['ProjData.' VarName '=squeeze(ProjData.' VarName '(iz,:,:));'])% select the z index iz 1397 end 1398 eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1399 eval(['ProjData.' AXProjName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates 1500 if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar 1501 ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; 1502 end 1503 eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']); 1504 end 1505 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 1506 ProjData.ListVarName=[ProjData.ListVarName 'FF']; 1507 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1508 ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; 1509 elseif ~testangle 1510 % unstructured z coordinate 1511 test_sup=(Coord{1}>=ObjectData.Coord(1,3)); 1512 iz_sup=find(test_sup); 1513 iz=iz_sup(1); 1514 if iz>=1 & iz<=npz 1515 %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; 1516 %ProjData.DimValue=[ProjData.DimValue npY npX]; 1517 for ivar=VarIndex 1518 VarName=FieldData.ListVarName{ivar}; 1519 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1520 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1521 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz 1522 %TODO : do a vertical average for a thick plane 1523 if test_interp(2) || test_interp(3) 1524 eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');']) 1525 end 1526 end 1527 end 1400 1528 else 1401 if NbDim==3 1402 DimCell(1)=[]; %suppress z variable 1403 DimValue(1)=[]; 1404 if test_direct(1) 1405 iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1; 1406 else 1407 iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1; 1408 end 1409 end 1410 max_indy=min(max_indy,DimValue(1));%introduce bounds in y and x indices 1411 max_indx=min(max_indx,DimValue(2)); 1412 for ivar=VarIndex% loop on non coordinate variables 1413 VarName=FieldData.ListVarName{ivar}; 1414 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1415 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1416 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar 1417 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; 1418 end 1419 if NbDim==3 1420 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']); 1421 else 1422 eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']); 1423 end 1424 end 1425 eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1426 eval(['ProjData.' AXProjName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates 1427 end 1428 end 1429 else % case with rotation and/or interpolation 1430 if NbDim==2 %2D case 1431 [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates 1432 XIMA=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-Y*sin(PlaneAngle(3));%corresponding coordinates in the original image 1433 YIMA=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+Y*cos(PlaneAngle(3)); 1434 XIMA=(XIMA-minAX)/DXinit+1;% image index along x 1435 YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y 1436 XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' 1437 YIMA=reshape(round(YIMA),1,npX*npY); 1438 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 1439 if isequal(ObjectData.ProjMode,'filter') 1440 npx_filter=ceil(abs(DX/DAX)); 1441 npy_filter=ceil(abs(DY/DAY)); 1442 Mfilter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter); 1443 test_filter=1; 1444 else 1445 test_filter=0; 1446 end 1447 eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates 1448 eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates 1449 for ivar=VarIndex 1450 VarName=FieldData.ListVarName{ivar}; 1451 if test_interp(1) || test_interp(2)%interpolate on a regular grid 1452 eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST 1453 end 1454 %filter the field (image) if option 'filter' is used 1455 if test_filter 1456 Aclass=class(FieldData.A); 1457 eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');']) 1458 if ~isequal(Aclass,'double') 1459 eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values 1460 end 1461 end 1462 eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line 1463 %ind_in=find(flagin); 1464 ind_out=find(~flagin); 1465 ICOMB=(XIMA-1)*DimValue(1)+YIMA; 1466 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 1467 vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); 1468 for icolor=1:nbcolor 1469 vec_B(ind_out,icolor)=zeros(size(ind_out)); 1470 end 1471 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1472 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1473 if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar 1474 ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; 1475 end 1476 eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']); 1477 end 1478 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 1479 ProjData.ListVarName=[ProjData.ListVarName 'FF']; 1480 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1481 ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; 1482 elseif ~testangle 1483 % unstructured z coordinate 1484 test_sup=(Coord{1}>=ObjectData.Coord(1,3)); 1485 iz_sup=find(test_sup); 1486 iz=iz_sup(1); 1487 if iz>=1 & iz<=npz 1488 %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; 1489 %ProjData.DimValue=[ProjData.DimValue npY npX]; 1490 for ivar=VarIndex 1491 VarName=FieldData.ListVarName{ivar}; 1492 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1493 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1494 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz 1495 %TODO : do a vertical average for a thick plane 1496 if test_interp(2) || test_interp(3) 1497 eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');']) 1498 end 1499 end 1500 end 1501 else 1502 errormsg='projection of structured coordinates on oblique plane not yet implemented'; 1503 %TODO: use interp3 1504 return 1505 end 1506 end 1529 errormsg='projection of structured coordinates on oblique plane not yet implemented'; 1530 %TODO: use interp3 1531 return 1532 end 1533 end 1507 1534 end 1508 1535
Note: See TracChangeset
for help on using the changeset viewer.