Scipy: how to convert KD-Tree distance from query to kilometers (Python/Pandas)

Question

This post builds upon this one.

I got a Pandas dataframe containing cities with their geo-coordinates (geodetic) as longitude and latitude.

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}]);

For each city I'm trying to find two other cities that are closest. Therefore I tried the scipy.spatial.KDTree. To do so, I had to convert the geodetic coordinates into 3D catesian coordinates (ECEF = earth-centered, earth-fixed):

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

This give me this:

With this I can create the KDTree:

coordinates = list(zip(df['x'], df['y'], df['z']))

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

Now I'm testing it with Berlin,

tree.query(coordinates[0], 2)

which correctly gives me Berlin (itself) and Potsdam as the two cities from my list that are closest to Berlin.

Question: But I wonder what to do with the distance from that query? It says 1501 - but how can I convert this to meters or kilometers? The real distance between Berlin and Potsdam is 27km and not 1501km.

Remark: I know I could get longitude/latitude for both cities and calculate the haversine-distance. But would be cool that use the output from KDTree instead.

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

Any help is appreciated.

Angarsk

Answer 1

The KDTree is computing the euclidean distance between the two points (cities). The two cities and the center of the earth form an isosceles triangle.

The German wikipedia entry contains a nice overview of the geometric properties which the English entry lacks. You can use this to compute the distance.

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

---

Update

After the comment about obtaining the opposite I re-read the question and realised that while it seems that one can use the proposed function above, the real problem lies somewhere else.

cos and sin in your function to_Cartesian expect the input to be in radians (documentation) whereas you are handing them the angles in degree. You can use the function deg2rad defined above to transform the latitude and longitude to radians. This should give you the distance in km directly from the KDTree.

Scutari

Can you give the opposite too? Convert let's say 400m to euclidian distance? – Angarsk

@Angarsk I updated my answer. Now that I thought more in depth about it, I am not sure why the angle between the two vectors (from the center of the earth to the two cities respectively) seems to be conserved despite handing degree instead of radians to your conversion function. Sorry that I didn't spot that earlier. – Scutari

Excellent, thanks for the hint. You're right, now the distance that is returned from the KD-Tree is fine too. – Angarsk

the deg2rad function is not necessary if you use the radians function from the math package (from math import radians) – Shae

@Scutari ... this function just returns the input converted to radians, nothing more! ( so it is IDENTICAL to np.deg2rad(x) ) – Marolda