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))
