r = 3
c = 3
lenr = a.shape[0] / r
lenc = a.shape[1] / c
np.array([a[i * r:(i + 1) * r, j * c:(j + 1) * c] for i, j in np.ndindex(
    lenr, lenc)]).reshape(lenr, lenc, r, c)
Индексы подматриц фиксированного размера массива numpy