识别大型阵列中由最大距离分隔的 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)))

执行时间将高度依赖于数据。对于 grins，我创建了您指定大小 (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)

识别大型阵列中由最大距离分隔的 Python 阵列单元对？

Identifying pairs of Python array cells separated by maximum distance in a large array?

python

arrays

numpy

scipy

distance