3

早上好,我正在使用 Cressman 过滤器在 Numpy 中进行距离加权平均。我使用 Ball Tree 实施(感谢 Jake VanderPlas)返回请求数组中每个点的位置列表。查询数组(q ) 是形状 [n,3] 并且在每个点都有 x,y,z 点我想要对存储在树中的点进行加权平均.. 包裹在树周围的代码返回一定距离内的点所以我得到一个可变长度数组的数组。我使用 where 来查找非空条目(即在影响半径内至少有一些点的位置)创建 isgood 数组...

然后我遍历所有查询点以返回值 self.z 的加权平均值(请注意,这可以是 dims=1 或 dims=2 以允许多个共同网格化)

所以使用 map 或其他更快的方法复杂化的事情是 self.distances 和 self.locations 中数组长度的不均匀性......我对 numpy/python 仍然相当陌生,但我想不出办法这个数组明智(即不恢复到循环)

self.locations, self.distances = self.tree.query_radius( q, r, return_distance=True)
t2=time()
if debug: print "Removing voids"
isgood=np.where( np.array([len(x) for x in self.locations])!=0)[0]
interpol = np.zeros( (len(self.locations),) + np.shape(self.z[0]) )
interpol.fill(np.nan)
for dist, ix, posn, roi in zip(self.distances[isgood], self.locations[isgood], isgood, r[isgood]):
    interpol[isgood[jinterpol]] = np.average(self.z[ix], weights=(roi**2-dist**2) / (roi**2 + dist**2), axis=0)
    jinterpol += 1

所以...有关如何加快循环的任何提示?...

对于将天气雷达数据从范围、方位角、仰角网格映射到笛卡尔网格的典型映射,我有 240x240x34 点和 4 个变量需要 99 秒来查询树(由 Jake 在 C 和 cython 中编写。这是艰难的一步,因为您需要搜索数据!)和 100 秒来进行计算......在我看来这很慢?我的开销在哪里?np.mean 是有效的还是被称为数百万次的,这里是否有加速?我会通过使用float32而不是default64来获得......甚至缩放到整数(这将很难避免在加权中回绕......任何提示都非常感激!

4

1 回答 1

0

您可以在以下位置找到有关 Cressman 方案与使用高斯权重函数的相对优点的讨论:

http://www.flame.org/~cdoswell/publications/radar_oa_00.pdf

关键是将平滑参数与数据匹配(我建议使用接近数据点之间平均间距的值)。一旦知道了平滑参数,就可以将“影响半径”设置为等于权重函数下降到 0.01(或其他值)的半径。

速度有多重要?如果您愿意,而不是调用指数函数来确定权重,您可以为一些固定数量的半径增量组成一个离散的权重表,这大大加快了计算速度。理想情况下,您应该拥有网格边界之外的数据,这些数据可用于映射网格点周围的值(甚至在网格的边界点上)。请注意,这不是一个真正的插值方案 - 它不会准确返回数据点的观察值。像 Cressman 方案一样,它是一个低通滤波器。

于 2011-06-03T15:37:15.887 回答