2d interpolation from a non-uniform grid

I have looked at the griddata function but am not clear if it applies to data on an irregular mesh. In my case, I have a 2d spatial dataset available on an [X,Y] grid where X and Y, which are the cartesian coordinates, are themselves 2d arrays. That is, both X and Y have sizes (m,n) but

X(:,1) ~= X(:,2) ~= X(:,3) and so on.

Similarly for the array Y. Visually, you can picturize the grid as follows, where X marks the location of the input data set:

X X X X

X X X X

X X X X

Along each row, the spacing between two Xs is uniform. But the spacing changes from row to row, leading to X(:,1), X(:,2), etc. all being different.

It looks like griddata will not accept such an input grid. Am I correct?

Thanks for your help.

[Akira Agata]

In that case, the data will not be treated as meshgrid. I belive there are 2 possible solutions:

Solution 1 (I think is better)

Using scatteredInterpolant function and conduct 2D interpolation.

Solution 2

Extract half of your data ( X(1:2:end, 1:2:end) and Y(1:2:end, 1:2:end) ) and apply interp2 function.

[oceanmod]

Thanks, I will try the first option and see how it goes.

I am not sure I understand the logic of the second option. Even if I extract half the data, I am still going to have a data set with irregular spacings, am I not?

[oceanmod]

Hmm... scatteredInterpolant is not working either as that is also expecting "a double array" but I am trying to feed in 2-d arrays for the first and second arguments (x and y below) in

F = scatteredInterpolant(x,y,V)

There must be a way to do this in Matlab:) I am sure I am not the first to have data on a grid where both the x and y coordinate vary with (i,j).

[Akira Agata]

Well, let me clarify my intention by the following small example.

I hope this will be somehow helpful for you to understand how these two functions works.

% Original Data (X and Y are non-uniform grid)
X = [...
  1 3 5 7;...
  2 4 6 8;...
  1 3 5 7;...
  2 4 6 8];

Y = [...
  1 1 1 1;
  2 2 2 2;
  3 3 3 3;
  4 4 4 4];

Z = X.^2 + Y.^2;

% Query points
[xGrid, yGrid] = meshgrid(1:0.2:8,1:0.2:4);

% [Solution 1]
F = scatteredInterpolant(X(:),Y(:),Z(:));
zGrid = F(xGrid,yGrid);

% [Solution 2]
X2 = X(1:2:end,:);
Y2 = Y(1:2:end,:);
Z2 = Z(1:2:end,:);

zGrid2 = interp2(X2,Y2,Z2,xGrid,yGrid);

Answer 1

[Bjorn Gustavsson]

No you can use griddata and scatteredInterpolant. You can do something like this:

Zi = griddata(X(:),Y(:),Z(:),Xi,Yi);

And you do the same thing with scatteredInterpolant - the (:) construct just unwraps an array into a 1-D column array.

HTH

[oceanmod]

Thanks for your response. Griddata does not work, even if you use the colon construct to expand out a 2-d array into a 1-d array. That is because if the size of the field to be interpolated is (m,n), griddata is expecting arrays of size m and size n for the first and second arguments in the call to griddata. In my case, the input coordinates themselves are 2d arrays and of the form x(i,j) and y(i,j). So, if I convert them to 1-d arrays, each will have size m x n, which is not what griddata is expecting.

[Bjorn Gustavsson]

Yes it does:

[X,Y] = meshgrid(1:12,1:11);
X(2:2:end,:) = X(2:2:end,:)+.5;
plot(X,Y,'+')
Z = peaks(12);
Z = Z(2:end,:);
subplot(1,2,1),
pcolor(X,Y,Z)
[Xi,Yi] = meshgrid(1:.1:12,1:.1:11);
Zi = griddata(X(:),Y(:),Z(:),Xi,Yi);
subplot(1,2,2)
pcolor(Xi,Yi,Zi),shading flat

You have to give the function the input it expects, for x, y and z. Then it demonstrably works.

[oceanmod]

Thanks for persisting! I understood what I was doing wrong. Once I use meshgrid to create [X,Y] at the beginning, I can then populate the values as I want, be they regularly spaced or not (as illustrated in your example). This is what I am doing now:

% xorig, yorig are the original 2d arrays describing the cartesian coordinates
[X,Y]=meshgrid(size(xorig,1),size(xorig,2));

for i=1:size(xorig,1);
  for j=1:size(xorig,2);
     X(j,i) = xorig(i,j);
     Y(j,i) = yorig(i,j);
  end
end

% x, y are the new locations where I want the interpolated values
[x,y] = meshgrid(-5:0.1:5, 5:0.1:5);

% The interpolation step (need to transpose the original field)
var_new = griddata(X,Y,var_old',x,y,'cubic');

I found that the mean value is changing at the second decimal after interpolation but that amount of interpolation error might be unavoidable.