[FIXED] numpy-based spatial reduction

Issue

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,

Solution

With the max (or other function supporting axis), you can just reshape the array:

a.reshape(a.shape[0]//2, 2, a.shape[1]//2, 2).max(axis=(1,3))

In general, you can reshape, swap the axis, flatten the 2x2 into a new axis 4, then work on that axis.

Answered By - Quang Hoang