【发布时间】:2016-02-05 15:20:19
【问题描述】:
我想检查我是否正确使用了 scipy 的 KD 树,因为它看起来比简单的暴力破解要慢。
我对此有三个问题:
第一季度。
如果我创建以下测试数据:
nplen = 1000000
# WGS84 lat/long
point = [51.349,-0.19]
# This contains WGS84 lat/long
points = np.ndarray.tolist(np.column_stack(
[np.round(np.random.randn(nplen)+51,5),
np.round(np.random.randn(nplen),5)]))
并创建三个函数:
def kd_test(points,point):
""" KD Tree"""
return points[spatial.KDTree(points).query(point)[1]]
def ckd_test(points,point):
""" C implementation of KD Tree"""
return points[spatial.cKDTree(points).query(point)[1]]
def closest_math(points,point):
""" Simple angle"""
return (min((hypot(x2-point[1],y2-point[0]),y2,x2) for y2,x2 in points))[1:3]
我希望 cKD 树是最快的,但是 - 运行这个:
print("Co-ordinate: ", f(points,point))
print("Index: ", points.index(list(f(points,point))))
%timeit f(points,point)
结果时间 - 简单的暴力破解方法更快:
closest_math: 1 loops, best of 3: 3.59 s per loop
ckd_test: 1 loops, best of 3: 13.5 s per loop
kd_test: 1 loops, best of 3: 30.9 s per loop
这是因为我用错了 - 不知何故?
第二季度。
我假设即使要获得最近点的排名(而不是距离),仍然需要投影数据。但是,投影点和未投影点似乎给了我相同的最近邻居:
def proj_list(points,
inproj = Proj(init='epsg:4326'),
outproj = Proj(init='epsg:27700')):
""" Projected geo coordinates"""
return [list(transform(inproj,outproj,x,y)) for y,x in points]
proj_points = proj_list(points)
proj_point = proj_list([point])[0]
这仅仅是因为我的点分布不够大而不会引入失真吗?我重新运行了几次,仍然从返回的投影和未投影列表中得到相同的索引。
第三季度。
与在(未投影的)纬度/经度上计算半正弦或文森特距离相比,投影点(如上)并计算斜边距离通常更快吗?还有哪个选项更准确?我做了一个小测试:
from math import *
def haversine(origin,
destination):
"""
Find distance between a pair of lat/lng coordinates
"""
lat1, lon1, lat2, lon2 = map(radians, [origin[0],origin[1],destination[0],destination[1]])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
c = 2 * asin(sqrt(a))
r = 6371000 # Metres
return (c * r)
def closest_math_unproj(points,point):
""" Haversine on unprojected """
return (min((haversine(point,pt),pt[0],pt[1]) for pt in points))
def closest_math_proj(points,point):
""" Simple angle since projected"""
return (min((hypot(x2-point[1],y2-point[0]),y2,x2) for y2,x2 in points))
结果:
所以这似乎是说先投影然后做距离比不做更快 - 但是,我不确定哪种方法会带来更准确的结果。
在online vincenty calculation 上进行测试似乎是预测坐标是可行的方法:
【问题讨论】:
-
一个几乎不相关的建议:使用
%timeit -n 10 f(points,point)可能比使用%timeit for x in range(10): f(points,point)更方便。 -
顺便说一下,github.com/storpipfugl/pykdtree 可能值得一看。与蛮力方法相比,这可能无法解决效率问题,但可能会比 scipy 的默认实现快一点。
标签: python scipy spatial kdtree map-projections