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,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.
  [1]: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html