Numpy aggregate across multiple axes

Let's say I have a 3d numpy array.shape of (27,27,27). I want to compress this to (9,9,9) by averaging every 3 elements across every axis simultaneously (e.g. make 3x3x3 pixels into 1x1x1).  The objective is to effectively compress by a single integer across all three axes simultaneously (with the assumption that any array will have a multiple of that integer for the shape of each axes).

My initial attempt was to use np.apply_over_axes, though I'm worried it is not getting the cubic mean of all 3 axes but instead averaging each sequentially.

def mean_over(arr, axis):

    np.average(arr.reshape(-1, 3), axis=axis)

the_array_small = np.apply_over_axes(mean_over, the_array, \[0,1,2\])

However this returns an error:

Traceback (most recent call last):

  File "\<stdin\>", line 1, in \<module\>

  File "\<\__array_function_\_ internals\>", line 180, in apply_over_axes

  File "/opt/homebrew/Caskroom/mambaforge/base/envs/seaborn/lib/python3.10/site-packages/numpy/lib/shape_base.py", line 496, in apply_over_axes

    if res.ndim == val.ndim:

AttributeError: 'NoneType' object has no attribute 'ndim'

I'm not convinced my apply_over_axes solution gets the aggregation I'm aiming for though. Ideally the mean of each (3,3,3) component is returned.

Just a first answer (but again, depending on your answer to my comment, it could be better to rely on some convolution)

arr=np.randon.rand(27,27,27)
the_array_small = arr.reshape(9,3,9,3,9,3).mean(axis=(1,3,5))

Another solution but take advantage of as_strided

a = np.arange(27**3).reshape(27, 27, 27)
tile_size = (3, 3, 3)
tile_shape = tuple(np.array(a.shape) // np.array(tile_size))
tile_strides = tuple(np.array(a.strides) * np.array(tile_size)) + tuple(a.strides)
tile_view = np.lib.stride_tricks.as_strided(
    a,
    shape=tile_shape + tile_size,
    strides=tile_strides,
    writeable=False,
)
result = np.mean(tile_view, axis=(-3, -2, -1))

Using the cubify function from this answer, you can break your array into cubes. With that result, you can then use apply_over_axes to get the averages and reshape for your desired result. Here I use an example of a 9x9x9 cube since it's easier to see the result that way.

import numpy as np

def cubify(arr, newshape):
    """https://stackoverflow.com/a/42298440/12131013"""
    oldshape = np.array(arr.shape)
    repeats = (oldshape / newshape).astype(int)
    tmpshape = np.column_stack([repeats, newshape]).ravel()
    order = np.arange(len(tmpshape))
    order = np.concatenate([order[::2], order[1::2]])
    # newshape must divide oldshape evenly or else ValueError will be raised
    return arr.reshape(tmpshape).transpose(order).reshape(-1, *newshape)

a = np.arange(9**3).reshape(9,9,9)
c = cubify(a, (3,3,3))
res = np.apply_over_axes(np.mean, c, [1,2,3]).reshape(3,3,3)

Result:

array([[[ 91.,  94.,  97.],
        [118., 121., 124.],
        [145., 148., 151.]],

       [[334., 337., 340.],
        [361., 364., 367.],
        [388., 391., 394.]],

       [[577., 580., 583.],
        [604., 607., 610.],
        [631., 634., 637.]]])