首先，我想说你想要的东西有点不明确。你想在一个底层表面上精确地绘制点，以使它们总是显示在图中，即。即 * 刚好在表面上方 *，而不明确地将它们向上移动。这已经是一个问题了，因为
1.浮点运算意味着点和曲面的精确坐标可以随着机器精度的变化而变化，因此试图依赖精确等式是行不通的。
1.即使这些数字精确到无限精确，曲面也是用近似平面的集合绘制的。这意味着你的 * 精确 * 数据点将 * 在函数为凸的近似曲面下 *。
然而，最大的问题是，matplotlib中的3d绘图在图形中绘制多个或复杂对象时是不可靠的。特别是，渲染器本质上是2d的，并且在试图找出对象的相对明显位置时经常遇到问题。要克服这个问题，你可以尝试hacking around the problem，或者切换到mayavi这样的3d渲染器。
不幸的是，zorder可选关键字参数通常被3d轴对象忽略。所以我在pyplot中唯一能想到的就是你几乎拥有的东西，注解了：使用ax.plot而不是ax.scatter。后者会生成第一张图中所示的图（无论视角如何，由于某种原因，每个散点都隐藏），前者会生成第二张图中所示的图（点可见）。通过从绘图样式中删除线条，我们几乎可以得到你想要的：

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

def f(x, y):                        
    return np.sin(2*x) * np.cos(2*y)

# data for the surface
x = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, x)
Z = f(X, Y)

# data for the scatter
xx = 4*np.random.rand(1000) - 2
yy = 4*np.random.rand(1000) - 2
zz = f(xx,yy)

fig = plt.figure(figsize=(8,6))
ax = plt.axes(projection='3d')
#ax.scatter(xx, yy, zz, c='r', marker='o')
ax.plot(xx, yy, zz, 'ro', alpha=0.5) # note the 'ro' (no '-') and the alpha

ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
                cmap='viridis', edgecolor='none')

但不完全是：很明显，在这种情况下，这些点总是可见的，即使它们应该隐藏在表面的一部分后面：

# note the change in points: generate only in the "back" quadrant
xx = 2*np.random.rand(1000) - 2
yy = 2*np.random.rand(1000)
zz = f(xx,yy)

fig = plt.figure(figsize=(8,6))
ax = plt.axes(projection='3d')
ax.plot(xx,yy,zz, 'ro', alpha=0.5)

ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
                cmap='viridis', edgecolor='none')

很容易看出，前面的凸起应该隐藏了背景中的一大块点，但这些点是可见的。这正是pyplot在复杂的3d可视化中遇到的问题。因此，我的底线是，你不能可靠地使用matplotlib做你想做的事情。无论如何，我不确定这样的情节有多容易理解。
为了更积极地结束，以下是如何使用mayavi（为此，您需要安装vtk，最好通过包管理器完成）：

import numpy as np
from mayavi import mlab
from matplotlib.cm import get_cmap # for viridis

def f(x, y):
    return np.sin(2*x) * np.cos(2*y)

# data for the surface
x = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, x)
Z = f(X, Y)

# data for the scatter
xx = 4*np.random.rand(1000) - 2
yy = 4*np.random.rand(1000) - 2
zz = f(xx,yy)

fig = mlab.figure(bgcolor=(1,1,1))
# note the transpose in surf due to different conventions compared to meshgrid
su = mlab.surf(X.T, Y.T, Z.T)
sc = mlab.points3d(xx, yy, zz, scale_factor=0.1, scale_mode='none',
                   opacity=1.0, resolution=20, color=(1,0,0))

# manually set viridis for the surface
cmap_name = 'viridis'
cdat = np.array(get_cmap(cmap_name,256).colors)
cdat = (cdat*255).astype(int)
su.module_manager.scalar_lut_manager.lut.table = cdat

mlab.show()

如您所见，结果是一个交互式3d图，其中表面上的数据点是适当的球体。人们可以尝试使用不透明度和球体比例设置，以获得令人满意的可视化效果。由于正确的3d渲染，您可以看到每个点的适当数量，而不管视角如何。

