将坐标元组信息转换为numpy数组

Question

我有一个有限元程序的结果，这些结果在三维空间中的网格点上给出了各种感兴趣的测量值（比如温度、密度、压力）。

这些值在每个坐标轴上是均匀分布的，但不同坐标轴之间的间距可能会有所不同。例如，

x1 = [0, 0.1, 0.2, ..., 1.0]      (a total of NX1 pts) 
x2 = [0, 0.5, 1.0, ..., 20]       (a total of NX2 pts) 
x3 = [0, 0.2, 0.4, ..., 15]       (a total of NX3 pts)

软件输出的结果格式如下：

x1_1, x2_1, x3_1, f_x, g_x, h_x
x1_1, x2_1, x3_2, f_x, g_x, h_x
x1_1, x2_1, x3_3, f_x, g_x, h_x
...
x1_1, x2_2, x3_1, f_x, g_x, h_x
x1_1, x2_2, x3_2, f_x, g_x, h_x
x1_1, x2_2, x3_3, f_x, g_x, h_x
...
x1_2, x2_1, x3_1, f_x, g_x, h_x
x1_2, x2_1, x3_2, f_x, g_x, h_x
x1_2, x2_1, x3_3, f_x, g_x, h_x
...

其中 f_x、g_x、h_x 是特定网格点上的感兴趣的测量值。

我想把上面的数据格式转换成 (NX1 x NX2 x NX3) 的 numpy 数组，分别用于 f、g 和 h。

有些结果集的规模相当大（80 x 120 x 100）。

有没有人能给我一些高效转换的建议？

Answer 1

无论如何，你都得遍历整个文件，那为什么不直接初始化一个数组，把值放进去呢？

这里有个难点，如果你想要一个形状为 (NX1,NX2,NX3) 的数组，但你的 x1,x2,x3 值是 float 类型的，那你就得想办法给数组索引。也许有现成的数据结构可以用，但你可以试试这样的方式：

def xyz_index((x,y,z),(n1,n2,n3)):
    """ return integer indices for x,y,z position
        given a constant step """
    return tuple(map(int,[x/n1,y/n2,z/n3]))

import numpy as np
NX1,NX2,NX3 =  (80, 120, 100)
ns = n1, n2, n3 =   (.1,.5,.2)
x1, x2, x3 = np.arange(0,1+n1,n1), np.arange(0,20+n2,n2), np.arange(0,15+n3,n3),

data = np.empty((NX1,NX2,NX3),dtype=[('f',float),('g',float),('h',float)])
with open(filename,'r') as f:
    for line in f:
        x,y,z,f,g,h = map(float,line.split(', '))
        data[xyz_index((x,y,z),ns)] = (f,g,h)

这样你就可以像下面这样访问你的数据：

要获取在点 x1,x2,x3 = .2, .5, 0. 的 h 值，可以使用：

data[xyz_index((.2,.5,0),ns)]['h']

如果没有 ['h']，这会返回一个包含 (f,g,h) 的元组，里面的数据类型是上面提到的。

如果不加索引，它会返回一个包含所有 h 值的 (NX1,NX2,NX3) 数组。

---

现在我想想，如果 n1, n2, n3 总是相同的，你可能想在你的 xyz_index 函数里面定义它们，这样每次就不用传 ns 进来了：

data[xyz_index(.2,.5,0)]['h']

Answer 2

假设你把整个数组读入内存，形成一个叫做 data 的数组，形状是 (Nx1 * Nx2 * Nx3, 6)。

data = np.loadtxt('data.txt', dtype=float, delimiter=',')

如果像你的例子所说，这些点是按字典顺序生成的，那么你只需要提取 f、g 和 h 这几列，然后重新调整它们的形状：

f = data[:, 3].reshape(Nx1, Nx2, Nx3)
g = data[:, 4].reshape(Nx1, Nx2, Nx3)
h = data[:, 5].reshape(Nx1, Nx2, Nx3)

如果你需要弄清楚 Nx1、Nx2 和 Nx3 的具体值，可以使用 np.unique 来帮助你：

Nx1 = np.unique(data[:, 0]).shape[0]
Nx2 = np.unique(data[:, 1]).shape[0]
Nx3 = np.unique(data[:, 2]).shape[0]

如果点的顺序不一定，想要更稳妥的解决方案，可以用 np.unique 来提取网格值的索引：

Nx1, idx1 = np.unique(data[:, 0], return_inverse=True)
Nx1 = Nx1.shape[0]
Nx2, idx2 = np.unique(data[:, 1], return_inverse=True)
Nx2 = Nx2.shape[1]
Nx3, idx3 = np.unique(data[:, 2], return_inverse=True)
Nx3 = Nx3.shape[0]

f = np.empty((Nx1, Nx2, Nx3))
f[idx1, idx2, idx3] = data[:, 3]
g = np.empty((Nx1, Nx2, Nx3))
g[idx1, idx2, idx3] = data[:, 4]
h = np.empty((Nx1, Nx2, Nx3))
h[idx1, idx2, idx3] = data[:, 5]

这样会为 f、g 和 h 创建新的数组，而不是直接在原来的 data 数组上操作，所以会占用更多的内存。

当然，别像我上面的代码那样丑陋地重复三次，应该用循环或者列表推导式来简化代码！