How to plot projections of a 3D surface onto coordinate planes?

edit

Consider the following code:

from sage.plot.plot3d.index_face_set import IndexFaceSet
from sage.plot.plot3d.texture import Texture

def project_xy_xz_yz(surf, position):
    oxy,oxz,oyz = position

    # Initialize surf as an IndexFaceSet
    surf.triangulate()

    # Vertices projected onto planes parallel to xy, xz and yz
    vxy_list = [(v[0],v[1],oxy) for v in surf.vertices()]
    vxz_list = [(v[0],oxz,v[2]) for v in surf.vertices()]
    vyz_list = [(oyz,v[1],v[2]) for v in surf.vertices()]

    # Compute projections. We have to consider whether or not
    # each face has its own color
    f_list = surf.index_faces()
    if surf.has_local_colors():
        f_data = surf.index_faces_with_colors()
        col_list = [Texture(color=item[1]) for item in f_data]
        pxy = IndexFaceSet(f_list,vxy_list,texture_list=col_list)
        pxz = IndexFaceSet(f_list,vxz_list,texture_list=col_list)
        pyz = IndexFaceSet(f_list,vyz_list,texture_list=col_list)
    else:
        texture = surf.get_texture()
        pxy = IndexFaceSet(f_list,vxy_list)
        pxy.set_texture(texture)
        pxz = IndexFaceSet(f_list,vxz_list)
        pxz.set_texture(texture)
        pyz = IndexFaceSet(f_list,vyz_list)
        pyz.set_texture(texture)

    # End function
    return [pxy, pxz, pyz]

This code defines a function, project_xy_xz_yz, which computes the projection of a surface, given by the first argument, onto planes parallel to the $xy-$, $xz-$ and $yz-$planes. Specifically, if the second argument is [a,b,c], we get the projections onto the planes $z=a$, $y=b$ and $x=c$. The functions used in the code are mostly documented here.

Let us test this function. Let us consider the surface $$ z=\frac{2x+3y^2}{(x^2+y^2+1)^2} $$ and plot it in the rectangle $[-2,2]\times[-2,2]$. First, the surface will have a uniform color and will be a bit transparent. We add the projections onto planes parallel to $xz$ and $yz$, which inherit the surface color. The orthographic projection will allow you to check that the projections match the surface.

var("x,y")
def f(x,y):
    return (2*x+3*y^2) / (x^2+y^2+1)^2

s1 = plot3d(f(x,y), (x,-2,2), (y,-2,2), color="lightgreen", opacity=0.6)
p1 = project_xy_xz_yz(s1,[-1,-3,-3])
show(s1+p1[1]+p1[2], aspect_ratio=[1,1,2], projection="orthographic")

Now, let us get a full color picture and the three projections. Note that you can freely choose where locate the projection planes.

z1, z2 = -0.6, 0.8
def col(x,y):
    return float((f(x,y)-z1)/(z2-z1))

s2 = plot3d(f(x,y), (x,-2,2), (y,-2,2), color=(col,colormaps.coolwarm), plot_points=80)
p2 = project_xy_xz_yz(s2,[-1,3,3])
show(s2+sum(p2), aspect_ratio=[1,1,2], frame=False)

I think that the function extrema are easily seen in the projections.

You can play this examples in this SageMath cell.