Andras Deak给出了一个非常全面的answer，讨论了matplotlibs 3D绘图功能在手头任务方面的问题/限制。在他回答的最后--以肯定的语气结束--他给出了一个使用替代库的解决方案。
我开始尝试在matplotlib中找到一个黑客/专业的解决方案。让我先谈谈原因。我想在2D表面上绘制轨迹，我被设置使用matplotlib。我使用它为我所有的二维绘图，并希望找到一个解决方案，为这个特定的三维绘图应用程序。matplotlibs 3D绘图的好处是它们是矢量化的，因为它们基本上只是通过将3D元素投影到相机平面上并覆盖它们（根据它们到相机的距离按顺序绘制它们）生成的2D绘图。可以单独控制图中每个元素的栅格化，而不会影响轴、标签等。使用光线跟踪的“真实的”3D绘图库通常不具有生成完全矢量化绘图的能力。我认为mayavi就是一个例子，我知道Mathematica在这方面也非常有限。
我的解决方案：我研究了plot_surface的代码，它最终基于Poly3DCollection，以了解matplotlib如何决定首先绘制表面上的哪些多边形/元素。Poly3DCollection的方法_do_3d_projection按照（原始3D对象的）到相机的距离对投影到2d相机平面上的多边形进行排序。首先绘制远离相机的元素，然后绘制靠近相机的元素。对于大多数情节来说，这足以创建一个合适的视角（但这种方法有局限性，见E。例如，mplot3d FAQ。然而，这个顺序是我解决方案的关键。给定一组点 pts 和一个表面 surf（必须用show或savefig绘制以设置其相机/投影变量）：
1.计算 surf 中的所有3D多边形到相机平面中的2D投影 segments_2d，包括它们基于到相机的距离的排序（存储在 segments_idxs 中）。
1.所有点都与三维曲面上的元素/多边形相关联。
1.计算3D点到相机平面中的2D投影。
1.为了确定一个点是否可见，我们检查它是否被一个多边形覆盖，该多边形在它所关联的多边形之后绘制（从步骤2。）。为此，我们使用matplotlib.path中的contains_points方法，另请参阅相关问题What's the fastest way of checking if a point is inside a polygon in python。
1.我包含了一个动态更新（改编自How to obscure a line behind a surface plot in matplotlib?）。警告：编码/绘制具有大量多边形的曲面可能会非常慢。
这里的必要的代码/最小的工作示例与表面给定的OP一个样本点位于单位圆。

import matplotlib.pyplot as plt
import numpy as np
import copy
import matplotlib.path as mpltPath
from mpl_toolkits.mplot3d import proj3d
from matplotlib import cm

def clip_on_surface(surf,pts):
    ## Get projection of 3d surface onto 2d camera plane
    ## [Code form [mpl_toolkits/mplot3d/art3d.py -  Poly3DCollection._do_3d_projection(self, renderer=None)] to ]
    txs, tys, tzs = proj3d._proj_transform_vec(surf._vec, surf.axes.M)
    xyzlist = [(txs[sl], tys[sl], tzs[sl]) for sl in surf._segslices]
    cface = surf._facecolor3d
    cedge = surf._edgecolor3d
    if len(cface) != len(xyzlist):
        cface = cface.repeat(len(xyzlist), axis=0)
    if len(cedge) != len(xyzlist):
        if len(cedge) == 0:
            cedge = cface
        else:
            cedge = cedge.repeat(len(xyzlist), axis=0)
                    
    if xyzlist:
        # sort by depth (furthest drawn first)
        z_segments_2d = sorted(
            ((surf._zsortfunc(zs), np.column_stack([xs, ys]), fc, ec, idx)
              for idx, ((xs, ys, zs), fc, ec)
              in enumerate(zip(xyzlist, cface, cedge))),
            key=lambda x: x[0], reverse=True)
        
    # z_segments_2d = sorted(z_segments_2d,key=lambda x:x[4])    
    segments_zorder, segments_2d, facecolors2d, edgecolors2d, segments_idxs =  zip(*z_segments_2d)
    segments_paths = [mpltPath.Path(x) for x in segments_2d]
    
    ## Get polygons in 3d space
    xs, ys, zs = surf._vec[0:3,:]
    xyzlist = [(xs[sl], ys[sl], zs[sl]) for sl in surf._segslices]
    segments_3d=[]
    segments_3d_centroid=[]
    for q in xyzlist:    
        vertices = np.transpose(np.array([q[0],q[1],q[2]]))
        segments_3d.append( vertices )
        segments_3d_centroid.append( sum(list(vertices))/len(list(vertices)) ) # centroid of polygon (mean of vertices)
    
    ## Process points
    pts_info = [[0,0,True] for x in range(len(pts))] 
        # 0: index of closest 3d polygon
        # 1: index of closest 3d polygon in segments_idxs: drawing order
        # 2: True if visible (not overlapped by polygons drawn after associated polygon), False else
    
    pts_visible = copy.copy(pts) # visible points (invisible set to np.nan)
    pts_invisible = copy.copy(pts) # invisible points (visible set to np.nan)
    
  
    # compute pts_info[:,0] and  pts_info[:,1] -- index of closest 3d polygon and its position in segments_idxs
    for i in range(len(pts)):
        # Associate by distance
        dist = np.inf
        x=[pts[i][0],pts[i][1],pts[i][2]]
        for j in range(len(segments_3d_centroid)):
            yc=segments_3d_centroid[j]
            dist_tmp = np.sqrt( (x[0]-yc[0])**2 + (x[1]-yc[1])**2 + (x[2]-yc[2])**2 )
            if dist_tmp<dist:
                dist=dist_tmp
                pts_info[i][0]=j        
        pts_info[i][1] = segments_idxs.index( pts_info[i][0] )
    
    # compute projection of 3d points into 2d camera plane
    pts_2d_x, pts_2d_y, pts_2d_z = proj3d._proj_transform_vec(np.transpose(np.array([[x[0],x[1],x[2],1.0] for x in pts])), surf.axes.M) 
  
    # decide visibility    
    for i in range(len(pts_info)):
        for j in range(pts_info[i][1]+1,len(segments_paths)):
            b=segments_paths[j].contains_points( [[pts_2d_x[i],pts_2d_y[i]]] )
            if b==True:
                pts_info[i][2]=False
                break
        if pts_info[i][2]:
            pts_invisible[i][0]=np.nan
            pts_invisible[i][1]=np.nan
            pts_invisible[i][2]=np.nan
        else:
            pts_visible[i][0]=np.nan
            pts_visible[i][1]=np.nan
            pts_visible[i][2]=np.nan
            
            
    return { 'pts_visible': pts_visible, 'pts_invisible':pts_invisible, 'pts_info':pts_info }

