python - Matplotlibplot_surface 伪面

Question

我在使用 matplotlib.plot_surface 时遇到了问题。当我重现this example时，我得到了我应该得到的，一切都很好:

但是，当我自己做某事时(绘制地球的等势 EGM96 geoid by NASA )，我会在图形上出现奇怪的面孔(蓝色区域):

生成我的图形的代码如下。我尝试更改 plot_surface 的 antialiased 和 shadow 参数，但无济于事。我已经不知道如何解决这个问题了，所以如果有人知道，甚至怀疑某些事情，我很乐意听到。

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot, matplotlib.cm, scipy.special, numpy, math

" Problem setup. "
GM = 3986004.415E8 # m**3/s**2, from EGM96.
a = 6378136.3 # m, from EGM96.
N_POINTS = 50 # Number of lattitudes and longitudes used to plot the geoid.
latitudes = numpy.linspace(0, 2*numpy.pi, N_POINTS) # Geocentric     latitudes and longitudes where the geoid will be visualised.
longitudes = numpy.linspace(0, 2*numpy.pi, N_POINTS)
radius = 6378136.3 # Radius at which the equipotential will be computed, m.
MAX_DEGREE = 2 # Maximum degree of the geopotential to visualise.

" EGM96 coefficients - minimal working example. "
Ccoeffs={2:[-0.000484165371736, -1.86987635955e-10, 2.43914352398e-06]}
Scoeffs={2:[0.0, 1.19528012031e-09, -1.40016683654e-06]}

" Compute the gravitational potential at the desired locations. "
gravitationalPotentials = numpy.ones( latitudes.shape ) # Gravitational potentials computed with the given geoid. Start with ones and just add the geoid corrections.

for degree in range(2, MAX_DEGREE+1): # Go through all the desired orders and compute the geoid corrections to the sphere.
    temp = 0. # Contribution to the potential from the current degree and all corresponding orders.
    legendreCoeffs = scipy.special.legendre(degree) # Legendre polynomial coefficients corresponding to the current degree.
    for order in range(degree): # Go through all the orders corresponding to the currently evaluated degree.
        temp += legendreCoeffs[order] * numpy.sin(latitudes) * (Ccoeffs[degree][order]*numpy.cos( order*longitudes ) + Scoeffs[degree][order]*numpy.sin( order*longitudes ))

    gravitationalPotentials = math.pow(a/radius, degree) * temp # Add the contribution from the current degree.

gravitationalPotentials *= GM/radius # Final correction.

" FIGURE THAT SHOWS THE SPHERICAL HARMONICS. "
fig = matplotlib.pyplot.figure(figsize=(12,8))
ax = Axes3D(fig)
ax.set_aspect("equal")
ax.view_init(elev=45., azim=45.)
ax.set_xlim([-1.5*radius, 1.5*radius])
ax.set_ylim([-1.5*radius, 1.5*radius])
ax.set_zlim([-1.5*radius, 1.5*radius])

# Make sure the shape of the potentials is the same as the points used to plot the sphere.
gravitationalPotentialsPlot = numpy.meshgrid( gravitationalPotentials, gravitationalPotentials )[0] # Don't need the second copy of colours returned by numpy.meshgrid
gravitationalPotentialsPlot /= gravitationalPotentialsPlot.max() # Normalise to [0 1]

" Plot a sphere. "
Xs = radius * numpy.outer(numpy.cos(latitudes), numpy.sin(longitudes))
Ys = radius * numpy.outer(numpy.sin(latitudes), numpy.sin(longitudes))
Zs = radius * numpy.outer(numpy.ones(latitudes.size), numpy.cos(longitudes))
equipotential = ax.plot_surface(Xs, Ys, Zs, facecolors=matplotlib.cm.jet(gravitationalPotentialsPlot), rstride=1, cstride=1, linewidth=0, antialiased=False, shade=False)

fig.show()

Answer 1

这些方程

