# Calculate the average of neighbor pixels for raster edge

I have a small raster (or matrix) and I would like to know what is the difference between each pixel and the average of its' 8 neighbors in the case where the kernel window is on the edges and no values exist. When the 8 neighbor pixels have "real" values I can calculate the average of the 8 for each pixel, by using `scipy.ndimage.convolve` with kernel window as `np.array([[0.125,0.125,0.125],[0.125,0,0.125],[0.125,0.125,0.125]])`.

The problem is on the edges.

For example, if this is my matrix:

``````array([[13, 21, 13,  8],
[ 5, 10, 22, 14],
[21, 33,  9,  0],
[ 0,  0,  0,  0]])
``````

I would like to calculate the average of the neighbors of the upper left cell (value 13). The desired calculation will be

(21+10+5)/3 = 12

For the pixel in the first column and second row (value 5) the desired calculation will be -

(13+21+10+5)/4 = 12.25

Note that I need to ignored the central pixel (value 13 or 5) and all the cells that are outside of the 3x3 kernel window and also that the number of neighboring pixels with "real" values is different from pixel to pixel.

I read about the `mode` that deals with edges here and in other places but could not find a solution.

The easiest way is to combine numpy's `NaN` functions (in this case `nanmean`) and `ndimage.generic_filter`, like so:

``````import numpy as np
from scipy import ndimage

m = np.array([[13, 21, 13,  8],
[ 5, 10, 22, 14],
[21, 33,  9,  0],
[ 0,  0,  0,  0]], dtype=np.float)

result = ndimage.generic_filter(a, np.nanmean, size=3, mode='constant', cval=np.NaN)
``````

At this point `result` comes out as:

``````array([[ 12.25      ,  14.        ,  14.66666667,  14.25      ],
[ 17.16666667,  16.33333333,  14.44444444,  11.        ],
[ 11.5       ,  11.11111111,   9.77777778,   7.5       ],
[ 13.5       ,  10.5       ,   7.        ,   2.25      ]])
``````

If you want to ignore the centre cell as well you can do so by supplying a footprint:

``````mask = np.ones((3, 3))

result = ndimage.generic_filter(a, np.nanmean, footprint=mask, mode='constant', cval=np.NaN)
``````

This returns the `result`

``````array([[ 12.        ,  12.6       ,  15.        ,  16.33333333],
[ 19.6       ,  17.125     ,  13.5       ,  10.4       ],
[  9.6       ,   8.375     ,   9.875     ,   9.        ],
[ 18.        ,  12.6       ,   8.4       ,   3.        ]])
``````

Which we can check for the top left cell as `(21 + 5 + 10) / 3 == 12` and the second row, first column (value 5) becomes `(13 + 21 + 10 + 21 + 33) / 5 == 19.6`, and the bottom right (value 0) becomres `(9 + 0 + 0) / 3 == 3`.

@om_henners gives a `generic_filter` method that works well for small arrays, which is the intended use case from the original question; however, this method can be slow for medium and large arrays. A similar approach using `convolve2d` will produce identical results and can provide substantial speed improvements, as demonstrated below. With a `(2048, 512)` array, I see a speedup of ~300 when using the `convolve2d` approach; even with a `(4, 4)` array, I see a speedup of ~20 when using the `convolve2d` approach. Cheers!

``````import numpy as np
from scipy.signal import convolve2d
from scipy.ndimage import generic_filter
import time

def eight_neighbor_average_generic_filter(x):
footprint = np.ones((3, 3))
footprint[1, 1] = 0

result = generic_filter(
x.astype('float'), np.nanmean, footprint=footprint,
mode='constant', cval=np.nan)

return result

def eight_neighbor_average_convolve2d(x):
kernel = np.ones((3, 3))
kernel[1, 1] = 0

neighbor_sum = convolve2d(
x, kernel, mode='same',
boundary='fill', fillvalue=0)

num_neighbor = convolve2d(
np.ones(x.shape), kernel, mode='same',
boundary='fill', fillvalue=0)

return neighbor_sum / num_neighbor

if __name__ == '__main__':
# Test 8-neighbor averaging on small example array
x = np.array([
[13, 21, 13,  8],
[ 5, 10, 22, 14],
[21, 33,  9,  0],
[ 0,  0,  0,  0]])

y1 = eight_neighbor_average_generic_filter(x)
y2 = eight_neighbor_average_convolve2d(x)
print('\nmethods equivalent: %r' % np.alltrue(y1 == y2))

# Test 8-neighbor averaging on "large" array
x_large = np.random.randn(2048, 512)

t0 = time.time()
y1_large = eight_neighbor_average_generic_filter(x_large)
t_generic_filter = time.time() - t0

t0 = time.time()
y2_large = eight_neighbor_average_convolve2d(x_large)
t_convolve2d = time.time() - t0

print('\nExecution times for "large" array:')
print('    generic filter: %f s' % t_generic_filter)
print('    convolve2d: %f s' % t_convolve2d)
``````