Draw line over surface plot

Question

I want to be able to see lines (and points) that are ON THE 3D SURFACE on top of the surface (second image), not behind (first image). This is my 3D function:

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

X,Y,Z for the 3D-surface:

x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

I've generated a vector of x points (xx) and of y points (yy), where 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, c= 'r')

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

As you can see, the points are behind the plot, the figure covers the points. I want to see the points over the plot. What should I do?

I want to be able to see points and lines like this:

EDIT: This is the way I've generated the points:

for i in range(1,2000):
    [a, b] =  np.random.rand(2,1)*np.sign(np.random.randn(2,1))*2
    xx = np.append(xx, a)
    yy = np.append(yy, b)

I've noticed that if I write zz = f(xx,yy) + epsilon I can see the points. If epsilon = 0, then the points are, mathematically, on the surface, and I can't see them clearly, like in the first image. If epsilon > 0.05, I can see the points, but this means to move the points upper. I don't really like this solution. If a point is on a surface, the surface has priority, she surface appears to be over the point. I want my graphic to be the other way around.

Hammonds

Since you know your points are on the surface you can get away with shifting your points for plotting by a non-visible amount... floating point errors make your goal somewhat impossible anyway. Finally, the surface is drawn with planes, so even in an exact scenario your points would be hidden where the function is convex. – Entresol

Yes, I know I can do that, this is why I edited my question. However, it's hard to choose a good epsilon that would make all the points that "are" on the surface visible. Sometimes, for a small non-zero epsilon, some poins are still "under" the surface. This is why I am looking for a more elegant solution – Hammonds

Answer 1

Let me start by saying that what you want is somewhat ill-defined. You want to plot points exactly on an underlying surface in a way that they are always shown in the figure, i.e. just above the surface, without explicitly shifting them upwards. This is already a problem because

1. Floating-point arithmetic implies that your exact coordinates for the points and the surface can vary on the order of machine precision, so trying to rely on exact equalities won't work.
2. Even if the numbers were exact to infinite precision, the surface is drawn with a collection of approximating planes. This implies that your exact data points will be under the approximating surface where the function is convex.

The largest problem, however, is that 3d plotting in matplotlib is known to be unreliable when it comes to plotting multiple or complex objects in a figure. In particular, the renderer is inherently 2d, and it often runs into problems when trying to figure out the relative apparent position of objects. To overcome this one can either try hacking around the problem, or switching to something like mayavi with a proper 3d renderer.

Unfortunately, the zorder optional keyword argument is usually ignored by 3d axes objects. So the only thing I could come up with within pyplot is what you almost had, commented out: using ax.plot rather than ax.scatter. While the latter produces a plot shown in your first figure (with every scatter point hidden for some reason, regardless of viewing angle), the former leads to a plot shown in your second figure (where the points are visible). By removing the lines from the plotting style we almost get what you want:

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

But not quite: it quickly becomes evident that in this case the points are always visible, even when they should be hidden behind part of the surface:

# 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')

It's easy to see that the bump in the front should hide a huge chunk of points in the background, however these points are visible. This is exactly the kind of problem that pyplot has with complex 3d visualization. My bottom line is thus that you can't reliably do what you want using matplotlib. For what it's worth, I'm not sure how easy such a plot would be to comprehend anyway.

---

Just to end on a more positive note, here's how you could do this using mayavi (for which you need to install vtk which is best done via your package manager):

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

As you can see, the result is an interactive 3d plot where the data points on the surface are proper spheres. One can play around with the opacity and sphere scale settings in order to get a satisfactory visualization. Due to proper 3d rendering you can see an appropriate amount of each point irrespective of viewing angle.

Gallonage

Your response is excelent, I will look forward to use mayavi in the next future. Thank you very much! I have only one question. I didn't understand what you mean by "an interactive 3D plot". Can I rotate the figure/view from another angle easily using mayavi? (This would be amazing) – Hammonds

@Hammonds exactly, just like in matplotlib. Well, not exactly, because in matplotlib the z axis is always vertical, but in mayavi (vtk) you get arbitrary camera angles based on yaw/pitch/roll :) – Entresol

Answer 2

Andras Deak gave a very comprehensive answer discussing the problems/limitations of matplotlibs 3D plotting capabilities when it comes to the task at hand. At the end of his answer -- ending on a positive note -- he gave a solution using an alternative library.

I set out to try to find a hacky/specialized solution in matplotlib. Let me first address why. I wanted to plot trajectories on a 2D surface and I was set on using matplotlib. I use it for all my 2D plots and wanted to find a solution for this particular 3D plot application. The nice thing about matplotlibs 3D plots is that they are vectorized, since they are basically just 2D plots generated from projecting 3D elements onto a camera plane and overlaying them (drawing them in an order based on their distance to the camera). Rasterisation can be controlled individually for each element of the plot without affecting axes, labels ect. 'Real' 3D plotting libraries which use ray tracing usually do not have the capability to produce completely vectorized plots. I think mayavi is one example and I know Mathematica is also very limited in this regard.

