python - 如何使用 haversine 距离矩阵返回满足距离阈值标准的所有点对

Question

我有 39,803 个不同的经纬度点(这里是前 5 个点)

lat lon
0   27.269987   -82.497004
1   27.537598   -82.508422
2   27.337793   -82.533753
3   27.497719   -82.570261
4   27.512062   -82.722158

我已经计算了 haversine 距离矩阵

array([[  0.        ,  29.77832527,   8.36846516, ...,  30.07072193,
      7.60700598,   7.63477669],
   [ 29.77832527,   0.        ,  22.35749757, ...,   2.6836159 ,
     22.17639199,  23.07090099],
   [  8.36846516,  22.35749757,   0.        , ...,  23.07825172,
      3.10333262,   0.75441483],
   ..., 
   [ 30.07072193,   2.6836159 ,  23.07825172, ...,   0.        ,
     22.53766911,  23.75965211],
   [  7.60700598,  22.17639199,   3.10333262, ...,  22.53766911,
      0.        ,   3.03795035],
   [  7.63477669,  23.07090099,   0.75441483, ...,  23.75965211,
      3.03795035,   0.        ]])

因此我的矩阵是 39,803 x 39,803。我想找出一种方法来找到所有距离小于 25 米的对。例如，如果我们考虑第一个数组

[  0.        ,  29.77832527,   8.36846516, ...,  30.07072193,
  7.60700598,   7.63477669]

一对

(lat[0],lon[0])-(lat[2],lon[2])=(27.269987,-82.497004)-(27.337793,-82.533753) = 8.36846516

满足这个条件但是

(lat[0],lon[0])-(lat[1],lon[1])=(27.269987,-82.497004)-(27.537598,-82.508422) = 29.77832527

没有。我想获得满足此标准的点的子集。这是我目前所拥有的:

X=df[['lat','lon']].dropna(axis=0) 
coors=np.array(X)

from math import radians, cos, sin, asin, sqrt

from scipy.spatial.distance import pdist, squareform
from sklearn.cluster import DBSCAN

import matplotlib.pyplot as plt
import pandas as pd


def haversine(lonlat1, lonlat2):

    lat1, lon1 = lonlat1
    lat2, lon2 = lonlat2
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r

distance_matrix = squareform(pdist(X, (lambda u,v: haversine(u,v))))

如果能帮助我计算这个，我将不胜感激

Answer 1

如何使用 numpy.argwhere()功能？这可用于查找矩阵中满足特定条件的条目的行/列索引。 在您的情况下，这可能类似于:

numpy.argwhere(distances < 25)

这将返回一个 Nx2 数组，给出 N 个点的行/列索引，这些点之间的距离小于 25m。 显然，您可能需要一些额外的过滤以避免与对称距离矩阵相关的重复，并消除每个点与其自身之间的零距离。也许您可以通过获取距离矩阵并执行以下转换来处理这个问题:

distances = distances + 26 * (numpy.tril(numpy.ones_like(distances), k=1))

这会强制距离矩阵的对角线上和上三角中的所有元素的值都高于 25m 阈值。

在您的特定情况下，因为您在所有点对上隐式使用 for 循环来计算您的距离矩阵，所以您可能不会从使用 numpy 例程来执行过滤的最后阶段。相反，您可能会考虑更像是:

indices = [ (i, j) for i in range(Npoints)
                     for j in range(i+1, Npoints)
                     if haversine(points[i], points[j]) < minDistance ]

这种方法虽然不如纯 numpy 方法有效，但在应用于大量点时应该不会占用大量内存。

混合方法可能涉及获取点列表的子集，在每对子集之间形成距离矩阵，并使用上述 numpy.argwhere() 方法过滤它们。