Interpolating to a known grid point

Hi,

point1 = (latitude, longitude) % in degrees
point2 = (latitude, longitude)
point3 = (latitude, longitude)
interpoint = (latitude, longitude)

point1,2,3 are stations where data is collected. point1,2,3 form a triangle and "interpoint" is located inside this triangle. "interpoint" is also a station.

P1 = [123, 345, 341; 986, 665, -645; 543, 86, -40; 0, 100, 0; .....cont.]
P2 = [-123, 381, 741; 986, 0, 645; 543, -86, 10; 12, 0, 54; ....cont.]
P3 = [723, 455, 224; -958, 0, 654; 743, 46, 11; 23, 0, 68; ....cont]

P4 are the values that I want to predict.

P1 are the values from point1 and P2 from point2 and so on.

As the output, I'm expecting the P4 matrix.

I experimented with interp2 and griddedInterpolant but the general error was assigning for e.g, point1 as the coordinates of P1 matrix (the whole thing) and expressing "interpoint" as the coordinates of the P4 values is to be interpolated.

Also with griddedInterpolant I get error, "Interpolation requires at least two sample points in each dimension" when running the following.

[X,Y] = ndgrid(lon1:lon3,lat3:lat1); % lon1 is the longitude of point1, etc..
V = P1;
F = griddedInterpolant(X,Y,V);
[Xq,Yq] = ndgrid(lat4,lon4);
Vq = F(Xq,Yq);
mesh(Xq,Yq,Vq);

scatteredInterpolant too ended in a disaster because I do not know how to call it to work. I'v been reading everything I can find on internet to learn how to get the output as a whole matrix of interpolated data, but no help.

Matlab keeps screaming that each cell needs a location but how can I tell it that all cells of P1 have the same location. To circumnavigate this issue I tried repmat location coordinates to the size of P1 but then the error was identical locations identified or something along those lines..

Anyways I think there IS a way but I can't seem to find it.

Alternatively , I followed a mathematical approach to calculate a single value at P4 and it's working.

Latitude and longitude converted to radians and then solve the following 3 equations (for a sample cell at row 8, column 4.

[S1, S2, S3] = solve(a0+(a1.*(lat1))+(a2.*(lon1)) == P1(8,4), a0+(a1.*(lat2))+(a2.*(lon2)) == P2(8,4), a0+(a1.*(lat3))+(a2.*(lon3)) == P3(8,4), a0, a1, a2);

SQ = S1+(S2.*(latq))+(S3.*(lonq)); % latq, lonq are the values from interpoint.

SQ is P4(8,4).

How do I get Matlab to run this in a loop and find every single value of SQ excluding where any of the input values are zero. i.e., if for a given cell the value in either P1, P2 or P3 is zero there should not be a result for that cell. Because 0 in my data set are artificial and thus the result would be wrong.

I don't like this method because the precision could get very low.

I really would like to compare this with Matlab's scatteredInterpolant or interp2 only if I know how to..

what's the best method to approach this?

Thanks!

[Jan]

It is not clear, how "time" comes into play. Is "interpoint" given or searched? What does this sentence mean:

 time for a particular point always remains the same P1 through P4 (e.g., P1(1,1) is at the same time as P2(1,1), P3(1,1) and P4(1,1).

[Geant Bepi]

time is the time at which each cell of data was generated.

yes all points are known including "interpoint".

so ideally when P4 line is created I would just have to match it up with the timescale.

[Geant Bepi]

I took time matrix off the question as it could probably have confused some of you as to what I wanted to be interpolated.

Answer 1

[Kelly Kearney]

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

[Geant Bepi]

Thank you Kelly, this is great! Interpolated values are off because "0" in the data set are artificial as mentioned above.

I replaced the 0 with NaN.

At what stage and how do I add the following;

if
isnan(P1) OR
isnan(P2) OR
isnan(P3) THEN
move to the next and *repeat*

Thanks heaps!

I tried the following;

out1 = isnan(P1);
out2 = isnan(P2);
out3 = isnan(P3);
outnan = or(out1,out2);
outlier = or(outnan,out3);
P = cat(3, P1, P2, P3);
P = P(~outlier);

but this is messing up the v1 matrix.

[Kelly Kearney]

Do you want to skip over just the individual time/variable combos that include a NaN in the points? Or skip an entire time slice? If the former, I'd just check during the loop:

P4 = nan(size(P1));

for it = 1:size(P,1)
    for ii = 1:size(P,2)
        if ~any(isnan(P(it,ii,:)))
            f.Values = permute(P(it,ii,:), [3 1 2]);
            P4(it,ii) = f(ltln(4,:));
        end
    end
end

Note that I also changed the preallocation of P4 to be NaNs instead of 0s, so entries that are skipped will remain NaNs in the result. The skipping isn't really necessary, though, since the interpolator will just return a NaN anyway if it finds one in the Values array; you might save a little time if there are a lot of NaNs in the dataset, but not much (since the more time-consuming part is the construction of the interpolant, not the evaluation of it).

[Geant Bepi]

did the trick! it's really strange how some things could be really simple yet so far from you!

Thanks heaps!

Answer 2

[Chad Greene]

You could fit a plane to the three points via threepoints2planez, then you'll be able to calculate P4 = z(x4,y4).

[Geant Bepi]

P4's coordinates are known. What I need is to interpolate to P4 using P1,P2,P3.