Overlapping Rasters as numpy arrays

Question

I have two raster files as numpy arrays. I'm combining the rasters into one file. The rasters overlap on several pixel but have different values for the the overlapping pixels (one raster has nodata value for pixels, the other has the actual elevation data point). How do I combine the two arrays such that I only take the actual elevation data point instead of the nodatavalue?

file one: [-999,-999,-999,-999]
file two: [7,4,3,2]

desired output for pixels after combination: [7,4,3,2]

My code that I've been using produces this image. https://i.stack.imgur.com/ejyX1.jpg

from osgeo import gdal
from gdalconst import GA_ReadOnly
f1 = r'C:\Users\User\Data\NHDPlusCA\NHDPlus18\NEDSnapshot\Ned18a\elev_cm'
f2 = r'C:\Users\User\Data\NHDPlusCA\NHDPlus18\NEDSnapshot\Ned18b\elev_cm'


xmins = []
ymins = []
xmaxs = []
ymaxs = []

files = [f1, f2]

for f in files:
    raster = gdal.Open(f, GA_ReadOnly)
    projection = raster.GetProjection()
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform()
    xmins.append(f_xmin)
    xmaxs.append(f_xmin + (f_pwidth*cols))
    ymaxs.append(f_ymax)
    ymins.append(f_ymax + (f_pheight*rows))
    del raster

#extents
xmax = max(xmaxs)
ymax = max(ymaxs)
xmin = min(xmins)
ymin = min(ymins)
newwidth = round((xmax - (xmin))/f_pwidth)
newheight = round((ymax -(ymin))/abs(f_pheight))
print('xmins', xmins)
print('ymins', ymins)
print('xmaxs', xmaxs)
print('ymaxs', ymaxs)
print('largest x', xmax)
print('smallest x', xmin)
print('largest y', ymax)
print('smallest y', ymin)
print('new raster width', newwidth)
print('new raster height', newheight)


driver = gdal.GetDriverByName('GTiff')
new_ds = driver.Create(
    r'C:\Users\User\Desktop\out.tif',
    newwidth,  # xsize
    newheight,  # ysize
    1,    # number of bands
    gdal.GDT_Int32,  # The datatype of the pixel values
    options=[  # Format-specific creation options.
        'TILED=YES',
        'BIGTIFF=IF_SAFER',
        'BLOCKXSIZE=256',  # must be a power of 2
        'BLOCKYSIZE=256',# also power of 2, need not match BLOCKXSIZE
        'COMPRESS=LZW'
    ])

# fill the new raster with nodata values
outband = new_ds.GetRasterBand(1)
outband.SetNoDataValue(-2147483647)
outband.Fill(-2147483647)
geotransform = [xmin, f_pwidth, 0, ymax, 0, f_pheight]
new_ds.SetGeoTransform(geotransform)
new_ds.SetProjection(projection)
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)
    outband.WriteArray(data,xoffset,yoffset)
    outband.FlushCache()


# When all references to new_ds are unset, the dataset is closed and flushed to disk
new_ds = None

Answer 1

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.