Problems With Contours Using Python's matplotlib 3D API

Question

I'm trying to do something similar to this 3D example from the docs, but with a point cloud instead of a smooth surface. The example projects 2D contours onto each of the three coordinate planes. This shows that I'm able to do that onto the xy-plane.

When I try doing the same onto the other two planes, I get either a weird contour collapsed down to a thin strip,

or an exception way beyond my reach in the bowels of matplotlib.

Traceback (most recent call last):
  File ".../matplotlib/backends/backend_qt5.py", line 519, in _draw_idle
    self.draw()
  File ".../matplotlib/backends/backend_agg.py", line 433, in draw
    self.figure.draw(self.renderer)
  File ".../matplotlib/artist.py", line 55, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File ".../matplotlib/figure.py", line 1475, in draw
    renderer, self, artists, self.suppressComposite)
  File ".../matplotlib/image.py", line 141, in _draw_list_compositing_images
    a.draw(renderer)
  File ".../mpl_toolkits/mplot3d/axes3d.py", line 281, in draw
    reverse=True)):
  File ".../mpl_toolkits/mplot3d/axes3d.py", line 280, in <lambda>
    key=lambda col: col.do_3d_projection(renderer),
  File ".../mpl_toolkits/mplot3d/art3d.py", line 226, in do_3d_projection
    self._segments3d]
  File ".../mpl_toolkits/mplot3d/art3d.py", line 225, in <listcomp>
    proj3d.proj_trans_points(points, renderer.M) for points in
  File ".../mpl_toolkits/mplot3d/proj3d.py", line 188, in proj_trans_points
    xs, ys, zs = zip(*points)
ValueError: not enough values to unpack (expected 3, got 0)

Here's an example of the problem. This version works. Uncomment one or both of the calls to ax.contour() in the plot() function to see the weird contours, or more likely, the exception.

import math
import sys
import matplotlib as mpl
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d
import numpy as np
import scipy.spatial

#np.random.seed(4)
numPts = 1000                                       # number of points in cloud
scatter = False                         # adds scatter plot to show point cloud


def main():
    (pts, f) = mkData()                                # create the point cloud
    tris = mkTris(pts)                                         # triangulate it
    plot(pts, tris, f)                                                # plot it


def plot(pts, tris, f):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")

    cmap = plt.get_cmap("coolwarm")

    collec = ax.plot_trisurf(tris, pts[:, 2], cmap=cmap)
    colors = np.mean(f[tris.triangles], axis=1)
    collec.set_array(colors)
    collec.autoscale()

    (xr, yr, zr, xMin, yMin, zMin, zMax) = resample(ax, tris, f)
    ax.set_zlim(zMin, zMax)      # default always includes zero for some reason

    #
    # Uncomment one or both of these lines to see the problem.
    #
    ax.contour(xr, yr, zr, 10, zdir="z", cmap=cmap, offset=zMin)
    #ax.contour(xr, yr, zr, 10, zdir="x", cmap=cmap, offset=xMin)
    #ax.contour(xr, yr, zr, 10, zdir="y", cmap=cmap, offset=yMin)

    if scatter:
        ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], alpha=0.1)

    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    fig.colorbar(collec, shrink=0.5, aspect=5)
    plt.show()


def mkData():
    """
    Create a random point cloud near a random plane, and define a function on
    the plane for colors and contours.
    """
    offset = 1                   # generate points near a unit square, xy-plane
    pts = 2 * np.random.rand(numPts, 3) - 1
    pts[:, 2] = offset * (2 * np.random.rand(numPts) - 1)

    x = 2 * np.ravel(pts[:, 0])
    y = 2 * np.ravel(pts[:, 1])
    f = x * np.exp(-x**2 - y**2)      # just some function for colors, contours

    width = 100                       # scale unit square to a larger rectangle
    height = 20
    pts[:, 0] *= width
    pts[:, 1] *= height

    (e1, e2, e3) =[2 * np.pi * np.random.rand() for _ in range(3)]
    (c1, s1) = (math.cos(e1), math.sin(e1))           # rotate scaled rectangle
    (c2, s2) = (math.cos(e2), math.sin(e2))
    (c3, s3) = (math.cos(e3), math.sin(e3))
    Ta2b = np.array((                     # do not have scipy.spatial.transform
        [               c1  *c2,       s2,              -s1 * c2],
        [s1 * s3 - c1 * s2 * c3,  c2 * c3, c1 *s3 + s1 * s2 * c3],
        [s1 * c3 + c1 * s2 * s3, -c2 * s3, c1 *c3 - s1 * s2 * s3]))

    pts = (Ta2b @ pts.T).T

    dist = 500                                 # translate away from the origin
    Ra2bNb = dist * (2 * np.random.rand(3, 1) - 1)

    pts += Ra2bNb.T

    return (pts, f)


def mkTris(pts):                                  # triangulate the point cloud
    u = np.ravel(pts[:, 0])
    v = np.ravel(pts[:, 1])

    try:
        return mpl.tri.Triangulation(u, v)
    except (ValueError, RuntimeError) as ex:
        sys.exit(f"Unable to compute triangulation: {ex}.")


