Extracting coordinates from meshgrid data
I have a cubic grid as shown in the image below.
I would like to list the vertices of each subcube, so I would end up nesting the subcube list with the corresponding vertex list.
My initial attempt was to use a generator,
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
dims = [9,9,9]
spacer = 3
subBoxCoords = np.array([(x, y, z) for x in range(0, dims[0], spacer) for y in range(0, dims[1], spacer) for z in range(0, dims[2], spacer)])
ax.scatter(subBoxCoords[:,0], subBoxCoords[:,1], subBoxCoords[:,2], c='k', marker='o')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
This gives me the shape I want, but the coordinates are ordered in such a way that removing vertices in subsamples is not straight forward. Also I would like to generalize this to blocks of arbitrary dimension, so spacing hardcoding is not a solution.
So then I decided to use meshgrid
,
nx,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') ax.scatter(xv, yv, zv, c='g', marker='^')
This seems to be a very powerful way to achieve what I want, but I am confused. Are there direct access tops in the meshgrid
way vertex(x,y,z)
? Or even a direct way to extract subcubes?
The solution seems painfully close to me, but I just can't figure it out!
source to share
meshgrid
is probably what you want, but the shape of the array returned meshgrid
is where it gets confusing. Meshgrid returns three coordinate arrays, all the same. The shape of each xv, yv, zv
is equal (len(x), len(y), len(z))
. So, to extract the coordinate in a corner (0, 2, 1)
, you have to writexv[0, 2, 1], yv[0, 2, 1], zv[0, 2, 1]
To extract all the coordinates of the corners of the subcubes, it's helpful to note that due to the ordering of the arrays returned meshgrid
, xv[:-1, :-1, :-1]
only returns the x-coordinates of the near-left corner, the lower corners of each subcube. Likewise, xv[1:, 1:, 1:]
returns the top-right corners of each subcube. The other six corners are given by the other six combinations of slices :-1
and 1:
( xv[:-1, 1:, :-1]
gives, for example, the upper left corner).
So, repeat all eight combinations of :-1
and 1:
to get eight parallel arrays from three parallel arrays of x, y, z coordinates for the eight corners of all subcubes len(x)-1 * len(y-1) * len(z-1)
. (If you need your angular subcube coordinate arrays in a specific shape or axis order, or if you want to use one index to indicate the subcube instead of three, use rollaxis
, swapaxis
and shape
as needed.)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import itertools
nx, 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 is always a way to get the indices in xv, yv, zv
instead of getting the values, since the values are duplicated corners
multiple times in the array . This will involve cutting the arrays of indices in xv, yv, zv
instead of cutting the arrays themselves. My head is already spinning after going this far in the ndarray voodoo, so I'll leave that as an exercise.
source to share