I have a DEM of a lake bottom, that contains a few holes (nodata) which I'd like to fill. For that purpose i'm using rasterio's fillnodata with a mask. Unfortunatly, I don't get the results I am expecting. Namely, that the holes within the lake are filled, but no changes at the 'lakeshore' are made. I'm creating a mask in the process, to indicate where the interpolation should happen and where it shouldn't but that doesn't seem to work. Any suggestions on what I might be doing wrong? Here's my code:
import os
import numpy as np
import rasterio as rio
from rasterio.fill import fillnodata
from rasterio.features import sieve
path_in = r'Z:\..\dhms\Lake_process_test\raw'
path_out = r'Z:\...\dhms\Lake_process_test\filled'
path_msk = r'Z:\...\dhms\Lake_process_test\msk'
with rio.open(os.path.join(path_in, 'Lake.tif'), "r") as src:
# collect metadata for output file
profile = src.profile
# create array which will be filled
arr = src.read(1)
# create mask indicating valid data area
msk = src.read_masks(1)
# fill holes in valid data area
sieved_msk = sieve(msk, size=2000)
# sieved_msk_inv = np.where((sieved_msk == 255), 0, 255)
# fill nodata areas
arr_filled = fillnodata(arr, mask=sieved_msk, max_search_distance = 1, smoothing_iterations=0)
# write filled file
with rio.open(os.path.join(path_out,'Lake_filled.tif'), "w", **profile) as filled:
filled.write(arr_filled, 1)
# write mask file
with rio.open(os.path.join(path_msk,'Lake_mask.tif'), "w", **profile) as mask:
mask.write(sieved_msk, 1)
Oh, and as you can see, I tried also inverting the mask, but that doesn't seem to be the problem.
EDIT: I've experimented a bit further: Inverting the mask leads to actually 'shrinking' the lake, whereas if I don't, it 'grows' but then the holes aren't filled.