def resample(ax, tris, f):         # resample triangulation onto a regular grid
    (xMin, xMax) = ax.get_xlim()
    (yMin, yMax) = ax.get_ylim()
    (zMin, zMax) = ax.get_zlim()

    x = np.linspace(xMin, xMax, 30)
    y = np.linspace(yMin, yMax, 30)

    (xm, ym)=np.meshgrid(x, y)

    interp = mpl.tri.triinterpolate.LinearTriInterpolator(tris, f)
    zm = interp(xm, ym)

    return (xm, ym, zm, xMin, yMin, zMin, zMax)

if __name__ == "__main__":
    main()

This is with matplotlib 2.2.2 and 3.1.1. Thanks for any help you can provide to get contours on all three planes, like the demo.

Jim

Answer 1

A matplotlib developer pointed out that the resampling was wrong. After fixing that, here's the corrected plot.

For coordinate planes that see the data edge-on, like the yz-plane in this case, the contours can look a little wonky. That's expected, since the data can approach being multi-valued. The xz-plane contours look pretty ragged too. I suspect both problems would improve by triangulating and contouring each plane individually, instead of favoring the xy-plane as done here.

Here's the fixed test script. The only important changes were in plot() and resample().

import math
import sys
import matplotlib as mpl
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d
import numpy as np
import scipy.spatial

#np.random.seed(4)
numPts = 1000                                       # number of points in cloud
numGrid = 120      # number of points in meshgrid used in resample for contours
scatter = False                         # adds scatter plot to show point cloud


def main():
    (pts, f) = mkData()                                # create the point cloud
    tris = mkTris(pts)                                         # triangulate it
    plot(pts, tris, f)                                                # plot it


def plot(pts, tris, f):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")

    cmap = plt.get_cmap("coolwarm")

    collec = ax.plot_trisurf(tris, pts[:, 2], cmap=cmap)
    colors = np.mean(f[tris.triangles], axis=1)
    collec.set_array(colors)
    collec.autoscale()

    (xr, yr, zr, fr, xMin, xMax, yMin, yMax, zMin, zMax) = resample(ax, tris,
       pts, f)

    ax.set_xlim(xMin, xMax)      # default always includes zero for some reason
    ax.set_ylim(yMin, yMax)
    ax.set_zlim(zMin, zMax)

    ax.contour(xr, yr, fr, 10, zdir="z", cmap=cmap, offset=zMin)
    ax.contour(fr, yr, zr, 10, zdir="x", cmap=cmap, offset=xMin)
    ax.contour(xr, fr, zr, 10, zdir="y", cmap=cmap, offset=yMax)

    if scatter:
       ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], alpha=0.1)

    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    fig.colorbar(collec, shrink=0.5, aspect=5)
    plt.show()


def mkData():
    """
    Create a random point cloud near a random plane, and define a function on
    the plane for colors and contours.
    """
    offset = 1                   # generate points near a unit square, xy-plane
    pts = 2 * np.random.rand(numPts, 3) - 1
    pts[:, 2] = offset * (2 * np.random.rand(numPts) - 1)

    x = 2 * np.ravel(pts[:, 0])
    y = 2 * np.ravel(pts[:, 1])
    f = x * np.exp(-x**2 - y**2)      # just some function for colors, contours

    width = 100                       # scale unit square to a larger rectangle
    height = 20
    pts[:, 0] *= width
    pts[:, 1] *= height

    (e1, e2, e3) =[2 * np.pi * np.random.rand() for _ in range(3)]
    (c1, s1) = (math.cos(e1), math.sin(e1))           # rotate scaled rectangle
    (c2, s2) = (math.cos(e2), math.sin(e2))
    (c3, s3) = (math.cos(e3), math.sin(e3))
    Ta2b = np.array((                     # do not have scipy.spatial.transform
        [               c1  *c2,       s2,              -s1 * c2],
        [s1 * s3 - c1 * s2 * c3,  c2 * c3, c1 *s3 + s1 * s2 * c3],
        [s1 * c3 + c1 * s2 * s3, -c2 * s3, c1 *c3 - s1 * s2 * s3]))

    pts = (Ta2b @ pts.T).T

    dist = 500                                 # translate away from the origin
    Ra2bNb = dist * (2 * np.random.rand(3, 1) - 1)

    pts += Ra2bNb.T

    return (pts, f)


def mkTris(pts):
    "Triangulate the point cloud."
    u = np.ravel(pts[:, 0])
    v = np.ravel(pts[:, 1])

    try:
        return mpl.tri.Triangulation(u, v)
    except (ValueError, RuntimeError) as ex:
        sys.exit(f"Unable to compute triangulation: {ex}.")


def resample(ax, tris, pts, f):
    "Resample the triangulation onto a regular grid for contours."
    (xMin, xMax) = ax.get_xlim()
    (yMin, yMax) = ax.get_ylim()
    (zMin, zMax) = ax.get_zlim()

    x = np.linspace(xMin, xMax, numGrid)
    y = np.linspace(yMin, yMax, numGrid)

    (xm, ym)=np.meshgrid(x, y)

    fInterp = mpl.tri.CubicTriInterpolator(tris, f)
    fm = fInterp(xm, ym)

    zInterp = mpl.tri.CubicTriInterpolator(tris, pts[:,2])
    zm = zInterp(xm, ym)

    return (xm, ym, zm, fm, xMin, xMax, yMin, yMax, zMin, zMax)

if __name__ == "__main__":
    main()