Trilinear Interpolation - Vectorising without Scipy

Question

I am looking to vectorize this piece of code but don't know where to begin. There has been another answer on this site answering a similar question to mine: 3D interpolation of NumPy arrays without SciPy , but I was confused as to how to implement it. I cannot use anything from SciPy, but I need my code to be fast as at the moment it is taking absolutely ages (as you would expect) for large arrays. I realise my code is slow because of a) using lists and b) using nested for loops.

I have a regular 3D grid of velocities (V) for latitude (lat)/longitude(lon)/time (t). I want to be able to get a single value back for a specific latitude/longitude/time.

My code is as follows:

def intrp_vel(lon_2_interp, lat_2_interp, t_2_interp, X, Y, T, V):
# lonlat should be a tuple of (lon, lat)
# t_2_interp should be a float of the current time step
# This function returns the surface velocity at any location in time.

Tlen = len(T)
Ylen = len(Y)-1
Xlen = len(X)-1

t_2_index = Tlen*(1-(T[-1]-t_2_interp)/(T[-1]-T[0]))
lat_2_index = Ylen*(1-(Y[-1]-lat_2_interp)/(Y[-1]-Y[0]))
lon_2_index = Xlen*(1-(X[-1]-lon_2_interp)/(X[-1]-X[0]))

time = np.linspace(0, Tlen, V.shape[0])
latitudes = np.linspace(0, Ylen, V.shape[1])
longitudes = np.linspace(0, Xlen, V.shape[2])

V1 = [] # To be brought down to 2D intp
V2 = [] # To be brought down to 1D intp
append1 = V1.append
for lats in latitudes:
    for lons in longitudes:
        append1(np.interp(t_2_index, time, V[:,lats,lons]))
V1 = np.reshape(V1, (len(Y),len(X)))
append2 = V2.append
for lons in longitudes:
    append2(np.interp(lat_2_index, latitudes, V1[:,lons]))

intrpvel = np.interp(lon_2_index, longitudes, V2)

return intrpvel

Any help would be much appreciated, or if you understand the code in the link above, and would be able to tell me how to use it in my own case.

score 0 · Answer 1 · answered Dec 20 '14 at 12:15

Here is some idea to remove the for loop call of np.interp().

Since t_2_index and time don't change in the loop, you can use np.interp() to calculate the linear mix parameter of V. Here is the code to confirm the idea:

y = np.sin(x)
x2 = np.linspace(x[0]-1, x[-1]+1, 100)
y2 = np.interp(x2, x, y)

calculate the mix parameters:

idx = np.arange(len(x))
idx_float = np.interp(x2, x, idx)
idx0 = idx_float.astype(int)
idx1 = np.clip(idx0 + 1, 0, len(x)-1)
a = idx_float - idx0
b = 1 - a

then calculate the interp result:

y3 = y[idx0] * b + y[idx1] * a
print np.allclose(y2, y3)

For your code, if you calculate the idx0, idx1, a, b from t_2_index and time, then you can calculate the interp results by:

V[idx0, :, :] * b + V[idx1, :, :] * a

Trilinear Interpolation - Vectorising without Scipy

1 Answers1