迭代一个numpy数组，然后索引另一个数组中的值

Question

我正在努力使这段代码正常工作我想迭代一个numpy数组，并根据结果，索引到另一个numpy数组中的一个值，然后基于该值将其保存到一个新的位置。在

     # Convert the sediment transport and the flow direction rasters into Numpy arrays
    sediment_transport_np = arcpy.RasterToNumPyArray(sediment_transport_convert, '#', '#', '#', -9999)
    flow_direction_np = arcpy.RasterToNumPyArray(flow_direction_convert, '#', '#', '#', -9999)

    [rows,cols]= sediment_transport_np.shape
    elevation_change = np.zeros((rows,cols), np.float)

   # Main body for calculating elevation change        

    # Attempt 1
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        if flow == 32:
            elevation_change[i, j]  = sediment_transport_np[i - 1, j - 1]
        elif flow == 16:
            elevation_change[i, j] = sediment_transport_np[i, j - 1]
        elif flow == 8:
            elevation_change[i, j] = sediment_transport_np[i + 1, j - 1]
        elif flow == 4:
            elevation_change[i, j] = sediment_transport_np[i + 1, j]
        elif flow == 64:
            elevation_change[i, j] = sediment_transport_np[i - 1, j]
        elif flow == 128:
            elevation_change[i, j] = sediment_transport_np[i - 1, j + 1]
        elif flow == 1:
            elevation_change[i, j] = sediment_transport_np[i, j + 1]
        elif flow == 2:
            elevation_change[i, j] = sediment_transport_np[i + 1, j + 1]

将Numpy数组转换回raster

^{pr2}$

我得到的错误是：

正在运行脚本提升更改。。。在

回溯（最近一次调用）：文件“”，第606行，在执行索引器错误：索引（655）超出范围（0<；=索引<；655）在维度0中

执行失败（提升更改）

Answer 1

问题的原因

错误是因为您试图索引到sediment_transport网格的边界之外（例如i+1和j+1部分）。现在，你试图得到一个当你在网格边界时不存在的值。而且，它不会引发错误，但是当你在i=0或j=0时，你正在抓取相反的边（由于i-1和j-1部分）。在

您提到您希望elevation_change的值在边界处为0（这显然是合理的）。另一个常见的边界条件是“包装”这些值并从另一个边获取一个值。在本例中，这可能没有多大意义，但我将在几个示例中展示它，因为使用某些方法很容易实现。在

诱人但不正确

很容易捕捉到异常并将值设置为0。例如：

for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 32:
            ...
        elif ...
            ...
    except IndexError:
        elevation_change[i, j] = 0

然而，这种方法实际上是不正确的。负索引有效，将返回网格的相对边。因此，这将基本上在网格的右边缘和下边缘实现“零”边界条件，在左侧和顶部边缘实现“环绕”边界条件。在

用零填充

在“零”边界条件的情况下，有一种非常简单的方法来避免索引问题：用零填充sediment_transport网格。这样，如果索引超出原始网格的边缘，我们将得到0。（或任何您想填充数组的常量值。）

旁注：这是使用numpy.pad的完美地方。但是，它是在v1.7中添加的。我将在这里跳过使用它，因为OP提到了ArcGIS，Arc并没有提供最新版本的numpy。

例如：

^{2}$

干（不要重复）

我们可以使用dict来编写这段代码：

elevation_change = np.zeros_like(sediment_transport)
nrows, ncols = flow_direction.shape
lookup = {32: (-1, -1),
          16:  (0, -1), 
          8:   (1, -1),
          4:   (1,  0),
          64: (-1,  0),
          128:(-1,  1),
          1:   (0,  1),
          2:   (1,  1)}

padded_transport = np.zeros((nrows + 2, ncols + 2), float)
padded_transport[1:-1, 1:-1] = sediment_transport  

for [i, j], flow in np.ndenumerate(flow_direction):
    # Need to take into account the offset in the "padded_transport"
    r, c = i + 1, j + 1
    # This also allows for flow_direction values not listed above...
    dr, dc = lookup.get(flow, (0,0))
    elevation_change[i,j] = padded_transport[r + dr, c + dc]

