Code snippet to apply convolutional filter to a raster

run_ntl_filter.py

from math import sqrt

import numpy as np
from scipy import signal
import rasterio

# Code from https://github.com/carderne/gridfinder/blob/master/gridfinder/prepare.py
# And https://github.com/carderne/gridfinder/blob/master/gridfinder/_util.py

def filter_func(i, j):
    d_rows = abs(i - 20)
    d_cols = abs(j - 20)
    d = sqrt(d_rows**2 + d_cols**2)
    if i == 20 and j == 20:
        return 0
    elif d <= 20: 
        return 1 / (1 + d/2)**3
    else:
        return 0.0

def create_filter():
    vec_filter_func = np.vectorize(filter_func)
    ntl_filter = np.fromfunction(vec_filter_func, (41, 41), dtype=float)
    ntl_filter = ntl_filter / ntl_filter.sum()
    return ntl_filter
  
def save_raster(file, raster, affine):
   os.makedirs(os.path.dirname(file))
   filtered_out = rasterio.open(file, 'w', driver='GTiff',
                                height=raster.shape[0], width=raster.shape[1],
                                count=1, dtype=raster.dtype,
                                crs='+proj=latlong', transform=affine)
   filtered_out.write(raster, 1)
   filtered_out.close()
  
ntl_rd = rasterio.open(ntl_in)  # ntl_in is file path
ntl = ntl_rd.read(1)
affine = ntl_rd.transform
ntl_filter = create_filter()  # create_filter() code from the repo
ntl_convolved = signal.convolve2d(ntl, ntl_filter, mode='same')
ntl_filtered = ntl - ntl_convolved
save_raster('path/to/file.tif', ntl_filtered, affine)