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