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