2

Is it possible to interpolate geotiff raster no value cells that are not connected to the geotiff border by other no value cells with QGIS or other freely available tools?

raster with gaps

The grey areas are gaps. I would like to have the small grey areas surrounded by RGB dots interpolated but leave the grey area in the upper part of the image (the sea) as it is. The image was created with the QGIS Raster -> Analysis -> Fill nodata with max distance=10. If I make the max distance larger, the dots at the coast start to change the sea, too.

QGIS Raster -> Analysis -> Fill nodata does not have such an option.

8
  • Have you tried interpolating then clip resulting raster with the shoreline?
    – BERA
    Dec 19 '20 at 18:36
  • Lots of manual work. This is just an example, usually the coastline is much longer and more complicated. Dec 19 '20 at 19:43
  • No it's not. One more proccessing tool, takes 10 s
    – BERA
    Dec 19 '20 at 20:12
  • Sorry, my answer was a bit sloppy. I meant to say that I did not have the coastline as a separate raster to use as a mask. Creating it by editing the geotiff is certainly possible but would take some time. Or have I grossly misunderstood anything? Dec 19 '20 at 21:09
  • I just made the mask tiff with PAINT.NET: gis.stackexchange.com/questions/115552/… Will update my QGIS and try again. Thanks for the idea! Dec 19 '20 at 21:54
4

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: enter image description here

After: enter image description here

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))
3
  • Great that you were able to help yourself. I am not so much in raster data, but I really wonder why there is no fitting tool for such task yet. The world is full of coast lines and no value areas. Your answer would be a lot more helpful for others , if there were at least minimal comments in the script, describing your approach.
    – Bernd V.
    Dec 24 '20 at 14:52
  • Just updated the code with comments. Merry ongoing Christmas! Dec 26 '20 at 12:38
  • Grrrrreat! Very kind of you! Also happy christmassing!
    – Bernd V.
    Dec 27 '20 at 15:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.