Xs = radius * np.outer(np.cos(latitudes), np.sin(longitudes))
Ys = radius * np.outer(np.sin(latitudes), np.sin(longitudes))
Zs = radius * np.outer(np.ones(latitudes.size), np.cos(longitudes))

正在计算给定半径、纬度和经度的球坐标的笛卡尔 X、Y、Z 坐标。但如果是这样，那么经度的范围应该是 0 到 pi，而不是 0 到 2pi。因此，改变

longitudes = np.linspace(0, 2*np.pi, N_POINTS)

到

longitudes = np.linspace(0, np.pi, N_POINTS)

<小时/>

import math
import numpy as np
import scipy.special as special
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

" Problem setup. "
# m**3/s**2, from EGM96.
GM = 3986004.415E8 
# m, from EGM96.
a = 6378136.3 
# Number of lattitudes and longitudes used to plot the geoid.
N_POINTS = 50 
# Geocentric     latitudes and longitudes where the geoid will be visualised.
latitudes = np.linspace(0, 2*np.pi, N_POINTS) 
longitudes = np.linspace(0, np.pi, N_POINTS)
# Radius at which the equipotential will be computed, m.
radius = 6378136.3 
# Maximum degree of the geopotential to visualise.
MAX_DEGREE = 2 

" EGM96 coefficients - minimal working example. "
Ccoeffs={2:[-0.000484165371736, -1.86987635955e-10, 2.43914352398e-06]}
Scoeffs={2:[0.0, 1.19528012031e-09, -1.40016683654e-06]}

" Compute the gravitational potential at the desired locations. "
# Gravitational potentials computed with the given geoid. Start with ones and
# just add the geoid corrections.
gravitationalPotentials = np.ones( latitudes.shape ) 

# Go through all the desired orders and compute the geoid corrections to the
# sphere.
for degree in range(2, MAX_DEGREE+1): 
    # Contribution to the potential from the current degree and all
    # corresponding orders.
    temp = 0. 
    # Legendre polynomial coefficients corresponding to the current degree.
    legendreCoeffs = special.legendre(degree) 
    # Go through all the orders corresponding to the currently evaluated degree.
    for order in range(degree): 
        temp += (legendreCoeffs[order] 
                 * np.sin(latitudes) 
                 * (Ccoeffs[degree][order]*np.cos( order*longitudes ) 
                    + Scoeffs[degree][order]*np.sin( order*longitudes )))

    # Add the contribution from the current degree.
    gravitationalPotentials = math.pow(a/radius, degree) * temp 

# Final correction.
gravitationalPotentials *= GM/radius 

" FIGURE THAT SHOWS THE SPHERICAL HARMONICS. "
fig = plt.figure(figsize=(12,8))
ax = Axes3D(fig)
ax.set_aspect("equal")
ax.view_init(elev=45., azim=45.)
ax.set_xlim([-1.5*radius, 1.5*radius])
ax.set_ylim([-1.5*radius, 1.5*radius])
ax.set_zlim([-1.5*radius, 1.5*radius])

# Make sure the shape of the potentials is the same as the points used to plot
# the sphere.

# Don't need the second copy of colours returned by np.meshgrid
gravitationalPotentialsPlot = np.meshgrid(
    gravitationalPotentials, gravitationalPotentials )[0] 
# Normalise to [0 1]
gravitationalPotentialsPlot /= gravitationalPotentialsPlot.max() 

" Plot a sphere. "
Xs = radius * np.outer(np.cos(latitudes), np.sin(longitudes))
Ys = radius * np.outer(np.sin(latitudes), np.sin(longitudes))
Zs = radius * np.outer(np.ones(latitudes.size), np.cos(longitudes))
equipotential = ax.plot_surface(
    Xs, Ys, Zs, facecolors=plt.get_cmap('jet')(gravitationalPotentialsPlot), 
    rstride=1, cstride=1, linewidth=0, antialiased=False, shade=False)

plt.show()

产量