Numpy Meshgrid in 3D

Numpy meshgrid in 3D

Here is the source code of meshgrid:

def meshgrid(x,y):
    """
    Return coordinate matrices from two coordinate vectors.

    Parameters
    ----------
    x, y : ndarray
        Two 1-D arrays representing the x and y coordinates of a grid.

    Returns
    -------
    X, Y : ndarray
        For vectors `x`, `y` with lengths ``Nx=len(x)`` and ``Ny=len(y)``,
        return `X`, `Y` where `X` and `Y` are ``(Ny, Nx)`` shaped arrays
        with the elements of `x` and y repeated to fill the matrix along
        the first dimension for `x`, the second for `y`.

    See Also
    --------
    index_tricks.mgrid : Construct a multi-dimensional "meshgrid"
                         using indexing notation.
    index_tricks.ogrid : Construct an open multi-dimensional "meshgrid"
                         using indexing notation.

    Examples
    --------
    >>> X, Y = np.meshgrid([1,2,3], [4,5,6,7])
    >>> X
    array([[1, 2, 3],
           [1, 2, 3],
           [1, 2, 3],
           [1, 2, 3]])
    >>> Y
    array([[4, 4, 4],
           [5, 5, 5],
           [6, 6, 6],
           [7, 7, 7]])

    `meshgrid` is very useful to evaluate functions on a grid.

    >>> x = np.arange(-5, 5, 0.1)
    >>> y = np.arange(-5, 5, 0.1)
    >>> xx, yy = np.meshgrid(x, y)
    >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2)

    """
    x = asarray(x)
    y = asarray(y)
    numRows, numCols = len(y), len(x)  # yes, reversed
    x = x.reshape(1,numCols)
    X = x.repeat(numRows, axis=0)

    y = y.reshape(numRows,1)
    Y = y.repeat(numCols, axis=1)
    return X, Y

It is fairly simple to understand. I extended the pattern to an arbitrary number of dimensions, but this code is by no means optimized (and not thoroughly error-checked either), but you get what you pay for. Hope it helps:

def meshgrid2(*arrs):
    arrs = tuple(reversed(arrs))  #edit
    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)

Efficient way to map 3D function to a meshgrid with NumPy

For each item of data you're scanning pixels of cuboid to check if it's inside. There is an option to skip this scan. You could calculate corresponding indices of these pixels by yourself, for example:

data = np.array([[1, 2, 3], #14 (corner1)
                 [4, 5, 6], #77 (corner2)
                 [2.5, 3.5, 4.5], #38.75 (duplicated pixel)
                 [2.9, 3.9, 4.9], #47.63 (duplicated pixel)
                 [1.5, 2, 3]]) #15.25 (one step up from [1, 2, 3])
step = 0.5                            
data_idx = ((data - data.min(axis=0))//step).astype(int)
M = np.zeros(np.max(data_idx, axis=0) + 1)
x, y, z = data_idx.T
M[x, y, z] = F

Note that only one value of duplicated pixels is being mapped to M.

Creating a numpy array of 3D coordinates from three 1D arrays

To use numpy mesh grid on the above example the following will work:

np.vstack(np.meshgrid(x_p,y_p,z_p)).reshape(3,-1).T

Numpy meshgrid for grids of more then two dimensions require numpy 1.7. To circumvent this and pulling the relevant data from the source code.

def ndmesh(*xi,**kwargs):
    if len(xi) < 2:
        msg = 'meshgrid() takes 2 or more arguments (%d given)' % int(len(xi) > 0)
        raise ValueError(msg)

    args = np.atleast_1d(*xi)
    ndim = len(args)
    copy_ = kwargs.get('copy', True)

    s0 = (1,) * ndim
    output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)]

    shape = [x.size for x in output]

    # Return the full N-D matrix (not only the 1-D vector)
    if copy_:
        mult_fact = np.ones(shape, dtype=int)
        return [x * mult_fact for x in output]
    else:
        return np.broadcast_arrays(*output)

Checking the result:

print np.vstack((ndmesh(x_p,y_p,z_p))).reshape(3,-1).T

[[ 1.  2.  8.]
 [ 1.  2.  9.]
 [ 1.  3.  8.]
 ....
 [ 5.  3.  9.]
 [ 5.  4.  8.]
 [ 5.  4.  9.]]

