Index: /issm/trunk/src/m/utils/DataProcessing/CreateDataBoundaries.m =================================================================== --- /issm/trunk/src/m/utils/DataProcessing/CreateDataBoundaries.m (revision 1079) +++ /issm/trunk/src/m/utils/DataProcessing/CreateDataBoundaries.m (revision 1079) @@ -0,0 +1,42 @@ +function [Xedge,Yedge,EdgeValues]=CreateDataBoundaries(x_m,y_m,x_matrix,y_matrix,datamatrix); +%CREATEDATABOUNDARIES - Create values on the edge of the matrix +% +% This function create vectors with coordinates and values to constrain a matrix +% on its edge from a given matrix tiwh data everywhere. +% x_m and y_m belongs to the matrix to be constrained +% x_matrix, y_matrix belongs to the matrix to be constrained with +% datamatrix is the matrix with values for the constraint +% +% Usage: +% [Xedge,Yedge,EdgeValues]=CreateDataBoundaries(x_m,y_m,x_matrix,y_matrix,datamatrix); +% +% See also: TRACKSTOMATRIX, CREATEDATAMATRIX + +%Check the consistency of the data matrix +if(length(x_matrix)~=(size(datamatrix,2)+1) | length(y_matrix)~=(size(datamatrix,1)+1)), + error('CreateDataBoundaries error message: size of matrix and vectors to constrained not consistent') +end + +nxglobal=length(x_m); +nyglobal=length(y_m); + +%Create the edgevalues and coordinates +xedge1=x_m(1)*ones(nyglobal,1); +xedge2=x_m(2:end-1); +xedge3=x_m(end)*ones(nyglobal,1); +xedge4=x_m(2:end-1); +xedge5=x_m(2)*ones(nyglobal-2,1); +xedge6=x_m(3:end-2); +xedge7=x_m(end-1)*ones(nyglobal-2,1); +xedge8=x_m(3:end-2); +yedge1=y_m; +yedge2=y_m(1)*ones(nxglobal-2,1); +yedge3=y_m; +yedge4=y_m(end)*ones(nxglobal-2,1); +yedge5=y_m(2:end-1); +yedge6=y_m(2)*ones(nxglobal-4,1); +yedge7=y_m(2:end-1); +yedge8=y_m(end-1)*ones(nxglobal-4,1); +Xedge=[xedge1;xedge2;xedge3;xedge4;xedge5;xedge6;xedge7;xedge8]; +Yedge=[yedge1;yedge2;yedge3;yedge4;yedge5;yedge6;yedge7;yedge8]; +EdgeValues=DataInterp(x_matrix,y_matrix,datamatrix,Xedge,Yedge); Index: /issm/trunk/src/m/utils/DataProcessing/CreateDataMatrix.m =================================================================== --- /issm/trunk/src/m/utils/DataProcessing/CreateDataMatrix.m (revision 1079) +++ /issm/trunk/src/m/utils/DataProcessing/CreateDataMatrix.m (revision 1079) @@ -0,0 +1,90 @@ +function [Mvalue Mx My]=CreateDataMatrix(x_m,y_m,track_coord,track_values), +%CREATEDATAMATRIX - Create a map with average values of map +% +% This routine creates a map with average values of tracks. +% x_m1 and y_m1 are two vectors containing the coordinates of the matrix +% trac_coord is an exp file containing the coordinates of the tracks (x and y) +% trav_values is a vector with the values along the track coordinates +% +% Usage: +% [Mvalue Mx My]=CreateDataMatrix(x_m,y_m,track_coord,track_values), +% +% Example: +% [Mvalue Mx My]=CreateDataMatrix(x_m,y_m,'trackcoord.exp',thickness_track) +% +% See also: CREATEDATABOUNDARIES, TRACKSTOMATRIX + +%Read the points of the tracks +stru=expread(track_coord,1); +nods=stru.nods; +xtracks=stru.x'; +ytracks=stru.y'; + +%First check that the parameters are ok: +if (size(track_values,1)~=nods) || (size(xtracks,2)~=nods) || (size(ytracks,2)~=nods), + error('CreateDataMatrix error message : track coordinates and track values must have the same size'); +end + +%Compute number of rows and columns +numrow=size(y_m,1)-1; +numcol=size(x_m,1)-1; + +%Remove useless points of the track +points=find(track_values==0); +track_values(points)=[]; +xtracks(points)=[]; +ytracks(points)=[]; +points=find(isnan(track_values)); +track_values(points)=[]; +xtracks(points)=[]; +ytracks(points)=[]; + +points=find(xtracksx_m(end) | ytracksy_m(end)); +track_values(points)=[]; +xtracks(points)=[]; +ytracks(points)=[]; + +%initialize some matrices +numpoints=zeros(numrow,numcol); +value=zeros(numrow,numcol); +coordx=zeros(numrow,numcol); +coordy=zeros(numrow,numcol); + +%Loop over the points of the track +nel=size(track_values,1); +fprintf('%s',' track processing progress: 0.00 %'); +for i=1:nel; + if mod(i,1000)==0, + fprintf('\b\b\b\b\b\b\b') + fprintf('%5.2f%s',i/nel*100,' %'); + end + + x=xtracks(i); + y=ytracks(i); + + %get indices for the matrix + indexx=max(find(x_m1000, + fprintf('\b\b\b\b\b\b\b\b') + fprintf('%4.2f%s\n',100,' %'); +end + +%Change the values of numoverlap to 1 if 0 since we are going to devide by this matrix +numpoints(find(~numpoints))=1; + +%Create the center of mass for coordiantes and values. +Mvalue=value./numpoints; +Mx=coordx./numpoints; +My=coordy./numpoints; Index: /issm/trunk/src/m/utils/DataProcessing/TracksToMatrix.m =================================================================== --- /issm/trunk/src/m/utils/DataProcessing/TracksToMatrix.m (revision 1079) +++ /issm/trunk/src/m/utils/DataProcessing/TracksToMatrix.m (revision 1079) @@ -0,0 +1,69 @@ +function [x_m y_m MatData MatTracks]=TracksToMatrix(xmin,ymin,nx,ny,posting,track_coord,track_values,theta,varargin), +%TRACKSTOMATRIX - Create a map from tracks +% +% This routine creates a map with average values of tracks using the DACE toolbox. +% xmin and ymin are two scalars are the extreme values of the data matrix +% trac_coord is an exp file containing the coordinates of the tracks (x and y) +% trav_values is a vector with the values along the track coordinates +% posting is a scalar giving the posting of the matrix in meters +% theta is a parameter of the correlation function +% +% Usage: +% [x_m y_m MatData]=TracksToMatrix(xmin,ymin,nx,ny,posting,track_coord,track_values), +% +% Example: +% [x_m y_m Thickness]=TracksToMatrix(-10^6,!0^6,200,180,1000,'trackcoord.exp',thickness_values), +% +% See also: CREATEDATABOUNDARIES, CREATEDATAMATRIX + +if nargin~=8 & nargin~=11, + error(' TracksToMatrix error message: wrong number of arguments') +end + +%First create the x_m and y_m fot the matrix +x_f=[xmin:posting:xmin+posting*nx]'; +y_f=[ymin:posting:ymin+posting*ny]'; + +%Now create a bigger map we will then truncate +x_m=[xmin-posting*nx/2:posting:xmin+posting*nx*3/2]'; +y_m=[ymin-posting*ny/2:posting:ymin+posting*ny*3/2]'; + +%Create DataMatrix with local averaging of tracks +[Mvalue Mx My]=CreateDataMatrix(x_m,y_m,track_coord,track_values); + +MatTracks=Mvalue; +%Create vector for these coordinates and values +Svalue=sparse(Mvalue); +Sx=sparse(Mx); +Sy=sparse(My); +[i,j,Values]=find(Svalue); +[i,j,X]=find(Sx); +[i,j,Y]=find(Sy); +Coord=[X,Y]; + +%Create boundaries for the track if specified +if nargin==11, + %Find out the elements to counstrain the border + x_matrix=varargin{1}; + y_matrix=varargin{2}; + datamatrix=varargin{3}; + + %Create values on the border + [Xedge,Yedge,EdgeValues]=CreateDataBoundaries(x_m,y_m,x_matrix,y_matrix,datamatrix); + + %Add these values to the track values + Values=[Values;EdgeValues]; + Coord=[X,Y;Xedge,Yedge]; +end + +%Create model for data +[dmodel,perf]=dacefit(Coord,Values,@regpoly1,@corrgauss,theta); + +%Create design site(points where we are looking for the data) +Points=gridsamp([x_f(1)+posting/2,y_f(1)+posting/2;x_f(end)-posting/2,y_f(end)-posting/2],[length(x_f)-1;length(y_f)-1]); + +%Compute data on these points +VecData=predictor(Points,dmodel); + +%Reshape to get a matrix +MatData=reshape(VecData,ny,nx);