使用 Python 正确插入 4D 数据（在网格上）
Correctly Interpolate 4D Data (on Grid) using Python
目标
我在特定 3D 坐标 x y z 处给出了值 v。数据存储为 pandas 数据帧：
x y z v
0 -68.5 68.50 -10.00 0.297845
1 -68.5 -23.29 61.10 0.148683
2 -68.5 -23.29 63.47 0.142325
3 -68.5 -23.29 65.84 0.135908
4 -68.5 -23.29 68.21 0.129365
... ... ... ...
91804 68.5 23.29 151.16 0.118460
91805 68.5 23.29 153.53 0.119462
91806 68.5 23.29 155.90 0.120386
91807 68.5 23.29 139.31 0.112257
91808 68.5 -68.50 227.00 0.127948
我想在不属于数据框的新坐标处找到值，因此我正在研究如何有效地插入数据。
我做了什么：
由于坐标在网格上，我可以使用 interpn:
import numpy as np
from scipy.interpolate import interpn
# Extract the list of coordinates (I know that they are on a grid)
xs = np.array(df["x"].to_list())
ys = np.array(df["y"].to_list())
zs = np.array(df["z"].to_list())
# Extract the associated values
vs = np.array(df["v"].to_list())
重塑数据以适应 scipy 函数：
points = (np.unique(xs), np.unique(ys), np.unique(zs))
values= vs.reshape(len(np.unique(xs)), len(np.unique(ys)), len(np.unique(zs)))
为了测试插值，我想看看我是否得到相同的值，如果我输入与原始点相同的点：
request = (xs,ys,zs)
output = interpn(points, values, request)
...但是
我在想，我做错了什么？？
其他：
数据集
请在此处找到完整的数据集：https://filebin.net/u10lrw956enqhg5i
可视化
from mayavi import mlab
# Create figure
fig = mlab.figure(1, fgcolor=(0, 0, 0), bgcolor=(0, 0, 0))
mlab.points3d(xs,ys,zs,output)
mlab.view(azimuth=270, elevation=90, roll=180, figure=fig)
# View plot
mlab.show()
我强烈怀疑您的数据虽然在网格上，但未按顺序排列以允许对值进行简单的重塑。您有两个可用的解决方案，都涉及以不同方式重新排序数据。
解决方案 1
由于您已经在使用 np.unique 提取网格，因此您可以使用 return_inverse 参数获得 vs 的正确顺序：
px, ix = np.unique(xs, return_inverse=True)
py, iy = np.unique(ys, return_inverse=True)
pz, iz = np.unique(zs, return_inverse=True)
points = (px, py, pz)
values = np.empty_like(vs, shape=(px.size, py.size, pz.size))
values[ix, iy, iz] = vs
return_inverse 有点神奇，主要是因为它太违反直觉了。在这种情况下，对于值的每个元素，它会告诉您它对应于哪个唯一的、排序的总位置。
顺便说一句，如果您缺少网格元素，您可能需要将 np.empty_like(vs, shape=(px.size, py.size, pz.size)) 替换为 np.zeros_like(vs, shape=(px.size, py.size, pz.size)) 或 np.empty_like(vs, np.nan, shape=(px.size, py.size, pz.size))。在后一种情况下，您可以先在网格中插入 nan。
解决方案 2
更明显的解决方案是重新排列索引，以便您可以按照您尝试的方式重塑 vs。这只有在您确定没有丢失的网格元素时才有效。最简单的方法是对整个数据帧进行排序，因为 pandas 方法比 np.lexsort (IMO):
更不烦人
df.sort_values(['x', 'y', 'z'], inplace=True, ignore_index=True)
提取时，高效地提取：
xs, ys, zs, vs = df.to_numpy().T
由于所有内容都已排序，您不再需要 np.unique 来识别网格。唯一 x 个值的数量是：
nx = np.count_nonzero(np.diff(xs)) + 1
唯一值是：
bx = xs.size // nx
ux = xs[::bx]
y 值每 bx 个元素经历一个完整的循环，所以
ny = np.count_nonzero(np.diff(ys[:bx])) + 1
by = bx // ny
uy = ys[:bx:by]
并且 z (bz == 1):
nz = by
uz = zs[:nz]
现在您可以构造您的原始数组了：
points = (ux, uy, uz)
values = vs.reshape(nx, ny, nz)
目标
我在特定 3D 坐标 x y z 处给出了值 v。数据存储为 pandas 数据帧：
x y z v
0 -68.5 68.50 -10.00 0.297845
1 -68.5 -23.29 61.10 0.148683
2 -68.5 -23.29 63.47 0.142325
3 -68.5 -23.29 65.84 0.135908
4 -68.5 -23.29 68.21 0.129365
... ... ... ...
91804 68.5 23.29 151.16 0.118460
91805 68.5 23.29 153.53 0.119462
91806 68.5 23.29 155.90 0.120386
91807 68.5 23.29 139.31 0.112257
91808 68.5 -68.50 227.00 0.127948
我想在不属于数据框的新坐标处找到值，因此我正在研究如何有效地插入数据。
我做了什么：
由于坐标在网格上，我可以使用 interpn:
import numpy as np
from scipy.interpolate import interpn
# Extract the list of coordinates (I know that they are on a grid)
xs = np.array(df["x"].to_list())
ys = np.array(df["y"].to_list())
zs = np.array(df["z"].to_list())
# Extract the associated values
vs = np.array(df["v"].to_list())
重塑数据以适应 scipy 函数：
points = (np.unique(xs), np.unique(ys), np.unique(zs))
values= vs.reshape(len(np.unique(xs)), len(np.unique(ys)), len(np.unique(zs)))
为了测试插值，我想看看我是否得到相同的值，如果我输入与原始点相同的点：
request = (xs,ys,zs)
output = interpn(points, values, request)
...但是
我在想，我做错了什么？？
其他：
数据集
请在此处找到完整的数据集：https://filebin.net/u10lrw956enqhg5i
可视化
from mayavi import mlab
# Create figure
fig = mlab.figure(1, fgcolor=(0, 0, 0), bgcolor=(0, 0, 0))
mlab.points3d(xs,ys,zs,output)
mlab.view(azimuth=270, elevation=90, roll=180, figure=fig)
# View plot
mlab.show()
我强烈怀疑您的数据虽然在网格上，但未按顺序排列以允许对值进行简单的重塑。您有两个可用的解决方案，都涉及以不同方式重新排序数据。
解决方案 1
由于您已经在使用 np.unique 提取网格，因此您可以使用 return_inverse 参数获得 vs 的正确顺序：
px, ix = np.unique(xs, return_inverse=True)
py, iy = np.unique(ys, return_inverse=True)
pz, iz = np.unique(zs, return_inverse=True)
points = (px, py, pz)
values = np.empty_like(vs, shape=(px.size, py.size, pz.size))
values[ix, iy, iz] = vs
return_inverse 有点神奇，主要是因为它太违反直觉了。在这种情况下，对于值的每个元素，它会告诉您它对应于哪个唯一的、排序的总位置。
顺便说一句，如果您缺少网格元素，您可能需要将 np.empty_like(vs, shape=(px.size, py.size, pz.size)) 替换为 np.zeros_like(vs, shape=(px.size, py.size, pz.size)) 或 np.empty_like(vs, np.nan, shape=(px.size, py.size, pz.size))。在后一种情况下，您可以先在网格中插入 nan。
解决方案 2
更明显的解决方案是重新排列索引，以便您可以按照您尝试的方式重塑 vs。这只有在您确定没有丢失的网格元素时才有效。最简单的方法是对整个数据帧进行排序，因为 pandas 方法比 np.lexsort (IMO):
df.sort_values(['x', 'y', 'z'], inplace=True, ignore_index=True)
提取时，高效地提取：
xs, ys, zs, vs = df.to_numpy().T
由于所有内容都已排序，您不再需要 np.unique 来识别网格。唯一 x 个值的数量是：
nx = np.count_nonzero(np.diff(xs)) + 1
唯一值是：
bx = xs.size // nx
ux = xs[::bx]
y 值每 bx 个元素经历一个完整的循环，所以
ny = np.count_nonzero(np.diff(ys[:bx])) + 1
by = bx // ny
uy = ys[:bx:by]
并且 z (bz == 1):
nz = by
uz = zs[:nz]
现在您可以构造您的原始数组了：
points = (ux, uy, uz)
values = vs.reshape(nx, ny, nz)