Plot solid 3D surface on top of another one

[Michele_Castellana]

Hello,
I want to plot the red and the green sphere as two solid, three-dimensional objects in this minimal working example code:

import matplotlib.pyplot as plt
import numpy as np

r = 0.05
fig = plt.figure(figsize = (12,8), constrained_layout=True)
ax = plt.axes(projection='3d')

# draw sphere 1
u, v = np.mgrid[0:2*np.pi:100j, 0:np.pi:100j]
x =  r * np.cos(u)*np.sin(v)
y = r * np.sin(u)*np.sin(v)
z =  -r + r * np.cos(v)
ax.plot_surface(x, y, z, color="r", zorder = 1, alpha = 1)

# draw sphere 2
u, v = np.mgrid[0:2*np.pi:100j, 0:np.pi:100j]
x =  r * np.cos(u)*np.sin(v)
y =  r * np.sin(u)*np.sin(v)
z =  -2*r + r * np.cos(v)
ax.plot_surface(x, y, z, color="g", zorder = 1, alpha = 1)

plt.show()

In the plot, I would like to see the boundary where the two spheres touch, but this is not what I get with the minimal working example: as I turn the image by changing the viewpoint, I see that sometimes the red sphere is on top, sometimes the green sphere is on top, see the images attached. How may I achieve this ?

I read here that the zorder optional keyword argument is usually ignored by 3d axes objects, so I don’t know how to do it.

Thanks

[QuLogic]

Unfortunately, I think you are running in to the limits of our 3D renderer, which is not a true 3D renderer. It cannot calculate proper occlusion or culling, and objects are fully drawn in one pass.

The only way to get what you intend is to calculate the intersections yourself, or use something that does true 3D rendering like https://s3dlib.org/ (I have not tried the latter, but just going by its examples.)

[Dr.Z]

Applying s3dlib for your example, the code would be:

import matplotlib.pyplot as plt
import s3dlib.surface as s3d

rez, r = 5, 0.05                      # rez of 5 to illustrate a smooth boundary.
red_surface = s3d.SphericalSurface(rez,color='r').transform(scale=r, translate=[0,0,-r]   )
grn_surface = s3d.SphericalSurface(rez,color='g').transform(scale=r, translate=[0,0,-2*r] )
surface = red_surface + grn_surface   #  <--- create ONE surface object to render.

fig = plt.figure(constrained_layout=True)
fig.text(0.95,0.95,str(surface),ha='right')   # document object faces & vertices.
ax = plt.axes(projection='3d')
ax.view_init(12)                      # tilt the surface to illustrate the boundary.
s3d.auto_scale(ax,surface)
ax.add_collection3d(surface.shade().hilite()) # shade/hilite to emphasize the boundary.

plt.show()

changing the view_init to -12 yields: