[gdal-dev] Clipping an image via Python API problems

[Graeme Merrall]

The smallest almost working code chunk I have right now is below.  This is
taaking a 40k by 40k and pulling out the top left quarter.  The original
has 4 bands - RGBA

from osgeo import gdal
import numpy
img = gdal.Open("original.tif")
geotransform = img.GetGeoTransform()
band1 = img.GetRasterBand(1)

driver = gdal.GetDriverByName("GTiff)
# coordinates will be calculated but hardcoded for now
dst = driver.Create('sliver.tif', 20000, 20000, 1, gdal.GDT_Byte)
dst.SetProjection(img.GetProjection())
scanline1 = band1.ReadAsArray(0, 0,  20000, 20000)
dst.GetRasterBand(1).WriteArray(scanline1)

Should I be looping through each band and reading/writing to the image? I'm
assuming yes

Also. band 1 of the original image is red but in the new image the band
color is grey. I've tried a few variations with
GDALRasterBand.SetColorInterpretation() with no luck.  What have I missed?

Cheers,
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20111128/1fe62e5d/attachment.html