5

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.

Improve this question
1
  • Not an expert on GDAL but it appears that CreateCopy is not meant for transformation of the "metadata" of a dataset, e.g. the number and order of bands can be considered part of the metadata of the file, which CreateCopy intends to simplify for you. So yes, I think you must use Create instead, which is a bit more involved. The source to gdal_merge.py could be useful to look at. Commented Oct 25, 2013 at 22:47

1 Answer 1

Reset to default
5

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
Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.