¿Cómo aplicar los ajustes de la "Banda" utilizando los enlaces gdal Python?

Question

Estoy intentando replicar el comando gdal_translate, pero usando Python puro:

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

Ya he mirado ¿Cómo llamar a gdal_translate desde código Python? pero estoy atascado con la forma de aplicar la información de banda utilizando los enlaces GDAL Python.

Hasta ahora puedo ver cómo hacer todo esto con el código Python, excepto aplicar la configuración de banda (-b).

#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

Después de leer el Tutorial de la API GDAL Parece que sólo puedo hacer esto con el método Create(), no con el método CreateCopy(). Si tengo que utilizar el método Create(), ¿cuál es la forma correcta de hacerlo? Parece que el uso de un VRT es probablemente la forma correcta, pero no puedo entender cómo escribir el código que "haría una copia" del archivo existente mientras se utilizan los métodos Create() o VRT.

Answer 1

Actualización:

A partir de GDAL 2.1, puede utilizar la función gdal.Translate() y pasar las bandas en la función bandList argumento ( lista de argumentos ).

from osgeo import gdal

infn='in.tif'
outfn='out.tif'

gdal.Translate(outfn, infn, bandList=[1,2,3], **other_kwargs)

Respuesta original:

gdal_translate lo hace construyendo un VRT desde cero, copiando todos los metadatos del conjunto de datos de origen (a menos que todo coincida, entonces sólo utiliza CreateCopy).

Puedes hacerlo en python, pero es un poco complicado. Una manera más fácil es utilizar el controlador VRT para crear una copia VRT, luego modificar el VRT XML para eliminar las bandas que no te interesan antes de escribir el VRT modificado en un nuevo archivo. Mira el ejemplo de abajo:

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