Upgrading from 2.1.1 to 3.3.2 Causes plot3 lines to be behind plot_surface objects

Bug report

Bug summary

Running the code below in matplotlib version 2.1.1 and 3.3.2 produces different results. In 2.1.1 the surface is behind the orbit and in 3.3.2 the surface is in front of the orbit. The code below is supposed to show the orbit around a planet in which case part of the orbit should be behind and some should be in front.

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

##Kerbin Parameters
G = 6.6742*10**-11; #%%Gravitational constant
Mkerbin = 5.2915158*10**22 #
#MMun = 9.7599066*10**20 #
muKerbin = G*Mkerbin
#muMun = G*MMun
Rkerbin = 600000 #meters
#RMun = 200000 #meters

##True Anamoly
nu = np.linspace(0,2*np.pi,100)
##Semi Major Axis of an 80 km parking orbit
alt_AGL = 80000
rp = Rkerbin + alt_AGL
#ra = 12000000
ra = rp
a = (ra+rp)/2.
#alt_AGL = 6000
#a = RMun + alt_AGL
##Eccentricity
e = (ra - rp)/(ra+rp)
##inclination
i = 0*98.0*np.pi/180.0 ##Drew's random satellite he wants a just slightly over polar retrograde orbit
###Longitude of the Ascending Node
W = 0*45*np.pi/180.0
#Argument of the periaps
w = 0.

##plot in the orbital plane
###Phat and Qhat
p = a*(1-e**2)
r = p/(1+e*np.cos(nu))
xp = r*np.cos(nu)
yq = r*np.sin(nu)

plt.plot(xp,yq,'r-')
plt.plot(xp[0],yq[0],'r*',markersize=20)
theta = np.linspace(0,2*np.pi,100)
xkerbin = Rkerbin*np.cos(theta)
ykerbin = Rkerbin*np.sin(theta)
plt.plot(xkerbin,ykerbin,'b-')
#xMun = RMun*np.cos(theta)
#yMun = RMun*np.sin(theta)
#plt.plot(xMun,yMun,'b-')
#plt.axis('equal')
plt.title('Orbital Plane')

###Rotate to Kerbin Centered Inertial Frame (KCI)
zr = 0*xp

TPI = np.asarray([[np.cos(W)*np.cos(w)-np.sin(W)*np.sin(w)*np.cos(i),-np.cos(W)*np.sin(w)-np.sin(W)*np.cos(w)*np.cos(i),np.sin(W)*np.sin(i)],
                 [np.sin(W)*np.cos(w)+np.cos(W)*np.sin(w)*np.cos(i),-np.sin(W)*np.sin(w)+np.cos(W)*np.cos(w)*np.cos(i),-np.cos(W)*np.sin(i)],
                 [np.sin(w)*np.sin(i),np.cos(w)*np.sin(i),np.cos(i)]])

xi = 0*xp
yj = 0*yq
zk = 0*zr
for x in range(0,len(xp)):
  xyzO = np.asarray([xp[x],yq[x],zr[x]]) ##3x1 vector
  xyzi = np.matmul((TPI),xyzO)
  xi[x] = xyzi[0]
  yj[x] = xyzi[1]
  zk[x] = xyzi[2]
  
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
xsph = np.cos(u)*np.sin(v)
ysph = np.sin(u)*np.sin(v)
zsph = np.cos(v)
ax.plot_surface(Rkerbin*xsph,Rkerbin*ysph,Rkerbin*zsph,color='blue')
ax.plot(xi,yj,zk,'r-')
ax.scatter(xi[0],yj[0],zk[0],'r*',s=20)
#ax.plot_wireframe(RMun*xsph,RMun*ysph,RMun*zsph,color='grey')
#ax.axis('square')

###Now let's plot velocity
vx = np.sqrt(muKerbin/p)*(-np.sin(nu))
vy = np.sqrt(muKerbin/p)*(e+np.cos(nu))
#vx = np.sqrt(muMun/p)*(-np.sin(nu))
#vy = np.sqrt(muMun/p)*(e+np.cos(nu))
v = np.sqrt(vx**2 + vy**2)

