Numpy在三维中的网格格点化

Question

Numpy的meshgrid函数非常有用，可以将两个向量转换为坐标网格。如何最容易地将此扩展到三维？因此，给定三个向量x、y和z，构造3个3D数组（而不是2个2D数组），可以用作坐标。

- astrofrog

---

截至2023年，它适用于任意维度，因此np.meshgrid(x, y, z)。 - Guimoute

Answer 1

据我所知，Numpy（版本在1.8及以上）现在支持使用meshgrid生成高于2D的位置网格。其中一个非常有用的新增功能是能够选择索引顺序（如笛卡尔或矩阵索引的xy或ij），我通过以下示例进行了验证：

import numpy as np

x_ = np.linspace(0., 1., 10)
y_ = np.linspace(1., 2., 20)
z_ = np.linspace(3., 4., 30)

x, y, z = np.meshgrid(x_, y_, z_, indexing='ij')

assert np.all(x[:,0,0] == x_)
assert np.all(y[0,:,0] == y_)
assert np.all(z[0,0,:] == z_)

- leifdenby

---

assert语句有何作用？为什么有必要使用它？ - Archie

assert 只是检查沿着 x、y 和 z 轴的索引是否等于预期的 *_ 变量。无需使用。 - Kwame

这应该是被接受的答案，以避免人们重新发明轮子。 - divenex

Answer 2

这里是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

这很简单易懂。我将该模式扩展到任意维度，但是这段代码并没有经过优化（也没有彻底检查错误），但你得到了你所付出的代价。希望对你有所帮助：

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)

- Paul

---

在3D网格的情况下，使用类似于numpy文档中提供的meshgrid示例，返回的顺序将是Z，Y，X而不是X，Y，Z。通过将返回语句替换为“return tuple(ans [::-1])”可以解决这个问题。 - levesque

@Paul 如果x或y数组的长度很长，则x.repeat()命令会崩溃并发送内存错误。有没有办法避免这个错误？ - Dalek

@ Dalek 数组有多长？你可能会用完内存吗？例如，如果有3个数组，每个数组有4096个条目，每个条目都包含一个double（即8字节），那么仅对于条目，我们需要（8 * 4 * 2 ** 10）** 3字节= 2 ** 45字节= 32 * 2 ** 40字节= 32 TB的内存，这显然是巨大的。我希望我没有犯错。 - aignas

Answer 3

你能展示一下你如何使用np.meshgrid吗？非常有可能你并不需要meshgrid，因为numpy广播可以在不生成重复数组的情况下完成相同的操作。

例如：

import numpy as np

x=np.arange(2)
y=np.arange(3)
[X,Y] = np.meshgrid(x,y)
S=X+Y

print(S.shape)
# (3, 2)
# Note that meshgrid associates y with the 0-axis, and x with the 1-axis.

print(S)
# [[0 1]
#  [1 2]
#  [2 3]]

s=np.empty((3,2))
print(s.shape)
# (3, 2)

# x.shape is (2,).
# y.shape is (3,).
# x's shape is broadcasted to (3,2)
# y varies along the 0-axis, so to get its shape broadcasted, we first upgrade it to
# have shape (3,1), using np.newaxis. Arrays of shape (3,1) can be broadcasted to
# arrays of shape (3,2).
s=x+y[:,np.newaxis]
print(s)
# [[0 1]
#  [1 2]
#  [2 3]]

重点是S=X+Y可以替换为s=x+y[:,np.newaxis]，因为后者不需要形成（可能很大的）重复数组。它还可以轻松地推广到更高维（更多轴）。只需在必要时添加np.newaxis以实现广播即可。

有关numpy广播的更多信息，请参见http://www.scipy.org/EricsBroadcastingDoc。

- unutbu

Answer 4

我认为你想要的是：

X, Y, Z = numpy.mgrid[-10:10:100j, -10:10:100j, -10:10:100j]

例如。

- Autoplectic

---

谢谢，但这不是我需要的 - meshgrid实际上使用向量的值来生成2D数组，并且这些值可能是不规则间隔的。 - astrofrog

Answer 5

不需要编写新的函数，numpy.ix_可以满足您的需求。

以下是文档中的示例：

>>> ixgrid = np.ix_([0,1], [2,4])
>>> ixgrid
(array([[0],
   [1]]), array([[2, 4]]))
>>> ixgrid[0].shape, ixgrid[1].shape
((2, 1), (1, 2))'

- j13r

Answer 6

这是我编写的meshgrid的多维版本:

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

请注意，返回的数组是原始数组数据的视图，因此更改原始数组将影响坐标数组。

- user545424

Answer 7

你可以通过改变顺序来实现：

import numpy as np
xx = np.array([1,2,3,4])
yy = np.array([5,6,7])
zz = np.array([9,10])
y, z, x = np.meshgrid(yy, zz, xx)

- Kemono