具有NaN值或掩码的大数组的双变量结构插值

Question

我正在尝试使用Scipy的RectBivariateSpline类来处理规则网格的风应力数据。在某些网格点上，输入数据包含无效的条目，这些条目被设置为NaN值。最开始，我参考了Scott的问题，使用了二维插值的解决方案。用我的数据进行插值时，结果返回的数组全是NaN。我还尝试了另一种方法，假设我的数据是非结构化的，使用了SmoothBivariateSpline类。显然，我的数据点太多，无法使用非结构化插值，因为数据数组的形状是(719 x 2880)。

为了说明我的问题，我创建了以下脚本：

from __future__ import division
import numpy
import pylab

from scipy import interpolate


# The signal and lots of noise
M, N = 20, 30  # The shape of the data array
y, x = numpy.mgrid[0:M+1, 0:N+1]
signal = -10 * numpy.cos(x / 50 + y / 10) / (y + 1)
noise = numpy.random.normal(size=(M+1, N+1))
z = signal + noise


# Some holes in my dataset
z[1:2, 0:2] = numpy.nan
z[1:2, 9:11] = numpy.nan
z[0:1, :12] = numpy.nan
z[10:12, 17:19] = numpy.nan


# Interpolation!
Y, X = numpy.mgrid[0.125:M:0.5, 0.125:N:0.5]
sp = interpolate.RectBivariateSpline(y[:, 0], x[0, :], z)
Z = sp(Y[:, 0], X[0, :])

sel = ~numpy.isnan(z)
esp = interpolate.SmoothBivariateSpline(y[sel], x[sel], z[sel], 0*z[sel]+5)
eZ = esp(Y[:, 0], X[0, :])


# Comparing the results
pylab.close('all')
pylab.ion()

bbox = dict(edgecolor='w', facecolor='w', alpha=0.9)
crange = numpy.arange(-15., 16., 1.)

fig = pylab.figure()
ax = fig.add_subplot(1, 3, 1)
ax.contourf(x, y, z, crange)
ax.set_title('Original')
ax.text(0.05, 0.98, 'a)', ha='left', va='top', transform=ax.transAxes, 
    bbox=bbox)

bx = fig.add_subplot(1, 3, 2, sharex=ax, sharey=ax)
bx.contourf(X, Y, Z, crange)
bx.set_title('Spline')
bx.text(0.05, 0.98, 'b)', ha='left', va='top', transform=bx.transAxes, 
    bbox=bbox)

cx = fig.add_subplot(1, 3, 3, sharex=ax, sharey=ax)
cx.contourf(X, Y, eZ, crange)
cx.set_title('Expected')
cx.text(0.05, 0.98, 'c)', ha='left', va='top', transform=cx.transAxes, 
    bbox=bbox)

这段代码的结果如下图所示：

图中展示了一个构建的数据图（a），以及使用Scipy的RectBivariateSpline（b）和SmoothBivariateSpline（c）类的结果。第一次插值的结果是一个全是NaN的数组。理想情况下，我希望得到的结果类似于第二次插值的结果，虽然那种方法计算量更大。我并不一定需要在领域区域之外进行数据外推。

Answer 1

你可以使用 griddata 这个工具，不过它在处理边缘部分的时候表现得不是很好。你可以尝试一些方法来改善这个问题，比如根据你的数据进行反射处理……下面是一个例子：

from __future__ import division
import numpy
import pylab
from scipy import interpolate


# The signal and lots of noise
M, N = 20, 30  # The shape of the data array
y, x = numpy.mgrid[0:M+1, 0:N+1]
signal = -10 * numpy.cos(x / 50 + y / 10) / (y + 1)
noise = numpy.random.normal(size=(M+1, N+1))
z = signal + noise


# Some holes in my dataset
z[1:2, 0:2] = numpy.nan
z[1:2, 9:11] = numpy.nan
z[0:1, :12] = numpy.nan
z[10:12, 17:19] = numpy.nan

zi = numpy.vstack((z[::-1,:],z))
zi = numpy.hstack((zi[:,::-1], zi))
y, x = numpy.mgrid[0:2*(M+1), 0:2*(N+1)]
y *= 5 # anisotropic interpolation if needed.

zi = interpolate.griddata((y[~numpy.isnan(zi)], x[~numpy.isnan(zi)]),
                zi[~numpy.isnan(zi)], (y, x), method='cubic')
zi = zi[:(M+1),:(N+1)][::-1,::-1]

pylab.subplot(1,2,1)
pylab.imshow(z, origin='lower')
pylab.subplot(1,2,2)
pylab.imshow(zi, origin='lower')
pylab.show()

如果你遇到内存不足的问题，可以考虑把你的数据分开处理，像这样：

def large_griddata(data_x, vals, grid, method='nearest'):
    x, y = data_x
    X, Y = grid
    try:
        return interpolate.griddata((x,y),vals,(X,Y),method=method)
    except MemoryError:
        pass

    N = (np.min(X)+np.max(X))/2.
    M = (np.min(Y)+np.max(Y))/2.

    masks = [(x<N) & (y<M),
             (x<N) & (y>=M),
             (x>=N) & (y<M),
             (x>=N) & (y>=M)]

    grid_mask = [(X<N) & (Y<M),
             (X<N) & (Y>=M),
             (X>=N) & (Y<M),
             (X>=N) & (Y>=M)]

    Z = np.zeros_like(X)
    for i in range(4):
        Z[grid_mask[i]] = large_griddata((x[masks[i]], y[masks[i]]),
                    vals[masks[i]], (X[grid_mask[i]], Y[grid_mask[i]]), method=method)

    return Z