CODEWITHSUNDEEP
pythonnumpymatplotlib
Tengis
Tengis

Reputation: 2809

Plot periodic trajectories

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.

Upvotes: 8

Views: 5242

Answers (4)

kho
kho

Reputation: 1291

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

Upvotes: 3

Thorsten Kranz
Thorsten Kranz

Reputation: 12765

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.

Upvotes: 17

tiago
tiago

Reputation: 23492

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

Upvotes: 2

Michael J. Barber
Michael J. Barber

Reputation: 25052

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.

Upvotes: 3

Related Questions