Interpolating 2 numpy arrays
<p>Is there any numpy or scipy or python function to interpolate between two 2D numpy array's? I have two 2D numpy arrays, and I want to apply changes to the first numpy array to make it similar to the second 2D array. The constraint is that I want the changes to be smooth. e.g., let the arrays be:</p>

<pre><code>A
[[1 1 1
  1 1 1
  1 1 1]]
</code></pre>

<p>and</p>

<pre><code>B
[[34  100 15
  62  17  87
  17  34  60]]
</code></pre>

<p>To make A similar to B, I could add 33 to the first grid cell of <code>A</code> and so on.. However, to make the changes smoother, I plan to compute a mean using a 2x2 window on array <code>B</code> and then apply the resulting changes to array <code>A</code>. Is there a built in numpy or scipy method to do this or follow this approach without using for loop.</p>

<p>You've just described a Kalman Filtering / data fusion problem. You have an initial state <strong>A</strong> that has some errors and you have some observations <strong>B</strong> that also have some noise. You want to improve your estimate of state <strong>A</strong> by injecting some information from <strong>B</strong>, all while accounting for spatially correlated errors in both datasets. We don't have any prior information about the errors in <strong>A</strong> and <strong>B</strong>, so we can just make it up. Here's an implementation:</p>

<pre><code>import numpy as np

# Make a matrix of the distances between points in an array
def dist(M):
    nx = M.shape[0]
    ny = M.shape[1]
    x = np.ravel(np.tile(np.arange(nx),(ny,1))).reshape((nx*ny,1))
    y = np.ravel(np.tile(np.arange(ny),(nx,1))).reshape((nx*ny,1))
    n,m = np.meshgrid(x,y)
    d = np.sqrt((n-n.T)**2+(m-m.T)**2)
    return d

# Turn a distance matrix into a covariance matrix. Here is a linear covariance matrix.
def covariance(d,scaling_factor):
    c = (-d/np.amax(d) + 1)*scaling_factor
    return c

A = np.array([[1,1,1],[1,1,1],[1,1,1]]) # background state
B = np.array([[34,100,15],[62,17,87],[17,34,60]]) # observations

x = np.ravel(A).reshape((9,1)) # vector representation
y = np.ravel(B).reshape((9,1)) # vector representation

P_a = np.eye(9)*50 # background error covariance matrix (set to diagonal here)
P_b = covariance(dist(B),2) # observation error covariance matrix (set to a function of distance here)

# Compute the Kalman gain matrix
K = P_a.dot(np.linalg.inv(P_a+P_b))

x_new = x + K.dot(y-x)
A_new = x_new.reshape(A.shape)

print(A)
print(B)
print(A_new)
</code></pre>

<p>Now, this method only works if your data are unbiased. So mean(A) must equal mean(B). But you'll still get okay results regardless. Also, you can play with the covariance matrices however you like. I'd recommend reading the Kalman filter wikipedia page for more details. </p>

<p>By the way, the example above yields:</p>

<pre><code>[[ 27.92920141  90.65490699   7.17920141]
 [ 55.92920141   7.65490699  79.17920141]
 [ 10.92920141  24.65490699  52.17920141]]
</code></pre>

python numpy