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

你无法计算每对城市之间的距离.相反,您需要将城市放在空间分区数据结构中,以便进行快速最近邻查询. SciPy带有kd -tree实现,scipy.spatial.KDTree适用于此应用程序.

这里有两个困难.首先,scipy.spatial.KDTree使用点之间的欧几里德距离,但是你想要沿着地球表面使用大圆距离.其次,经度环绕,因此最近的邻居可能有相差360°的经度.如果采用以下方法,这两个问题都可以解决:

1. 将您的位置从大地坐标(纬度,经度)转换为ECEF(以地球为中心,地球固定)坐标(x,y,z).

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

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

4. 打电话scipy.spatial.KDTree.query_ball_point让城市在范围内.

下面是一些示例代码来说明这种方法.该功能geodetic2ecef来自David Parunakian的PySatel,并根据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公里范围内有三个城市,即原始列表中索引为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. 这种方法似乎足够快.在这里,我们发现前1000个城市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中的一个大圆距离函数来过滤结果点.