Write multi-layer TIFF using GDAL

Question

I am working with 4 band GeoTIFF images. I used GDAL to open these images. Each layer needs pre-processing and after completion of processing, I need to write back into a TIFF file. I used the below code, but I couldn't get all the bands together. When I check the written file, it contains only a single band.

for data_file in os.listdir(destSAR):
    data = data_file.split('_')
    joinStr = [data[6],data[7],data[8],data[9]]
    joinedName = '_'.join(joinStr)
    imgSAR = joinSAR+joinedName
    imgEO = joinEO+joinedName
    SARt = gdal.Open(imgSAR)
    opref = cv2.imread(imgEO)
    SAR=SARt.ReadAsArray()

    output=getSARmask(SAR,opref) #processing the SAR images
    os.chdir(outputPath) 
    filename = 'post_process_'+joinedName
#     print(len(output)) # have 4 bands in output
    for i in range(len(output)):
        ds = SARt #original SAR image with 4 bands, using for Geo information
        arr = output[i]
#         arr = band.ReadAsArray()
        [cols, rows] = arr.shape
        arr_min = arr.min()
        arr_max = arr.max()
        arr_mean = int(arr.mean())
        arr_out = np.where((arr < arr_mean), 10000, arr)
        driver = gdal.GetDriverByName("GTiff")
        outdata = driver.Create(filename, rows, cols, 1, gdal.GDT_UInt16)
        outdata.SetGeoTransform(ds.GetGeoTransform())##sets same geotransform as input
        outdata.SetProjection(ds.GetProjection())##sets same projection as input
        outdata.GetRasterBand(1).WriteArray(arr_out)
        outdata.GetRasterBand(1).SetNoDataValue(10000)##if you want these values transparent
        outdata.FlushCache() ##saves to disk!!
        outdata = None
        band=None
        ds=None

Answer 1

1. You're creating a single band output raster every iteration of the inner loop, move the following out of the inner loop (and change the band count to len(output)):

   driver = gdal.GetDriverByName("GTiff")
   outdata = driver.Create(filename, cols, rows, len(output), gdal.GDT_UInt16)
   outdata.SetGeoTransform(ds.GetGeoTransform())##sets same geotransform as input
   outdata.SetProjection(ds.GetProjection())##sets same projection as input
   

2. You're writing the 1st band every loop - outdata.GetRasterBand(1). You need:

   band = outdata.GetRasterBand(i+1)
   

3. This may be incorrect - [cols, rows] = arr.shape, numpy raster arrays are usually row, col (or band, row, col) order, you may need:

   rows, cols = arr.shape
   

   As is this:

   outdata = driver.Create(filename, rows, cols, 1, gdal.GDT_UInt16)
   

   Change it to cols, rows

4. Move the following out of the inner loop:

   outdata.FlushCache() ##saves to disk!!
   outdata = None
   band=None
   ds=None
   

You'll get something like (completely untested):

for data_file in os.listdir(destSAR):
    data = data_file.split('_')
    joinStr = [data[6],data[7],data[8],data[9]]
    joinedName = '_'.join(joinStr)
    imgSAR = joinSAR+joinedName
    imgEO = joinEO+joinedName
    SARt = gdal.Open(imgSAR)
    opref = cv2.imread(imgEO)
    SAR=SARt.ReadAsArray()

    output=getSARmask(SAR,opref) #processing the SAR images
    os.chdir(outputPath) 
    filename = 'post_process_'+joinedName

    driver = gdal.GetDriverByName("GTiff")
    outdata = driver.Create(filename, cols, rows, len(output), gdal.GDT_UInt16)
    outdata.SetGeoTransform(ds.GetGeoTransform())##sets same geotransform as input 
    outdata.SetProjection(ds.GetProjection())##sets same projection as 

    for i in range(len(output)):
        ds = SARt #original SAR image with 4 bands, using for Geo information
        arr = output[i]
        rows, cols = arr.shape
        arr_min = arr.min()
        arr_max = arr.max()
        arr_mean = int(arr.mean())
        arr_out = np.where((arr < arr_mean), 10000, arr)

        band = outdata.GetRasterBand(i+1)
        band.WriteArray(arr_out)
        band.SetNoDataValue(10000)##if you want these values transparent

    # outdata.FlushCache()  # Not required!!!
    band = None
    outdata = None
    ds=None  ##saves to disk!!