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

cmontalvo251 · 2020-09-30T16:40:16Z

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

The text was updated successfully, but these errors were encountered:

dopplershift · 2020-09-30T17:03:19Z

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.

jklymak · 2020-09-30T18:23:22Z

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)

cmontalvo251 · 2020-10-02T21:39:43Z

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

dstansby added the topic: mplot3d label Sep 30, 2020

timhoffm added the status: won't or can't fix label Oct 3, 2020

timhoffm closed this as completed Oct 3, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

cmontalvo251 commented Sep 30, 2020 •

edited

Loading

dopplershift commented Sep 30, 2020

jklymak commented Sep 30, 2020

cmontalvo251 commented Oct 2, 2020 •

edited

Loading

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

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

Comments

cmontalvo251 commented Sep 30, 2020 • edited Loading

Bug report

dopplershift commented Sep 30, 2020

jklymak commented Sep 30, 2020

cmontalvo251 commented Oct 2, 2020 • edited Loading

cmontalvo251 commented Sep 30, 2020 •

edited

Loading

cmontalvo251 commented Oct 2, 2020 •

edited

Loading