Python - 如何加快计算城市之间的距离

Question

我的数据库中有 55249 个城市。每个人都有纬度经度值。对于每个城市，我想计算与其他每个城市的距离，并存储不超过 30 公里的距离。这是我的算法:

# distance function
from math import sin, cos, sqrt, atan2, radians

def distance(obj1, obj2):
    lat1 = radians(obj1.latitude)
    lon1 = radians(obj1.longitude)
    lat2 = radians(obj2.latitude)
    lon2 = radians(obj2.longitude)
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    return round(6373.0 * c, 2)

def distances():
    cities = City.objects.all()  # I am using Django ORM
    for city in cities:
        closest = list()
        for tested_city in cities:
            distance = distance(city, tested_city)
            if distance <= 30. and distance != 0.:
                closest.append(tested_city)
        city.closest_cities.add(*closest)  # again, Django thing
        city.save()  # Django

这可行，但需要花费大量时间。需要数周才能完成。有什么办法可以加快速度吗？

Answer 1

您无法计算每对城市之间的距离。相反，您需要将您的城市放在 space-partitioning data structure 中。您可以为此进行快速的最近邻查询。 SciPy带有 kd-tree实现，scipy.spatial.KDTree , 适用于此应用。

这里有两个难点。一、scipy.spatial.KDTree使用点之间的欧氏距离，但您想使用地球表面的大圆距离。其次，经度环绕，因此最近的邻居的经度可能相差 360°。如果您采用以下方法，这两个问题都可以得到解决:

1. 从 geodetic coordinates 转换您的位置(纬度，经度)到ECEF (地心、地固)坐标(x、y、z)。

2. 将这些 ECEF 坐标放入 scipy.spatial.KDTree .

3. 将您的大圆距离(例如 30 公里)转换为欧氏距离。

4. 调用scipy.spatial.KDTree.query_ball_point获取范围内的城市。

下面是一些示例代码来说明这种方法。函数 geodetic2ecef 来自 PySatel by David Parunakian并根据 GPL 获得许可。

from math import radians, cos, sin, sqrt

# Constants defined by the World Geodetic System 1984 (WGS84)
A = 6378.137
B = 6356.7523142
ESQ = 6.69437999014 * 0.001

def geodetic2ecef(lat, lon, alt=0):
    """Convert geodetic coordinates to ECEF."""
    lat, lon = radians(lat), radians(lon)
    xi = sqrt(1 - ESQ * sin(lat))
    x = (A / xi + alt) * cos(lat) * cos(lon)
    y = (A / xi + alt) * cos(lat) * sin(lon)
    z = (A / xi * (1 - ESQ) + alt) * sin(lat)
    return x, y, z

def euclidean_distance(distance):
    """Return the approximate Euclidean distance corresponding to the
    given great circle distance (in km).

    """
    return 2 * A * sin(distance / (2 * B))

让我们组成五万个随机城市位置并将它们转换为 ECEF 坐标:

>>> from random import uniform
>>> cities = [(uniform(-90, 90), uniform(0, 360)) for _ in range(50000)]
>>> ecef_cities = [geodetic2ecef(lat, lon) for lat, lon in cities]

将它们放入 scipy.spatial.KDTree :

>>> import numpy
>>> from scipy.spatial import KDTree
>>> tree = KDTree(numpy.array(ecef_cities))

查找距离伦敦大约 100 公里范围内的所有城市:

>>> london = geodetic2ecef(51, 0)
>>> tree.query_ball_point([london], r=euclidean_distance(100))
array([[37810, 15755, 16276]], dtype=object)

对于您查询的每个点，该数组包含距离 r 内的城市数组。每个邻居都作为其在您传递给 KDTree 的原始数组中的索引给出。所以距离伦敦100公里左右的范围内有3个城市，即原列表中索引为37810、15755、16276的城市:

>>> from pprint import pprint
>>> pprint([cities[i] for i in [37810, 15755, 16276]])
[(51.7186871990946, 359.8043453670437),
 (50.82734317063884, 1.1422052710187103),
 (50.95466110717763, 0.8956257749604779)]

注意事项:

1. 您可以从示例输出中看到，正确发现了经度相差大约 360° 的邻居。

2. 该方法似乎足够快。在这里，我们为前千个城市找到 30 公里范围内的邻居，大约需要 5 秒:

   >>> from timeit import timeit
   >>> timeit(lambda:tree.query_ball_point(ecef_cities[:1000], r=euclidean_distance(30)), number=1)
   5.013611573027447
   

   推断，我们预计在大约四分钟内为所有 50,000 个城市找到 30 公里范围内的邻居。

3. 我的 euclidean_distance 函数高估了对应于给定大圆距离的欧几里得距离(以免遗漏任何城市)。对于某些应用程序来说，这可能已经足够好了——毕竟，城市不是点对象——但如果你需要比这更高的精度，那么你可以使用来自 geopy 的大圆距离函数之一来过滤结果点。 .