Gridding and interpolate data

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]

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!

Answer 1

[Dr. Seis]

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

[Robert]

I believe there is an error in initializing gg. Correct me if I am wrong.

for i = 1 : size(ZI,1) for j = 1 : size(ZI,2) gg = zeros(size(m)); 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); end end ZI(i,j) = sum(gg.*m); end end

Also, I found for a 3D face, gg(k)= magx; worked whereas gg(k) = (magx^2) * (log(magx)-1); did not.

[Neil]

Hi Elige,

   I tried your implementation of the biharmonic spline interpolation. Isn't a spline interpolation supposed to exactly reproduce the data? biharmonic_spline_interp2() seems to create a surface that approximates the data but doesn't go through all of them. Another thing I noticed is that MATLAB's own Curve Fitting Tool, when used with interpolant fitting method of biharmonicinterp, produces a surface that differs from biharmonic_spline_interp2. Do you know why?

Thanks,

Neil

[Dr. Seis]

The interpolation should exactly reproduce the data... in my implementation it did not because log of zero produced a NaN. I made an edit to the code above so that if it is trying to interpolate to an XI/YI location that is the same as an X/Y location, the point will be weighted more heavily and should now (to machine precision) reproduce that data point.

I do not know how/why Matlab's version differs from Sandwell's implementation.

[Neil]

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]

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]

Valid point... and something I should implement !

Answer 2

[joo tan]

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]

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.

Answer 3

[joo tan]

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]

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

Answer 4

[Neil]

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]

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]

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

Answer 5

[Christopher Gordon]

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?