[Gdal-dev] issues merging 7 bands of Landsat scene into 1 image

[Maciej Sieczka]

Hi,

I'm having issues trying to merge 7 bands of Landsat scene into a
multiband image. Please could somebody help me out with this? I have
put a problematic sample (183 KB) on http://kufaya.googlepages.com/tm.zip.

I'd prefer a VRT output, but the following command fails:

$ gdal_merge.py -of VRT -separate -o merge.vrt nn*.tif

ERROR 1: Writing through VRTSourcedRasterBand is not supported.
Traceback (most recent call last):
  File "/usr/local/bin/gdal_merge.py", line 505, in ?
    fi.copy_into( t_fh, 1, t_band, nodata )
  File "/usr/local/bin/gdal_merge.py", line 303, in copy_into
    nodata_arg )
  File "/usr/local/bin/gdal_merge.py", line 138, in raster_copy
    data, t_xsize, t_ysize, t_band.DataType )
  File "/usr/lib/python2.4/site-packages/gdal.py", line 863, in WriteRaster
    buf_string, buf_xsize, buf_ysize,buf_type)
TypeError: Writing through VRTSourcedRasterBand is not supported.

Does that mean VRT output is not possible?

So I tried GeoTIFF:

$ gdal_merge.py -of GTIFF -separate -o merge.tif nn*.tif

which went ok, BUT colortables of bands 2-6 are "Undefined" instead of
"Gray":

$ gdalinfo merge.tif | grep ColorInterp
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray
Band 2 Block=254x32 Type=Byte, ColorInterp=Undefined
Band 3 Block=254x32 Type=Byte, ColorInterp=Undefined
Band 4 Block=254x32 Type=Byte, ColorInterp=Undefined
Band 5 Block=254x32 Type=Byte, ColorInterp=Undefined
Band 6 Block=254x32 Type=Byte, ColorInterp=Undefined
Band 7 Block=254x32 Type=Byte, ColorInterp=Undefined

while the for each input image the ColorInterp is "Gray":

$ for i in `ls | grep nn.\.tif`; do gdalinfo $i | grep ColorInterp; done
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray
Band 1 Block=254x32 Type=Byte, ColorInterp=Gray

In that case I tried the -pct switch, but it makes gdal_merge.py crash:

$ gdal_merge.py -pct -of GTIFF -separate -o merge_pct.tif nn*.tif
Traceback (most recent call last):
  File "/usr/local/bin/gdal_merge.py", line 479, in ?
    t_fh.GetRasterBand(1).SetRasterColorTable(file_infos[0].ct)
  File "/usr/lib/python2.4/site-packages/gdal.py", line 892, in
SetRasterColorTable
    return _gdal.GDALSetRasterColorTable( self._o, ct._o )
AttributeError: 'NoneType' object has no attribute '_o



After some browsing I found this nice script:
http://www.vso.cape.com/~nhv/files/gdal/gdal_vrtmerge.py, but to my bad
luck, it fails as well:

$ gdal_vrtmerge.py -separate -o merge_2nd.vrt nn*.tif
Traceback (most recent call last):
  File "/home/shoofi/bin/gdal_vrtmerge.py", line 268, in ?
    t_fh.write('\t\t<ColorInterp>%s</ColorInterp>\n' %
IndexError: list index out of range

I will be very gratefull for hints on how to output a VRT multiband
image, or a GeoTIFF as a last resource, from such seven 8bit bands,
having ColorInterp set to "Gray" for each.

Best,
Maciek