plt.figure()
plt.plot(nu,vx,label='Vx')
plt.plot(nu,vy,label='Vy')
plt.plot(nu,v,label='V')
plt.grid()
plt.xlabel('True Anomaly (rad)')
plt.ylabel('Velocity (m/s)')
plt.legend()
    
plt.show()

Expected Outcome

See attached screenshot. The left is python 2.7.1 and the right is 3.6.1
Notice the graphs produce different outputs.

Matplotlib version

2.2.1 and 3.3.2

Operating system: Ubuntu 18.04
Matplotlib version: 2.2.1 and 3.3.2
Matplotlib backend (print(matplotlib.get_backend())): Qt5Agg
Python version: 2.7.1 for matplotlib 2.2.1 and 3.6.1 for matplotlib 3.3.2
Jupyter version (if applicable): N/A
Other libraries: numpy version 1.19.2

I used pip3 and pip to install

If you add zorder to your call to plot on the 3D axis it works fine--like this:

ax.plot(xi,yj,zk,'r-', zorder=10)

Given the complexities of ordering, drawing in 3D, etc., and the simple work-around, I'm not sure if this is worth figuring out what changed or not. (master shows the same draw order as 3.3.2 for this.) It's also interesting to me that #18216 doesn't change it, though I would have expected it to.

The code below is supposed to show the orbit around a planet in which case part of the orbit should be behind and some should be in front.

You really need a proper 3-D rendering library for this. I recently used mayavi for exactly this problem (well, not exactly, but pretty similar) and it was quite easy.

import numpy as np
from mayavi import mlab
import scipy.spatial.transform as strans

def make_scene(time=0):
    dphi, dtheta = np.pi/250.0, np.pi/250.0
    [phi,theta] = np.mgrid[0:np.pi+dphi*1.5:dphi,0:2*np.pi+dtheta*1.5:dtheta]

    phi0 = np.arange(0, np.pi+dphi*1.5, dphi) # azimuth
    theta0 = np.arange(0, 2*np.pi+dtheta*1.5, dtheta) #

    orbit_radius = 5
    decl = 20 * np.pi / 180

    periodmoon = 27.3*24
    periodearth = 24


    moonrot = strans.Rotation.from_euler('zyx',
        [np.pi/2 + time / periodmoon * np.pi * 2,
        -decl, 0], degrees=False)

    def apply_surf(x, y, z, rot):
        shape = np.shape(x)
        X = np.vstack([x.flat, y.flat, z.flat])
        Xnew = rot.apply(X.T).T

        x = Xnew[0, :].reshape(shape)
        y = Xnew[1, :].reshape(shape)
        z = Xnew[2, :].reshape(shape)
        return x, y, z


    # Moon orbit


    xorb = orbit_radius * np.cos(theta0)
    yorb = orbit_radius * np.sin(theta0)
    zorb = 0 * np.cos(theta0)

    X = rot.apply(np.vstack([xorb, yorb, zorb]).T).T

    s = mlab.plot3d(X[0, :], X[1, :], X[2, :])

    # Moon

    r = 0.4
    x = r*np.cos(phi)*np.sin(theta) + orbit_radius
    y = r*np.sin(phi)*np.sin(theta)
    z = r*np.cos(theta)

    x,y,z = apply_surf(x, y, z, moonrot)

    s = mlab.mesh(x, y, z, color=(0.7, 0.7, 0.7))

    # Earth
    r = 1
    x = r*np.cos(phi)*np.sin(theta)
    y = r*np.sin(phi)*np.sin(theta)
    z = r*np.cos(theta)

    s = mlab.mesh(x, y, z, opacity=1, color=(0., 0.9, 0.0))

    # Equator
    x = r*np.cos(theta0)
    y = r*np.sin(theta0)
    z = r*np.cos(theta0) * 0

    s = mlab.plot3d(x, y, z, )

    # Axis
    x = np.zeros(10)
    y = np.zeros(10)
    z = np.linspace(-r*1.3, r*1.3, 10)
    s = mlab.plot3d(x, y, z, color=(0, 0, 0))

    # Victoria
    ph = time * np.pi * 2 / periodearth
    x = r * np.cos(ph) * np.cos(49*np.pi/180)
    y = r * np.sin(ph) * np.cos(49*np.pi/180)
    z = r* np.sin(49*np.pi/180)

    s = mlab.points3d(x, y, z, mode='sphere', scale_factor=.4)

    # parallel at Victoria

    x = r * np.cos(theta0) * np.cos(49*np.pi/180)
    y = r * np.sin(theta0) * np.cos(49*np.pi/180)
    z = r* np.sin(49*np.pi/180) + 0 * theta0

    s = mlab.plot3d(x, y, z)

    # Tidal Bulge

    x = 1.7 * r * np.cos(phi)*np.sin(theta)
    y = 1.2 * r * np.sin(phi)*np.sin(theta)
    z = 1.2 * r * np.cos(theta)

    x, y, z = apply_surf(x, y, z, moonrot)

    s = mlab.mesh(x, y, z, opacity=0.5,
                  color=(0., 0.5, 1.0))


    X,Y = np.mgrid[-10:10:0.2,-10:10:0.2]

    # Bulge Equator

    x = 1.7 * r * np.cos(theta0)
    y = 1.2*r*np.sin(theta0)
    z = 0 * r*np.cos(theta0)