此时，填充原始数组有点多余。如果使用numpy.pad，通过填充实现不同的边界条件非常容易，但我们可以直接写出逻辑：

elevation_change = np.zeros_like(sediment_transport)
nrows, ncols = flow_direction.shape
lookup = {32: (-1, -1),
          16:  (0, -1), 
          8:   (1, -1),
          4:   (1,  0),
          64: (-1,  0),
          128:(-1,  1),
          1:   (0,  1),
          2:   (1,  1)}

for [i, j], flow in np.ndenumerate(flow_direction):
    dr, dc = lookup.get(flow, (0,0))
    r, c = i + dr, j + dc
    if not ((0 <= r < nrows) & (0 <= c < ncols)):
        elevation_change[i,j] = 0
    else:
        elevation_change[i,j] = sediment_transport[r, c]

“矢量化”计算

在python中迭代numpy数组是相当缓慢的，原因是我将不在这里深入研究。因此，在numpy中有更有效的方法来实现这一点。诀窍是将numpy.roll与布尔索引一起使用。在

对于“环绕”边界条件，它很简单：

elevation_change = np.zeros_like(sediment_transport)
nrows, ncols = flow_direction.shape
lookup = {32: (-1, -1),
          16:  (0, -1), 
          8:   (1, -1),
          4:   (1,  0),
          64: (-1,  0),
          128:(-1,  1),
          1:   (0,  1),
          2:   (1,  1)}

for value, (row, col) in lookup.iteritems():
    mask = flow_direction == value
    shifted = np.roll(mask, row, 0)
    shifted = np.roll(shifted, col, 1)
    elevation_change[mask] = sediment_transport[shifted]

return elevation_change

如果你不熟悉numpy，这可能看起来有点像希腊语。这有两个部分。第一种是使用布尔索引。作为一个简单的例子：

In [1]: import numpy as np

In [2]: x = np.arange(9).reshape(3,3)

In [3]: x
Out[3]: 
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

In [4]: mask = np.array([[False, False, True],
...                      [True, False, False],
...                      [True, False, False]])


In [5]: x[mask]
Out[5]: array([2, 3, 6])

如您所见，如果我们用相同形状的布尔网格索引数组，则返回其为真的值。类似地，您可以这样设置值。在

下一个技巧是numpy.roll。这将使值在一个方向上移动给定的量。它们会在边缘“缠绕”。在

In [1]: import numpy as np

In [2]: np.array([[0,0,0],[0,1,0],[0,0,0]])
Out[2]: 
array([[0, 0, 0],
       [0, 1, 0],
       [0, 0, 0]])

In [3]: x = _

In [4]: np.roll(x, 1, axis=0)
Out[4]: 
array([[0, 0, 0],
       [0, 0, 0],
       [0, 1, 0]])

In [5]: np.roll(x, 1, axis=1)
Out[5]: 
array([[0, 0, 0],
       [0, 0, 1],
       [0, 0, 0]])

希望这有点道理，无论如何。在

要实现“零”边界条件（或使用numpy.pad的任意边界条件），我们可以这样做：

def vectorized(flow_direction, sediment_transport):
    elevation_change = np.zeros_like(sediment_transport)
    nrows, ncols = flow_direction.shape
    lookup = {32: (-1, -1),
              16:  (0, -1), 
              8:   (1, -1),
              4:   (1,  0),
              64: (-1,  0),
              128:(-1,  1),
              1:   (0,  1),
              2:   (1,  1)}

    # Initialize an array for the "shifted" mask
    shifted = np.zeros((nrows+2, ncols+2), dtype=bool)

    # Pad "sediment_transport" with zeros
    # Again, `np.pad` would be better and more flexible here, as it would
    # easily allow lots of different boundary conditions...
    tmp = np.zeros((nrows+2, ncols+2), sediment_transport.dtype)
    tmp[1:-1, 1:-1] = sediment_transport
    sediment_transport = tmp

    for value, (row, col) in lookup.iteritems():
        mask = flow_direction == value

        # Reset the "shifted" mask
        shifted.fill(False)
        shifted[1:-1, 1:-1] = mask

        # Shift the mask by the right amount for the given value
        shifted = np.roll(shifted, row, 0)
        shifted = np.roll(shifted, col, 1)

        # Set the values in elevation change to the offset value in sed_trans
        elevation_change[mask] = sediment_transport[shifted]

    return elevation_change

