Python中将不规则间隔的数据重新采样到正则网格上

Question

我需要将2D数据重新采样到一个规则网格中。

这是我的代码：

import matplotlib.mlab as ml
import numpy as np

y = np.zeros((512,115))
x = np.zeros((512,115))

# Just random data for this test:
data = np.random.randn(512,115)

# filling the grid coordinates:    
for i in range(512):
    y[i,:]=np.arange(380,380+4*115,4)

for i in range(115):
    x[:,i] = np.linspace(-8,8,512)
    y[:,i] -=  np.linspace(-0.1,0.2,512)

# Defining the regular grid
y_i = np.arange(380,380+4*115,4)
x_i = np.linspace(-8,8,512)

resampled_data = ml.griddata(x,y,data,x_i,y_i)

(512,115) 是这个二维数据的形状，并且我已经安装了 mpl_toolkits.natgrid。

我的问题是，返回的是掩码数组，其中大多数条目都是 nan，而不是由大多数常规条目组成，仅在边界处为 nan 的数组。

有人能指出我做错了什么吗？

谢谢！

- Dzz

---

作为第一步，我只是尝试运行您的示例代码，但它无法正常工作。我得到了一个ValueError: x,y must be equal length 1-D arrays的错误。我正在从Enthought发行版中运行matplotlib v0.99.3。您能否修复示例以重现NaN，并在可能的情况下包括有关计算所需输入/输出的更多信息，也许可以用图表说明？ - dtlussier

我认为你必须安装mpl_toolkits中的natgrid才能使示例正常工作。 - Dzz

Answer 1

比较你的代码示例和问题标题，我觉得你有点困惑...

在你的示例代码中，你创建了一个规则网格的随机数据，然后将其重新采样到另一个规则网格上。你的示例中没有任何不规则数据...

(另外，这段代码不能直接运行，你应该研究一下meshgrid，而不是通过循环生成x和y网格。)

如果你想要对已经规则采样的网格进行重新采样，就像你在示例中所做的那样，有比griddata更高效的方法，或者说比我下面要描述的任何方法都更高效。（在这种情况下，scipy.ndimage.map_coordinates非常适合解决你的问题。）

根据你的问题，听起来你有一些不规则间距的数据，你想要插值到一个规则网格上。

在这种情况下，你可能有一些像这样的点：

import numpy as np
import matplotlib.pyplot as plt
#import matplotlib.mlab as mlab # 2023 use instead:
from scipy.interpolate import griddata

# Bounds and number of the randomly generated data points
ndata = 20
xmin, xmax = -8, 8
ymin, ymax = 380, 2428

# Generate random data
x = np.random.randint(xmin, xmax, ndata)
y = np.random.randint(ymin, ymax, ndata)
z = np.random.random(ndata)

# Plot the random data points
plt.scatter(x,y,c=z)
plt.axis([xmin, xmax, ymin, ymax])
plt.colorbar()
plt.show()

你可以像之前一样对数据进行插值处理...（继续上面的代码片段...）

# Size of regular grid
ny, nx = 512, 115

# Generate a regular grid to interpolate the data.
xi = np.linspace(xmin, xmax, nx)
yi = np.linspace(ymin, ymax, ny)
xi, yi = np.meshgrid(xi, yi)

# Interpolate using delaunay triangularization 
#zi = mlab.griddata(x,y,z,xi,yi) # 2023 use instead:
zi = griddata( (x,y), z, (xi,yi) )

# Plot the results
plt.figure()
plt.pcolormesh(xi,yi,zi)
plt.scatter(x,y,c=z)
plt.colorbar()
plt.axis([xmin, xmax, ymin, ymax])
plt.show()

然而，你会注意到在网格中出现了很多伪影。这是因为你的x坐标范围从-8到8，而y坐标范围从约300到约2500。插值算法试图使事物各向同性，而你可能希望进行高度各向异性的插值（这样在绘制网格时看起来各向同性）。

为了纠正这个问题，你需要创建一个新的坐标系统来进行插值。没有一种正确的方法来做这件事。下面我使用的方法可以工作，但是“最好”的方式取决于你的数据实际代表什么。

换句话说，根据你对数据所测量系统的了解来决定如何处理它。这在插值中始终是真实的！除非你知道结果应该是什么样子，并且对插值算法足够熟悉，能够利用先验信息获得优势，否则不应进行插值！默认情况下，griddata使用的Delaunay三角插值算法也有比它更灵活的插值算法，但对于简单的示例来说，它已经足够了...

无论如何，一种方法是重新调整x和y坐标的比例，使它们大致处于相同的数量级。在这种情况下，我们将把它们重新调整为从0到1的范围内...（请原谅这段代码看起来有点混乱...我只是想用这个作为一个例子...）

# (Continued from examples above...)
# Normalize coordinate system
def normalize_x(data):
    data = data.astype(np.float)
    return (data - xmin) / (xmax - xmin)

def normalize_y(data):
    data = data.astype(np.float)
    return (data - ymin) / (ymax - ymin)

x_new, xi_new = normalize_x(x), normalize_x(xi)
y_new, yi_new = normalize_y(y), normalize_y(yi)

# Interpolate using delaunay triangularization 
#zi = mlab.griddata(x_new, y_new, z, xi_new, yi_new) # 2023 use instead:
zi = griddata( (x_new, y_new), z, (xi_new, yi_new) )

# Plot the results
plt.figure()
plt.pcolormesh(xi,yi,zi)
plt.scatter(x,y,c=z)
plt.colorbar()
plt.axis([xmin, xmax, ymin, ymax])
plt.show()

希望这能有所帮助，无论如何... 对回答的长度感到抱歉！

- Joe Kington

---

我只想说谢谢你的回答。对于确实需要处理不规则网格数据的人来说，它非常有帮助。然而，我有一个问题。你能否指出任何参考资料，以便我更好地了解一些可用的插值方法，以便我可以选择最佳方法？ - Vorticity

真正出色的解释。远远超过了mlab.griddata文档，相比之下几乎是晦涩难懂的。 - hobs

这是一个很好的答案，你有没有想过用Java怎么实现？ - River

@JoeKington 谢谢你的回答，这真的帮了我很多（即使是在4年之后...）！ - Mathias711

“mlab.griddata” 使用“scipy.interpolate.griddata”作为其基础，而后者的文档更为详尽。 - Cuadue