Scipy:如何将KD-Tree查询距离转换为千米(Python/Pandas)

5

本文建立在这篇文章的基础上。

我有一个Pandas数据框,其中包含城市及其地理坐标(大地)的经度和纬度。

import pandas as pd

df = pd.DataFrame([{'city':"Berlin", 'lat':52.5243700, 'lng':13.4105300},
                   {'city':"Potsdam", 'lat':52.3988600, 'lng':13.0656600},
                   {'city':"Hamburg", 'lat':53.5753200, 'lng':10.0153400}]);

我正在尝试为每个城市找到最接近的两个其他城市。因此,我尝试了scipy.spatial.KDTree。为此,我必须将地理坐标转换为三维笛卡尔坐标(ECEF = 以地球为中心,固定于地球):

from math import *

def to_Cartesian(lat, lng):
    R = 6367 # radius of the Earth in kilometers

    x = R * cos(lat) * cos(lng)
    y = R * cos(lat) * sin(lng)
    z = R * sin(lat)
    return x, y, z

df['x'], df['y'], df['z'] = zip(*map(to_Cartesian, df['lat'], df['lng']))
df

这给我带来了这个结果:

python output

有了这个,我可以创建 KDTree:
coordinates = list(zip(df['x'], df['y'], df['z']))

from scipy import spatial
tree = spatial.KDTree(coordinates)
tree.data

现在我正在用Berlin进行测试。
tree.query(coordinates[0], 2)

从我的列表中正确地给出了柏林(本身)和波茨坦作为离柏林最近的两个城市。

问题:但我想知道如何处理查询结果中的距离?它显示1501,但我如何将其转换为米或千米?柏林和波茨坦之间的实际距离是27公里,而不是1501公里。

备注:我知道我可以获取这两个城市的经度/纬度并计算哈弗辛距离。但是使用KDTree的输出会更好。

(array([ 0. , 1501.59637685]), array([0, 1]))

任何帮助都将不胜感激。

1个回答

4
KDTree正在计算两个点(城市)之间的欧几里得距离。这两个城市和地球的中心形成一个等腰三角形。德语维基百科条目包含了英语条目缺少的几何属性的良好概述。您可以使用此内容来计算距离。
import numpy as np

def deg2rad(degree):
    rad = degree * 2*np.pi / 360
    return(rad)

def distToKM(x):
    R = 6367 # earth radius
    gamma = 2*np.arcsin(deg2rad(x/(2*R))) # compute the angle of the isosceles triangle
    dist = 2*R*sin(gamma/2) # compute the side of the triangle
    return(dist)

distToKM(1501.59637685)
# 26.207800812050056

更新

在获取相反意见的评论后,我重新阅读了问题,并意识到虽然似乎可以使用上面提出的函数,但真正的问题在于其他地方。

to_Cartesian 函数中的 cossin 期望输入以弧度为单位(documentation),而你将角度传递给它们。你可以使用上面定义的 deg2rad 函数将纬度和经度转换为弧度。这应该直接从 KDTree 中给出距离,单位是千米。


你能给出相反的吗?比如将400米转换为欧几里得距离? - Matthias
2
@Matthias 我更新了我的回答。现在我更深入地思考了一下,我不确定为什么两个向量之间的角度(从地球中心到两个城市分别)似乎保持不变,尽管将度数而不是弧度传递给您的转换函数。很抱歉我之前没有注意到这一点。 - Niklas
太好了,谢谢你的提示。你说得对,现在从KD-Tree返回的距离也很好。 - Matthias
如果您使用 math 包中的 radians 函数(from math import radians),则不需要使用 deg2rad 函数。 - Oren Pinsky
@Niklas …这个函数只是将输入转换为弧度并返回,没有其他作用!(因此与np.deg2rad(x)完全相同) - raphael

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接