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