Recherche de quartier pour un ensemble de points (coordonnées géographiques) sur une sphère

Bonjour. Pour un ensemble de (coordonnées géographiques) points sur une sphère $ \ {(\ lambda_i, , \ phi_i), , i = 0,1,2, ..., n-1 \} $ Envisagez de rechercher d'autres points contenus dans un cercle avec un rayon spécifié centré sur un point. J'ai calculé le problème de la création d'un tel sous-ensemble pour tous les points. Un exemple de résultat (un sous-ensemble de $ i $) est

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

Recommended Posts

Recherche de quartier pour un ensemble de points (coordonnées géographiques) sur une sphère
Remarques sur l'installation de Chainer 1.5 pour GPU sous Windows
Pour ceux d'entre vous qui ne savent pas comment définir un mot de passe avec Jupyter sur Docker