Search code examples
pythonnumpyscipyinterpolation

NumPy Interpolate Between Two 2-D Arrays At Various Timesteps

I have a pair of two-dimensional arrays from a gridded dataset (in GeoTIFF format), both with the exact same resolution and number of rows/columns.

Suppose Array #1 is at timestep +0 hours and Array #2 is a timestep +3 hours. I am looking to interpolate and create additional arrays at intervals [1,2] in a linear interpolation.

I have perused this link about leveraging scipy to achieve this with stacked 1-D arrays, but would prefer to execute this inside NumPy if possible.

What is the best method or starting point to generate these additional interpolated 2-D arrays that are at the pre-determined intervals?

Solution

  • I believe you are just asking how to do this:

    import numpy as np
t0 = 0
a0 = np.array([[0,0,0],[1,1,1],[2,2,2]], dtype=np.float64)
t3 = 3
a3 = np.array([[3,3,3],[0,5,0],[-2,0,2]], dtype=np.float64)
aDiff = a3 - a0
t1 = 1
a1 = a0 + aDiff * ((t1 - t0) / (t3 - t0))
t2 = 2
a2 = a0 + aDiff * ((t2 - t0) / (t3 - t0))
print(a0)
print(a1)
print(a2)
print(a3)

    Output:

    [[0. 0. 0.]
 [1. 1. 1.]
 [2. 2. 2.]]
[[1.         1.         1.        ]
 [0.66666667 2.33333333 0.66666667]
 [0.66666667 1.33333333 2.        ]]
[[ 2.          2.          2.        ]
 [ 0.33333333  3.66666667  0.33333333]
 [-0.66666667  0.66666667  2.        ]]
[[ 3.  3.  3.]
 [ 0.  5.  0.]
 [-2.  0.  2.]]

    If you want a more generic solution, you can do this:

    • reshape the difference between the 2D input arrays to be 1D
    • repeat this row n times for the n required interpolations
    • multiply the i'th row by the respective interpolation factor (t[i] - tFirst) / (tNext - tFirst) to make each row contain the necessary interpolation steps for that interpolated time
    • add a repeated 1D row version of the first 2D array to these interpolation steps to get the results in row-by-row 1D shape
    • reshape into a 3D array containing 2D results at each interpolated time

    Sample code:

    import numpy as np
tFirst, aFirst = 0, np.array([[0,0,0],[1,1,1],[2,2,2]], dtype=np.float64)
tNext, aNext = 3, np.array([[3,3,3],[0,5,0],[-2,0,2]], dtype=np.float64)
tInterp = np.array([1, 2])

aDiff1D = np.reshape(aNext - aFirst, (1, np.size(aFirst)))
aDiffRepeated = np.repeat(aDiff1D, np.size(tInterp), axis=0)
aStep = aDiffRepeated * ((tInterp[:, None] - tFirst) / (tNext - tFirst))
aInterp = np.repeat(np.reshape(aFirst, (1, np.size(aFirst))), np.size(tInterp), axis = 0) + aStep
aInterp = np.reshape(aInterp, (np.size(tInterp), aFirst.shape[0], aFirst.shape[1]))
a1 = aInterp[0]
a2 = aInterp[1]