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
I finally came up with a quite nice result calling GDAL expressions from Python using "Call ()" method and specifying the path with "os.chdir ()" method. Here is a script example:
import os
from subprocess import call
call(["ls", "-l"])
path= 'C:/users/unimi/Desktop/'
os.chdir(path)
cmd= 'gdalwarp -dstnodata 0 -q -cutline \
study_area_foscagno.shp -crop_to_cutline -of GTiff sept2000_B5_refl.tif clip_sept2000_B5_refl.tif'
call (cmd)
The "Call ()" method is here used to do all the GDAL stuff without even importing GDAL in Python ("Call ()" method is even better than "os.system ()" method, see "https://stackoverflow.com/questions/89228/calling-an-external-command-in-python").
"os.chdir ()" method than have let me set the path where the program should look for my data. Hope this could help.
Best Answer
You have now a better way to do.
Since RFC 59.1 : GDAL/OGR utilities as a library, you can use
gdalwarpfrom Python directly without using any call to the command line utility but using really the function from Python.This solution is a bit "on the edge" as you need at the moment to use the latest GDAL version (version 2.1, in fact the master/trunk version)