我在使用matplotlib.plot_surface时遇到了问题。当我重现this example时，我得到了应有的一切：
但是，当我自己做某事（绘制地球的EGM96 geoid by NASA的等势线）时，图中出现了怪异的面孔（蓝色区域）：
生成我的图形的代码如下。我尝试更改antialiased的shadow和plot_surface参数，但无济于事。我没有解决该问题的方法，所以，如果有人知道，甚至怀疑某事，我很高兴听到。
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()
最佳答案
这些方程式
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()
产量
关于python - Matplotlib plot_surface虚假面，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/28787069/