The scatteredInterpolant function will handle this; with a linear interpolation scheme, it basically does the plane-fitting calculation that Chad suggested.
I'm assuming that the column dimension of your P data represents separate variables to be interpolated. So start with your data (I just chose four random points, with the 4th one inside the triangle of the first 3):
ltln = [...
49.511 -120.29
46.292 -124.23
45.261 -116.42
47.54 -120.3];
P1 = [123, 345, 341; 986, 665, -645; 543, 86, -40; 0, 100, 0];
P2 = [-123, 381, 741; 986, 0, 645; 543, -86, 10; 12, 0, 54];
P3 = [723, 455, 224; -958, 0, 654; 743, 46, 11; 23, 0, 68];
P = cat(3, P1, P2, P3);
Step 1 is to set up the interpolant geometry using a single column-variable and a single point in time across the three locations:
v1 = permute(P(1,1,:), [3 1 2]);
f = scatteredInterpolant(ltln(1:3,:), v1);
Now you can reuse that interpolant, just switching out the Values property for each relevant time and variable:
P4 = zeros(size(P1));
for it = 1:size(P,1)
for ii = 1:size(P,2)
f.Values = permute(P(it,ii,:), [3 1 2]);
P4(it,ii) = f(ltln(4,:));
end
end
Best Answer