[PYTHON] Nachbarschaftssuche nach einer Reihe von (geografischen Koordinaten) Punkten auf einer Kugel

Hallo. Für eine Menge von (geografischen Koordinaten) Punkten auf einer Kugel $ \ {(\ lambda_i, , \ phi_i), , i = 0,1,2, ..., n-1 \} $ Suchen Sie nach anderen Punkten in einem Kreis mit einem bestimmten Radius, der auf einem Punkt zentriert ist. Ich habe das Problem berechnet, eine solche Teilmenge für alle Punkte zu erstellen. Ein Beispielergebnis (eine Teilmenge von $ i $) ist

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