How to interpolate 3d data based on changes in 3rd dimension?

I have a 3d matrix A of gridded data with dimensions lat, lon, and time (size 180x360x135). I want to create a new matrix B from this using linear interpolation where the size of the 3rd dimension (time) has increased (new size 180x360x147). Also, the time from B is entirely within the range of A and also contains the first and last values of A. So, in essence I want to interpolate B from A based on changes in the 3rd dimension (time). Does anyone know of a quick and dirty way to do this?

Answer 1

[Chad Greene]

The dirty way which is not quick but everyone seems to do it, is to perform 1D interpolation for each lat/lon grid box. That would look like this:

% Given this data: 
t = 1:135; 
A = rand(180,360,135); 

% You want Ai at times ti:
ti = linspace(1,135,147); 

% Preallocate Ai to speed up the loop: 
Ai = NaN(180,360,147); 

% Perform 1D interpolation at each lat/lon grid box:
for lati = 1:180
    for loni = 1:360; 
        Ai(lati,loni,:) = interp1(t,squeeze(A(lati,loni,:)),ti,'linear'); 
    end
end

However, nested loops are dreadful. The loops above call interp1 more than sixty four thousand times. That makes computation slow, and debugging nested loops tends to be a headache. A better way to do this is 3D interpolation. Grid your source data, grid the data points to which you want to interpolate, then call interpn just once:

[LON,LAT,T] = meshgrid(-179.5:179.5,-89.5:89.5,t); 
[LONI,LATI,TI] = meshgrid(-179.5:179.5,-89.5:89.5,ti); 

Ai2 = interpn(LAT,LON,T,A,LATI,LONI,TI,'linear');

[Image Analyst]

64 thousand iterations is not that many. I can do tens of millions of for loop iterations in less than a second. If it's slow, it would be the memory access rather than the looping per se. So, to avoid that, you could use permute() to turn the array "on it's side" and then call imresize() just once for every column (which is now a slice in the Z direction). imresize() requires the Image Processing Toolbox.

[Kurt]

Thanks for your responses. It looks like this might work too, which also uses permute.

X = permute(A,[3 1 2]);

B = interp1(time,X,time_new);

Where time is identical to the 3rd dimension found in A and time_new is the new time dimension I'm interpolating to.