python - 在 matplotlib 中分层轮廓图和表面图

Question

我正在努力分层和zorder在 python 中。我正在使用带有三个相关元素的 matplotlib 制作 3D 绘图:A 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))

代码生成这样的图像:



灰色阴影应该位于行星前面的环上。但是，它不会显示在行星前面，因此实际上只出现了行星右侧的一小片阴影。阴影在所有情况下都能正确显示，除非它需要在行星前面。

对此有任何修复吗？

Answer 1

我目前正在努力解决这段代码，但与此同时，至少到目前为止，这似乎是 matplotlib3d 的一个已知问题。

正如@TheImportanceOfBeingErnest 很久以前指出的那样，这个问题出现在 mpl3d faq 中。

My 3D plot doesn’t look right at certain viewing angles

This is probably the most commonly reported issue with mplot3d. The problem is that – from some viewing angles – a 3D object would appear in front of another object, even though it is physically behind it. This can result in plots that do not look “physically correct.”

Unfortunately, while some work is being done to reduce the occurrence of this artifact, it is currently an intractable problem, and can not be fully solved until matplotlib supports 3D graphics rendering at its core.

The problem occurs due to the reduction of 3D data down to 2D + z-order scalar. A single value represents the 3rd dimension for all parts of 3D objects in a collection. Therefore, when the bounding boxes of two collections intersect, it becomes possible for this artifact to occur. Furthermore, the intersection of two 3D objects (such as polygons or patches) can not be rendered properly in matplotlib’s 2D rendering engine.

This problem will likely not be solved until OpenGL support is added to all of the backends (patches are greatly welcomed). Until then, if you need complex 3D scenes, we recommend using MayaVi.