Interpolate only no value cells not connected to the egde of the raster

As I did not find a good solution, I set out to write it by myself.

I started with a routine to detect the outer cells to exclude them from the interpolation mask. Then I discovered that I could not get rasterio to fill the no data cells for reasons beyond my comprehension. I therefore included a custom interpolation as well.

It will be slow for large datasets but it works.

Before:

After:

The Python code:

import rasterio
from rasterio.fill import fillnodata
from os.path import join
from sys import argv
from time import time
from numpy import ndarray
from scipy.spatial.kdtree import KDTree

try:
    fpath = argv[1]
    fname = argv[2]
    fname_out = argv[3]
except IndexError:
    print("Usage: path, input_filename, output_filename, nodata if not in tif")
    exit(1)

dataset = rasterio.open(join(fpath, fname))
metadata = dataset.meta

if dataset.nodata is None:
    try:
        nodata = int(argv[4])
        #dataset.nodata = int(nodata)
    except IndexError:
        print("Nodata value missing!")
        exit(1)
else:
    nodata = dataset.nodata

print("Nodata value:", nodata)

red, green, blue, alpha = dataset.read(1), dataset.read(2), dataset.read(3), dataset.read(4)

image = dataset.read()

# input for kdtree during interpolation
point_array = []

# remaining no value cells updated every iteration
no_values = []

# size of the image
y_max = len(red)
x_max = len(red[0])

# temporary mask to store outer no value cells
mask = ndarray(shape=(y_max,x_max), dtype=int)
for i in range(0, y_max):
    for j in range(0, x_max):
        mask[i][j] = 1
        if (red[i][j] == nodata) and (green[i][j] == nodata) and (blue[i][j] == nodata):
            no_values.append((i, j))
        else:
            point_array.append((i,j))

# for interpolation
kdt = KDTree(point_array)

# find no value cells adjacent to the borders of the image or already found no value cells.
# return updated temporary mask & list of remaining no value cells
def search(novalues, m, x_max, y_max):

    result_mask = m.copy()
    result_novalues = []

    for y, x in novalues:
        if (x == 0) or (y == 0) or (x == x_max) or (y == y_max):
            result_mask[y][x] = 0

        elif (m[y-1][x] == 0) or (m[y+1][x] == 0) or (m[y][x-1] == 0) or (m[y][x+1] == 0):
            result_mask[y][x] = 0

        else:
            result_novalues.append((y,x))

    return result_mask, result_novalues

print("Searching for outside cells...")

# repeat the search until the temporary mask remains the same
count = 0
while True:

    start_time = time()

    mask1, no_values = search(no_values, mask, x_max-1, y_max-1)
    equal = True
    for i in range(0, y_max):
        if (mask[i] == mask1[i]).all():
            pass
        else:
            equal = False
            break
    if equal is True:
        break
    mask = mask1
    count += 1

    print(time()-start_time)

print(count, "iterations")
print("Remaining no value cells:", len(no_values))

# construct the mask actually used for interpolation
mask = ndarray(shape=(y_max,x_max), dtype=bool)
for i in range(0, y_max):
    for j in range(0, x_max):
        mask[i][j] = 1

# fill it with remaining (inland) no value cells
for y,x in no_values:
    mask[y][x] = 0

#image = fillnodata(image=image, mask=mask) # Does not work for no apparent reason :(

start_time = time()
print("Interpolating...")

# straightforward inverse nearest neighbour interpolation
for i in range(0, y_max):
    for j in range(0, x_max):

        if mask[i][j] == 0:

            image[3][i][j] = 255 # Alpha

            dist, ind = kdt.query((i, j), k=4)

            r, g, b = [], [], []

            print('n', i, j)
            for p_i in ind:
                print(point_array[p_i][0], point_array[p_i][1])
                y = point_array[p_i][0]
                x = point_array[p_i][1]
                r.append(red[y, x])
                g.append(green[y, x])
                b.append(blue[y, x])

            d_sum = sum(dist)
            new_r, new_g, new_b = 0, 0, 0
            for l, d in enumerate(dist):
                w = d / d_sum
                new_r, new_g, new_b = new_r + r[l]*w, new_g + g[l]*w, new_b + b[l]*w

            new_r, new_g, new_b = int(new_r), int(new_g), int(new_b)

            image[0][i][j], image[1][i][j], image[2][i][j] = new_r, new_g, new_b

            print(new_r, new_g, new_b)

print(time()-start_time)

# save the output
metadata.update({'driver': 'GTiff'})
rgb_dataset = rasterio.open(join(fpath, fname_out), 'w', **metadata)
rgb_dataset.write(image)
rgb_dataset.close()

print("Saved", join(fpath, fname_out))