KDTree 正在返回半径之外的点答案

【问题标题】：KDTree is Returning Points Outside of RadiusKDTree 正在返回半径之外的点
【发布时间】：2018-08-22 07:17:32
【问题描述】：

我有一个经纬度坐标数组，我正在尝试使用 KDTree 和 scipy 的 query_ball_point 返回指定纬度和经度 1 英里半径内的所有数据点。

问题在于query_ball_point 返回的点超出了指定的 1 英里半径。这是我的代码：

import pandas as pd
import scipy as sp
import geocoder
import pysal as psl


search_list = df['coordinates'].tolist()
tree = psl.cg.KDTree(search_list, distance_metric='Arc', radius=psl.cg.RADIUS_EARTH_MILES)
latlong = (39.698840000000004, -104.975916)
index = tree.query_ball_point(latlong,r=1)

结果是一个坐标数组，如下所示：

+---------------------------------------+
|              coordinates              |
+---------------------------------------+
| (39.676973877551, -104.966231826172)  |
| (39.6777407534644, -104.988982458831) |
| ...                                   |
+---------------------------------------+

当我尝试使用半正弦公式验证这些结果时，我看到第一个坐标是 1.6 英里

from haversine import haversine
haversine((39.676973877551, -104.966231826172),
         (39.698840000000004, -104.975916),miles=True)

1.5961362762187963

【问题讨论】：

标签： python scipy geolocation geospatial latitude-longitude

【解决方案1】：

Pysal 不使用 hasrsine 函数来计算 query_ball_point 方法的距离。它使用了 pysal.cg.sphere.arcdist 函数，这是不同的。

import pysal
from pysal.cg.kdtree import KDTree    

locations = [(40.702566, -73.816859),
         (40.70546, -73.810708),
         (40.709179, -73.820574),
         (40.700486, -73.807969),
         (40.694624, -73.820593),
         (40.695132, -73.820841),
         (40.694095, -73.821334),
         (40.694165, -73.822368),
         (40.695077, -73.822817),
         (40.6747769261, -73.8092618174)] 
tree = KDTree(locations, distance_metric='Arc', radius=pysal.cg.RADIUS_EARTH_MILES)
current_point = (40.709523, -73.802472)
# get all points within X miles of 'current_point'
indices = tree.query_ball_point(current_point, 1)
for i in indices:
    print(locations[i])

1英里内有3个点

(40.70546, -73.810708)
(40.700486, -73.807969)
(40.6747769261, -73.8092618174)

根据半正弦公式，并非所有这些点都在 1 英里范围内：

from haversine import haversine
for i in indices:
    print(haversine(current_points, locations[i], miles = True))

0.5146716729994124
0.6875825817591269
2.4269297885659022

但是根据 pysal 的 arcdist 公式（半径为 3958.756 英里），它们在 1 英里范围内：

from pysal.cg.sphere import arcdist
for i in indices:
    print(arcdist(current_points, locations[i], 3958.756))

0.5744128196875283
0.4178272122350164
0.8175408580090955

【讨论】：

【解决方案2】：

PySAL 期望输入为 (longitude, latitude)（即 x,y），而 hasrsine python 包期望输入为 (latitude, longitude)。否则 arcdist 和 hasrsine 应该返回几乎相同的结果。

from libpysal.cg.sphere import arcdist, RADIUS_EARTH_MILES
from haversine import haversine

locations = [(40.702566, -73.816859),
         (40.70546, -73.810708),
         (40.709179, -73.820574),
         (40.700486, -73.807969),
         (40.694624, -73.820593),
         (40.695132, -73.820841),
         (40.694095, -73.821334),
         (40.694165, -73.822368),
         (40.695077, -73.822817),
         (40.6747769261, -73.8092618174)] 
current_point = (40.709523, -73.802472)

H = [haversine(current_point, loc, unit='mi') for loc in locations]
print(', '.join(["%0.5f"%dist for dist in H]))
A = [arcdist(current_point[::-1], loc[::-1], radius=RADIUS_EARTH_MILES) for loc in locations]
print(', '.join(["%0.5f"%dist for dist in A]))
print(', '.join(['%0.8f'%(h-a) for h,a in zip(H,A)]))

输出：

0.89381, 0.51467, 0.94839, 0.68758, 1.40024, 1.38364, 1.45343, 1.48732, 1.46011, 2.42693
0.89381, 0.51467, 0.94838, 0.68758, 1.40024, 1.38364, 1.45343, 1.48732, 1.46011, 2.42693
0.00000123, 0.00000071, 0.00000131, 0.00000095, 0.00000193, 0.00000191, 0.00000201, 0.00000205, 0.00000202, 0.00000335

【讨论】：