2

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
0
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!!
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.