Python -如何加快城市间距离的计算
我的数据库里有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这是可行的，但花费了大量的时间。要花几个星期才能完成。我能加快速度吗？
回答 5
Stack Overflow用户
发布于 2013-12-18 14:26:12
你负担不起每一对城市之间的距离。相反，您需要将您的城市放在一个space-partitioning data structure中，您可以对其进行快速的近邻查询。SciPy附带了一个适用于此应用程序的-tree实现scipy.spatial.KDTree。
这里有两个困难。首先，scipy.spatial.KDTree使用点之间的欧几里德距离，但你想使用地球表面的大圆距离。其次，经度环绕，所以近邻的经度可能相差360°。如果采取以下方法，这两个问题都可以解决：
- 将您的位置从geodetic coordinates (纬度、经度)转换为ECEF (以地球为中心，地球固定)坐标(x，y，z)。
- 将这些ECEF坐标放入
scipy.spatial.KDTree中。 - 将你的大圆距离(例如，30公里)转换成欧几里得距离。
- 打电话给
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))让我们将5万个随机城市位置转换为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)]备注：
- 从这个例子的输出可以看出，正确地发现了与经度相差约360°的邻域。
- 这个方法似乎足够快了。在这里，我们在前1000个城市的30公里内找到邻居，花费了大约5秒：从时间导入>>> >>>r=euclidean_distance(30)，number=1) 5.013611573027447 根据推断，我们预计在大约4分钟内，所有5万座城市都能在30公里以内找到邻居。
- 我的
euclidean_distance函数高估了对应于给定大圆距离的欧几里德距离(以免错过任何城市)。对于某些应用程序来说，这可能已经足够好了--毕竟，城市不是点对象--但是如果您需要更高的精度，那么您可以使用(比方说)从geopy的一个大圆距离函数来过滤产生的点。
Stack Overflow用户
发布于 2013-12-18 11:57:34
如果你知道城市之间的距离超过30公里，你可以通过不输入复杂的三角公式来加速你的距离计算，因为它们的纬度差异相当于一个30公里以上的弧度。长度为a= 30公里的弧对应于a/r = 0.00470736的角度，因此：
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
if dlat > 0.00471:
return 32
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)半径32只是一个虚拟值，表明城市之间的距离超过30公里。您应该对经度应用类似的逻辑，对此您必须考虑最大的绝对纬度：
if cos(lat1) * dlon > 0.00471 and cos(lat2) * dlon > 0.00471:
return 32如果你知道你的城市是在一个固定的纬度范围，你可以调整常量限制到最坏的情况。例如，如果你所有的城市都在毗连的美国，它们应该在北纬49°以下，你的极限是0.00471 / cos(49°) = 0.00718。
if dlon > 0.00718:
return 32这个更简单的标准意味着，在得克萨斯州或佛罗里达州的太多城市中，你将进入精确的计算范围。你也可以连锁这些标准。首先使用近似极限，然后是基于最大绝对纬度的精确极限，然后计算所有剩余候选人的精确距离。
你可以用你最大的绝对纬度预先计算这个极限。正如RemcoGerlich所建议的那样，这种启发也可以帮助您将城市放入固定的经度和纬度的桶中。他的方法应该大大加快你的进程，只考虑合理的对城市之前。
看到上面的代码没有检查这个限制的绝对值，我感到有点惭愧。无论如何，这里真正的教训是，无论你如何加快距离计算，大数据集的真正好处来自于选择一种智能搜索机制，比如其他评论者建议的桶搜索或kd树，可能还会加上一些回忆录，以消除双重检查。
Stack Overflow用户
发布于 2013-12-18 11:49:01
我首先要创建“扇区”，每个扇区都有两个纬度相距X公里，两个经度相距X公里。X应该尽可能大，但有一个限制:一个部门内的所有城市都不超过30公里。
这些扇区可以存储在一个数组中：
Sector[][] sectors;在这个数组中，很容易识别包含特定坐标的扇区。确定某一特定部门的相邻部门也很容易。
然后：
(1)每个城市都有自己的分区。每个部门都有一个城市清单。
(2)对于每个城市，查找其所在区域内的所有城市。那些立即达到30公里的标准。
(3)对于每个城市C，查找所有8个相邻区域的所有城市C‘。对于每个C'，检查距离C‘和输出C’，如果它是< 30公里.
这个算法仍然是O(n^2)，但是它应该要快得多，因为对于每个城市，您只检查整个集合的一个小子集。
https://stackoverflow.com/questions/20654918
复制相似问题
腾讯云开发者