矢量化的优势

“矢量化”版本速度更快，但将使用更多的RAM。在

对于1000 x 1000网格：

In [79]: %timeit vectorized(flow_direction, sediment_transport)
10 loops, best of 3: 170 ms per loop

In [80]: %timeit iterate(flow_direction, sediment_transport)
1 loops, best of 3: 5.36 s per loop

In [81]: %timeit lookup(flow_direction, sediment_transport)
1 loops, best of 3: 3.4 s per loop

这些结果是通过将以下实现与随机生成的数据进行比较得出的：

import numpy as np

def main():
    # Generate some random flow_direction and sediment_transport data...
    nrows, ncols = 1000, 1000
    flow_direction = 2 ** np.random.randint(0, 8, (nrows, ncols))
    sediment_transport = np.random.random((nrows, ncols))

    # Make sure all of the results return the same thing...
    test1 = vectorized(flow_direction, sediment_transport)
    test2 = iterate(flow_direction, sediment_transport)
    test3 = lookup(flow_direction, sediment_transport)
    assert np.allclose(test1, test2)
    assert np.allclose(test2, test3)


def vectorized(flow_direction, sediment_transport):
    elevation_change = np.zeros_like(sediment_transport)
    sediment_transport = np.pad(sediment_transport, 1, mode='constant')
    lookup = {32: (-1, -1),
              16:  (0, -1), 
              8:   (1, -1),
              4:   (1,  0),
              64: (-1,  0),
              128:(-1,  1),
              1:   (0,  1),
              2:   (1,  1)}

    for value, (row, col) in lookup.iteritems():
        mask = flow_direction == value
        shifted = np.pad(mask, 1, mode='constant')
        shifted = np.roll(shifted, row, 0)
        shifted = np.roll(shifted, col, 1)
        elevation_change[mask] = sediment_transport[shifted]

    return elevation_change

def iterate(flow_direction, sediment_transport):
    elevation_change = np.zeros_like(sediment_transport)
    padded_transport = np.pad(sediment_transport, 1, mode='constant')  

    for [i, j], flow in np.ndenumerate(flow_direction):
        r, c = i + 1, j + 1
        if flow == 32:
            elevation_change[i, j] = padded_transport[r - 1, c - 1]
        elif flow == 16:
            elevation_change[i, j] = padded_transport[r, c - 1]
        elif flow == 8:
            elevation_change[i, j] = padded_transport[r + 1, c - 1]
        elif flow == 4:
            elevation_change[i, j] = padded_transport[r + 1, c]
        elif flow == 64:
            elevation_change[i, j] = padded_transport[r - 1, c]
        elif flow == 128:
            elevation_change[i, j] = padded_transport[r - 1, c + 1]
        elif flow == 1:
            elevation_change[i, j] = padded_transport[r, c + 1]
        elif flow == 2:
            elevation_change[i, j] = padded_transport[r + 1, c + 1]
    return elevation_change

def lookup(flow_direction, sediment_transport):
    elevation_change = np.zeros_like(sediment_transport)
    nrows, ncols = flow_direction.shape
    lookup = {32: (-1, -1),
              16:  (0, -1), 
              8:   (1, -1),
              4:   (1,  0),
              64: (-1,  0),
              128:(-1,  1),
              1:   (0,  1),
              2:   (1,  1)}

    for [i, j], flow in np.ndenumerate(flow_direction):
        dr, dc = lookup.get(flow, (0,0))
        r, c = i + dr, j + dc
        if not ((0 <= r < nrows) & (0 <= c < ncols)):
            elevation_change[i,j] = 0
        else:
            elevation_change[i,j] = sediment_transport[r, c]

    return elevation_change

if __name__ == '__main__':
    main()