Matplotlib 3D plot - parametric curve “wraparound” from certain perspectives

[dvida]

I have already posted a question on StackOverflow about this issue, and the people there pointed me that this is probably an issue inside matplotlib. Here is the original question: LINK
One comment suggested that this is probably because a wrong culling index while rendering the 3D scene.

EDIT: I am also reposting my Stack Overflow post here, as it was suggested by WeatherGod. Thanks for the suggestion!

I have been writing a Python script (LINK) for visualizing asteroid/comet/meteoroid orbits. The script also plots the position of planets and their orbits.

It works just right for orbits with small semi-major axis (i.e. "smaller" orbits). But when I have an orbit that goes way beyond Neptune (e.g. of a Halley-type comet), and from certain perspectives, there is a weird "wraparound" (for lack of a better word) effect.

Let me show you what I mean:

1. This image shows the plot from a perspective where it does not break.
2. When you rotate the same plot a bit to the right, it is as if the orbit folded in half and reversed its direction!
3. And if you look at the plot from a great distance, you can see that the elipse is plotted as it should be.

And here is a minimal version of the code with which the issue can be reproduced. The "wraparound" occurs only when the perspective of the camera is closely parallel with the large orbit.

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

def orbitalElements2Cartesian(a, e, I, peri, node, E):
    """ Convert orbital elements to Cartesian coordinates in the Solar System.

    Args: 
        a (float): semi-major axis (AU)
        e (float): eccentricity
        I (float): inclination (degrees)
        peri (float): longitude of perihelion (degrees)
        node (float): longitude of ascending node (degrees)
        E (float): eccentric anomaly (radians)

    """

    # The source of equations used:
    # http://farside.ph.utexas.edu/teaching/celestial/Celestialhtml/node34.html

    # Check if the orbit is parabolic or hyperbolic
    if e >=1:
        e = 0.99999999

    # Convert degrees to radians
    I, peri, node = map(np.radians, [I, peri, node])

    # True anomaly
    theta = 2*np.arctan(np.sqrt((1.0 + e)/(1.0 - e))*np.tan(E/2.0))

    # Distance from the Sun to the poin on orbit
    r = a*(1.0 - e*np.cos(E))

    # Cartesian coordinates
    x = r*(np.cos(node)*np.cos(peri + theta) - np.sin(node)*np.sin(peri + theta)*np.cos(I))
    y = r*(np.sin(node)*np.cos(peri + theta) + np.cos(node)*np.sin(peri + theta)*np.cos(I))
    z = r*np.sin(peri + theta)*np.sin(I)

    return x, y, z


if __name__ == '__main__':

    # Example orbital elements
    # a, e, incl, peri, node
    orb_elements = np.array([
        [2.363, 0.515, 4.0, 205.0, 346.1],
        [0.989, 0.089, 3.1, 55.6, 21.2],
        [0.898, 0.460, 1.3, 77.1, 331.2],
        [104.585332285, 0.994914, 89.3950, 130.8767, 282.4633]
        ])

    # Setup the plot
    fig = plt.figure()
    ax = fig.gca(projection='3d')


    # Eccentric anomaly (full range)
    E = np.linspace(-np.pi, np.pi, 100)

    # Plot the given orbits
    for i, orbit in enumerate(orb_elements):
        a, e, I, peri, node = orbit

        # Take extra steps in E if the orbit is very large
        if a > 50:
            E = np.linspace(-np.pi, np.pi, (a/20.0)*100)

        # Get the orbit in the cartesian space
        x, y, z = orbitalElements2Cartesian(a, e, I, peri, node, E)

        # Plot orbits
        ax.plot(x, y, z, c='#32CD32')

    # Add limits (in AU)
    ax.set_xlim3d(-5,5)
    ax.set_ylim3d(-5,5)
    ax.set_zlim3d(-5,5)

    plt.tight_layout()
    plt.show()

More info about my system and environment (according to the guidelines):