for time in np.arange(0, 49, 0.5) + 0*24:
    make_scene(time=time)
    mlab.view(azimuth=-0, elevation=0, distance=20)
    mlab.orientation_axes()
    mlab.savefig(f'moonearth/Above{int(time*10):04d}.png')
    mlab.clf()

mlab.show()

Result is the top frame (I then added the line plot in Matplotlib)

Thanks! Mayavi was definitely the way to go. I made a simple script here just to make sure I can do it my self. I'd say you can go ahead and close this bug report.

import numpy as np
from mayavi import mlab
import matplotlib.pyplot as plt

###Alright Let's plot an orbit around Kerbin using this new Mayavi Library
##Kerbin Parameters
G = 6.6742*10**-11; #%%Gravitational constant
Mkerbin = 5.2915158*10**22 #
#MMun = 9.7599066*10**20 #
muKerbin = G*Mkerbin
#muMun = G*MMun
Rkerbin = 600000 #meters
#RMun = 200000 #meters

##True Anamoly
nu = np.linspace(0,2*np.pi,100)
##Semi Major Axis of an 80 km parking orbit
alt_AGL = 80000
rp = Rkerbin + alt_AGL
#ra = 12000000
ra = rp
a = (ra+rp)/2.
#alt_AGL = 6000
#a = RMun + alt_AGL
##Eccentricity
e = (ra - rp)/(ra+rp)
##inclination
i = 45.0*np.pi/180.0 ##Drew's random satellite he wants a just slightly over polar retrograde orbit
###Longitude of the Ascending Node
W = 0*45*np.pi/180.0
#Argument of the periaps
w = 0.

##plot in the orbital plane
###Phat and Qhat
p = a*(1-e**2)
r = p/(1+e*np.cos(nu))
xp = r*np.cos(nu)
yq = r*np.sin(nu)

plt.plot(xp,yq,'r-')
plt.plot(xp[0],yq[0],'r*',markersize=20)
theta = np.linspace(0,2*np.pi,100)
xkerbin = Rkerbin*np.cos(theta)
ykerbin = Rkerbin*np.sin(theta)
plt.plot(xkerbin,ykerbin,'b-')
#xMun = RMun*np.cos(theta)
#yMun = RMun*np.sin(theta)
#plt.plot(xMun,yMun,'b-')
#plt.axis('equal')
plt.title('Orbital Plane')

###Rotate to Kerbin Centered Inertial Frame (KCI)
zr = 0*xp

