0

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.

1 Answer 1

0

So i solved it by just clipping back to the original extent with the mask created in the beginning. I'm still confused though, because the mask in the 'fillnodata' step doesn't do what i'd expect. So further hint's on whats going on are still welcome ;).

with rio.open(os.path.join(path_in, 'Lake.tif'), "r") as src:
    
    # collect metadata for output file
    profile = src.profile
    
    # read band which should be filled as array
    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)
    
    # fill nodata areas
    arr_filled = fillnodata(arr, mask=src.read_masks(1), max_search_distance = 50, smoothing_iterations=0)
    
    # use mask to clip back to original extent
    arr_clipped = np.where((sieved_msk == 255), arr_filled, np.nan)  
    
    with rio.open(os.path.join(path_out, 'Lake_filled.tif'), "w", **profile) as filled:
        filled.write(arr_clipped, 1)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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