MATLAB: Gridding and interpolate data

griddinginterpolatekrigging

Dear friends,
I have a lot of data and i need to interpolate and grid the data..I plan to use krigging method to grid the data.but my problem, the data are not uniform..Can someone suggest to me or example matlab program to solve my problem.
Thanks

Best Answer

You could use the formulation found in the following reference:
Sandwell, D. T. (1987), Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data, Geophysical Research Letters, Vol. 2, p. 139 – 142.
The function below can take and interpolate data collected on an irregularly spaced grid and output the result on a regularly spaced grid. It is setup similarly to "interp2" except the input X, Y, and Z points are in column vectors. The XI and YI define the desired regular grid spacing and can be constructed using "meshgrid" before running.
To see an example using the "Peaks" data set, run the code from the command line without any input arguments.
function ZI = biharmonic_spline_interp2(X,Y,Z,XI,YI)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   2D Biharmonic spline interpolation implemented from:
%
%   Sandwell, D. T. (1987), Biharmonic spline interpolation of GEOS-3 and
%       SEASAT altimeter data, Geophysical Research Letters, Vol. 2, 
%       p. 139 – 142.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Run an example if no input arguments are found
if nargin ~= 5
    fprintf('Running Peaks Example \n');
    X = rand(100,1)*6 - 3;
    Y = rand(100,1)*6 - 3;
    Z = peaks(X,Y);
    [XI,YI] = meshgrid(-3:.25:3);
end
%   Check to make sure sizes of input arguments are correct/consistent
if size(Z,1) < size(Z,2)
    error('X, Y and Z should all be column vectors !');
end
if (size(Z,1) ~= size(X,1)) || (size(Z,1) ~= size(Y,1))
    error('Length of X, Y, and Z must be equal !');
end
if (size(XI,1) ~= size(YI,1)) || (size(XI,2) ~= size(YI,2))
    error('Size of XI and YI must be equal !');
end
%   Initialize output
ZI = zeros(size(XI));
%   Compute GG matrix for GG*m = d inversion problem
GG = zeros(length(Z),length(Z));
for i = 1 : length(Z)
    for j = 1 : length(Z)
        if i ~= j
            magx = sqrt((X(i)-X(j))^2 + (Y(i)-Y(j))^2);
            if magx >= 1e-7
                GG(i,j) = (magx^2) * (log(magx)-1);
            end
        end
    end
end
%   Compute model "m" where data "d" is equal to "Z"
m = GG\Z;
%   Find 2D interpolated surface through irregular/regular X, Y grid points 
gg = zeros(size(m));
for i = 1 : size(ZI,1)
    for j = 1 : size(ZI,2)
        for k = 1 : length(Z)
            magx = sqrt((XI(i,j)-X(k))^2 + (YI(i,j)-Y(k))^2);
            if magx >= 1e-7
                gg(k) = (magx^2) * (log(magx)-1);
            else
                gg(k) = (magx^2) *(-100);
            end 
        end
        ZI(i,j) = sum(gg.*m);
    end
end
%   Plot result if running example or if no output arguments are found
if nargin ~= 5 || nargout ~= 1
    figure;
    subplot(3,1,1);
        [XE,YE] = peaks(-3:.25:3);
        ZE = peaks(XE,YE);
        mesh(XE,YE,ZE); hold on;
        scatter3(X,Y,Z,'filled'); hold off;
        title('Peaks');
        axis([-3 3 -3 3 -5 5]);
        caxis([-5 5]); colorbar;
    subplot(3,1,2);
        mesh(XI,YI,ZI); hold on;
        scatter3(X,Y,Z,'filled'); hold off;
        title('Peaks Interpolated'); 
        axis([-3 3 -3 3 -5 5]);   
        caxis([-5 5]); colorbar;
    subplot(3,1,3);
        mesh(XI,YI,ZI-ZE); hold on;
        scatter3(X,Y,Z-Z,'filled'); hold off;
        title('Peaks Interpolated Difference'); 
        axis([-3 3 -3 3 -5 5]);  
        caxis([-5 5]); colorbar; 
end
Related Question