numpy-based spatial reduction

Question

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,

score 2 · Accepted Answer · answered Aug 29 '22 at 13:28

2

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 Aug 29 '22 at 13:28

Quang Hoang

That works! thanks for the quick response. I guess you're 'in general' statement is meant to apply to functions without an axis kwarg. For those as slow as me, the answer code is stacking each of the spatial blocks/windows into axis 1 and 3.. – cefect Aug 29 '22 at 14:16

1 Answers1