How to create an rgb GeoTIFF file raster from bands using the GDAL python module?

Question

I'm trying to create a new raster from three bands (rgb). I got the script mostly from here. I can create a single band raster with no problems, but I'm having trouble inputting a tree band numpy array. I'm greeted with this error message when I attempt it.

ValueError: expected array of dim 2

If I try to create a three band GeoTIFF image I end up with a grey,alpha,alpha image instead of an rgb image, nor am I able to populate these with my three-dimensional numpy array. The following script works for creating a greyscale raster. with three bands as grey,alpha,alpha.

from osgeo import gdal, osr
import os
import io
import numpy as np

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,format,spatialRef):

    bands = array.shape[0]
    cols = array.shape[1]
    rows = array.shape[2]

    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName(format)
    outRaster = driver.Create(newRasterfn, cols, rows, bands,gdal.GDT_UInt16,options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF',])
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

    for band in range(bands):
        outRaster.GetRasterBand(band+1).WriteArray(array[band,:,:])
        outRaster.GetRasterBand(band+1).FlushCache()
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(spatialRef)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())


def get_rgb_bands():
    source = gdal.Open("C:\\RGBtest.tif")

    red = source.GetRasterBand(1).ReadAsArray()
    green = source.GetRasterBand(2).ReadAsArray()
    blue = source.GetRasterBand(3).ReadAsArray()

    #rgbOutput = source.ReadAsArray() #Easier method
    rgbOutput = np.zeros((3,16,16), 'uint16')

    rgbOutput[0,...] = red
    rgbOutput[1,...] = green
    rgbOutput[2,...] = blue

    #Clear so file isn't locked
    source = None
    return rgbOutput

exampleRGB = get_rgb_bands()

exampleOutput = "C:\\ExampleOutput.tif"
rasterOrigin=[0,0]

array2raster(exampleOutput,rasterOrigin,1, 1,exampleRGB,'GTiff',3879) 

createdFile = gdal.Open("C:\\ExampleOutput.tif")

print 'old file info:'+gdal.Info(source)
print 'new file info:'+gdal.Info(createdFile)

createdFile = None

exampleOutput = None

EDIT: Luke answered this question perfectly I have updated the code to reflect a working example to output a 16-bit rgb raster.

Answer 1

You can set the colour interpretation for most rasters using Band.SetColorInterpretation, i.e

outband.SetColorInterpretation(gdal.GCI_RedBand)

However, this is only understood by GDAL based software. To specifically create a GeoTIFF with an RGB baseline PhotometricInterpretation tag (this is not a GDAL tag, but part of the TIFF spec) pass a PHOTOMETRIC=RGB option to the GeoTIFF driver when creating the raster:

driver = gdal.GetDriverByName('GeoTIFF')
options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']
outRaster = driver.Create(filepath, cols, rows, bands, datatype, options=options)

You're also trying to write 3 bands to a single band:

(outRaster.GetRasterBand(1).WriteArray(etc...))

Get used to having your arrays in [bands, rows, cols] order not [cols, rows, bands] as this is how gdal arranges them when reading into a numpy array. Then you need to loop through the array "bands" and write each of them to the raster bands:

bands = array.shape[0]
for band in range(bands):
    outRaster.GetRasterBand(band+1).WriteArray( array[band, :, :] )

Try something like:

from osgeo import gdal, osr
import os
import io
import numpy as np

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,format='GTiff'):

    bands = array.shape[0]
    rows = array.shape[1]
    cols = array.shape[2]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    driver = gdal.GetDriverByName(format)
    #Here is where I assign three bands to the raster with int 3
    options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']
    outRaster = driver.Create(newRasterfn, cols, rows, bands, gdal.GDT_UInt16, options=options)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    #outband = outRaster.GetRasterBand(1)
    #outband.WriteArray(array)
    for band in range(bands):
        outRaster.GetRasterBand(band+1).WriteArray( array[band, :, :] )

    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(6962)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())

    return outRaster


def example_rgb_creation():
    source = gdal.Open("C:\\RGBtest.tif")
    rgb_array = source.ReadAsArray()

    return rgb_array

exampleRGB = example_rgb_creation()
exampleOutput = "C:\\Temp\\ExampleOutput.tif"
rasterOrigin=[0,0]

outRaster = array2raster(exampleOutput,rasterOrigin,1, 1,exampleRGB)

outRaster.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
outRaster.GetRasterBand(2).SetColorInterpretation(gdal.GCI_GreenBand)
outRaster.GetRasterBand(3).SetColorInterpretation(gdal.GCI_BlueBand)

del outRaster