绘制周期轨迹

Question

我有一些物体在闭合边界的走廊内移动的数据。 将轨迹绘制出来会得到一个锯齿形轨迹。

我想知道如何防止plot()连接粒子回到起点的点，就像图片上面部分那样，但没有"."

我最初的想法是找到numpy数组a[:-1]-a[1:]变为正的索引，然后从0到该索引进行绘图。 但是我怎样才能获得a[:-1]-a[1:]的第一次正值元素的索引呢？ 可能还有其他的想法。

- Tengis

---

仅仅不连接数据点就有很多值得探讨的地方。 - Michael J. Barber

Answer 1

我会选择不同的方法。首先，我不会通过查看导数的符号来确定跳跃点，因为运动可能向上或向下，甚至可能具有一些周期性。我会查看那些具有最大导数的点。

其次，一个优雅的方法在绘制线条时具有突破点是在每个跳跃点掩盖一个值。然后Matplotlib将自动创建线段。我的代码如下：

import pylab as plt
import numpy as np

xs = np.linspace(0., 100., 1000.)
data = (xs*0.03 + np.sin(xs) * 0.1) % 1

plt.subplot(2,1,1)
plt.plot(xs, data, "r-")

#Make a masked array with jump points masked
abs_d_data = np.abs(np.diff(data))
mask = np.hstack([ abs_d_data > abs_d_data.mean()+3*abs_d_data.std(), [False]])
masked_data = np.ma.MaskedArray(data, mask)
plt.subplot(2,1,2)
plt.plot(xs, masked_data, "b-")

plt.show()

结果如下图所示：

当然，缺点是每次断开连接都会损失一个点数-但是根据你所拥有的采样率，我想你可以用更简单的代码来交换这一点。

- Thorsten Kranz

---

+1) -- 我学到了一个非常优雅的用例，使用matplotlib中的掩码数组。 - tzelleke

我发现这种解决方案失败的情况。 这种情况是轨迹没有中断的情况。 我的修复方法只是检查标准差是否大于某个最小值（ad-hoc选择）。 - Tengis

Answer 2

要找出粒子穿过上边界的位置，你可以这样做：

>>> import numpy as np
>>> a = np.linspace(0, 10, 50) % 5
>>> a = np.linspace(0, 10, 50) % 5 # some sample data
>>> np.nonzero(np.diff(a) < 0)[0] + 1
array([25, 49])
>>> a[24:27]
array([ 4.89795918,  0.10204082,  0.30612245])
>>> a[48:]
array([ 4.79591837,  0.        ])
>>> 

np.diff(a) 函数计算数组 a 的离散差分，而 np.nonzero 函数则用于找出条件 np.diff(a) < 0 不成立的位置，即粒子向下移动的位置。

- Michael J. Barber

Answer 3

基于Thorsten Kranz的回答，这个版本在'y'穿过周期时向原始数据添加点。如果数据点的密度不是很高，例如np.linspace(0., 100., 100)与原始np.linspace(0., 100., 1000)相比，则这一点非常重要。曲线的x位置进行线性插值。封装成一个函数如下：

import numpy as np
def periodic2plot(x, y, period=np.pi*2.):
    indexes = np.argwhere(np.abs(np.diff(y))>.5*period).flatten()
    index_shift = 0
    for i in indexes:
        i += index_shift
        index_shift += 3  # in every loop it adds 3 elements

        if y[i] > .5*period:
            x_transit = np.interp(period, np.unwrap(y[i:i+2], period=period), x[i:i+2])
            add = np.ma.array([ period, 0., 0.], mask=[0,1,0])
        else:
            # interpolate needs sorted xp = np.unwrap(y[i:i+2], period=period)
            x_transit = np.interp(0, np.unwrap(y[i:i+2], period=period)[::-1], x[i:i+2][::-1])
            add = np.ma.array([ 0., 0., period], mask=[0,1,0])
        x_add = np.ma.array([x_transit]*3, mask=[0,1,0])

        x = np.ma.hstack((x[:i+1], x_add, x[i+1:]))
        y = np.ma.hstack((y[:i+1], add, y[i+1:]))
    return x, y

将Thorsten Kranz的原始答案与较低数据点密度的比较代码。

import matplotlib.pyplot as plt

x = np.linspace(0., 100., 100)
y = (x*0.03 + np.sin(x) * 0.1) % 1

#Thorsten Kranz: Make a masked array with jump points masked
abs_d_data = np.abs(np.diff(y))
mask = np.hstack([np.abs(np.diff(y))>.5, [False]])
masked_y = np.ma.MaskedArray(y, mask)

# Plot
plt.figure()
plt.plot(*periodic2plot(x, y, period=1), label='This answer')
plt.plot(x, masked_y, label='Thorsten Kranz')

plt.autoscale(enable=True, axis='both', tight=True)
plt.legend(loc=1)
plt.tight_layout()

- Andrew

Answer 4

为避免连接线，您需要按段绘制。以下是在导数a变号时按段绘制的快速方法：

import numpy as np
a = np.linspace(0, 20, 50) % 5  # similar to Micheal's sample data
x = np.arange(50)  # x scale

indices = np.where(np.diff(a) < 0)[0] + 1  # the same as Micheal's np.nonzero
for n, i in enumerate(indices):
    if n == 0:
        plot(x[:i], a[:i], 'b-')
    else:
        plot(x[indices[n - 1]:i], a[indices[n - 1]:i], 'b-')

- tiago