A 3-D grid of regularly spaced points

I want to create a list containing the 3-D coords of a grid of regularly spaced points, each as a 3-element tuple. I'm looking for advice on the most efficient way to do this.

In C++ for instance, I simply loop over three nested loops, one for each coordinate. In Matlab, I would probably use the meshgrid function (which would do it in one command). I've r开发者_Go百科ead about meshgrid and mgrid in Python, and I've also read that using numpy's broadcasting rules is more efficient. It seems to me that using the zip function in combination with the numpy broadcast rules might be the most efficient way, but zip doesn't seem to be overloaded in numpy.

Use ndindex:

import numpy as np
ind=np.ndindex(3,3,2)
for i in ind:
    print(i)

# (0, 0, 0)
# (0, 0, 1)
# (0, 1, 0)
# (0, 1, 1)
# (0, 2, 0)
# (0, 2, 1)
# (1, 0, 0)
# (1, 0, 1)
# (1, 1, 0)
# (1, 1, 1)
# (1, 2, 0)
# (1, 2, 1)
# (2, 0, 0)
# (2, 0, 1)
# (2, 1, 0)
# (2, 1, 1)
# (2, 2, 0)
# (2, 2, 1)

Instead of meshgrid and mgrid, you can use ogrid, which is a "sparse" version of mgrid. That is, only the dimension along which the values change are filled in. The others are simply broadcast. This uses much less memory for large grids than the non-sparse alternatives.

For example:

>>> import numpy as np
>>> x, y = np.ogrid[-1:2, -2:3]
>>> x
array([[-1],
      [ 0],
      [ 1]])
>>> y
array([[-2, -1,  0,  1,  2]])
>>> x**2 + y**2
array([[5, 2, 1, 2, 5],
      [4, 1, 0, 1, 4],
      [5, 2, 1, 2, 5]])

I would say go with meshgrid or mgrid, in particular if you need non-integer coordinates. I'm surprised that Numpy's broadcasting rules would be more efficient, as meshgrid was designed especially for the problem that you want to solve.

for multi-d (greater than 2) meshgrids, use numpy.lib.index_tricks.nd_grid like so:

import numpy
grid = numpy.lib.index_tricks.nd_grid()
g1 = grid[:3,:3,:3]
g2 = grid[0:1:0.5, 0:1, 0:2]
g3 = grid[0:1:3j, 0:1:2j, 0:2:2j]

where g1 has x values of [0,1,2] and g2 has x values of [0,.5], and g3 has x values of [0.0,0.5,1.0] (the 3j defining the step count instead of the step increment. see the documentation for more details.

Here's an efficient option similar to your C++ solution, which I've used for exactly the same purpose:

import numpy, itertools, collections
def grid(xmin, xmax, xstep, ymin, ymax, ystep, zmin, zmax, zstep):
    "return nested tuples of grid-sampled coordinates that include maxima"
    return collections.deque( itertools.product( 
        numpy.arange(xmin, xmax+xstep, xstep).tolist(),
        numpy.arange(ymin, ymax+ystep, ystep).tolist(),
        numpy.arange(zmin, zmax+zstep, zstep).tolist() ) )

Performance is best (in my tests) when using a.tolist(), as shown above, but you can use a.flat instead and drop the deque() to get an iterator that will sip memory. Of course, you can also use a plain old tuple() or list() instead of deque() for a slight performance penalty (again, in my tests).