Python：加速地理比较

Question

Python：加速地理比较

pythonoptimizationdistancegeography

13

我写了一些代码，其中包含一个嵌套循环，内部循环执行大约150万次。在这个循环中，我有一个函数需要进行优化。我已经做了一些工作并得到了一些结果，但我需要一些反馈来检查我的做法是否合理。

一些背景：

我有两个地理点（纬度，经度）的集合，一个相对较小的集合和一个相对非常大的集合。对于小集合中的每个点，我需要找到大集合中最接近的点。

显然的方法是使用Haversine公式。这里的好处是距离肯定是准确的。

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

def haversine(point1, point2):
    """Gives the distance between two points on earth.
    """
    earth_radius_miles = 3956
    lat1, lon1 = (radians(coord) for coord in point1)
    lat2, lon2 = (radians(coord) for coord in point2)
    dlat, dlon = (lat2 - lat1, lon2 - lon1)
    a = sin(dlat/2.0)**2 + cos(lat1) * cos(lat2) * sin(dlon/2.0)**2
    great_circle_distance = 2 * asin(min(1,sqrt(a)))
    d = earth_radius_miles * great_circle_distance
    return d

然而，在我的机器上运行1.5百万次需要大约9秒钟（根据timeit）。由于准确距离并不重要，我只需要找到最接近的点，所以我决定尝试一些其他函数。

对勾股定理的简单实现让我加速了约30％。认为我可以做得更好，我写了以下内容：

def dumb(point1, point2):
    lat1, lon1 = point1
    lat2, lon2 = point2
    d = abs((lat2 - lat1) + (lon2 - lon1))

我使用一种新的haversine函数，可以让我的代码运行速度提升10倍。然而，现在我担心这个函数不能保持三角不等式。

所以，我的最终问题有两个方面：我想要一个像那样快速运行的函数，但仍然正确。可以工作吗？如果不行，你有任何改进haversine函数的建议吗？

- Wilduck

6个回答

7

你可以考虑使用一种图形哈希技术，即快速查找最接近的点，然后在这些点上进行计算。例如，你可以创建一个均匀的网格，并将大集合中的点分配到由网格创建的桶中。

现在，如果有一个来自小集合的点，你只需要处理少量的点（即仅限于相关桶中的点）。

- Drakosha

虽然这完全回避了我所问的问题，但我必须接受它，因为这就是我的最终做法。感谢您提供这个视角。 - Wilduck

4

算法优化始终是最佳答案-网格建议非常好，是四叉树（或八叉树）空间划分方案的一个特例，而且相对容易实现。 - jheriko

2

我曾经遇到过类似的问题，于是决定编写一个Cython函数。在我的2008年款MacBook Pro上，它可以每秒处理约120万次迭代。去除类型检查可以进一步提高25%的速度。毫无疑问，还可以进行进一步的优化（但会降低代码清晰度）。

你可能还想看看scipy.spatial.distance.cdist函数。

from libc.math cimport sin, cos, acos

def distance(float lat1, float lng1, float lat2, float lng2):
    if lat1 is None or lat2 is None or lng1 is None or lng2 is None: return None
    cdef float phi1
    cdef float phi2
    cdef float theta1
    cdef float theta2
    cdef float c
    cdef float arc

    phi1 = (90.0 - lat1)*0.0174532925
    phi2 = (90.0 - lat2)*0.0174532925
    theta1 = lng1*0.0174532925
    theta2 = lng2*0.0174532925

    c = (sin(phi1)*sin(phi2)*cos(theta1 - theta2) + cos(phi1)*cos(phi2))
    arc = acos( c )
    return arc*6371

- Will

2

您所写的公式（d=abs(lat2-lat1)+(lon2-lon1)）并不能保持三角不等式：如果您找到lat、lon，使得d最小，您并没有找到最近的点，而是找到了两条对角线相交的直线最接近的点！

我认为您应该按照纬度和经度对大量的点进行排序（这意味着：（1,1），（1,2），（1,3）...（2,1），（2,2）等等。然后使用Gunner方法来查找在纬度和经度方面一些最接近的点（这应该非常快，它将花费CPU时间与ln2（n）成比例，其中n是点数）。您可以很容易地做到这一点，例如：选择围绕要检查的点的10x10正方形中的所有点，这意味着：找到所有在lat从-10到+10之间的点（Gunner方法），再次找到在lon从-10到+10之间的点（Gunner方法）。现在您有一个非常少的数据需要处理，应该非常快！

- Matteo Monti

1

这不是原帖作者写的公式，我认为在他的版本下三角不等式实际上是成立的。 - user3850

