在给定点周围计算网格值平均值的最快方法

Question

我有一个二维的numpy数组，里面的每个元素代表一个网格点。这个网格的每个方块边长是13公里。我需要计算一个特定网格点周围50英里内所有点的平均值。

我现在的做法是先确定一个边界框，然后通过索引来引用这个框内的数组元素，但这样在numpy中速度比较慢。我想找一个更快的解决方案。

目前的解决方案：

num_x = 400        #horizontal dimension of the 2D array
num_y = 300        #vertical dimension of the 2D array
num_dx = 6         #maximum number of horizontal grid points that fit within 50 miles
num_dy = 6         #same as above but for vertical (square grid)
radius_m = 80467.2 #50 miles expressed in meters
values = []        # stores the extracted values

for ix in range(-num_dx,num_dx+1):
    for jy in range(-num_dy,num_dy+1):
        # Determine distance to this point
        dist = ((ix*dx)**2+(jy*dy)**2)**0.5
        if dist <= radius_m:
            # Ensure this grid point actually exists within the grid
            if (j+jy) < num_y and (i+ix) < num_x:
                value = myarray[i+ix,j+jy]
                    if value is not masked and value >= 0:
                        values.append(float(value))

average = sum(values) / float(len(values))

这个方法比较慢（大约需要1.5秒），因为我需要访问myarray超过100次才能提取一个元素的值。有没有更好的向量化方法可以用在这里？我似乎无法用掩码来解决这个问题，因为条件是基于网格点相对于另一个点的位置，而不是元素本身的值。

Answer 1

对于图像内部的点（也就是半径不会超出图像的地方），你可以只计算一个掩码，这个掩码可以用于任何内部点。首先，创建一个全是零的数组：

mask = np.zeros((2 * num_dx + 1, 2 * num_dy + 1), dtype=np.int)

假设你关注的点在这个数组的中心位置，那么把半径范围内的每个元素都设置为1（这里没有显示这个过程）。然后，

indices = np.argwhere(mask.ravel() == 1)

对于myarray中任何一个内部元素(i, j)，你可以这样获取半径内的值：

values = myarray[i-num_dx: i+num_dx+1, j-num_dy: j+num_dy+1].ravel()[indices]

对于靠近边缘的点，你需要先复制一份mask，然后把图像外的行和列设置为零，再设置indices。

Answer 2

你的代码无法运行，而且似乎在 i < num_dx 或 j < num_dy 的时候会出现一个错误（这时它会绕到数组的另一边）。不过根据你变量的命名，我会这样做：

# First make sure we stay in the grid
i1, i2 = max(i-num_dx, 0), min(i+num_dx+1, num_x)
j1, j2 = max(j-num_dy, 0), min(j+num_dy+1, num_y)

# Get the radius in blocks, grid should be homogeneous
radius_i = radius_m / 13000.0

# Calc distances per element by broadcasting
DX = np.arange(i1, i2) - i
DY = np.arange(j1, j2)[:, None] - j
mask = DX*DX + DY*DY <= radius_i*radius_i

# Get block of interest and apply mask
values = myarray[i1:i2, j1:j2][mask]