gdal.Translate deletes GTiff instead of overwriting?

Question

Something is wrong with the following attempt to add a band to a raster.

I have a GTiff in which I stored numerical values from a field of an ESRI Shapefile. The shapefiles have another field which I want to add as a second band of the same GTiff, if need be. To this end, I defined the following function:

def add_wband(raster, shapefile):
    '''
    adds a new raster band with either bad data areas (add = "bad_data"),
    polygon boundaries (add = "boundaries"), or both (add = "both")
    '''
    # create new temporary raster
    src_ds = gdal.Open(raster, gdalconst.GA_ReadOnly)
    xres, yres, geo_tf = src_ds.RasterXSize, src_ds.RasterYSize, \
        src_ds.GetGeoTransform()
    tmp = gdal.GetDriverByName("GTiff"). \
        Create("/vsimem/tmp.tif", xres, yres, n_bands, gdal.GDT_Byte)
    tmp.SetGeoTransform(geo_tf)
    tmp.SetProjection(src_ds.GetProjection())
    src_ds = None
    tmp.GetRasterBand(1).Fill(1)
    tmp.FlushCache()
    # write values to temporary raster
    drv = ogr.GetDriverByName("ESRI Shapefile")
    shp = drv.Open(shapefile)
    layer = shp.GetLayer()
    band1 = tmp.GetRasterBand(1)
    gdal.RasterizeLayer(tmp, [1], layer, \
                            options = ["ATTRIBUTE=Bad_data"])
    band1.FlushCache()
    band1 = None
    # stack new layers to original raster
    tifs = [raster, "/vsimem/tmp.tif"]
    virtual_stack = gdal.BuildVRT("/vsimem/stacked.vrt", tifs, separate=True)
    out_ds = gdal.Translate(raster, virtual_stack, format = "GTiff")
    del out_ds

I would expect that the GTiff is exchanged for a new one, which contains the additional band from the virtual layer which should have the same properties and contain the values from the "Bad_data" field of the shapefile.

However, when I run

add_wband(raster = 'G:\\test\\B1_6_0003_MASK.tif', \
          shapefile = 'G:\\test\\B1_6_0003\\B1_6_0003_UTM.shp')

there is no additional band but the entire GTiff file just disappears.

There are no error messages, so I have no clue how this happens. When I run the lines of the function separately, the disappearing happens at the last step, when I run gdal.Translate.

How can I fix this?

Answer 1

It might be that the core of your problem is that you are translating a file onto itself, in this case your raster variable.

out_ds = gdal.Translate(raster, virtual_stack, format = "GTiff")

See the thread here. Its related to gdal warp, but I think it shares a lot of the internals with gdal translate. https://github.com/OSGeo/gdal/issues/2867

I would suggest you create a new file as your final output (not raster) and then do some renaming and deleting with python's os library if you only want to keep 1 file.

And if you want to see some error messages form gdal, you can turn on the gdal exceptions module with gdal.UseExceptions()

Let me know if this works, otherwise I can provide other ways.