How to average slices of a 3D Matrix maintaining its shape

Question

I got this working code snippet:

import numpy as np
from matplotlib import pyplot as plt

in_raster = np.random.randn(36, 3, 2151)
matrix = np.reshape(in_raster, [(np.shape(in_raster)[0] * np.shape(in_raster)[1]), np.shape(in_raster)[2]])
#  reshaping the matrix to prepare loop
out_raster = np.empty([np.shape(in_raster)[0]/3, np.shape(in_raster)[1]/3, np.shape(in_raster)[2]])
#  creating empty output matrix

i = 0
j = 0
while i <= len(in_raster)-9 or j < len(out_raster):
    if i % 9 == 0:
        avg_in_raster = np.nanmean(matrix[i:i+9, :], axis=0)
        out_raster[j] = avg_in_raster
    i += 9
    j += 1
out_raster = np.reshape(out_raster, [np.shape(out_raster)[0], np.shape(in_raster)[1]/3, np.shape(in_raster)[2]])

# plot example
low = 0
high = 50

for row in range(0, 3):
    for col in range(np.shape(in_raster)[1]):
        plt.plot(range(low,high), (in_raster[row, col, low:high]))
plt.plot(range(low,high), (out_raster[0,0,low:high]), 'k')
plt.show()

The program averages (aggregates) 3x3 slices of the input matrix (a raster image) and sets up a new one maintainig the dimensionality of the original matrix.

Now I got the feeling that there must be an easier way to achieve this. Does somebody have an idea how to obtain the same result in a more pythonic way?

Thank you!

Answer 1

To my knowledge, there is no easier or quicker way to perform blockwise averaging. Your code might look big, but most of it is just preparation of arrays and resizing or plotting stuff. Your main function is a well-placed while-loop and the averaging itself you leave to numpy which is already a shortcut and should run quickly.

I don't see any reason to further shorten this, without losing readability.