For the above example:

%timeit sol2()
10000 loops, best of 3: 56.1 us per loop

%timeit np.vstack((ndmesh(x_p,y_p,z_p))).reshape(3,-1).T
10000 loops, best of 3: 55.1 us per loop

For when each dimension is 100:

%timeit sol2()
1 loops, best of 3: 655 ms per loop
In [10]:

%timeit points = np.vstack((ndmesh(x_p,y_p,z_p))).reshape(3,-1).T
10 loops, best of 3: 21.8 ms per loop

Depending on what you want to do with the data, you can return a view:

%timeit np.vstack((ndmesh(x_p,y_p,z_p,copy=False))).reshape(3,-1).T
100 loops, best of 3: 8.16 ms per loop

Evenly sampled 3D meshgrid

In [117]: x = np.linspace(-1,1,100)
     ...: xx, yy, zz = np.meshgrid(x, x, x)
In [118]: xx.shape
Out[118]: (100, 100, 100)

1000 equally spaced points in xx, similarly for all other grids:

In [119]: xx[::10,::10,::10].shape
Out[119]: (10, 10, 10)

Or with advanced indexing (making a copy)

In [123]: i=np.arange(0,100,10)
In [124]: xx[np.ix_(i,i,i)].shape
Out[124]: (10, 10, 10)

I think we could use np.ravel_multi_index to get an array of flattened indices. We'd have to generate 1000 tuples of indices to do that!

I don't see how we could get a 10,000 points. ::5 would give 8000 points.

numpy 3D meshgrid only as a view

>>> xi, yi, zi = [np.arange(3) for i in range(3)]
>>> xx, yy, zz = np.broadcast_arrays(xi,yi[:,np.newaxis],zi[:,np.newaxis,np.newaxis])
>>> xx.shape
(3, 3, 3)
>>> xx.strides
(0, 0, 8)

You can see it didn't create new copies since the strides are 0 in the first two dimensions.

I wrote a n dimensional version of this also:

def ndmesh(*args):
   args = map(np.asarray,args)
   return np.broadcast_arrays(*[x[(slice(None),)+(None,)*i] for i, x in enumerate(args)])

How to extract a 2D plane from a 3D numpy meshgrid

Remove projection='3d' and index X and Y:

fig = plt.figure()
ax = fig.gca()
zplane = round(N / 2)
ax.quiver(X[:, :, zplane], Y[:, :, zplane], sumEx[:, :, zplane], sumEy[:, :, zplane])
plt.show()

If you select a specific zplane your plot is no longer a 3D-plot.

Plotting a 3D Meshgrid:

Please try to use the search function; you can find working code here among many other threads on plotting surfaces here on Overflow.

---

I think you might have some terminology mix up here. You have a matrix of values which correspond to the values of a function of two dimensions; this is not a meshgrid. In matplotlib you can visualize a three-dimensional matrix in a number of ways, typically as a surface, a wireframe, or an image (as a heatmap).

The function np.meshgrid() gives you N-dimensional indices. You need this for your n and m vectors, not your matrix U, to turn them into multidimensional arrays the same shape as your matrix. Luckily for you, matplotlib doesn't care whether U is a matrix or a vector.

For example with np.meshgrid():

>>> t = np.linspace(0,60,6000)
>>> x = np.linspace(0,2*np.pi,3600)
>>> T, X = np.meshgrid(t, x)
>>> t.shape
(6000,)
>>> x.shape
(3600,)
>>> T.shape
(3600, 6000)
>>> X.shape
(3600, 6000)

And now to create a function U(x,t):

>>> U = np.matrix(x[:, np.newaxis] * t) # broadcast
>>> U.shape
(3600, 6000)

And to plot the surface, you can use the surf function in Axes3D from matplotlib for example, but you can also use the wireframe or image method linked above.

>>> import matplotlib.pyplot as plt
>>> import pylab
>>> from mpl_toolkits.mplot3d import Axes3D

>>> fig = plt.figure()
>>> ax = fig.gca(projection='3d')
>>> surf = ax.plot_surface(X,T,U)
>>> plt.show()

---

Related Topics