在图像的浮点高度图上，我迭代数组中的每个2x2正方形子矩阵，执行计算，对结果求和。这是缓慢的，因为海拔Map很大（16K x 16K），并且循环比NumPy或SciPy慢。如何遍历多维NumPy数组块？
这个3x3 NumPy数组可以是NxM：

[[0.0, 1.0, 2.0],
 [3.0, 4.0, 5.0],
 [6.0, 7.0, 8.0]]

我想要一个快速迭代器，它产生：

> [0.0, 1.0, 3.0, 4.0]
> [1.0, 2.0, 4.0, 5.0]
> [3.0, 4.0, 6.0, 7.0]
> [4.0, 5.0, 7.0, 8.0]

子矩阵中的值的顺序并不重要，只要它们是一致的（逆时针、顺时针、Z字形等）。我的代码不使用NumPy：

shape_dem_data = shape_dem.getdata() # shape_dem is a PIL image

for x in range(width - 1):
    for y in range(height - 1):
        i = y * width + x
        z1 = shape_dem_data[i]
        z2 = shape_dem_data[i + 1]
        z3 = shape_dem_data[i + width + 1]
        z4 = shape_dem_data[i + width]
        # Create a bit-mask indicating the available elevation data
        mask = (z1 != NULL_HEIGHT) << 3 |\
               (z2 != NULL_HEIGHT) << 2 |\
               (z3 != NULL_HEIGHT) << 1 |\
               (z4 != NULL_HEIGHT) << 0
        if mask == 0b1111:
            # All data available.
            surface_area += area_of_triangle(((0, 0, z1), (gsd, 0, z2), (gsd, gsd, z3)))
            surface_area += area_of_triangle(((0, 0, z1), (gsd, gsd, z3), (0, gsd, z4)))
            pass
        elif mask == 0b1101:
            # Top left triangle
            surface_area += area_of_triangle(((0, 0, z1), (gsd, 0, z2), (0, gsd, z4)))
        elif mask == 0b0111:
            # Bottom right triangle
            surface_area += area_of_triangle(((gsd, 0, z2), (gsd, gsd, z3), (0, gsd, z4)))
        elif mask == 0b1011:
            # Bottom left triangle
            surface_area += area_of_triangle(((0, 0, z1), (gsd, gsd, z3), (0, gsd, z4)))
        elif mask == 0b1110:
            # Top right triangle
            surface_area += area_of_triangle(((0, 0, z1), (gsd, 0, z2), (gsd, gsd, z3)))
return surface_area

目的是计算高度数组的表面积（给定像素之间的固定采样距离）。位屏蔽检查哪些像素组合不是“空”高度（并相应地调整计算）。