[PYTHON] Neighborhood search for a set of (geographical coordinates) points on a sphere

Hello. For a set of (geographical coordinates) points on a sphere $ \ {(\ lambda_i, , \ phi_i), , i = 0,1,2, ..., n-1 \} $ Consider searching for other points contained within a circle with a specified radius centered on one point. I calculated the problem of making such a subset for all points. An example result (a subset of $ i $) is

$ ./neighbors.py -n 10000
[[0, 69, 1720, 7675, 8172, 9198], [1, 1167, 4385, 5063, 7716, 9184], [2, 2129, 4849, 7180], [3, 374, 753, 1039, 1506, 3230, 7188, 7433], [4, 2836, 4777, 7088, 9563], [5, 2936], [6, 9214, 9839], [7, 3288, 3762, 8710]]

neighbors.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import print_function
import numpy as np
import scipy
from scipy import spatial

DEGREE = scipy.pi / 180
RE = 6378137.0  # GRS80

def mercator_projection(lon, lat):
    lon, lat = [x * DEGREE for x in [lon, lat]]
    return lon, scipy.log(scipy.tan(lat/2+scipy.pi/4))

def neighbors(points, radius):
    r = radius / RE
    points = [mercator_projection(*p) for p in points]
    tree = spatial.cKDTree(points)
    neigh = [tree.query_ball_point([x, y], r*scipy.cosh(y)) for x, y in points]
    return neigh

def labelsorted(labels, i):
    def key(x):
        if x == i: return -1
        return x
    return sorted(labels, key=key)

import sys
n = int(sys.argv[2]) if len(sys.argv) > 2 and sys.argv[1]=='-n' else 10000
radius = 1e3  # 1 km
n_print = 8
d = 4e-6 * radius * scipy.sqrt(n)
points = np.random.uniform(low=-1,high=1,size=(n,2))*d + [135, 35]  # lon, lat
neigh = neighbors(points, radius)
print([labelsorted(x, i) for i, x in enumerate(neigh[:n_print])])