python - 识别大数组中按最大距离分隔的成对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)，我怀疑其中任何一个都是值得的(可能的异常(exception)是在数组中而不是子列表中创建 numpy View )。

已更新 - 代码修复了一些错误 - (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)