You need to create a raster that has spatial extents that contain both of your previous rasters (see Create Raster Dataset). I suggest setting the value for this new raster to be 0.
Then, you can add the existing rasters to the new blank raster. To fill in the gap, you can use Spatial Analyst, choosing the option Neighborhood Statistics to interpolate or smooth the gap. Depending on how big it is you may need to adjust the radius you are choosing to smooth.
You may also want to run the smooth part, and then recombine your rasters so that you do not loose resolution on your entire dataset.
Because you are limited to writing contiguous blocks of data to GDAL rasters, the best way around this is to create the output array first, then write the array to the output raster. Using this method, you can use numpy's boolean indexing to dictate what is written. Your last for loop would work as follows using this logic:
import numpy
no_data = -2147483647
output_array = numpy.full((newheight, newwidth), no_data, 'int32')
for f in files:
raster = gdal.Open(f, GA_ReadOnly)
f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform()
cols = raster.RasterXSize
rows = raster.RasterYSize
xoffset = int((f_xmin - xmin) / f_pwidth)
yoffset = int((f_ymax - ymax) / f_pheight) # Assumes grids align
band = raster.GetRasterBand(1)
data = band.ReadAsArray()
# Create a boolean array that allows data where they exist
mask = (output_array == no_data) & (data != band.GetNoDataValue()) # Note: this assumes you prefer data in order based on the files list
# Allocate data to output array
output_array[yoffset:yoffset + rows, xoffset:xoffset + cols][mask] = data[mask]
# Write data to the output raster
outband.WriteArray(output_array)
outband.FlushCache()
However, you will need to fully load the raster in to memory prior to flushing it to the disk. If memory-management is an issue, you can either use numpy.memmap, or read the overlapping block from the output raster for each iteration, which would look like this:
for f in files:
raster = gdal.Open(f, GA_ReadOnly)
f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform() #geotransform for the old files being iterated over
cols = raster.RasterXSize
rows = raster.RasterYSize
xoffset = int((f_xmin - xmin) / f_pwidth)
yoffset = int((f_ymax - ymax) / f_pheight)
band = raster.GetRasterBand(1)
data = band.ReadAsArray(0, 0, cols, rows)
data_update = outband.ReadAsArray(xoffset, yoffset, cols, rows)
mask = (data_update == no_data) & (data != band.GetNoDataValue())
data_update[mask] = data[mask]
outband.WriteArray(data_update, xoffset, yoffset)
outband.FlushCache()
One last general note of caution: your method assumes that raster grids and data types align.
Best Answer
You can use
gdal_translate -a_nodata <your_no_data_value> input.tif output.tifto get rid of the no_data values. In QGIS Raster>Conversions>Translate. After that you can build a vrt (Raster>Miscellaneous>Build Virtual Raster;gdalbuildvrt) or you can merge them. With the help of a vrt layer you can save any supported raster file format you want and you can use map extends (right-klick on layer and Save As).E.g. using a Linux machine you can write a bash script when you have a lot of raster files.