球面上の（地理座標）点集合に対する近傍探索

こんにちは。
球面上の（地理座標）点の集合 $\{(\lambda_i,\,\phi_i),\, i=0,1,2,...,n-1\}$に対して、そのうちの一点を中心とする指定半径の円内に含まれる他の点を探索することを考えます。これを、全ての点に対してそのような部分集合を作る問題を計算してみました。結果例（ $i$ の部分集合）は、

$ ./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])])