def f(x, y):                        
    return np.sin(2*x) * np.cos(2*y)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev=30., azim=55.)

# Generate surface plot (surf)
xs = np.linspace(-2, 2, 25)
ys = np.linspace(-2, 2, 25)
Xs, Ys = np.meshgrid(xs, ys)
zs = np.array(f(np.ravel(Xs), np.ravel(Ys)))
Zs = zs.reshape(Xs.shape)

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

surf = ax.plot_surface(Xs, Ys, Zs, rstride=1, cstride=1, 
                        cmap=cm.get_cmap('viridis'),linewidth=0.0,edgecolor='black',
                        antialiased=True,rasterized=False)

# Generate pts on surf
t = np.linspace(0, 1, 200)
xp = np.sin(t*2*np.pi)
yp = np.cos(t*2*np.pi)
zp = f(xp,yp)
pts=np.transpose(np.array([xp,yp,zp]))

def rotate(event):
    if event.inaxes == ax:
        surf_pts=clip_on_surface(surf,pts)
        ax.plot(surf_pts['pts_visible'][:,0],surf_pts['pts_visible'][:,1],surf_pts['pts_visible'][:,2],'.', zorder=10,c='red',markersize=2)
        ax.plot(surf_pts['pts_invisible'][:,0],surf_pts['pts_invisible'][:,1],surf_pts['pts_invisible'][:,2],'.', zorder=10,c='green',markersize=2)
        
        fig.canvas.draw_idle()
        
c1 = fig.canvas.mpl_connect('motion_notify_event', rotate)
plt.show()

代码仍然有点混乱，它仍然不能完美地工作，但这里有一些结果，表面上有25*25=625个四边形，单位圆上有200个点。x1c 0d1xx 1c 1d 1xx 1c 2d 1x红色的点是可见的点，绿色的点是不可见的点（这里绘制用于说明目的，但为了最初的问题/问题，人们会从图中省略它们）。一些点应该清晰可见，但被检测为不可见。我还不确定那里出了什么问题，但对我来说，这种有限的错过检测不是太大的问题，因为我最终想绘制许多（任意密集）点的线/轨迹。如果错过不集群，我可以生活与一些失踪的。
另一个固有问题/限制是当前方法不具有点是在表面上方还是下方的真实的概念，这意味着当从表面下方观看时，表面上方/上的点是可见的。下面是这种行为的示例：

这与Andras Deak已经提出的观点相联系，即目前的问题在没有额外限制的情况下有些定义不清，或者至少是模棱两可的。例如，可以要求将所有点放置在指向相机的表面上。在本方法中实现这一点是困难的。在几何方面，当前的实现方式将有限大小的球放置在无限小的多边形上，使得它们从两侧可见（这在某些用例中实际上可能是可行的/期望的）。
代码仍在进行中，如果我发现重大改进，我可能会更新这个答案。欢迎对一般方法和/或执行情况提出意见。我绝不是一个pythonMaven（我几乎只使用它来绘图和相关的非常轻的数据处理），所以它们在代码性能和范围方面可能有很大的改进空间。

从你的图中，看起来你愿意显示非线性优化器局部解的路径，所以我认为你应该考虑在等高线图上绘制一个线框：
...
ax.scatter（xx，yy，zz，c='r'，marker='o'）###这将只绘制点，而不是线
ax.plot_surface（X，Y，Z，rstride=1，cstride= 1，cmap='viridis'，edgecolor='none'）
ax.plot_wireframe（xx，yy，zz）###这将绘制线和可能的点。
...
were（xx，yy，zz）包含了到达局部最大值的路径，正如我所假设的，它是从非线性回归方法得到的。