I'm looking for a flexible fast method for computing a custom reduction on an np.array using a square non-overlapping window. e.g.,

array([[4, 7, 2, 0],
       [4, 9, 4, 2],
       [2, 8, 8, 8],
       [6, 3, 5, 8]])

let's say I want the np.max , (on a 2x2 window in this case) I'd like to get:

array([[9, 4],
       [8, 8]])

I've built a slow function using for loops, but ultimately I need to apply this to large raster arrays.
scipy.ndimage.generic_filter is close, but this uses sliding windows (with overlap), giving a result with the same dimensions (no reduction).
numpy.lib.stride_tricks.as_strided combined with a reducing function doesn't seem to handle relationships between rows (i.e., 2D spatial awareness).
rasterio has some nice resampling methods built on GDAL, but these don't allow for custom reductions.
skimage.transform.downscale_local_mean does not support custom functions on the blocks.
I'm sure there's something out there for custom spatial anti-aliasing , but I can't seem to find a solution and am feeling dumb.
Any help is greatly appreciated,