识别大数组中由最大距离分隔的 Python 数组单元对？

Question

我有包含已转换为二维 numpy 数组的空间生态栖息地数据的栅格。在此数组中，值 1 = 数据，0 = 无数据。 从这些数据中，我想生成一个包含所有数据单元格对的数组，其中每个单元格之间的距离小于最大欧几里得截止距离（即相隔 2 个单元格）。

我发现this answer 很有用，但那里的答案似乎首先测量所有成对距离，然后通过最大截止值对结果进行阈值处理。我的数据集很大（13500*12000 数组中有超过 100 万个数据单元），因此任何试图计算 所有 对单元格之间的距离的成对距离测量都将失败：我需要一个以某种方式停止的解决方案在某个搜索半径（或类似的东西）之外寻找可能的邻居。

我已经尝试过scipy.spatial.distance.pdist，但到目前为止还没有运气将它应用于我的二维数据，或者找到一种方法来阻止pdist 计算即使是相距较远的单元格对之间的距离。我附上了一个示例数组和一个最大欧几里德截止距离 = 2 个单元的所需输出数组：

import numpy as np
import matplotlib.pyplot as plt

# Example 2-D habitat array (1 = data)
example_array = np.array([[0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
                          [0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1],
                          [0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1],
                          [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
                          [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
                          [1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1],
                          [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0],
                          [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                          [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                          [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Plot example array
plt.imshow(example_array, cmap="spectral", interpolation='nearest')

Answer 1

我不得不承认我的 numpy 很弱——也许有办法直接做到这一点。尽管如此，这个问题在纯 Python 中并不难。以下代码将输出匹配数据的 x/y 坐标对。有很多潜在的优化可能会掩盖代码并使其运行得更快，但考虑到数据集的大小和示例半径的大小（2.0），我怀疑其中任何一个都是值得的（除了可能的例外在数组而不是子列表中创建 numpy 视图）。

更新 -- 代码已经修复了几个错误 -- (1) 它在起始点下方的行上向左看太远了，并且 (2) 它是没有在左边缘附近做正确的事情。该函数的调用现在使用 2.5 的半径来显示如何拾取额外的对。

example_array = [[0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
                [0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1],
                [0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1],
                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
                [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
                [1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1],
                [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1],
                [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
                [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0],
                [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

def findpairs(mylist, radius = 2.0):
    """
    Find pairs with data within a given radius.
    If we work from the top of the array down, we never
    need to look up (because we already would have found
    those, and we never need to look left on the same line.
    """

    # Create the parameters of a half circle, which is
    # the relative beginning and ending X coordinates to
    # search for each Y line starting at this one and
    # working down.  To avoid duplicates and extra work,
    # not only do we not look up, we never look left on
    # the same line as what we are matching, but we do
    # on subsequent lines.

    semicircle = []
    x = 1
    while x:
        y = len(semicircle)
        x = int(max(0, (radius ** 2 - y ** 2)) ** 0.5)
        # Don't look back on same line...
        semicircle.append((-x if y else 1, x + 1))

    # The maximum number of y lines we will search
    # at a time.
    max_y = len(semicircle)

    for y_start in range(len(mylist)):
        sublists = enumerate(mylist[y_start:y_start + max_y], y_start)
        sublists = zip(semicircle, sublists)
        check = (x for (x, value) in enumerate(mylist[y_start]) if value)
        for x_start in check:
            for (x_lo, x_hi), (y, ylist) in sublists:
                # Deal with left edge problem
                x_lo = max(0, x_lo + x_start)
                xlist = ylist[x_lo: x_start + x_hi]
                for x, value in enumerate(xlist, x_lo):
                    if value:
                        yield (x_start, y_start), (x, y)

print(list(findpairs(example_array, 2.5)))

执行时间将高度依赖数据。对于咧嘴笑，我创建了您指定的大小 (13500 x 12000) 的数组来测试时间。我使用了更大的半径（3.0 而不是 2.0）并尝试了两种情况：没有匹配，以及每次匹配。为了避免一遍又一遍地重新分配列表，我只是运行迭代器并扔掉结果。执行此操作的代码如下。对于最佳情况（空）阵列，它在我的机器上运行了 7 秒；最坏情况（全为 1）阵列的时间约为 12 分钟。

def dummy(val):
    onelist = 13500 * [val]
    listolists = 12000 * [onelist]

    for i in findpairs(listolists, 3.0):
      pass

dummy(0)
dummy(1)