将坐标元组信息转换为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是特定网格点上的感兴趣的度量。

我希望将上述数据格式转换并获取f、g和h的（NX1 x NX2 x NX3）numpy数组。

其中一些结果集相当大（80 x 120 x 100）。

有人知道如何以高效的方式进行此转换吗？

- cytochrome

---

你能否提供一些Python代码的小样本数据来说明你的问题？我不确定我理解你的输出格式以及当你说想要"f"、"g"和"h"数组时的意思。 - YXD

Answer 1

假设您将整个数组作为形状为 (Nx1 * Nx2 * Nx3, 6) 的数组 data 读入内存中。

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 的视图，因此会使用更多内存。

当然，你应该使用循环或列表推导式，而不是像上面那样重复三次丑陋的代码！

- Jaime

---

所有对 data[n] 的引用都不是指第n列，而是第n行吧？我认为大部分应该替换为 data[:,n]。 - askewchan

@askewchan 我已经编辑过了，但是为了简单起见，在调用np.loadtxt时设置unpack=True可能更好，这相当于保持代码不变并执行data = data.T。 - Jaime

谢谢您提供的好主意！现在我已经有了一个可行且灵活的解决方案来解决我的问题。 - cytochrome

Answer 2

无论如何，您都必须循环遍历整个文件，那么为什么不初始化数组并输入值呢？

棘手的部分是，如果您想要一个形状为(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']，这将返回一个带有上述dtype的(f,g,h)元组。

如果没有索引，它将返回一个(NX1,NX2,NX3)数组，其中包含所有h值。

---

现在我看了一下，如果n1，n2，n3总是相同的，您可能希望在xyz_index函数内定义它们，这样您就不必每次都传递ns：

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

- askewchan