Расслоение контурного графика и Surface_plot в Matplotlib - PullRequest
0 голосов
/ 30 августа 2018

Я борюсь со слоями и zorder в питоне. Я делаю 3D-график, используя matplotlib с тремя соответствующими элементами: surface_plot планеты, surface_plot колец вокруг этой планеты и contourf изображение, которое показывает тень планеты, отбрасываемую на кольца.

Я хочу, чтобы графика точно показывала, как этот сценарий будет выглядеть в реальной жизни, когда кольца движутся вокруг планеты, а тень находится на кольцах в соответствующем месте. Если тень находится за планетой для данного POV, я хочу, чтобы тень была заблокирована планетой, и наоборот, если тень находится перед планетой для данного POV.

Чтобы было ясно, это ТОЛЬКО многоуровневая проблема. У меня есть планета, кольца и тень, которые все правильно строят. Тем не менее, тень никогда не будет отображаться перед планетой. Он действует так, как будто планета «блокирует» тень, хотя планета должна находиться под тенью с точки зрения слоистости.

Я перепробовал все, что могу придумать в терминах zorder и перестановки, в каком порядке различные элементы сюжета должны быть нарисованы. Кольца правильно отображаются перед планетой, но тень не будет.

Мой настоящий код очень длинный. Вот соответствующие части:

Настройка графика:

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)

Planet:

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.)

Shadow

def draw_shadowboundary(ax3, sundir):
    #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)

Объединить графику:

def plot_results(datafile, plot_names):
    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])

Код создает изображение, подобное этому:

enter image description here

Предполагается, что серая тень находится на кольцах перед планетой. Тем не менее, он не будет отображаться перед планетой, так что на самом деле появляется только маленький кусочек тени справа от планеты. Тень отображается правильно во всех сценариях, кроме случаев, когда она должна идти впереди планеты.

Есть какие-нибудь исправления для этого?

...