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.

  • 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. – blah238 Oct 25 '13 at 22:47
2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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