[GIS] Interpolate nodata in rasters using Python

gdalinterpolationpythonrasterrasterio

I have a raster of the shape (1000,1000) and some areas having no data values. I would like to fill the data gaps by interpolating or tinning (does not matter) over the surrounding areas, however I fail to do that using Python.

I have searched tried some procedures already discussed at stackexchange, but failed to succeed:

Attempt 1: using gdal.FillNodata

import gdal, gdalconst
file = "C:\\Foo\\bar.tif"


raster = gdal.Open(file, gdalconst.GA_Update) 
fill = raster.GetRasterBand(1)
f = gdal.FillNodata(targetBand = fill, maskBand = None, maxSearchDist = 50, smoothingIterations = 0)

fill= fi.ReadAsArray()
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create("D:\\Foo\\out.tif", cols, rows, 1, gdal.GDT_Float32)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(fir)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()


### does not work - it takes ages, even with a smaller example. And if it works, there are no changes compared to the input file

Attempt 2: Using rasterio.fill — However, this returns an AttributeError.

Attempt 3: Using Scipy

import tifffile
from scipy.interpolate import griddata

raster = tifffile.imread("D:\\Foo\\bar.tif")

grid_x, grid_y = np.mgrid[0:1000, 0:1000]

nans = np.argwhere(np.isnan(raster))
coordnotnan = np.argwhere(~np.isnan(raster))
valuenotnan = excr[coordnotnan]
vanona = valuenotnan[:,0,0]

grid_z1 = griddata(coordnotnan,vanona,nans, method=nearest) #for sake of testing

grT = grid_z1.T
rR= list(raster.ravel())
new = []
i = 0
for na in rR:
    if np.isnan(na) == True:
        new.append(grT[i])
        i += 1
    else:
        new.append(na)

arrayed = np.array(new, dtype=np.float32).reshape((1000,1000))
#...
## Besides being inefficient (loops), it does not seem to take the values of the neighboring pixels

Do you spot the mistakes I've made? Is there another way to get this job done?

Best Answer

rasterio "fill" function is worked for me. using rasterio 1.1.0

import rasterio
from rasterio.fill import fillnodata

tif_file = r"D:\dem.tif"
with rasterio.open(tif_file) as src:
    profile = src.profile
    arr = src.read(1)
    arr_filled = fillnodata(arr, mask=src.read_masks(1), max_search_distance=10, smoothing_iterations=0)

newtif_file = r"D:\dem_filled.tif"   
with rasterio.open(newtif_file, 'w', **profile) as dest:
    dest.write_band(1, arr_filled)
Related Question