surfergriddata.m

将matlab数据转化为surfer网格格式文件

function [Xi,Yi,Zi]=surfergriddata(X,Y,Z,varargin)
% Just like griddata but uses surfer instead. Default kringing
%
% Usage: 
% [|[Xi,Yi],Zi]=surfergriddata(X,Y,Z,|[Xi,Yi],|[method]) 
%
% Example:
%   X=randn(100,1);Y=randn(100,1);
%   [Xi,Yi,Zi]=surfergriddata(X,Y,peaks(X,Y));
%   surf(Xi,Yi,Zi)
%
% methods:
%   InverseDistance, Kriging(default), MinCurvature, NaturalNeighbor
%   NearestNeighbor, RadialBasis, Regression, Shepards, Triangulation
%   
%
% See also GRIDDATA, MESHGRID

% Aslak Grinsted 2006
% revised by degang 2006

if nargin<=4
    Xi=[];    Yi=[];
    if nargin==4,method=varargin{1};end
elseif nargin>=5
    Xi=varargin{1};Yi=varargin{2};
    if nargin==6,method=varargin{3};end
end

if isempty(Xi)
    NumCols=[];         NumRows=[];
    xMin=[];xMax=[];    yMin=[];yMax=[];
else
    if length(Xi(:))~=length(Yi(:)),[Xi,Yi]=meshgrid(Xi,Yi);end
    [msg,X,Y,Z,Xi,Yi] = xyzchk(X,Y,Z,Xi,Yi);
    if ~isempty(msg), error(msg); end    
    NumCols=(size(Xi,2));    NumRows=(size(Yi,1));
    xMin=min(Xi(:));         xMax=max(Xi(:));
    yMin=min(Yi(:));         yMax=max(Yi(:));
end

%methods:
methods={{'invdist','v4','inversedistance'},'kriging',{'mincurv','mincurvature'},...% 1 2 3
    'shepards',{'natural','naturalneighbor'},{'nearest','nearestneighbor'},...% 4 5 6
    'regression','radialbasis',{'linear','triangulation'},'movingaverage',...% 7 8 9 10
    {'local','localpolynomial'},'datametric'};% 11 12

if exist('method','var')
    method=optionchoose(method,methods);
else
    method=2;% srfKriging;
end

%some enum values:
srfGridFmtxyz = 4;
srfDupAvg = [];% Z value is equal to the midrange of the duplicate Z values

srf=actxserver('surfer.application');
try 
    temp=[X,Y,Z];
    save('~surferdata_temp.dat','temp','-ascii');
    srf.DefaultFilePath = cd;    

   b= srf.GridData('~surferdata_temp.dat', [], [], [], [], ...%  DataFile, xCol, yCol, zCol, ExclusionFilter, 
       srfDupAvg, [], [], ... %DupMethod, xDupTol, yDupTol, 
       NumCols, NumRows, xMin, xMax, yMin, yMax,...%NumCols, NumRows, xMin, xMax, yMin, yMax,
       method, 0, [], [], ... % Algorithm, ShowReport, SearchEnable, SearchNumSectors, 
       [], [], [], [], [], [], ... %SearchRad1, SearchRad2, SearchAngle, SearchMinData, SearchDataPerSect, SearchMaxEmpty,
       [], [], [], [],... % FaultFileName, BreakFileName, AnisotropyRatio, AnisotropyAngle, 
       [], [], [], [], [], [], ... %IDPower, IDSmoothing, KrigType, KrigDriftType, KrigStdDevGrid, KrigVariogram, 
       [], [], [], [], [],... %MCMaxResidual, MCMaxIterations, MCInternalTension, MCBoundaryTension, MCRelaxationFactor,
       [], [], [], [], [], ... %ShepSmoothFactor, ShepQuadraticNeighbors, ShepWeightingNeighbors, ShepRange1, ShepRange2, 
       [], [], [], [], [], ... %RegrMaxXOrder, RegrMaxYOrder, RegrMaxTotalOrder, RBBasisType, RBRSquared, 
       '~surfergrid_temp.dat', srfGridFmtxyz, ...    %OutGrid, OutFmt )
       [], [], [], [], [],[]);% SearchMaxData, KrigStdDevFormat, DataMetric, LocalPolyOrder, LocalPolyPower, TriangleFileName
    if ~b
        warning('something went wrong in srf.griddata')
    end    

    delete('~surferdata_temp.dat');
    srf.Quit;
    srf.delete;
catch
    err=lasterror;
    srf.Quit;
    srf.delete;
    rethrow(err);
end

temp=load('~surfergrid_temp.dat');
if isempty(NumCols),
    NumRows=sum(temp(:,1)-temp(1)==0);NumCols=size(temp,1)/NumRows;
end
Xi=reshape(temp(:,1),[NumCols, NumRows])';
Yi=reshape(temp(:,2),[NumCols, NumRows])';
Zi=reshape(temp(:,3),[NumCols, NumRows])';
if nargout==1
    Xi=Zi;    clear Zi;
end