Spatial Filters

def linear_filter(map, mask, range_, type = "mean", sigma = None):

    if type == "mean":
        n = (2*range_ + 1)**2
        filter = [[1/n for _ in range(2*range_+1)] for k in range(2*range_+1)]

    elif type == "gaussian":
        n = 2*range_ + 1
        if sigma is None:
            sigma = range_
        filter = gkern(n, sigma)
    else:
        raise ValueError("Cannot identify method.")

    pmap = np.pad(map, (range_, range_), mode = "constant")
    pmask = np.pad(mask, (range_, range_), mode = "constant")

    meanres = scipy.signal.convolve2d(np.multiply(pmap,pmask), np.asarray(filter), boundary = "symm", mode = "same")
    maskres = scipy.signal.convolve2d(pmask, np.asarray(filter), boundary = "symm", mode = "same")

    r = np.divide(np.multiply(meanres, pmask), maskres)
    r[~pmask.astype("bool")] = np.nan
    return r[range_:-range_, range_:-range_]

def median_filter(map, mask, range_):
    pmap = np.pad(map, (range_, range_), mode = "constant")
    pmask = np.pad(mask, (range_, range_), mode = "constant")

    medianres = np.zeros_like(map)
    k = 2*range_ + 1
    n,m = map.shape
    for i in range(n):
        for j in range(m):

            region = pmap[i:i+k, j:j+k][pmask[i:i+k, j:j+k].astype("bool")]
            if len(region > 0):

                medianres[i, j] += np.median(region)

    medianres[~mask] = np.nan

    return medianres