Extracting coordinates from meshgrid data

meshgrid is probably what you need, but the shape of the array returned by meshgrid is where it gets confusing. Meshgrid returns three coordinate arrays, all the same shape. The shape of each of xv, yv, zv is (len(x), len(y), len(z)). So, to extract the coordinate at the corner (0, 2, 1), you would write xv[0, 2, 1], yv[0, 2, 1], zv[0, 2, 1]

To extract all of the subcubes' corners' coordinates, it helps to observe that, because of the way the arrays returned by meshgrid are ordered sequentially, xv[:-1, :-1, :-1] returns only the x-coordinates of the near-left-bottom corners of each subcube. Likewise, xv[1:, 1:, 1:] returns the far-right-top corners of each subcube. The other six corners are given by the other six combinations of the slices :-1 and 1: (xv[:-1, 1:, :-1] gives the far-left-top corner, for example).

So, iterate through all eight combinations of :-1 and 1: to get eight parallel arrays of three parallel arrays of x, y, z coordinates for the eight corners of all len(x)-1 * len(y-1) * len(z-1) subcubes. (If you need your subcube corner coordinate arrays to be in a particular shape or axis order, or if you want to use a single index to specify the subcube rather than three, use rollaxis, swapaxis and shape as needed.)

import numpy as npimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dimport itertoolsnx, ny, nz = (3,3,3)x = np.linspace(0, 10, nx)y = np.linspace(0, 10, ny)z = np.linspace(0, 10, nz)xv, yv, zv = np.meshgrid(x, y, z, indexing='xy')slices = slice(None, -1), slice(1, None)cornerSlices = list(itertools.product(slices, slices, slices))corners = np.array([(xv[s], yv[s], zv[s]) for s in cornerSlices])# The shape of `corners` is `(len(cornerSlices), 3, len(x-1), len(y-1), len(z-1)`# The axes of `corners` represent, in the same order: the corner index; the cartesian # coordinate axis (the index into [x, y, z]); the x, y, and z indexes of the subcube. # Plot the first subcube (subcube 0, 0, 0)fig = plt.figure()ax = fig.add_subplot(111, projection='3d')subcube = corners[:, :, 0, 0, 0]subcubeX = subcube [:, 0]subcubeY = subcube [:, 1]subcubeZ = subcube [:, 2]ax.scatter(subcubeX , subcubeY , subcubeZ , c='g', marker='^')

There's invariably a way to get the indexes into xv, yv, zv instead of getting the values, since the values are duplicated quite a few times in the corners array. It would involve slicing arrays of indexes into xv, yv, zv instead of slicing the arrays themselves. My head is already spinning after getting this far into the ndarray voodoo, so I'll leave that as an exercise.