How to apply "Band" settings using gdal Python bindings?

Question

I am trying to replicate the gdal_translate command, but using pure Python:

gdal_translate infile.tif outfile.tif -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co TFW=YES -b 1 -b 2 -b 3

I have already looked at How to call gdal_translate from Python code? but am stuck with how to apply band information using the GDAL Python bindings.

So far I can see how to do all of this with the Python code, except apply the band (-b) settings.

#Import gdal
from osgeo import gdal

src_filename = r"infile.tif"
dst_filename = r"outfile.tif"

#Open existing dataset
src_ds = gdal.Open( src_filename )

#Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )

#Output to new format
dst_ds = driver.CreateCopy( dst_filename, src_ds, 0, ['COMPRESS=JPEG','PHOTOMETRIC=YCBCR','TFW=YES'])

#Properly close the datasets to flush to disk
dst_ds = None
src_ds = None

After reading the GDAL API Tutorial, it appears that I can only do this with the Create() method, not the CreateCopy() method. If I have to go the route of using the Create() method, what is the right approach to doing this? It looks like using a VRT is likely the right way, but I cannot wrap my head around how to write the code that would "make a copy" of the existing file while using either the Create() or VRT methods.

Answer 1

Update:

As of GDAL 2.1, you can use the gdal.Translate() function and pass the bands in the bandList argument (list of arguments).

from osgeo import gdal
infn='in.tif'
outfn='out.tif'
gdal.Translate(outfn, infn, bandList=[1,2,3], **other_kwargs)

Original answer:

gdal_translate does this by building a VRT from scratch, copying all the metadata from the source dataset (unless everything matches, then it just uses CreateCopy).

You can do this in python but it's a bit fiddly. An easier way is use the VRT driver to create a VRT copy, then modify the VRT XML to remove the bands you aren't interested in before writing the modified VRT out to new file. See example below:

from osgeo import gdal
from lxml import etree
import os,sys
def read_vsimem(fn):
    '''Read GDAL vsimem files'''
    vsifile = gdal.VSIFOpenL(fn,'r')
    gdal.VSIFSeekL(vsifile, 0, 2)
    vsileng = gdal.VSIFTellL(vsifile)
    gdal.VSIFSeekL(vsifile, 0, 0)
    return gdal.VSIFReadL(1, vsileng, vsifile)
def write_vsimem(fn,data):
    '''Write GDAL vsimem files'''
    vsifile = gdal.VSIFOpenL(fn,'w')
    size = len(data)
    gdal.VSIFWriteL(data, 1, size, vsifile)
    return gdal.VSIFCloseL(vsifile)
infn='test.tif'
outfn='test2.tif'
#List of bands to retain, output bands will be reordered to match
bands=[4,3,2]
ds = gdal.Open(os.path.abspath(infn)) #Ensure path is absolute not relative path
vrtdrv=gdal.GetDriverByName('VRT')
tifdrv=gdal.GetDriverByName('GTIFF')
#Create an in memory VRT copy of the input raster
vfn='/vsimem/test.vrt'
vrtds=vrtdrv.CreateCopy(vfn,ds)
#Read the XML from memory and parse it
#Could also write the VRT to a temp file on disk instead of /vsimem
#and used etree.parse(vrtfilepath) instead
#i.e.
#vfn='some/path/on/disk/blah.vrt'
#vrtds=vrtdrv.CreateCopy(vfn,ds)
#tree=etree.parse(vfn)
#os.unlink(vfn)
vrtds=vrtdrv.CreateCopy(vfn,ds)
vrtxml=read_vsimem(vfn)
tree=etree.fromstring(vrtxml)
#Ignore bands not in band list
#And put bands we want to keep in the correct order
for band in tree.findall('VRTRasterBand'):
    try:
        i=bands.index(int(band.attrib['band']))
        band.attrib['band']=str(i+1)
        bands[i]=band
    except ValueError:pass
    finally:
        tree.remove(band)
for band in bands:
    tree.append(band)
#Get the XML string from the tree
vrtxml=etree.tostring(tree, pretty_print=True)
#Write it to the in memory file
#Could also just write to a file on disk
write_vsimem(vfn,vrtxml)
#Open the VRT containing only the selected bands
vrtds=gdal.Open(vfn)
#Write to a new raster
tifds=tifdrv.CreateCopy(outfn,vrtds, 0, ['COMPRESS=JPEG','PHOTOMETRIC=YCBCR','TFW=YES'])
#Close dataset to ensure cache is properly flushed
del tifds