Index: /issm/trunk-jpl/src/m/utils/DataProcessing/gamv.m =================================================================== --- /issm/trunk-jpl/src/m/utils/DataProcessing/gamv.m (revision 12302) +++ /issm/trunk-jpl/src/m/utils/DataProcessing/gamv.m (revision 12302) @@ -0,0 +1,42 @@ +function output = gamv(x,y,data,varargin); +%GAMV - use gslib for Kriging +% +% Usage: +% output = gamv(x,y,data,varargin) + +%Write data file +fid=fopen('cluster.dat','w'); +fprintf(fid,'%s\n','Data file'); +fprintf(fid,'%i\n',3); +fprintf(fid,'%s\n','Xlocation'); +fprintf(fid,'%s\n','Ylocation'); +fprintf(fid,'%s\n','Data'); +fprintf(fid,'%g %g %g\n',[x y data]'); +fclose(fid); + +%Write parameter file +fid=fopen('gamv.par','w'); +fprintf(fid,'\t\t\t\t%s\n','Parameters for GAMV'); +fprintf(fid,'\t\t\t\t%s\n','*******************'); +fprintf(fid,'\n'); +fprintf(fid,'%s\n','START OF PARAMETERS:'); +fprintf(fid,'%-30s %s\n','./cluster.dat' ,'\file with data'); +fprintf(fid,'%-30s %s\n','1 2 0' ,'\columns for X, Y, Z coordinates'); +fprintf(fid,'%-30s %s\n','1 3 ' ,'\number of variables, column number'); +fprintf(fid,'%-30s %s\n','-1.0e21 1.0e21' ,'\trimming limits'); +fprintf(fid,'%-30s %s\n','gamv.out' ,'\file for variogram output'); +fprintf(fid,'%-30s %s\n','20' ,'\number of lags'); +fprintf(fid,'%-30s %s\n','5.0' ,'\lag separation distance'); +fprintf(fid,'%-30s %s\n','3.0' ,'\lag tolerance'); +fprintf(fid,'%-30s %s\n','3' ,'\number of directions'); +fprintf(fid,'%-30s %s\n','0.0 90.0 50.0 0.0 90.0 50.0','\azm, atol, bandh, dip, dtol, bandv'); +fprintf(fid,'%-30s %s\n','0.0 22.5 25.0 0.0 22.5 25.0','\azm, atol, bandh, dip, dtol, bandv'); +fprintf(fid,'%-30s %s\n','90. 22.5 25.0 0.0 22.5 25.0','\azm, atol, bandh, dip, dtol, bandv'); +fprintf(fid,'%-30s %s\n','0' ,'\standardize sill? (0=no, 1=yes)'); +fprintf(fid,'%-30s %s\n','2' ,'\number of variograms'); +fprintf(fid,'%-30s %s\n','1 1 1' ,'\tail var., head vars., variogram type'); +fprintf(fid,'%-30s %s\n','1 1 3' ,'\tail var., head vars., variogram type'); +fclose(fid); + +%Call gamv +system([issmdir() '/externalpackages/gslib/install/gamv gamv.par']); Index: /issm/trunk-jpl/src/m/utils/DataProcessing/gslib.m =================================================================== --- /issm/trunk-jpl/src/m/utils/DataProcessing/gslib.m (revision 12301) +++ /issm/trunk-jpl/src/m/utils/DataProcessing/gslib.m (revision 12302) @@ -1,25 +1,7 @@ -x=[]; -y=[]; -thickness=[]; -for i=1:5, - switch(i) - case 1, filename='2002_Antarctica_P3chile'; - case 2, filename='2004_Antarctica_P3chile'; - case 3, filename='2009_Antarctica_DC8'; - case 4, filename='2009_Antarctica_TO'; - case 5, filename='2010_Antarctica_DC8'; - otherwise, error('not supported'); - end - disp(['Loading ' filename '...']); - A=load(filename); - [tempx tempy]=ll2xy(A.lat(:),A.lon(:),-1,0,71); - x=[x;tempx]; - y=[y;tempy]; - thickness=[thickness;A.thickness]; -end -[xl yl]=basinzoom(pairoptions('basin','Pine Island Glacier')); -x=x(1:100:end); -y=y(1:100:end); -data=thickness(1:100:end); +function output = gslib(x,y,data,varargin); +%GSLIB - use gslib for Kriging +% +% Usage: +% output = gslib(x,y,data,varargin) %Output Matrix @@ -44,5 +26,4 @@ c = (sill-nugget); a = sqrt(3)*range; - %Write data file @@ -108,4 +89,5 @@ tic;system([issmdir() '/externalpackages/gslib/install/kb2d kb2d.par']);toc; + delete('kb2d.par'); %Read output Index: /issm/trunk-jpl/src/m/utils/DataProcessing/varmap.m =================================================================== --- /issm/trunk-jpl/src/m/utils/DataProcessing/varmap.m (revision 12302) +++ /issm/trunk-jpl/src/m/utils/DataProcessing/varmap.m (revision 12302) @@ -0,0 +1,55 @@ +function output = varmap(x,y,data,varargin); +%VARMAP - use gslib for Kriging +% +% Usage: +% output = varmap(x,y,data,varargin) + +options=pairoptions(varargin{:}); + +nxlag = getfieldvalue(options,'nxlag', 20); +nylag = getfieldvalue(options,'nylag', 20); +dxlag = getfieldvalue(options,'dxlag', 1000); +dylag = getfieldvalue(options,'dylag', 1000); + +%Write data file +fid=fopen('cluster.dat','w'); +fprintf(fid,'%s\n','Data file'); +fprintf(fid,'%i\n',3); +fprintf(fid,'%s\n','Xlocation'); +fprintf(fid,'%s\n','Ylocation'); +fprintf(fid,'%s\n','Data'); +fprintf(fid,'%g %g %g\n',[x y data]'); +fclose(fid); + +%Write parameter file +fid=fopen('varmap.par','w'); +fprintf(fid,'\t\t\t\t%s\n','Parameters for GAMV'); +fprintf(fid,'\t\t\t\t%s\n','*******************'); +fprintf(fid,'\n'); +fprintf(fid,'%s\n','START OF PARAMETERS:'); +fprintf(fid,'%-30s %s\n','./cluster.dat' ,'\file with data'); +fprintf(fid,'%-30s %s\n','1 3 ' ,'\number of variables, column number'); +fprintf(fid,'%-30s %s\n','-1.0e21 1.0e21' ,'\trimming limits'); +fprintf(fid,'%-30s %s\n','0 ' ,'\1=regular grid, 0=scattered values'); +fprintf(fid,'%-30s %s\n','50 50 1' ,'\if =1: nx, ny, nz'); +fprintf(fid,'%-30s %s\n','1.0 1.0 1.0' ,'\ xsiz, ysiz, zsiz if igrid=1'); +fprintf(fid,'%-30s %s\n','1 2 0' ,'\if =0: columns for x, y and z coordinates'); +fprintf(fid,'%-30s %s\n','varmap.out' ,'\file for variogram output'); +fprintf(fid,'%-30s %s\n',num2str([nxlag nylag 0],'%i '),'\nxlag, nylag, nzlag'); +fprintf(fid,'%-30s %s\n',num2str([dxlag dylag 1],'%g %g %i'),'\dxlag, dylag, dzlag'); +fprintf(fid,'%-30s %s\n','5' ,'\minimum number of pairs'); +fprintf(fid,'%-30s %s\n','0' ,'\standardize sill? (0=no, 1=yes)'); +fprintf(fid,'%-30s %s\n','1' ,'\number of variograms'); +fprintf(fid,'%-30s %s\n','1 1 1' ,'\tail, head, variogram type'); +fclose(fid); + +%Call varmap +system([issmdir() '/externalpackages/gslib/install/varmap varmap.par']); +delete('varmap.par'); + +%Read output +fid=fopen('varmap.out','r'); +A = textscan(fid,'%f %f %f %f %f %f','headerlines',8); +fclose(fid); +delete('varmap.out') +output = reshape(A{1},[2*nxlag+1 2*nylag+1]);