在一个大数组中识别以最大距离分隔的Python数组单元对？

Question

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

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

我曾经尝试过scipy.spatial.distance.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

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

更新了代码已经修复了几个错误（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)))

执行时间将高度依赖于数据。对于grins，我创建了您指定的大小（13500 x 12000）的数组来测试计时。我使用了一个更大的半径（3.0而不是2.0），并尝试了两种情况：没有匹配，每个匹配。为了避免一次又一次地重新分配列表，我只需运行迭代器并抛出结果。代码如下。对于一个best case（空）数组，它在我的机器上运行了7秒；对于最坏情况（all 1s）数组，时间大约是12分钟。在

^{pr2}$