I'm trying to create a new raster from three bands (rgb). I got the script mostly from [here.][1] 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)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    
    for band in range(bands):
        if band is 0:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
        if band is 1:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_GreenBand)
        if band is 2:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
        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: I edited the script to reflect Lukes recomendations. The file is written but I think there is a bug in GDAL where it doesn't actually write a 16-bit rgb image to the RGD image channels. The output still remains as grey,alpha,alpha.

Gdal.info says that the image is formatted correctly but that is of little consolement. [1]: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html