How to convert the output of meshgrid to the corresponding array of points?

Question

I want to create a list of points that would correspond to a grid. So if I want to create a grid of the region from (0, 0) to (1, 1), it would contain the points (0, 0), (0, 1), (1, 0) and (1, 0).

I know that that this can be done with the following code:

g = np.meshgrid([0,1],[0,1])
np.append(g[0].reshape(-1,1),g[1].reshape(-1,1),axis=1)

Yielding the result:

array([[0, 0],
       [1, 0],
       [0, 1],
       [1, 1]])

My question is twofold:

1. Is there a better way of doing this?
2. Is there a way of generalizing this to higher dimensions?

Answer 1

I just noticed that the documentation in numpy provides an even faster way to do this:

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

This can easily be generalized to more dimensions using the linked meshgrid2 function and mapping 'ravel' to the resulting grid.

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

The result is about 35 times faster than the zip method for a 3D array with 1000 ticks on each axis.

Source: http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html#scipy.stats.gaussian_kde

To compare the two methods consider the following sections of code:

Create the proverbial tick marks that will help to create the grid.

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)

Define the function that mgilson linked to for the meshgrid:

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)

Create the grid and time the two functions.

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

Answer 2

Are your gridpoints always integral? If so, you could use numpy.ndindex

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

Higher dimensions:

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

---

Unfortunately, this does not meet the requirements of the OP since the integral assumption (starting with 0) is not met. I'll leave this answer only in case someone else is looking for the same thing where those assumptions are true.

---

Another way to do this relies on zip:

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

This portion scales nicely to arbitrary dimensions. Unfortunately, np.meshgrid doesn't scale well to multiple dimensions, so that part will need to be worked out, or (assuming it works), you could use this SO answer to create your own ndmeshgrid function.

Answer 3

Yet another way to do it is:

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

Which can be generalized to higher dimensions, e.g.:

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]])

Answer 4

A simple example in 3D (can be extended to N-dimensions I guess, but beware of the final dimension and RAM usage):

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)

in: 2.76 s ± 185 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

which yields:

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

Answer 5

To get the coordinates of a grid from 0 to 1, a reshape can do the work. Here are examples for 2D and 3D. Also works with floats.

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

Answer 6

I am using the following to convert meshgrid to M X 2 array. Also changes the list of vectors to iterators can make it really fast.

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)

for N-D array using generator

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)

Answer 7

This is the n-d case.

The idea is that we combine the grids into an array to make something n x dim1 x .. dim_n; we then move this to the end (with moveaxis - a generalization of transpose) and flatten everything apart from the last dim.

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)