TPI = np.asarray([[np.cos(W)*np.cos(w)-np.sin(W)*np.sin(w)*np.cos(i),-np.cos(W)*np.sin(w)-np.sin(W)*np.cos(w)*np.cos(i),np.sin(W)*np.sin(i)],
                 [np.sin(W)*np.cos(w)+np.cos(W)*np.sin(w)*np.cos(i),-np.sin(W)*np.sin(w)+np.cos(W)*np.cos(w)*np.cos(i),-np.cos(W)*np.sin(i)],
                 [np.sin(w)*np.sin(i),np.cos(w)*np.sin(i),np.cos(i)]])

xi = 0*xp
yj = 0*yq
zk = 0*zr
for x in range(0,len(xp)):
  xyzO = np.asarray([xp[x],yq[x],zr[x]]) ##3x1 vector
  xyzi = np.matmul((TPI),xyzO)
  xi[x] = xyzi[0]
  yj[x] = xyzi[1]
  zk[x] = xyzi[2]


###Let's plot velocity
vx = np.sqrt(muKerbin/p)*(-np.sin(nu))
vy = np.sqrt(muKerbin/p)*(e+np.cos(nu))
#vx = np.sqrt(muMun/p)*(-np.sin(nu))
#vy = np.sqrt(muMun/p)*(e+np.cos(nu))
v = np.sqrt(vx**2 + vy**2)

plt.figure()
plt.plot(nu,vx,label='Vx')
plt.plot(nu,vy,label='Vy')
plt.plot(nu,v,label='V')
plt.grid()
plt.xlabel('True Anomaly (rad)')
plt.ylabel('Velocity (m/s)')
plt.legend()

###From here we can plot the orbit in 3D using matplotlib.
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
xsph = np.cos(u)*np.sin(v)
ysph = np.sin(u)*np.sin(v)
zsph = np.cos(v)
ax.plot_surface(Rkerbin*xsph,Rkerbin*ysph,Rkerbin*zsph,color='blue',zorder=1)
ax.plot(xi,yj,zk,'r-',zorder=0)
ax.scatter(xi[0],yj[0],zk[0],'r*',s=20)
#ax.plot_wireframe(RMun*xsph,RMun*ysph,RMun*zsph,color='grey')
#ax.axis('square')

#plt.show()

##Unfortunately When you use matplotlib the orbit does not plot properly.
##So instead we will use this mayavi module
##First let's plot the orbit of the satellite
#dphi, dtheta = np.pi/250.0, np.pi/250.0
#theta0 = np.arange(0, 2*np.pi+dtheta*1.5, dtheta) #
#orbit_radius = 5
#xorb = orbit_radius * np.cos(theta0)
#yorb = orbit_radius * np.sin(theta0)
#zorb = 0 * np.cos(theta0)
#s = mlab.plot3d(xorb,yorb,zorb)
f = 100000
s=mlab.plot3d(xi/f,yj/f,zk/f)

##Then let's plot a sphere the size of the satellite
rsat = 0.1*Rkerbin # km?
s=mlab.mesh((rsat*xsph+xi[0])/f,(rsat*ysph+yj[0])/f,(rsat*zsph+zk[0])/f,color=(0.7,0.7,0.7))

##Finally let's plot kerbin in blue at the center of the orbit
s=mlab.mesh(Rkerbin*xsph/f,Rkerbin*ysph/f,Rkerbin*zsph/f,color=(0.0,0.9,0.0))

##Then we set the proper view port
#mlab.view(azimuth=-0, elevation=45, distance=Rkerbin*10)
mlab.view(azimuth=-0, elevation=45, distance=5*ra/f)
mlab.orientation_axes()
#mlab.savefig(f'moonearth/Above{int(time*10):04d}.png')
mlab.show()'

Navigation Menu

Repository navigation

Upgrading from 2.1.1 to 3.3.2 Causes plot3 lines to be behind plot_surface objects #18612

Bug report

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Upgrading from 2.1.1 to 3.3.2 Causes plot3 lines to be behind plot_surface objects #18612

Description

Bug report

Activity

dopplershift commented on Sep 30, 2020

jklymak commented on Sep 30, 2020

cmontalvo251 commented on Oct 2, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions