我使用 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 除外)。因此，对于每个时间步长和每个粒子，我们都会有大量类似这样的线条。
我想绘制粒子位置的时间演化图。我该怎么做？如果这太难了，我会对所有粒子位置的不同时间步长的静态图感到满意。
非常感谢。
Edit1:我试图做的是这个，但我之前没有展示它，因为它很不合适:
from matplotlib import numpy
import matplotlib.pyplot as plot
plot.plot(*np.loadtxt('data.dat',unpack=True), linewidth=2.0)
或者这个:
plot.plotfile('data.dat', delimiter=' ', cols=(0, 1), names=('col1', 'col2'), marker='o')
请您参考如下方法:
我会使用 numpy.loadtxt 来读取输入，但这只是因为后处理也需要 numpy。您可以将所有数据读入内存，然后找到分隔线，然后 reshape 其余数据以适合您的粒子数。以下假设没有一个粒子能准确地 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)
上面的代码产生了一个 Nt 元素数组 times 和一个 position 形状的数组 (Nt,Nparticles,2) 每个粒子在每个时间步长的 2d 位置。通过计算粒子的数量，我们可以让 numpy 确定第一维的大小(这就是 reshape() 中的 -1 索引的作用)。
对于绘图，您只需切入您的 positions 数组以提取您真正需要的内容。对于 2d x 数据和 2d y 数据，matplotlib.pyplot.plot() 将自动尝试绘制输入数组的列作为彼此的函数。以下是如何使用实际输入数据进行可视化的示例:
import matplotlib.pyplot as plt
t_indices = slice(None,None,500) # every 500th time step
particle_indices = slice(None) # every particle
#particle_indices = slice(None,5) # first 5 particles
#particle_indices = slice(-5,None) # last 5 particles
plt.figure()
_ = plt.plot(times[myslice],positions[myslice,particle_indices,0])
plt.xlabel('t')
plt.ylabel('x')
plt.figure()
_ = plt.plot(times[myslice],positions[myslice,particle_indices,1])
plt.xlabel('t')
plt.ylabel('z')
plt.figure()
_ = plt.plot(positions[myslice,particle_indices,0],positions[myslice,particle_indices,1])
plt.xlabel('x')
plt.ylabel('z')
plt.show()
每条线对应一个粒子。前两个图分别显示了 x 和 z 分量的时间演化，第三个图显示了 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]，其中两个切片都是根据...定义的, 切片。您应该查看文档，但简而言之，arr[slice(from,to,stride)] 等同于 arr[from:to:stride]，如果任何变量为 None，则对应的索引为空:arr[slice(-5,None)] 等同于 arr[-5 :]，即它将对数组的最后 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()
