7

I have some data of a particle moving in a corridor with closed boundary conditions. Plotting the trajectory leads to a zig-zag trajectory.

enter image description here

I would like to know how to prevent plot() from connecting the points where the particle comes back to the start. Some thing like in the upper part of the pic, but without "."

The first idea I had was to find the index where the numpy array a[:-1]-a[1:] becomes positive and then plot from 0 to that index. But how would I get the index of the first occurrence of a positive element of a[:-1]-a[1:]? Maybe there are some other ideas.

tdy
  • 36,675
  • 19
  • 86
  • 83
Tengis
  • 2,721
  • 10
  • 36
  • 58

4 Answers4

16

I'd go a different approach. First, I'd determine the jump points not by looking at the sign of the derivative, as probably the movement might go up or down, or even have some periodicity in it. I'd look at those points with the biggest derivative.

Second, an elegant approach to have breaks in a plot line is to mask one value on each jump. Then matplotlib will make segments automatically. My code is:

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

And gives us as result: enter image description here

The disadvantage of course is that you lose one point at each break - but with the sampling rate you seem to have I guess you can trade this in for simpler code.

Michael J. Barber
  • 24,518
  • 9
  • 68
  • 88
Thorsten Kranz
  • 12,492
  • 2
  • 39
  • 56
  • 1
    +1) -- I learned a very elegant use case for masked array in matplotlib – tzelleke Jan 18 '13 at 13:54
  • I found a case where this solution fails. It is the case where the trajectory has no breaks in it. My fix was just to check the std if it is larger than some minimum (ad-hoc chosed). – Tengis Mar 01 '13 at 14:06
3

To find where the particle has crossed the upper boundary, you can do something like this:

>>> 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) calculates the discrete difference of a, while np.nonzero finds where the condition np.diff(a) < 0 is negative, i.e., the particle has moved downward.

Michael J. Barber
  • 24,518
  • 9
  • 68
  • 88
2

To avoid the connecting line you will have to plot by segments.

Here's a quick way to plot by segments when the derivative of a changes sign:

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

enter image description here

tiago
  • 22,602
  • 12
  • 72
  • 88
2

Based on Thorsten Kranz answer a version which adds points to the original data when the 'y' crosses the period. This is important if the density of data-points isn't very high, e.g. np.linspace(0., 100., 100) vs. the original np.linspace(0., 100., 1000). The x position of the curve transitions are linear interpolated. Wrapped up in a function its:

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

The code for comparison to the original answer of Thorsten Kranz with lower data-points density.

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
  • 817
  • 4
  • 9