我刚刚注意到numpy中的文档提供了一种更快的方法来做到这一点：

X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])

使用链接的meshgrid2函数并将“ravel”Map到结果网格，可以轻松地将其推广到更多维度。

g = meshgrid2(x, y, z)
positions = np.vstack(map(np.ravel, g))

对于每个轴上有1000个刻度的3D数组，结果比zip方法快约35倍。
来源：http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html#scipy.stats.gaussian_kde
要比较这两种方法，请考虑以下代码部分：
创建众所周知的刻度线，这将有助于创建网格。

In [23]: import numpy as np

In [34]: from numpy import asarray

In [35]: x = np.random.rand(100,1)

In [36]: y = np.random.rand(100,1)

In [37]: z = np.random.rand(100,1)

为meshgrid定义mgilson链接到的函数：

In [38]: def meshgrid2(*arrs):
   ....:     arrs = tuple(reversed(arrs))
   ....:     lens = map(len, arrs)
   ....:     dim = len(arrs)
   ....:     sz = 1
   ....:     for s in lens:
   ....:        sz *= s
   ....:     ans = []
   ....:     for i, arr in enumerate(arrs):
   ....:         slc = [1]*dim
   ....:         slc[i] = lens[i]
   ....:         arr2 = asarray(arr).reshape(slc)
   ....:         for j, sz in enumerate(lens):
   ....:             if j != i:
   ....:                 arr2 = arr2.repeat(sz, axis=j)
   ....:         ans.append(arr2)
   ....:     return tuple(ans)

创建网格并计算两个函数的时间。

In [39]: g = meshgrid2(x, y, z)

In [40]: %timeit pos = np.vstack(map(np.ravel, g)).T
100 loops, best of 3: 7.26 ms per loop

In [41]: %timeit zip(*(x.flat for x in g))
1 loops, best of 3: 264 ms per loop

你的网格点总是整数吗？如果是，你可以使用numpy.ndindex

print list(np.ndindex(2,2))

更高维度：

print list(np.ndindex(2,2,2))

不幸的是，这不符合OP的要求，因为积分假设（从0开始）不满足。我只会留下这个答案，以防其他人在寻找同样的东西，这些假设是真的。
另一种方法是使用zip：

g = np.meshgrid([0,1],[0,1])
zip(*(x.flat for x in g))

这部分可以很好地扩展到任意维度。不幸的是，np.meshgrid不能很好地扩展到多维，因此这部分需要解决，或者（假设它工作），您可以使用此SO answer创建自己的ndmeshgrid函数。

另一种方法是：

np.indices((2,2)).T.reshape(-1,2)

其可以推广到更高的维度，例如：

In [60]: np.indices((2,2,2)).T.reshape(-1,3)
Out[60]:
array([[0, 0, 0],
       [1, 0, 0],
       [0, 1, 0],
       [1, 1, 0],
       [0, 0, 1],
       [1, 0, 1],
       [0, 1, 1],
       [1, 1, 1]])

要获得从0到1的网格坐标，整形可以完成这项工作。这里有2D和3D的例子。也适用于浮动。

grid_2D = np.mgrid[0:2:1, 0:2:1]
points_2D = grid_2D.reshape(2, -1).T

grid_3D = np.mgrid[0:2:1, 0:2:1, 0:2:1]
points_3D = grid_3D.reshape(3, -1).T

一个简单的3D示例（我猜可以扩展到N维，但要注意最终的维度和RAM使用情况）：

import numpy as np

ndim = 3
xmin = 0.
ymin = 0.
zmin = 0.
length_x = 1000.
length_y = 1000.
length_z = 50.
step_x = 1.
step_y = 1.
step_z = 1.

x = np.arange(xmin, length_x, step_x)
y = np.arange(ymin, length_y, step_y)
z = np.arange(zmin, length_z, step_z)
%timeit xyz = np.array(np.meshgrid(x, y, z)).T.reshape(-1, ndim)

在：2.76 s ± 185 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
其产生：

In [2]: xyx
Out[2]: 
array([[  0.,   0.,   0.],
       [  0.,   1.,   0.],
       [  0.,   2.,   0.],
       ...,
       [999., 997.,  49.],
       [999., 998.,  49.],
       [999., 999.,  49.]])

In [4]: xyz.shape
Out[4]: (50000000, 3)

Python 3.6.9
Numpy：1.19.5

我用下面的方法将meshgrid转换成MX2数组。同时将向量列表改为迭代器可以使它变得非常快。

import numpy as np

# Without iterators
x_vecs = [np.linspace(0,1,1000), np.linspace(0,1,1000)]

%timeit np.reshape(np.meshgrid(*x_vecs),(2,-1)).T
6.85 ms ± 93.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# With iterators
x_vecs = iter([np.linspace(0,1,1000), np.linspace(0,1,1000)])

%timeit np.reshape(np.meshgrid(*x_vecs),(2,-1)).T
5.78 µs ± 172 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

对于使用发生器的N-D阵列

vec_dim = 3
res = 100

# Without iterators
x_vecs = [np.linspace(0,1,res) for i in range(vec_dim)]

>>> %timeit np.reshape(np.meshgrid(*x_vecs),(vec_dim,-1)).T
11 ms ± 124 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# With iterators
x_vecs = (np.linspace(0,1,res) for i in range(vec_dim))

>>> %timeit np.reshape(np.meshgrid(*x_vecs),(vec_dim,-1)).T
5.54 µs ± 32.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

这是N-D的情况。
我们的想法是，我们将网格组合成一个数组，使某个东西n x dim 1 x.. dimn;然后，我们将其移动到末尾（使用moveaxis -转置的概括），并使除了最后一个dims之外的所有内容变平。

grids = np.meshgrid([1, 2, 3], [10, 20, 30], [100, 200, 300])
 
def grids_to_points(grids):
    return np.moveaxis(np.array(grids), 0, grids[0].ndim).reshape(-1, len(grids))

grids_to_points(grids)