Python求和地理半径内的纬度/经度点并求和到网格

Question

基本上，我的实验程序试图找出在给定时间落在有效点半径（例如 50 公里）内的点数。我的数据是结构化的（但如果需要，我可以重组）三个单独的数组，例如：

1_LAT,1_LON,1_TIM

其中 1_LAT、1_LON、1_TIM 都包含大约 250 个值，分别对应于纬度、经度（十进制度）和时间。

我有 20 组这些数组（即 1_LAT、1_LON、1_TIM...20_LAT、20_LON、20_TIM）。

这是我想要完成的：

1) 计算出落在每个集合的特定半径内的纬度/经度集合的数量。例如，在1_TIM的有效时间，在1_LAT,1_LON的50km半径范围内有多少点在其他19组点中。然后我想遍历每个有效时间，以找出每个特定点和有效时间的有效半径中的点数。

我在下面附上了一张图片，以帮助直观地描述。

黑色方块代表 LAT_1/LON_1 数组中的所有点。 蓝色方块代表 LAT_n/LAT_n 数组中的所有点。

我想为每组纬度/经度数组计算每个有效时间每个半径中的点数。最终显示将是地理底图图像上每个网格点的密度求和栅格或网格网格（即计数数 / 20）。

我觉得 KDEtree 可能是完成此任务的最佳方式，但我对此几乎没有/没有经验。任何想法或建议将不胜感激。

Answer 1

您可以执行以下操作...首先，将每个组的 (x, y) 坐标分组到一个 points_x 数组中：

points_1 = np.column_stack((LAT_1, LON_1))
...
points_n = np.column_stack((LAT_n, LON_n))

将它们存储在数组列表中可能是个好主意：

points = [point_1, points_2, ..., points_n]

现在，用每组点制作一个 kdTree：

import scipy.spatial as spsp
kdtrees = [spsp.cKdTree(p) for p in point]

你准备好了。如果您现在运行以下代码：

r = whatever_your_threshold_value_is
points_within_r = np.zeros((len(kdtrees), len(kdtrees)), dtype=np.int)
for j in xrange(len(kdtrees)):
    for k in xrange(j+1, len(kdtrees)):
        points_within_r[j, k] = kdtrees[j].count_neighbors(kdtrees[k], r, 2)
points_within_r = points_within_r + points_within_r.T

您现在应该发现points_within_r[j, k] 包含points_j 中的多少点在points_k 中一个点的半径r 内。

请记住，这里的距离是坐标的欧几里得距离，忽略它们测量的是球面角这一事实。