Search code examples
pythonnumpyscipyrasteriondimage

Apply Boxcar average to geospatial image

Assuming that the following array A is the result of reading a GeoTIFF image, for example with rasterio where nodata values are masked which is the array B.

I would like to apply a boxcar average smoothing over a square neighbourhood. The first problem is that I am not sure which scipy function represents a boxcar average?

I thought it might be the ndimage.uniform_filter. However, in contrast to scipy.signal, ndimage is not applicable to masked arrays.

from scipy.signal import medfilt
from scipy.ndimage import uniform_filter
import numpy as np

A = np.array([[-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999],
    [-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999],
    [-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999],
    [-9999, -9999, -9999, 0, 300, 400, 200, -9999],
    [-9999, -9999, -9999, -9999, 200, 0, 400, -9999],
    [-9999, -9999, -9999, 300, 0, 0, -9999, -9999],
    [-9999, -9999, -9999, 300, 0, -9999, -9999, -9999],
    [-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999]])

B = np.ma.masked_array(A, mask=(A == -9999))
print(B)


filtered = medfilt(B, 3).astype('int')
result = np.ma.masked_array(filtered, mask=(filtered == -9999))
print(result)

boxcar = ndimage.uniform_filter(B)
print(boxcar)

So, how can I apply a boxcar average that accounts for nodata values such as scipy.signal.medfilt?

Solution

  • This seems to be a good solution:

    import numpy as np
from scipy.signal import fftconvolve

def boxcar(A, nodata, window_size=3):

    mask = (A==nodata)
    K = np.ones((window_size, window_size),dtype=int)

    out = np.round(fftconvolve(np.where(mask,0,A), K, mode="same")/fftconvolve(~mask,K, mode="same"), 2)
    out[mask] = nodata

    return np.ma.masked_array(out, mask=(out == nodata))

A = np.array([[100, 100, 100, 100, 100, 100, 100, 100],
              [100, 100, 100, 100, 100, 100, 100, 100],
              [100, 100, 100, 100, 100, 100, 100, 100],
              [100, 100, 100, 100, 1  , 0  , 1  , 100],
              [100, 100, 100, 1  , 0  , 1  , 0  , 100],
              [100, 100, 100, 0  , 1  , 0  , 1  , 100],
              [100, 100, 100, 100, 100, 100, 100, 100]])

print(boxcar(A, 100))

    Would be great to get some feedback, in particular on improvements!