• Running on Ubuntu 15.04, Python 2.7.10, matplotlib 1.5.1.
• I have installed matplotlib using pip.
• A short code example which reproduces the issue is given above. The dependencies are only matplotlib, numpy and Python standard libraries.
• I have not checked if the issue is present on other Operating Systems or matplotlib versions.

[WeatherGod]

Interesting. I wonder if the 3d->2d calculation is getting messed up
somewhere. One possibility is that the orbit is getting "behind" the
viewer, and as such, some of the calculations are getting wrapped around to
the back of the scene.

While mplot3d isn't great for 3d scenes, it should be able to handle this
situation. Could you please bring the example code and images into this
issue as SO links tend to rot?

On Fri, Apr 15, 2016 at 3:56 AM, Denis Vida notifications@github.com
wrote:

I have already posted a question on StackOverflow about this issue, and
the people there pointed me that this is probably an issue inside
matplotlib. Here is the original question: LINK
http://stackoverflow.com/questions/36608197/matplotlib-3d-plot-parametric-curve-wraparound-from-certain-perspectives
One comment suggested that this is probably because a wrong culling index
while rendering the 3D scene.

Here is more info (according to the guidelines):

• Running on Ubuntu 15.04, Python 2.7.10, matplotlib 1.5.1.
• I have installed matplotlib using pip.
• A short code example which reproduces the issue is given in the
  StackOverflow post. The dependencies are only matplotlib and Python
  standard libraries.
• I have not checked if the issue is present on other Operating
  Systems or matplotlib versions.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
#6305

[github-actions]

This issue has been marked "inactive" because it has been 365 days since the last comment. If this issue is still present in recent Matplotlib releases, or the feature request is still wanted, please leave a comment and this label will be removed. If there are no updates in another 30 days, this issue will be automatically closed, but you are free to re-open or create a new issue if needed. We value issue reports, and this procedure is meant to help us resurface and prioritize issues that have not been addressed yet, not make them disappear. Thanks for your help!

[ksunden]

This does persist in current main.

Slightly more minimal example (took out unnecessary knobs from the ellipse calc, used fewer sig fig numbers, renamed variables, and adjusted the viewport to include the problem without interaction):

import numpy as np
import matplotlib.pyplot as plt

def ellipse(semiaxis, ecc, pts):
    """ Calculate points of an ellipse

    Args: 
        semiaxis (float): semi-major axis
        ecc (float): eccentricity
        pts (array, in radians): Angle from center of elipse to compute the points
    """
    # Only allow ellipses
    if ecc >=1:
        ecc = 0.99999999

    # Compute polar coords
    theta = 2*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(pts/2.0))
    r = semiaxis*(1.0 - e*np.cos(pts))

    # Cartesian coords
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    z = 0

    return x, y, z


if __name__ == '__main__':
    orb_elements = [(2., 0.5), (1., 0.1), (1., 0.5), (100.0, 0.999)]

    fig = plt.figure()
    ax = fig.subplots(subplot_kw={'projection':'3d'})

    pts = np.linspace(-np.pi, np.pi, 100)

    for i, (a, e) in enumerate(orb_elements):
        x, y, z = ellipse(a, e, pts)
        ax.plot(x, y, z)

    ax.set_xlim(-5,5)
    ax.set_ylim(-5,5)
    ax.set_zlim(-5,5)

    ax.view_init(5, 180, 0)

    plt.show()

Produces:

I will also note that it does not display this behavior for orthographic projections, only for perspective projection. As such I think the theory of it being because it is "behind" the "camera" is probably quite right.

[scottshambaugh]

I looked into this and the suggestions above were right - it's because we're not properly culling points behind the camera for 3D lines. Here's an even more minimal reproducible example:

import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
xs = np.linspace(0, -50, 100)
ys = np.ones(100)
zs = np.zeros(100)

ax.plot(xs, ys, zs, 'k')
ax.set(xlim=(-5, 5), ylim=(-5, 5), zlim=(-5, 5))
ax.view_init(5, 180, 0)

Which produces this:

But should produce this:

Fix incoming.

[scottshambaugh]

With #27385, the plot correctly outputs this: