84 views (last 30 days)

Show older comments

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.

Thanks in advance!

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

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.

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

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

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?

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!