KDTree 正在返回半径之外的点

Question

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

问题是 query_ball_point 是 return 指定 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

Answer 1

Pysal 不使用 haversine 函数来计算 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)

根据 haversine 公式，并非所有这些点都在 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

Answer 2

PySAL 期望输入为 (longitude, latitude)（即 x,y），而半正弦 python 包期望（latitude,经度)。否则 arcdist 和 haversine 应该 return 接近相同的结果。

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

KDTree 正在返回半径之外的点

KDTree is Returning Points Outside of Radius

python

geolocation

geospatial

latitude-longitude

scipy