Python: 加速地理比较

Question

我写了一段代码，其中有一个嵌套循环，里面的内层循环大约要执行150万次。我在这个循环中有一个函数想要进行优化。我已经做了一些工作，得到了些结果，但我需要一点建议来确认我做的是否合理。

一些背景信息：

我有两组地理坐标（纬度和经度），一组相对较小，另一组则相对较大。对于小组中的每一个点，我需要找到大组中离它最近的点。

显而易见的方法是使用哈弗辛公式。这样做的好处是，计算出的距离是非常准确的。

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

不过，在我的机器上运行这个公式150万次大约需要9秒（根据timeit的结果）。因为准确的距离并不是很重要，我只需要找到最近的点，所以我决定尝试一些其他的函数。

使用简单的勾股定理实现让我速度提高了大约30%。我觉得可以做得更好，于是我写了以下代码：

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

这让我提高了10倍的速度。但是，现在我有点担心这样做会不符合三角不等式。

所以，我最后的问题有两个方面：我希望有一个运行速度和dumb一样快的函数，但仍然要正确。dumb能行吗？如果不行，有没有建议可以改进我的哈弗辛函数？

Answer 1

abs(lat2 - lat1) + abs(lon2 - lon1) 这个公式表示的是一种叫做“曼哈顿距离”的计算方法。简单来说，它就是计算两个点之间的水平和垂直距离之和。这个方法遵循一个叫做“三角不等式”的规则，也就是说，直接从一个点到另一个点的距离，永远不会比先到一个中间点再到目标点的距离要短。

Answer 2

这种计算是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函数组合成一个语句来进一步提高速度，但那样会变得很长且难以阅读。

Answer 3

你可以考虑一种图形哈希的方法，也就是快速找到最近的点，然后再对这些点进行计算。比如，你可以创建一个均匀的网格，然后把大量的点分配到这个网格划分出来的小区域里。

现在，当你有一个来自小集合的点时，你只需要处理更少的点（也就是那些在相关小区域里的点）