Python 中 3D 数据数组中缺失值的插值

Question

我在以下问题上挣扎了相当长一段时间，但尚未找到有效的解决方案。

我想在我的数据数组中插入缺失值（nan）[在Python中]。该数组的维度为 lon、lat、time – 它是随时间变化的栅格数据集。

不幸的是，在某些时间步长中，所有值均缺失，并且无法使用 lon 和 lat (2D) 进行插值。这就是为什么我提出在时间轴上进行插值的想法。我希望一个时间步长的缺失值用完全相同的像素处的前后时间步长的值进行插值。

你有什么想法如何做到这一点吗？

我当前的尝试是：

” def arr_interp(数组): arrN=np.array(数组,copy=False) arrN[np.isnan(arrN)]=interpolate.interp2d(my_array.lat, my_array.lon, my_array.time, fill_value="nan")

arr_interp(my_array)”

Answer 1

问题在于 NaN 数据可能会形成块，您无法从邻居中进行插值。

解决方案是进行高斯赛德尔插值求解拉普拉斯方程（这会创建最小化函数曲率的数据）。

此代码查找 NaN 值并进行 3D 插值。我无权访问您的数据，因此它是通过合成数据完成的。

import numpy as np
import matplotlib.pyplot as plt


# create data
print("Creating data...")
size = 10  # 3D matrix of size: size³
# create x,y,z grid
x, y, z = np.meshgrid(np.arange(0, size), np.arange(
    0, size), np.arange(0, size))


def f(x, y, z):
    """function to create synthetic data"""
    return np.sin((x+y+z)/2)


data = np.zeros((size, size, size))
data[x, y, z] = f(x, y, z)

# create corrupted data
sizeCorruptedData = int(data.size*.2)  # 20% of data is corrupted
# create random x,y,z index for NaN values
xc, yc, zc = np.random.randint(0, size, (3, sizeCorruptedData))

corruptedData = data.copy()
corruptedData[xc, yc, zc] = np.nan

# Interpolate on NaN values
print("Interpolating...")

# get index of nan in corrupted data
nanIndex = np.isnan(corruptedData).nonzero()

interpolatedData = data.copy()
# make an initial guess for the interpolated data using the mean of the non NaN values
interpolatedData[nanIndex] = np.nanmean(corruptedData)


def sign(x):
    """returns the sign of the neighbor to be averaged for boundary elements"""
    if x == 0:
        return [1, 1]
    elif x == size-1:
        return [-1, -1]
    else:
        return [-1, 1]

#calculate kernels for the averages on boundaries/non boundary elements
for i in range(len(nanIndex)):
    nanIndex = *nanIndex, np.array([sign(x) for x in nanIndex[i]])

# gauss seidel iteration to interpolate Nan values with neighbors
# https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method
for _ in range(100):
    for x, y, z, dx, dy, dz in zip(*nanIndex):
        interpolatedData[x, y, z] = (
            (interpolatedData[x+dx[0], y, z] + interpolatedData[x+dx[1], y, z] +
             interpolatedData[x, y+dy[0], z] + interpolatedData[x, y+dy[1], z] +
             interpolatedData[x, y, z+dz[0]] + interpolatedData[x, y, z+dz[1]]) / 6)


# plot results
f, axarr = plt.subplots(2, 2)
axarr[0, 0].imshow(data[:, :, 1])
axarr[0, 0].title.set_text('Original Data')
axarr[0, 1].imshow(corruptedData[:, :, 1])
axarr[0, 1].title.set_text('Corrupted Data')
axarr[1, 0].imshow(interpolatedData[:, :, 1])
axarr[1, 0].title.set_text('Fixed Data')
axarr[1, 1].imshow(data[:, :, 1]-interpolatedData[:, :, 1])
axarr[1, 1].title.set_text('Error = Original-Fixed')
f.tight_layout()
plt.show()

Answer 2

我还没有足够的声誉来发表评论，但是添加到@Colim的答案中，用矢量化版本替换for循环将显着提高性能

for _ in tqdm(range(100)):
    x, y, z, dx, dy, dz = nanIndex
    interpolatedData[x, y, z] = (
            (interpolatedData[x+dx[:,0], y, z] + interpolatedData[x+dx[:,1], y, z] +
             interpolatedData[x, y+dy[:,0], z] + interpolatedData[x, y+dy[:,1], z] +
             interpolatedData[x, y, z+dz[:,0]] + interpolatedData[x, y, z+dz[:,1]]) / 6)