python中的加速地理定位算法

Question

我有一组 100k 的地理位置（纬度/经度）和一个六边形网格（4k 多边形）。我的目标是计算位于每个多边形内的点总数。

我当前的算法使用 2 个 for 循环来遍历所有地理点和所有多边形，如果我增加多边形的数量，这真的很慢......你将如何加速算法？我上传了一个最小的示例，它创建了 100k 个随机地理点并在网格中使用了 561 个单元格...

我还看到读取 geo json 文件（带有 4k 多边形）需要一些时间，也许我应该将多边形导出到 csv 中？

hexagon_grid.geojson 文件： https://gist.github.com/Arnold1/9e41454e6eea910a4f6cd68ff1901db1

最小的python示例： https://gist.github.com/Arnold1/ee37a2e4b2dfbfdca9bfae7c7c3a3755

Answer 1

您无需显式测试每个六边形以查看给定点是否位于其中。

暂时让我们假设您的所有点都落在六边形网格的范围内。因为你的六边形形成了一个规则的格子，你只需要知道哪个六边形中心离每个点最近。

这可以使用scipy.spatial.cKDTree 非常有效地计算：

import numpy as np
from scipy.spatial import cKDTree
import json

with open('/tmp/grid.geojson', 'r') as f:
    data = json.load(f)

verts = []
centroids = []

for hexagon in data['features']:

    # a (7, 2) array of xy coordinates specifying the vertices of the hexagon.
    # we ignore the last vertex since it's equal to the first
    xy = np.array(hexagon['geometry']['coordinates'][0][:6])
    verts.append(xy)

    # compute the centroid by taking the average of the vertex coordinates
    centroids.append(xy.mean(0))

verts = np.array(verts)
centroids = np.array(centroids)

# construct a k-D tree from the centroid coordinates of the hexagons
tree = cKDTree(centroids)

# generate 10000 normally distributed xy coordinates
sigma = 0.5 * centroids.std(0, keepdims=True)
mu = centroids.mean(0, keepdims=True)
gen = np.random.RandomState(0)
xy = (gen.randn(10000, 2) * sigma) + mu

# query the k-D tree to find which hexagon centroid is nearest to each point
distance, idx = tree.query(xy, 1)

# count the number of points that are closest to each hexagon centroid
counts = np.bincount(idx, minlength=centroids.shape[0])

绘制输出：

from matplotlib import pyplot as plt

fig, ax = plt.subplots(1, 1, subplot_kw={'aspect': 'equal'})
ax.hold(True)
ax.scatter(xy[:, 0], xy[:, 1], 10, c='b', alpha=0.25, edgecolors='none')
ax.scatter(centroids[:, 0], centroids[:, 1], marker='h', s=(counts + 5),
           c=counts, cmap='Reds')
ax.margins(0.01)

根据您需要的准确度，我可以想出几种不同的方法来处理网格之外的点：

• 您可以排除落在六边形顶点外边界矩形之外的点（即x < xmin、x > xmax 等）。但是，这将无法排除位于网格边缘“间隙”内的点。

• 另一个直接的选择是根据六边形中心的间距在distance 上设置一个截止值，这相当于对外六边形使用圆形近似值。

• 如果准确性至关重要，那么您可以定义一个与六边形网格的外部顶点相对应的matplotlib.path.Path，然后使用它的.contains_points() method 来测试您的点是否包含在其中。与其他两种方法相比，这可能会更慢且更繁琐。