Mosaicing GeoTIFFs with gdal creates gaps due to alpha layers

Question

I have a series of GeoTIFFs which I'm trying to mosaic together, however when they are combined their individual alpha layers overlap sections of the other images and create gaps in the final mosaic:

My approach is first to generate a .vrt from all the GeoTIFFs:

gdalbuildvrt -hidenodata mosaic.vrt geo_pict20140910_131*

And then to convert that .vrt to a GeoTIFF with gdal_translate:

gdal_translate mosaic.vrt mosaic.tif

Is there a way I can mosaic only the non-alpha portions of these GeoTIFFs to avoid these gaps?

Here's a link to the original GeoTIFFs: https://www.dropbox.com/s/udjq9bter9fmgxz/gtiffs.zip?dl=0

Here's a link to the .vrt and final .tif mosaics: https://www.dropbox.com/s/oov9vatz69f31zg/mosaic.zip?dl=0

Here's the gdalinfo for the final mosaic:

Driver: GTiff/GeoTIFF
Files: mosaic.tif
Size is 4782, 2837
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-121.460884208999985,38.583432355500015)
Pixel Size = (0.000000591354474,-0.000000591354474)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-121.4608842,  38.5834324) (121d27'39.18"W, 38d35' 0.36"N)
Lower Left  (-121.4608842,  38.5817547) (121d27'39.18"W, 38d34'54.32"N)
Upper Right (-121.4580564,  38.5834324) (121d27'29.00"W, 38d35' 0.36"N)
Lower Right (-121.4580564,  38.5817547) (121d27'29.00"W, 38d34'54.32"N)
Center      (-121.4594703,  38.5825935) (121d27'34.09"W, 38d34'57.34"N)
Band 1 Block=4782x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=4782x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=4782x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=4782x1 Type=Byte, ColorInterp=Alpha

Answer 1

Found a quick solution - replacing "-hidenodata" with "-srcnodata 0" in the .vrt build:

gdalbuildvrt -srcnodata 0 mosaic.vrt geo_pict20140910_131*

Answer 2

Jeremy's answer is nice but I found an issue with when there is an RGB band with a value of 0

This causes, for example, a pixel in a dark area in the shadow of a tree to become nodata, because there was an RGB value there of 0

What has been working so far for me is setting nodata to all 4 bands of the image in areas where the alpha band was originally 0

The below code has been helping me so far

Curious to know if there are any more concise solutions than this

import glob, os, processing


#this is where we set the processing folders
inputFolder = "D:/Imagery 2022/QGIS Desperation Band Splitting/In/"
bandHoldFolder = "E:/Imagery 2022/QGIS Desperation Band Splitting/BandsP2/"
outputFolder = 'E:/Imagery 2022/QGIS Desperation Band Splitting/OutP2/'
os.chdir(inputFolder)



#now we use a for to count through the images
for fname in glob.glob("*.tif"):
    
    fullPath = inputFolder + fname
    print(fullPath)
    
    if (fname[-4])!= '.':
        print('Your raster extension needs to be 3 characters (I think)')
    
    
    #reset the list of the split out bands
    filelistBand = []
    
    
    #split out the first band and set any pixels with an alpha value less than 127 to be null through a division of zero
    output_raster = bandHoldFolder + fname[:-4] + "Band1.tif"
    proc = processing.run("qgis:rastercalculator", {'EXPRESSION':' \"' +  fname[:-4] + '@1\" * (1/(\"' +  fname[:-4] + '@4\" > 127))','LAYERS':[fullPath],'CELLSIZE':0,'EXTENT':None,'CRS':None,'OUTPUT':'TEMPORARY_OUTPUT'})
    tempOutput = proc['OUTPUT']
    #now we need to compress these bands because they're big as, note that I'm exporting to an Australian CRS
    processing.runAndLoadResults("gdal:translate", {'INPUT':tempOutput,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:7855'),'NODATA':None,'COPY_SUBDATASETS':False,'OPTIONS':'COMPRESS=LZW|PREDICTOR=2|ZLEVEL=9|NUM_THREADS=8|BIGTIFF=IF_SAFER','EXTRA':'','DATA_TYPE':0,'OUTPUT':output_raster})
    #add this to the list of bands to put in the virtual raster
    filelistBand.append(output_raster)
    os.remove(tempOutput)

    
    #split out the second band
    output_raster = bandHoldFolder + fname[:-4] + "Band2.tif"
    proc = processing.run("qgis:rastercalculator", {'EXPRESSION':' \"' +  fname[:-4] + '@2\" * (1/(\"' +  fname[:-4] + '@4\" > 127))','LAYERS':[fullPath],'CELLSIZE':0,'EXTENT':None,'CRS':None,'OUTPUT':'TEMPORARY_OUTPUT'})
    tempOutput = proc['OUTPUT']
    processing.runAndLoadResults("gdal:translate", {'INPUT':tempOutput,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:7855'),'NODATA':None,'COPY_SUBDATASETS':False,'OPTIONS':'COMPRESS=LZW|PREDICTOR=2|ZLEVEL=9|NUM_THREADS=8|BIGTIFF=IF_SAFER','EXTRA':'','DATA_TYPE':0,'OUTPUT':output_raster})
    filelistBand.append(output_raster)
    os.remove(tempOutput)

    
    #split out the third band
    output_raster = bandHoldFolder + fname[:-4] + "Band3.tif"
    proc = processing.run("qgis:rastercalculator", {'EXPRESSION':' \"' +  fname[:-4] + '@3\" * (1/(\"' +  fname[:-4] + '@4\" > 127))','LAYERS':[fullPath],'CELLSIZE':0,'EXTENT':None,'CRS':None,'OUTPUT':'TEMPORARY_OUTPUT'})
    tempOutput = proc['OUTPUT']
    processing.runAndLoadResults("gdal:translate", {'INPUT':tempOutput,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:7855'),'NODATA':None,'COPY_SUBDATASETS':False,'OPTIONS':'COMPRESS=LZW|PREDICTOR=2|ZLEVEL=9|NUM_THREADS=8|BIGTIFF=IF_SAFER','EXTRA':'','DATA_TYPE':0,'OUTPUT':output_raster})
    filelistBand.append(output_raster)
    os.remove(tempOutput)
    
    
    fullAlphaPath = bandHoldFolder + fname[:-4] + 'Alpha.tif'
    #we need to reclassify the alpha band to make sure it's either nodata or 255
    processing.run("native:reclassifybytable", {'INPUT_RASTER':fullPath,'RASTER_BAND':4,'TABLE':[128,256,255],'NO_DATA':-1,'RANGE_BOUNDARIES':0,'NODATA_FOR_MISSING':True,'DATA_TYPE':1,'OUTPUT':fullAlphaPath})
    filelistBand.append(fullAlphaPath)                                  
  

    #this is where we will send the output
    print(fname[:-4])
    outputVrt = outputFolder + fname[:-4] + '.vrt'

    #we make a virtual raster using the list of bands that we have been appending onto
    processing.runAndLoadResults("gdal:buildvirtualraster", {'INPUT':filelistBand,'RESOLUTION':0,'SEPARATE':True,'PROJ_DIFFERENCE':False,'ADD_ALPHA':False,'ASSIGN_CRS':None,'RESAMPLING':0,'SRC_NODATA':'','EXTRA':'','OUTPUT':outputVrt})

After this you can manually make a virtual raster of all of the constituent virtual rasters, before saving it out as a rendered image

The virtual raster of virtual rasters should show the blank areas as all nodata, with no gaps in between the original images

With the final rendered image showing those areas as 0 in all bands