我在python中尝试分层和zorder
。我正在使用具有三个相关元素的matplotlib制作3D图:一个行星的surface_plot
,该行星周围的环的surface_plot
和一个contourf
图像,该图像显示了转换到环上的行星阴影。
我希望这些图形准确地显示这种情况在现实生活中的样子,环在地球上绕转,阴影在适当的位置穿过环。如果阴影位于给定POV的行星后面,则我希望阴影被行星遮挡,反之亦然,如果阴影位于给定POV的行星前面,则相反。
需要明确的是,这只是一个分层问题。我已经正确绘制了行星,环和阴影。但是,阴影永远不会显示在地球前面。它的作用就好像行星正在“遮挡”阴影,即使在分层方面应该将行星置于阴影之下。
我已经尝试过用zorder
想到的每件事,并重新排列各种绘图元素的绘制顺序。环确实可以正确显示在行星的前面,但是阴影不会。
我的实际代码很长。以下是相关部分:
绘图设置:
def orthographic_proj(zfront, zback):
a = (zfront+zback)/(zfront-zback)
b = -2*(zfront*zback)/(zfront-zback)
return np.array([[1,0,0,0],
[0,1,0,0],
[0,0,a,b],
[0,0,0,zback]])
def setup_saturn_plot(ax3, elev, azim, drawz, drawxy,view):
#ax3.set_aspect('equal','box')
ax3.view_init(elev=elev, azim=azim)
if(view=="top" or view == "Top" or view == "TOP"):
ax3.dist = 5.5
if(view=="star" or view == "Star" or view == "STAR"):
ax3.dist = 5.0 #4.5 is best value
proj3d.persp_transformation = orthographic_proj
# hide grid and background
ax3.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax3.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax3.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax3.grid(False)
# hide z axis in orthographic top view, xy axes in star view
if (drawz == False):
ax3.w_zaxis.line.set_lw(0.)
ax3.set_zticks([])
if (drawz == True):
ax3.set_zlabel('Z (1000 km)',fontsize=12)
if (drawxy == False):
ax3.w_xaxis.line.set_lw(0.)
ax3.set_xticks([])
ax3.w_yaxis.line.set_lw(0.)
ax3.set_yticks([])
if (drawxy == True):
ax3.set_xlabel('X (1000 km)',fontsize=12)
ax3.set_ylabel('Y (1000 km)',fontsize=12)
行星:
def draw_saturn(ax3, elev, azim):
# Saturn dimensions
radius = 60268. / 1000.
radius_pole = 54364. / 1000.
# draw Saturn
phi, theta = np.mgrid[0.0:np.pi:100j, 0.0:2.0*np.pi:100j]
x = radius*np.sin(phi)*np.cos(theta)
y = radius*np.sin(phi)*np.sin(theta)
z = radius_pole*np.cos(phi)
line3 = ax3.plot_surface(x, y, z, color="w", edgecolor='b', rstride = 8, cstride=5, shade=False, lw=0.25)
#line3 = ax3.plot_wireframe(x, y, z, color="w", edgecolor='b', rstride = 5, cstride=5, lw=0.25)
ax3.tick_params(labelsize=10)
戒指:
def draw_rings(ax3, elev, azim, draw_mode):
# Saturn dimensions
radius = 60268. / 1000.
# Saturn rings
dringmin = 1.110 * radius
dringmax = 1.236 * radius
cringmin = 1.239 * radius
titanringlet = 1.292 * radius
maxwellgap = 1.452 * radius
cringmax = 1.526 * radius
bringmin = 1.526 * radius
bringmax = 1.950 * radius
aringmin = 2.030 * radius
enckegap = 2.214 * radius
keelergap = 2.265 * radius
aringmax = 2.270 * radius
fringmin = 2.320 * radius
gringmin = 2.754 * radius
gringmax = 2.874 * radius
eringmin = 2.987 * radius
eringmax = 7.964 * radius
if (draw_mode == 'back'):
offset = -azim*np.pi/180. - 0.5*np.pi
if (draw_mode == 'front'):
offset = -azim*np.pi/180. + 0.5*np.pi
rad, theta = np.mgrid[dringmin:dringmax:4j, 0.0-offset:1.0*np.pi-offset:100j]
x = rad * np.cos(theta)
y = rad * np.sin(theta)
z = 0. * rad
line1 = ax3.plot_surface(x, y, z, color="w", edgecolor='b', rstride = 8, cstride=25, shade=False, lw=0.25,alpha=0.)
rad, theta = np.mgrid[cringmin:cringmax:4j, 0.0-offset:1.0*np.pi-offset:100j]
x = rad * np.cos(theta)
y = rad * np.sin(theta)
z = 0. * rad
line2 = ax3.plot_surface(x, y, z, color="w", edgecolor='b', rstride = 8, cstride=25, shade=False, lw=0.25,alpha=0.)
rad, theta = np.mgrid[bringmin:bringmax:4j, 0.0-offset:1.0*np.pi-offset:100j]
x = rad * np.cos(theta)
y = rad * np.sin(theta)
z = 0. * rad
line3 = ax3.plot_surface(x, y, z, color="w", edgecolor='b', rstride = 8, cstride=25, shade=False, lw=0.25,alpha=0.)
rad, theta = np.mgrid[aringmin:aringmax:4j, 0.0-offset:1.0*np.pi-offset:100j]
x = rad * np.cos(theta)
y = rad * np.sin(theta)
z = 0. * rad
line4 = ax3.plot_surface(x, y, z, color="w", edgecolor='b', rstride = 8, cstride=25, shade=False, lw=0.25,alpha=0.)
rad, theta = np.mgrid[fringmin:1.005*fringmin:2j, 0.0-offset:1.0*np.pi-offset:100j]
x = rad * np.cos(theta)
y = rad * np.sin(theta)
z = 0. * rad
line7 = ax3.plot_surface(x, y, z, color="w", edgecolor='b', rstride = 8, cstride=25, shade=False, lw=0.1,alpha=0.)
阴影:
def draw_shadowboundary(ax3, sundir):
sqrt = np.sqrt
#azimuthal angle between x direction and direction of sun
alpha = np.arctan2(sundir[1],sundir[0])
#adjustments to keep -pi/2 < alpha < pi/2
alphaadj = 0.*np.pi/180.
if (alpha<0.):
alpha += 2.*np.pi
if ((alpha >= np.pi/2.) & (alpha <= np.pi)):
alpha += np.pi
alphaadj = np.pi
if ((alpha > np.pi) & (alpha <= 3.*np.pi/2.)):
alpha -= np.pi
alphaadj = np.pi
if (alpha>3.*np.pi/2.):
alpha-=2*np.pi
#azimuthal angle between x direction and northern summer -- found using VIMS_2005_14_OMICET and VIMS_2017_053_ALPORI to define eq. of plane of Sun's annual path in chosen coordinate system: -0.193318*x + 0.1963755*y + 0.5471502*z = 0
beta = 44.5505*np.pi/180.
#Saturn's obliquity -- from NASA fact sheet
psi = 26.73*np.pi/180.
#Saturn's oblateness -- from NASA fact sheet
obl = 0.09796
#helpful definitions for optimization
cpsic = np.cos(psi*np.cos(alpha+beta))
spsic = np.sin(psi*np.cos(alpha+beta))
calpha = np.cos(alpha)
salpha = np.sin(alpha)
#Saturn's projected shorter planetary axis as seen by the sun & ring inner edge
req = 60268. / 1000.
b = req*sqrt((1.-obl)*(1.-obl)*cpsic*cpsic + spsic*spsic)
ringstart = 1.239 * req
ringend = 2.270 * req
#shadow boundary of Saturn's rings -- can approximate using a=inf and cancelling terms
a = 9.582*1.496*10.**5
shadowline = lambda x,y : (1/a)*sqrt((req*salpha*(-a+x*calpha*cpsic+y*salpha)*(y*calpha-x*cpsic*salpha)/sqrt((y*calpha-x*cpsic*salpha)**2 + (x*spsic)**2) + calpha*(a*cpsic*(x*calpha*cpsic+y*salpha) + b*x*(a-x*calpha*cpsic-y*salpha)*spsic*spsic/sqrt((y*calpha-x*cpsic*salpha)**2 + (x*spsic)**2)))**2 + (req*calpha*(a-x*calpha*cpsic-y*salpha)*(y*calpha-x*cpsic*salpha)/sqrt((y*calpha-x*cpsic*salpha)**2 + (x*spsic)**2) + salpha*(a*cpsic*(x*calpha*cpsic+y*salpha)+b*x*(a-x*calpha*cpsic-y*salpha)*spsic*spsic/sqrt((y*calpha-x*cpsic*salpha)**2 + (x*spsic)**2)))**2)
#azimuthal radius & antisolar angle for inequalities
radius = lambda x,y : np.sqrt(x**2+y**2)
anti = lambda x,y : abs(np.arctan2(y,x)-(alpha-alphaadj))
#properties of shadow
samples=1200
d = np.linspace(-3*req,3*req,samples)
x,y = np.meshgrid(d,d)
#z = ((radius(x,y)<=shadowline(x,y)) & (ringstart<=radius(x,y)) & (np.pi/2<=anti(x,y)) & (anti(x,y)<=3.*np.pi/2)).astype(int)
z = ((radius(x,y)<=shadowline(x,y)) & (ringstart<=radius(x,y)) & (radius(x,y)<=ringend) & (np.pi/2<=anti(x,y)) & (anti(x,y)<=3.*np.pi/2)).astype(int)
cmap = matplotlib.colors.ListedColormap(["k","k"])
#add shadow to plot
ax3.contourf(x,y,z, [0.5,1.50001], cmap=cmap,alpha=0.5)
合并图形:
import matplotlib
import numpy
from math import *
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D # <--- This is important for 3d plotting
from mpl_toolkits.mplot3d import proj3d
def plot_results(phi, theta, sundir=[0.5, 0.5]):
#plot_names.append("occultation_track_" + starname)
fig2 = plt.figure(figsize=(9,9))
ax3 = fig2.add_subplot(111, projection='3d')
setup_saturn_plot(ax3, phi, theta, False, False, "star")
draw_saturn(ax3, phi, theta)
draw_rings(ax3, phi, theta, 'back')
draw_rings(ax3, phi, theta, 'front')
draw_shadowboundary(ax3,sundir)
ax3.set_xlim([-200, 200])
ax3.set_ylim([-200, 200])
ax3.set_zlim([-200, 200])
plot_results(phi=40, theta=50, sundir = (30,60))
该代码生成如下图像:
灰色阴影应该位于行星前方的环上。但是,它不会显示在行星的前面,因此实际上只出现了行星右侧的一小部分阴影。阴影在所有情况下均能正确显示,除非需要在地球前方。
有什么解决办法吗?
最佳答案
目前,我正在努力弄清这段代码,但与此同时,至少到目前为止,这似乎是matplotlib3d的已知问题。
正如@TheImportanceOfBeingErnest很久以前指出的那样,此问题出现在mpl3d faq中
关于python - 在matplotlib中分层contourf图和surface_plot,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/52087571/