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.

2 Answers 2

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

I actually had the same problem for a while and I would do a clip after, which is fine for most solutions, but sometimes you actually only want certain holes filled and others left alone in a complex fashion.

My problem is that for the following code:

mask_area = "../mask_area.shp"
input_image = "../some_image.tif"

gdf = gpd.read_file(mask_area)
geometries = [mapping(geom) for geom in gdf.geometry]
geom_mask = rasterio.features.geometry_mask(geometries, transform=src.transform, invert=True, out_shape=(src.height, src.width))

with rasterio.open(input_image) as src:
    arr = src.read(1)
    threshold = -9998
    arr[arr < threshold] = np.nan
    valid_data_mask = np.isnan(arr)
    fill_mask = ~valid_data_mask & geom_mask
        
    arr = fillnodata(arr, mask=fill_mask, max_search_distance=search_distance, smoothing_iterations=0)

I was creating the wrong mask via:

fill_mask = valid_data_mask & geom_mask

or this:

fill_mask = ~valid_data_mask & geom_mask

When I should have been doing this:

fill_mask = ~(valid_data_mask & geom_mask)

Your Answer

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

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