# Gridding and interpolate data

50 views (last 30 days)
joo tan on 8 Oct 2011
Answered: Christopher Gordon on 19 Feb 2021
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
Taindra Neupane on 1 Oct 2020
Hello Friends,
I have a 2d uniformly gridded data (like an image or dose with specific voxel size), I am looking for interpolation funtion to divide into +ve interger as well as non-integer values like 1, 1.5,2....! I am using interpn (A,k) or interpn2 (A,k) but its interpolating only into even +ve intergers points/values like k=0, no interpolation, k=1, two values, k=2, 4 values as it divide into 2^k -1 values. So, I am wondering if I want to divide 3mmx3mm pixel into 1mmx1mm pixel (3 sections). Please help me on this.

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
Dr. Seis on 7 Jan 2015
Valid point... and something I should implement !

joo tan 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 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 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 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 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
##### 2 CommentsShowHide 1 older comment
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

Christopher Gordon on 19 Feb 2021
Good afternoon!
It seems as though you can't run this function unless all the initial dimensions (X,Y,Z) are equal. I have a 2D matrix of spectra, (30 spectra, each being 1024 data points long) and I'm trying to interpolate it to be a 1024 x 1024 matrix. Is there a way to do that with this, or do I have to use another approach entirely?