Coming to my solution: I looked into the code for plot_surface which is ultimately based on Poly3DCollection to find out how matplotlib decides which polygons/elements on the surface to draw first. The method _do_3d_projection of Poly3DCollection orders the polygons projected onto the 2d camera plane by distance (of the original 3D object) to the camera. Elements far away from the camera are drawn first, elements close to the camera a drawn later. This works well enough for most plots to create a proper perspective (but the approach has limitations, see, e.g., mplot3d FAQ. This ordering however is the key to my solution. Given a set of points pts and a surface surf (which has to be plotted with show or savefig to have its camera/projection variables set):

1. The 2D projections segments_2d into the camera plane of all the 3D polygons in surf are computed including their ordering based on distance to the camera (stored in segments_idxs).
2. All points get associated to an element/polygon on the 3D surface.
3. The 2D projection into the camera plane of the 3D points is computed.
4. To decide wether a point is visible or not we check if it gets overlayed by a polygon drawn after the one it is associated to (from step 2.). To do so we use the contains_points method from matplotlib.path, see also the related question What's the fastest way of checking if a point is inside a polygon in python.
5. I included a dynamic updating (adapted from How to obscure a line behind a surface plot in matplotlib?). Warning: the code/plotting of surfaces with a lot of polygons can get quite slow.

Here the necessary code/minimal working example with the surface given by OP an sample points located at the unit circle.

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

The code is still a bit messy and it still does not work perfectly but here some results with 25*25=625 quads in the surface and 200 points on the unit circle. The points in red are the visible point and the points in green are the invisible ones (here plotted for illustrative purposes but for the sake of the initial question/problem one would of curse omit them from the plot). Some points should clearly be visible but are detected as invisible. I am not sure yet what goes wrong there but for me this limited miss detections are not too problematic since I ultimately want to plot lines/trajectories of a lot of (arbitrary dense) points. If the misses do not cluster I can live with a few missing ones.

Another inherent problem/limitation is that the current approach has no real notion of wether points are above or below a surface, meaning that when looking from below the surface the points above/on the surface are visible. Here an example for this behavior:

This ties in to the point already raised by Andras Deak that the present problem is somewhat ill-defined or at least ambiguous without additional restrictions. One could for example demand, that all points are placed on the surface face pointing towards the camera. Implementing this in the present approach is difficult. In geometrical terms the current implementation places finite-sized balls on infinitesimal polygons making them visible from both sides (which might actually be viable/desired in certain use cases).

The code is still a work in progress and if I find significant improvements I might update this answer. Comments on the general approach and/or the implementation are of curse very welcome. I am by no means a python expert (I use it almost exclusively for plotting and related very light data processing) so their might be significant room for improvement in terms of code performance and scope.

Indigested

Good job trying to find a matplotlib solution. Do I understand correctly that you basically have to reimplement part of the renderer? I've actually switched from mayavi to pyvista, I should probably update my half/answer with that :) – Entresol

I did not really reimplement the renderer but I use parts of it. I extract the order in which the elements on the surface are rendered and based on this order I decide wether a point on the surface should be visible or not by checking if it gets overlayed by elements of the surface rendered later. A way to generalize this approach to other composite 3D plots would be to split all plot elements in the smallest possible sub-elements (polygons, points,...) and then render all those elements in the order determined by distance to camera. – Indigested

I see, thanks for the explanation. Yeah, that makes sense: the rendering bug comes from complex surfaces being either behind or in front of each other. If you break down everything to small primitives, the problem goes away. I wonder about the memory overhead about that... but it does sound like a plausible workaround. – Entresol

Answer 3

From your graphs, it looks like you are willing to show the path to a non linear optimizer local solution, so that I think you should consider plotting a wireframe over the contour plot:

...

ax.scatter(xx, yy, zz, c='r', marker='o') ### this will plot only the points, not the lines

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

ax.plot_wireframe(xx,yy,zz) ### this will plot the lines and possibly the points.

...

were (xx,yy,zz) contains the path to the local maximum, obtained from the non linear regretion method, as I'm supposing.

Navarrette

Answer 4

I think you can solve this problem by manually controlling zorder with Matplotlib 3.5.0 by:

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d',computed_zorder=False)
ax.plot_surface(X, Y, Z, cmap='jet', rstride=1, cstride=1, alpha=0.5)
ax.set_proj_type('ortho')
ax.plot(x1,y1,z1,'-k',linewidth=3,zorder=1)

The result will be that the line will always show on top of the surface plot, even if it would be obscured by part of the surface. Additional lines can also be added with their own zorder. This solution is from: 3D parametric curve in Matplotlib does not respect zorder. Workaround?

Dogtrot

I'm not in front of a laptop now, but I suspect this will always show the points on top. See my second code block for a case where it can be obvious if this can be an issue. – Entresol

Yes, this is correct and is what the questioner asked for. To vary the point location(s), the zorder can also be specified for the surface plot and any other plotted lines. – Dogtrot

Ah, I think I see where you're coming from. My interpretation was that the asker only wants the line not to be obstructed by the surface piece it's situated on. To me the intuitive behaviour would be that if other, unrelated parts of the surface obscure the line then it should be hidden. I would recommend making your assumption on this front explicit and clear in your answer, because than it's clear how and why it differs from the existing answers, and to whom it might be useful. – Entresol

Thanks, I have added the clarification you suggested. – Dogtrot