## Gridding and interpolate data

### joo tan (view profile)

on 8 Oct 2011
Latest activity Commented on by Dr. Seis

on 7 Jan 2015

### Dr. Seis (view profile)

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 (view profile)

on 8 Oct 2011
Edited by Dr. Seis

### Dr. Seis (view profile)

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

### Neil (view profile)

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

### Neil (view profile)

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.
Dr. Seis

### Dr. Seis (view profile)

on 7 Jan 2015
Valid point... and something I should implement !

### joo tan (view profile)

on 8 Oct 2011

i try to run but still curious how to input my data..i have non uniform coordinate (x,y) with height like (x,y,h)..so i need to grid the data..when i run, the matlab progrm always use the sample data

Dr. Seis

### Dr. Seis (view profile)

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.

### joo tan (view profile)

on 9 Oct 2011

sorry friends..what you meant the range? you meant the spacing of data? can i know, the value 0.5=dx and 1=dy is for what? the spacing? hbelow is example of my data
Latitude= 2.013534 2.071721 2.129907 2.188094 2.246266 2.304453 2.362639 2.420812 2.478998 2.537185 2.595372
Longitude= 124.715424 124.702625 124.689824 124.677023 124.664224 124.651421 124.638618 124.625817 124.613012 124.600206 124.5874
zonal= -138 -186 -191 -92 -20 64 32767 32767 32767 32767 32767
meridional= -227 -293 -335 -308 -295 -295 32767 32767 32767 32767 32767
so, i need to grid them to grid spacing 0.25 for latitude and longitude .Really need your help because i already try a few times

Dr. Seis

### Dr. Seis (view profile)

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));

### Neil (view profile)

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.
Best regards,
Neil

John D'Errico

### John D'Errico (view profile)

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

### Neil (view profile)

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