python - 查找多边形点的二维数组中最近的点

Question

我有一个二维地理坐标数组，如下所示

coords = np.array(
[[[54.496163, 21.770491],
  [54.495438, 21.755107],
  [54.494713, 21.739723],
  [54.493988, 21.724339],
  [54.493263, 21.708955]],
 [[54.504881, 21.769271],
  [54.504157, 21.753884],
  [54.503432, 21.738497],
  [54.502707, 21.72311 ],
  [54.501983, 21.707723]],
 [[54.5136, 21.768052],
  [54.512875, 21.752661],
  [54.512151, 21.737271],
  [54.511426, 21.72188 ],
  [54.510702, 21.70649 ]],
 [[54.522318, 21.766832],
  [54.521594, 21.751439],
  [54.52087, 21.736045],
  [54.520145, 21.720651],
  [54.519421, 21.705257]],
 [[54.531037, 21.765613],
  [54.530312, 21.750216],
  [54.529588, 21.734819],
  [54.528864, 21.719421],
  [54.52814, 21.704024]]]
)

在空间中它定义了一个多边形

我想找到某个点的coords中最近点的索引，例如pt = [54.5, 21.7]

coords 在这里可能看起来像一个平行四边形，但实际上它是一个形状 (1200, 1500, 2) 的多边形。出于显而易见的原因，我在这里显示了坐标[0:5,0:5]。 多边形的真实形状可以在这个question中找到.

现在我正在计算整个coords数组相对于点pt的欧氏距离，以找到最近的点，位于[r1,c1]

flidx = ((coords - pt) ** 2).sum(2).argmin()
r1 = int(flidx / coords.shape[1])
c1 = flidx % coords.shape[1]

但这需要太多时间。

我正在考虑在多边形中实现二分搜索，我可以将其分为 4 个部分，检查点在哪个部分内部，然后循环直到我有一个相对较小的点数组，例如 16 x 16 。然后应用欧氏距离法。

问题是我不知道如何检查一个点是否在多边形内部。矩形相当简单，但这不是一个。

任何有关此方法或任何其他方法来查找最近点的帮助将不胜感激。

谢谢

Answer 1

首先注意，数据不是完美的网格，但它是“类似网格的”

from netCDF4 import Dataset
import numpy as np
from matplotlib import pyplot as plt

group = Dataset('./coords.nc', 'r', format='NETCDF4')

# reverse the input so that the bottom left point is at [0, 0]
lat = np.array(group['latitude_in'])[::-1]
lon = np.array(group['longitude_in'])[::-1]

# plot a sub-grid
slat = np.array([arr[::100] for arr in lat[::100]]).flatten()
slon = np.array([arr[::100] for arr in lon[::100]]).flatten()
plt.scatter(slat, slon)
plt.show()

要查找集合中距某个目标点最近的点的坐标，您可以通过执行“更改基础”来获得合理的近似值(搜索的初始猜测)。 IE。如果从左下角到右下角的向量是 x 方向，左下角到左上角是 y 方向向量，则应用基础矩阵的更改会将点映射到单位正方形(不完美)。然后就可以算出相对坐标了。

最后，你可以沿着网格(从最初的猜测开始)向目标点的方向行走(即移动到最近的邻居)

import itertools

class NearestIndex:
    def __init__(self, points):
        self.points = points 
        self.size = np.array(self.points.shape[:2]) - 1  # 1199 x 1499

        self.origin = points[0][0]  # origin must be at [0, 0]
        dX = points[-1, 0] - self.origin # the X-direction
        dY = points[0, -1] - self.origin # the Y-direction
        self.M = np.linalg.inv(np.array([dX, dY])) # change of basis matrix

    def guess(self, target):
        """ guess the initial coordinates by transforming points to the unit square """
        p = map(int, self.size * np.matmul(target - self.origin, self.M))
        return np.clip(p, 0, self.size)  # ensure the initial guess is inside the grid

    def in_grid(self, index):
        return (index == np.clip(index, 0, self.size)).all()

    def distance_to_target(self, index):
        return np.linalg.norm(self.points[index] - self.target)

    def neighbour_distances(self, index):
        i, j = index
        min_dist = np.inf
        min_index = None       
        for di, dj in itertools.product((-1, 0, 1), repeat=2):
            neighbour = (i + di, j + dj)
            if not (di == dj == 0) and self.in_grid(neighbour):
                dist = self.distance_to_target(neighbour)
                if dist < min_dist:
                    min_dist, min_index = dist, neighbour

        return min_index, min_dist

    def find_nearest(self, target):
        self.target = target
        index = self.guess(target)  # make an initial guess
        min_dist = self.distance_to_target(index)  # distance to initial guess
        while True:
            # check the distance to the target from each neighbour of index
            neighbour, dist = self.neighbour_distances(index)
            if dist < min_dist:
                index, min_dist = neighbour, dist
            else:
                return index, min_dist

像这样使用

points = np.dstack([lat, lon])
indexer = NearestIndex(points)
index, dist = indexer.find_nearest(np.array([46, 15])) 

print(index, coords[index], dist)  # (546, 556) [46.004955 14.999708] 0.004963596377623203

它已经相当快了，但还有很大的优化空间。您可以记住函数distance_to_target，或者在走向该点的过程中使用不同的步长。