1

I am trying to find an efficient way to calculate the distance to the nearest neighbour for a set of coordinates of form (lat, lon):

[[51.51045038114607, -0.1393407528617875],
[51.5084300350736, -0.1261805976142865],
[51.37912856172232, -0.1038613174724213]]

I previously had a working (i thought!) piece of code which used sklearn's NearestNeighbors to reduce the algorithmic complexity of this task:

from sklearn.neighbors import NearestNeighbors
from sklearn.metrics.pairwise import haversine_distances
from math import sin, cos, sqrt, atan2, radians

# coordinates
coords = [[51.51045038114607, -0.1393407528617875],
          [51.5084300350736, -0.1261805976142865],
          [51.37912856172232, -0.1038613174724213]]

# tree method that reduces algorithmic complexity from O(n^2) to O(Nlog(N))
nbrs = NearestNeighbors(n_neighbors=2,
                        metric=_haversine_distance
                        ).fit(coords)

distances, indices = nbrs.kneighbors(coords)

# the outputted distances
result = distances[:, 1]

The output is as follows:

array([ 1.48095104,  1.48095104, 14.59484348])

Which used my own version of the haversine distance as the distance metric

def _haversine_distance(p1, p2):
"""
p1: array of two floats, the first point
p2: array of two floats, the second point

return: Returns a float value, the haversine distance

"""
lon1, lat1 = p1
lon2, lat2 = p2

# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

# get the deltas
dlon = lon2 - lon1
dlat = lat2 - lat1

# haversine formula
a = np.sin(dlat/2)**2 + (np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2)
c = 2 * np.arcsin(np.sqrt(a))

# approximate radius of earth in km
R = 6373.0

# convert to km distance
distance = R * c

return distance

These distances are wrong, my first question is, why is this? Is there any way I can correct this while retaining the algorithmic simplicity of the NearestNeighbors method?

I then discovered I can get the correct answer by using the geopy.distance method, however this does not come with in-build techniques to reduce the complexity and therefore computation time

import geopy.distance

coords_1 = (51.51045038, -0.13934075)
coords_2 = (51.50843004, -0.1261806)

geopy.distance.geodesic(coords_1, coords_2).km

My second question is then, are there any implementations of this method that reduce the complexity, otherwise I will be forced to use nested for loops to check the distance between every point and all others.

Any help appreciated!

Related Question Vectorised Haversine formula with a pandas dataframe

user17033672
  • 15
  • 5

1 Answers1

0

Depending how big your dataset is, it might be more efficient to convert your data to a Pandas dataframe.

import pandas as pd
import numpy as np

def haversine(lon1, lat1, lon2, lat2):
    lon1, lat1, lon2, lat2 = np.radians([lon1, lat1, lon2, lat2])
    dlon = lon2 - lon1
    dlat = lat2 - lat1

    haver_formula = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2

    r = 6371 #6371 for distance in KM for miles use 3958.756
    dist = 2 * r * np.arcsin(np.sqrt(haver_formula))
    return pd.Series(dist)

#added random id number
df = pd.DataFrame({'id':[123,456,789],
                   'lat':[51.51045038114607, 51.5084300350736, 51.37912856172232],
                   'lon':[-0.1393407528617875, -0.1261805976142865, -0.1038613174724213]})

>>> df
    id        lat       lon
0  123  51.510450 -0.139341
1  456  51.508430 -0.126181
2  789  51.379129 -0.103861

#self merging (faster than iterating through rows in larger datasets)
df2 = pd.merge(df.assign(key=1),df.assign(key=1), on='key', suffixes=('', '_2')).drop('key', axis=1)


>>> df2
    id        lat       lon  id_2      lat_2     lon_2
0  123  51.510450 -0.139341   123  51.510450 -0.139341
1  123  51.510450 -0.139341   456  51.508430 -0.126181
2  123  51.510450 -0.139341   789  51.379129 -0.103861
3  456  51.508430 -0.126181   123  51.510450 -0.139341
4  456  51.508430 -0.126181   456  51.508430 -0.126181
5  456  51.508430 -0.126181   789  51.379129 -0.103861
6  789  51.379129 -0.103861   123  51.510450 -0.139341
7  789  51.379129 -0.103861   456  51.508430 -0.126181
8  789  51.379129 -0.103861   789  51.379129 -0.103861

#drop duplicates
df2 = df2[df2['id']!=df2['id_2']].reset_index(drop=True)

>>> df2
    id        lat       lon  id_2      lat_2     lon_2
0  123  51.510450 -0.139341   456  51.508430 -0.126181
1  123  51.510450 -0.139341   789  51.379129 -0.103861
2  456  51.508430 -0.126181   123  51.510450 -0.139341
3  456  51.508430 -0.126181   789  51.379129 -0.103861
4  789  51.379129 -0.103861   123  51.510450 -0.139341
5  789  51.379129 -0.103861   456  51.508430 -0.126181

#find distance
df2['dist'] = haversine(df2['lon'], df2['lat'], df2['lon_2'], df2['lat_2'])

>>> df2
    id        lat       lon  id_2      lat_2     lon_2       dist
0  123  51.510450 -0.139341   456  51.508430 -0.126181   0.938061
1  123  51.510450 -0.139341   789  51.379129 -0.103861  14.807897
2  456  51.508430 -0.126181   123  51.510450 -0.139341   0.938061
3  456  51.508430 -0.126181   789  51.379129 -0.103861  14.460639
4  789  51.379129 -0.103861   123  51.510450 -0.139341  14.807897
5  789  51.379129 -0.103861   456  51.508430 -0.126181  14.460639

#closest neighbor
>>> df2[['id', 'lat', 'lon', 'dist']].sort_values(['id', 'dist']).groupby('id').first().reset_index()
    id        lat       lon       dist
0  123  51.510450 -0.139341   0.938061
1  456  51.508430 -0.126181   0.938061
2  789  51.379129 -0.103861  14.460639
amance
  • 883
  • 4
  • 14