[Doc]: Gallery example showing 3D slice planes

[scottshambaugh]

@timhoffm suggested in #3919 (comment) that we make a gallery example for plotting X Y Z planar slices through a volume.

His result and its code:

import matplotlib.pyplot as plt
import numpy as np

def plot_quadrants(ax, array, fixed_coord, cmap):
    """For a given 3d *array* plot a plane with *fixed_coord*, using four individual quadrants."""
    nx, ny, nz = array.shape
    index = {
        'x': (nx//2, slice(None), slice(None)),
        'y': (slice(None), ny//2, slice(None)),
        'z': (slice(None), slice(None), nz//2),
    }[fixed_coord]
    plane_data = array[index]

    n0, n1 = plane_data.shape
    quadrants = [
        plane_data[:n0//2, :n1//2],
        plane_data[:n0//2, n1//2:],
        plane_data[n0//2:, :n1//2],
        plane_data[n0//2:, n1//2:]
    ]

    min_val = array.min()
    max_val = array.max()
    
    cmap = plt.get_cmap(cmap)
    
    for i, quadrant in enumerate(quadrants):
        facecolors = cmap((quadrant - min_val) / (max_val-min_val))
        if fixed_coord == 'x':
            Y, Z = np.mgrid[0:ny//2, 0:nz//2]
            X = nx//2 * np.ones_like(Y)
            Y_offset = (i // 2) * ny//2
            Z_offset = (i % 2) * nz//2
            ax.plot_surface(X, Y + Y_offset, Z + Z_offset, rstride=1, cstride=1, facecolors=facecolors, shade=False)
        elif fixed_coord == 'y':
            X, Z = np.mgrid[0:nx//2, 0:nz//2]
            Y = ny//2 * np.ones_like(X)
            X_offset = (i // 2) * nx//2
            Z_offset = (i % 2) * nz//2
            ax.plot_surface(X + X_offset, Y, Z + Z_offset, rstride=1, cstride=1, facecolors=facecolors, shade=False)
        elif fixed_coord == 'z':
            X, Y = np.mgrid[0:nx//2, 0:ny//2]
            Z = nz//2 * np.ones_like(X)
            X_offset = (i // 2) * nx//2
            Y_offset = (i % 2) * ny//2
            ax.plot_surface(X + X_offset, Y + Y_offset, Z, rstride=1, cstride=1, facecolors=facecolors, shade=False)


def figure_3D_array_slices(array, cmap=None):
    """Plot a 3d array using three intersecting centered planes."""
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.set_box_aspect(array.shape)
    plot_quadrants(ax, array, 'x', cmap=cmap)
    plot_quadrants(ax, array, 'y', cmap=cmap)
    plot_quadrants(ax, array, 'z', cmap=cmap)
    return fig, ax


nx, ny, nz = 70, 100, 50
r_square = (np.mgrid[-1:1:1j*nx, -1:1:1j*ny, -1:1:1j*nz]**2).sum(0)

figure_3D_array_slices(r_square, cmap='viridis_r')
plt.show()

[story645]

Um can it be split up into three subplots, one for each axes?

Just thinking through that while pretty, it might be a bit hard to follow from a "what do I do to make this slice for one axes"?

[timhoffm]

No. The point here is to illustrate multiple intersecting planes in one subplot, which requires the quadrant subdivision for correct drawing order.

This is somewhat inherently complex (with current API (?)) due to the involved indexing. Note also that indexing will get even worse if you don’t split in the center (which I thus have refained from to generalize to). I rather see this as an illustration/proof of concept how you can work around limitations of the 3d rendering.

Separate planes (well, one would be enough then) is a simpler example that should be done separately. Also note that drawing a single color-coded plane is currently already quite a hassle because you have to use plot_surface. Some quadmesh like function to draw axes-aligned planes from 2D arrays would be reasonable. But that’s even a bigger topic.

[story645]

This is somewhat inherently complex (with current API (?)) due to the involved indexing.

My concern more is that in this case the complexity/way that it's written may make it hard to unravel and tweak for someone who wanted to make something similar but not identical.

I wonder if this would work better as a tutorial where these stages were also included:

• Separate planes (well, one would be enough then) is a simpler example that should be done separately.
• Also note that drawing a single color-coded plane is currently already quite a hassle because you have to use plot_surface

And build up to the figure above.

[timhoffm]

My concern more is that in this case the complexity/way that it's written may make it hard to unravel and tweak for someone who wanted to make something similar but not identical.

It is difficult, but better than nothing at all. We don't have the capacity to pamper users for every edge case.

I'm -0.5 on a tutorial. First, I'd like to have both the separate plane and the intersecting plane show up in the gallery overview. Second, going into an extensive tutorial for working around limitations and insufficiencies of our API feels the wrong way. It would be better to invest that time into improving the API, i.e. the quadmesh-like parallel-axes plane.

But as always, if you are convinced that a tutorial is worth the effort, you are free to create one.

[story645]

We don't have the capacity to pamper users for every edge case.

Maybe the solution is more code comments and not wrapping the figure code in figure_3D_array_slices, but my concern is that the example as written might be hard to modify and cross apply to a slightly different use case, which I think is the core purpose of examples and not an edge case. And a very consistent complaint about the examples is that they tend to do too much in a way where it's hard to see what code is responsible for which aspect of the visualization.

Second, going into an extensive tutorial for working around limitations and insufficiencies of our API feels the wrong way. It would be better to invest that time into improving the API, i.e. the quadmesh-like parallel-axes plane.

I think of it more as a tutorial showing how to use the existing API to build something that the API doesn't provide out of the box, but yeah if a new function makes more sense than I'm not opposed to that.