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

【问题标题】：Converting coordinate tuple information to numpy arrays将坐标元组信息转换为 numpy 数组
【发布时间】：2013-02-27 14:38:13
【问题描述】：

我有一个有限元程序的结果，该程序在三维空间中的规则间隔网格位置给出了各种感兴趣的测量值（例如温度、密度、压力）。

值沿每个坐标等间距分布，但不同坐标的间距可能不同。例如，

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)。

有没有人有任何提示可以有效地进行这种转换？

【问题讨论】：

有没有办法可以提供一些 Python 代码形式的小样本数据来说明你的问题？我不确定我是否理解你的输出格式以及你说你想要数组“for”f、g 和 h 时的意思。

标签： python numpy scipy

【解决方案1】：

假设您将整个数组作为数组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 数组的视图，因此会占用更多内存。

当然，您应该使用循环或列表推导，而不是我上面的丑陋代码重复三遍！

【讨论】：

不是所有对data[n] 的引用都给出nth 行，而不是列？我应该认为他们中的大多数应该替换为data[:,n]
@askewchan 我已经编辑过了，但是为了简单起见，在对np.loadtxt的调用中简单设置unpack=True可能会更好，相当于做data = data.T，然后离开代码原样。
感谢您的好主意！我现在有一个可行且灵活的解决方案来解决我的问题。

【解决方案2】：

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

棘手的部分是，如果你想要一个形状为(NX1,NX2,NX3) 的数组，但如果你的x1,x2,x3 值是floats，那么你必须以某种方式索引你的数组。也许为此存在数据结构，但您可以使用类似

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) 元组和上面的dtype。

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

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

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

【讨论】：