3

基本上,我的实验程序试图找到在给定时间落在有效点(例如,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 可能是实现这一目标的最佳方式,但我对此几乎没有/没有经验。任何想法或建议将不胜感激。

4

1 回答 1

2

您将执行以下操作...首先,将(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_jrpoints_k

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

于 2013-07-30T22:07:55.260 回答