2

abs(lat2 - lat1) + abs(lon2 - lon1)是1-范数或曼哈顿度量，因此三角不等式成立。

- user3850

这是很好的事情。实际上我没有在我的代码中写这个，但显然这就是我想要的。感谢您的帮助。 - Wilduck

1

@Wilduck：实际上，在我回答的第一个草稿中，我还写到了你的公式并不完全是1-范数的再创造，但我为了讽刺而将其删减了...无论如何，我认为你的公式实际上是一种伪度量：三角不等式和对称性成立，但两个不同的点可以有0的距离。然而，我懒得去检查。 - user3850

1

最快的方法是避免为每对点计算一个函数，假设您相对较小的集合不是非常微小。

有一些数据库具有您可以使用的地理索引（mysql、oracle、mongodb..），或者自己实现。

您可以使用python-geohash。对于较小集合中的每个文档，您需要快速找到与具有匹配项的最长哈希大小的geohash.neighbors共享哈希的较大集合中的文档集。您将需要使用适当的数据结构进行查找，否则这将很慢。

要查找点之间的距离，简单方法的误差随着点之间的距离增加而增加，并且还取决于纬度。例如，请参见http://www.movable-type.co.uk/scripts/gis-faq-5.1.html。

- Shane Evans

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jterrace · Accepted Answer

这是numpy擅长的计算类型。不必循环整个坐标集，只需进行一次计算即可计算单个点与整个数据集之间的距离。通过下面的测试，您可以获得一个数量级的速度提升。

这里有一些使用您的haversine方法、dumb方法（不太清楚其功能）和我的numpy haversine方法进行计时测试的结果。它计算了两个点之间的距离——一个在弗吉尼亚州，另一个在加利福尼亚州，相距2293英里。

from math import radians, sin, cos, asin, sqrt, pi, atan2
import numpy as np
import itertools

earth_radius_miles = 3956.0

def haversine(point1, point2):
    """Gives the distance between two points on earth.
    """
    lat1, lon1 = (radians(coord) for coord in point1)
    lat2, lon2 = (radians(coord) for coord in point2)
    dlat, dlon = (lat2 - lat1, lon2 - lon1)
    a = sin(dlat/2.0)**2 + cos(lat1) * cos(lat2) * sin(dlon/2.0)**2
    great_circle_distance = 2 * asin(min(1,sqrt(a)))
    d = earth_radius_miles * great_circle_distance
    return d

def dumb(point1, point2):
    lat1, lon1 = point1
    lat2, lon2 = point2
    d = abs((lat2 - lat1) + (lon2 - lon1))
    return d
    
def get_shortest_in(needle, haystack):
    """needle is a single (lat,long) tuple.
        haystack is a numpy array to find the point in
        that has the shortest distance to needle
    """
    dlat = np.radians(haystack[:,0]) - radians(needle[0])
    dlon = np.radians(haystack[:,1]) - radians(needle[1])
    a = np.square(np.sin(dlat/2.0)) + cos(radians(needle[0])) * np.cos(np.radians(haystack[:,0])) * np.square(np.sin(dlon/2.0))
    great_circle_distance = 2 * np.arcsin(np.minimum(np.sqrt(a), np.repeat(1, len(a))))
    d = earth_radius_miles * great_circle_distance
    return np.min(d)
    
    
x = (37.160316546736745, -78.75)
y = (39.095962936305476, -121.2890625)

def dohaversine():
    for i in xrange(100000):
        haversine(x,y)
        
def dodumb():
    for i in xrange(100000):
        dumb(x,y)
        
lots = np.array(list(itertools.repeat(y, 100000)))
def donumpy():
    get_shortest_in(x, lots)

from timeit import Timer
print 'haversine distance =', haversine(x,y), 'time =',
print Timer("dohaversine()", "from __main__ import dohaversine").timeit(100)
print 'dumb distance =', dumb(x,y), 'time =',
print Timer("dodumb()", "from __main__ import dodumb").timeit(100)
print 'numpy distance =', get_shortest_in(x, lots), 'time =',
print Timer("donumpy()", "from __main__ import donumpy").timeit(100)

这是它打印的内容：

haversine distance = 2293.13242188 time = 44.2363960743
dumb distance = 40.6034161104 time = 5.58199882507
numpy distance = 2293.13242188 time = 1.54996609688

使用numpy方法计算与使用您的函数方法计算相同数量的距离计算所需的时间不到 1.55 秒，而使用您的函数方法需要 44.24 秒。您可能可以通过将一些 numpy 函数组合成单个语句来获得更多的加速，但这会变成一行很长、难以阅读的代码。