使用python读取.dat文件

Question

我用MITgcm模拟了内波和水中移动的粒子。对于每个时间步，其输出如下所示：

   -9999 0.0000000000000000000  #Time step (0.00000000000000 seconds)
 1308.2021183321899       -14.999709364517091 # Particle 1 (X,Z)
 1308.2020142528656       -24.999521595698688 # Particle 2 (X,Z)
 1308.2018600072618       -34.999345597877536 # .
 1308.2016593336587       -44.999185870669805 # .
 1308.2014165588744       -54.999046508237896 # .
 1308.2011370083103       -64.998931076248894
 1308.2008269116873       -74.998842490305705
 1308.2004933548124       -84.998782925797485
 1308.2001441978532       -94.998753764086956
 1308.1997879652938       -104.99875557384759
 1308.1994336881464       -114.99878812280582
 1308.1990906721119       -124.99885041328211
 1308.1987681881285       -134.99894073461562
 1308.1984750963150       -144.99905672694641
 1308.1982194336249       -154.99919545294702
 1308.1980080134056       -164.99935347476733
 1308.1978461242272       -174.99952693694112
 1308.1977378137256       -184.99971163492469
 1308.2000000000000       -195.00000000000000
 5232.8000000000002       -15.000038916290352
 5232.8000000000002       -25.000064153684303
 5232.8000000000002       -35.000089286157163
 5232.8000000000002       -45.000114270293523
 5232.8000000000002       -55.000139061712051 # Particle 57

其中-9999#number是时间步长（以秒为单位），左栏是X位置，右栏是Z位置（以米为单位）；每行都是不同的粒子（除了-9999）。所以我们会有大量的线，每一个时间步，每一个粒子都有这样的东西。在

我想画出粒子位置的时间演化。我该怎么做？如果这太难了，我会很乐意用所有粒子位置的不同时间步的静态图。在

非常感谢你。在

编辑1：我想做的是这个，但我之前没有展示过，因为它太不恰当了：

^{pr2}$

或者这个：

 plot.plotfile('data.dat', delimiter=' ', cols=(0, 1), names=('col1', 'col2'), marker='o')

Answer 1

我将使用numpy.loadtxt读取输入，但这仅仅是因为后处理也需要numpy。你可以把所有的数据读到内存中，然后找到分隔线，然后改变其余数据的形状以适应你的粒子数。下面假设没有一个粒子能精确到达的x=-9999，这应该是一个合理的（虽然不是万无一失的）假设。在

import numpy as np
filename = 'input.dat'
indata = np.loadtxt(filename, usecols=(0,1)) # make sure the rest is ignored
tlines_bool = indata[:,0]==-9999
Nparticles = np.diff(np.where(tlines_bool)[0][:2])[0] - 1
# TODO: error handling: diff(np.where(tlines_bool)) should be constant
times = indata[tlines_bool,1]
positions = indata[np.logical_not(tlines_bool),:].reshape(-1,Nparticles,2)

上面的代码为每个粒子在每个时间步的2d位置生成一个Nt-元素数组times和形状为(Nt,Nparticles,2)的数组position。通过计算粒子数，我们可以让numpy确定第一维的大小（这就是reshape()中的-1索引的用途）。在

要进行绘图，您只需将positions数组切片，以提取您真正需要的内容。如果是2dx数据和2dy数据，matplotlib.pyplot.plot()将自动尝试将输入数组的列作为彼此的函数来绘制。在这里，您可以使用一个实际的输入来可视化：

^{pr2}$

每条线对应一个粒子。前两个图分别显示了x和{}分量的时间演化，第三个图显示了z(x)轨迹。请注意，数据中有许多粒子根本不移动：

>>> sum([~np.diff(positions[:,k,:],axis=0).any() for k in range(positions.shape[1])])
15

（这将逐个计算每个粒子的两个坐标的时间定向差，并计算两个维度上的每个差均为0的粒子数，即粒子不移动）。这解释了前两个图中所有的水平线；这些静止粒子在第三个图中根本没有出现（因为它们的轨迹是一个单点）。在

我特意引入了一个有点花哨的索引，这样可以更轻松地处理数据。如您所见，索引如下所示：times[myslice]，positions[myslice,particle_indices,0]，其中两个片段都是根据slice定义的。您应该看一下文档，但是短篇故事是arr[slice(from,to,stride)]等价于arr[from:to:stride]，如果其中任何一个变量是None，那么相应的索引是空的：arr[slice(-5,None)]相当于{}，即它将对数组的最后5个元素进行切片。在

因此，如果您使用较少的轨迹进行打印（因为57是很多），您可以考虑添加一个图例（只要matplotlib的默认颜色循环允许您区分粒子，这才有意义，否则您必须设置手动颜色或更改轴的默认颜色循环）。为此，您必须保留从plot返回的句柄：

particle_indices = slice(None,5)   # first 5 particles
plt.figure()
lines = plt.plot(positions[myslice,particle_indices,0],positions[myslice,particle_indices,1])
plt.xlabel('x')
plt.ylabel('z')
plt.legend(lines,['particle {}'.format(k) for k in range(len(t))])
plt.show()