[gdal-dev] python CreateCopy just one band

[Even Rouault]

Le Wednesday 17 June 2009 23:17:42 Gong, Shawn (Contractor), vous avez écrit :
> hi list,
>
> I use these codes for CreateCopy()
>     src_ds = gdal.Open( src_filename )
>     dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 )
>
>
> If the src_ds has two bands and my dst_ds has only one band, is there a
> way to CreateCopy just one band?

No, not directly. Not even in C/C++ GDAL

> Or is there a way to remove bands in 
> Gdal Python?

No, not even in C/C++ GDAL

>
>

But... you could create a VRT file that just references one band of your 
source raster and use CreateCopy from it. Below I've written a sample script 
that should do what you want to do. The vrt_xml has been deduced by observing 
the output of "gdal_translate -of VRT -b 1 rgbsmall.tif rgbsmall_red.vrt" and 
generalizing it a bit. It does the basic stuff (copy georeferencing and pixel 
values), but it should be improved if you want to copy metadata, color table, 
etc

#!/usr/bin/python

import gdal

src_filename = "rgbsmall.tif"
dst_filename = "rgbsmall_red.tif"
band_list = [1]

src_ds = gdal.Open( src_filename )
xsize = src_ds.RasterXSize
ysize = src_ds.RasterYSize
gt = src_ds.GetGeoTransform()
srs = src_ds.GetProjectionRef()

vrt_xml = '''<VRTDataset rasterXSize="%d" rasterYSize="%d">
    <SRS>%s</SRS>
    <GeoTransform>%f, %f, %f, %f, %f, %f</GeoTransform>''' % \
     (xsize, ysize, gdal.EscapeString(srs, scheme = gdal.CPLES_XML), \
      gt[0], gt[1], gt[2], gt[3], gt[4], gt[5])

dst_band = 1
for src_band in band_list:
    data_type = src_ds.GetRasterBand(src_band).DataType
    data_type_name = gdal.GetDataTypeName(data_type)

    vrt_xml = vrt_xml + '''
    <VRTRasterBand dataType="%s" band="%d">
        <SimpleSource>
        <SourceFilename relativeToVRT="1">%s</SourceFilename>
        <SourceBand>%d</SourceBand>
        <SrcRect xOff="0" yOff="0" xSize="%d" ySize="%d"/>
        <DstRect xOff="0" yOff="0" xSize="%d" ySize="%d"/>
        </SimpleSource>
    </VRTRasterBand>''' % \
    (data_type_name, dst_band, src_filename, src_band, \
     xsize, ysize, xsize, ysize)

    dst_band = dst_band + 1

vrt_xml = vrt_xml + '''
</VRTDataset>'''

vrt_ds = gdal.Open(vrt_xml)
driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.CreateCopy( dst_filename, vrt_ds )

> thanks,
> Shawn