101 views (last 30 days)

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

Dr. Seis
on 8 Oct 2011

Edited: Dr. Seis
on 30 Dec 2014

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

Neil
on 30 Dec 2014

Hi Elige. Thanks for the fix. Comparing your implementation of Sandwell's algorithm with matlab's biharmonic spline in Curve Fitting Tool (aka 'v4' in griddata), your implementation bridges over large gaps in the data more than matlab's; I guess the surface looks like it has more tension. GMT's implementation of biharmonic spline interpolation has a tension parameter, because it's based on Wessel's extension to Sandwell's idea, but neither your nor matlab's implementation uses tension, so the difference is unexpected. Both biharmonic spline implementations bridge over gaps in the data more than gridfit does, however, which is apparently also to be expected from biharmonic interpolations.

Regards,

Neil

Neil
on 6 Jan 2015

Hi Elige. I found out the reason for the discrepancy. When I was using biharmonic spline in the Curve Fitting Tool, I had checked the Center and Scale parameter (same as doing fit(...,'biharmonicinterp','Normalize','on'). This makes a difference in the case of my data. When I unchecked the Center and Scale parameter, the fit was practically the same as the one done by your implementation in biharmonic_spline_interp2().

I should point out that at least in the case of my data, centering and scaling the data seems to have been necessary prior to calling the fitting algorithm. I was already log transforming the data prior to passing it to the fitting functions, but that didn't remove the need for centering and scaling. Anyone using biharmonic_spline_interp2() should be aware that it doesn't center and scale the data, and consider doing it themselves if necessary.

Sign in to comment.

joo tan
on 8 Oct 2011

Dr. Seis
on 8 Oct 2011

The program needs five input arguments. The first three represent the non-uniform points that you know, the last two represent the uniform x, y grid you want to interpolate to. Let's say you have x data with range 300 to 400 and y data with range 600 to 950. Then you would construct the last two input arguments for the program by running, say:

dx = 0.5;

dy = 1;

[XI, YI] = meshgrid(300:dx:400, 600:dy:950);

Then run the program:

HI = biharmonic_spline_interp2(x,y,h,XI,YI);

HI will be a matrix similar to what you would get out of "interp2" program.

Sign in to comment.

joo tan
on 9 Oct 2011

Dr. Seis
on 9 Oct 2011

This program isn't really meant for plotting on a spherical surface, unless you are looking at a small enough region. 1. Your latitude range would be from min(Latitude) to max(Latitude) and your longitude range would be from min(Longitude) to max(Longitude). 2. If you want to create a latitude grid increment (dLat) and a longitude grid increment (dLon), then you decide how many grid points you want. If you want to divide your Latitude into 101 points and your Longitude into 201 points, then you would:

dLon = (max(Longitude)-min(Longitude))/200;

dLat = (max(Latitude)-min(Latitude))/100;

[XI, YI] = meshgrid(min(Longitude):dLon:max(Longitude) , min(Latitude):dLat:max(Latitude));

Sign in to comment.

Neil
on 13 Aug 2014

Thanks for the implementation of the Sandwell algorithm. Is there a way to adapt this implementation to get it to work on larger sets of input data? I have a 600x600 pixel image, with many holes in it, and I tried to fill in the holes using biharmonic spline interpolation. However the part of the code where it uses

gg = zeros(length(Z),length(Z));

blows up because length(Z) is 124504 (in the example it's 100). I could break up the image into smaller parts and process each individually but I'm not sure what would happen at the boundaries when I put the reconstructed sub-images back together.

I would appreciate your feedback.

Best regards,

Neil

John D'Errico
on 26 Dec 2014

I would point out that you got no response, because you added an ANSWER! How can you possibly expect that someone will know that you had a question and not a real answer?

In the future, if you have a question, then ask it, as a question.

Neil
on 26 Dec 2014

You're right, John. I was trying to ask a question by commenting on E. Grant's original answer, but used the wrong link. (I think I did that correctly last night when I asked another question about biharmonic_spline_interp2.)

For what it's worth, since I asked this question, I've realized that biharmonic spline interpolation is not the best approach for fitting a surface to a very large number of data, unless one would like to incorporate slope data, which I don't; also, when there are error bars, forcing the fit to exactly reproduce the data is not a good idea. (It makes more sense to use something like gridfit in that case.) But I'm curious to know why biharmonic_spline_interp2 behaves differently from MATLAB's own cftool when Model='biharmonicinterp' or griddata when Method='v4'.

Neil

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.