0

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?

1 Answer 1

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.

1
  • Thanks, I noticed the issue only occurs when I write the file onto itself. I hoped to be able to simply update existing files instead of creating new ones, deleting the old ones and maybe renaming the new ones. But, alas, I can do that. At least it is just a cosmetic problem. Nov 4, 2021 at 19:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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