0

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?

1 Answer 1

3

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)

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.