# script to merge several individual Tiffs (bands) into one multiband GeoTiff
# and potentially do calculations like NDVI, albedo, etc...

# import the GDAL and numpy libraries
from osgeo import gdal
from numpy import *

# ***************************************************************
# ok lets load in the first 4 bands of Landsat imagery into their own numpy arrays
# the numpy arrays are named band1, band2, etc.
# I'm using the Boulder images that I provided you, but -
# alternatively you could change these file names to be Landast imagery that you downloaded

g = gdal.Open("boulderBand1.tif") # blue band
band1 = g.ReadAsArray()
g = gdal.Open("boulderBand2.tif") # green band
band2 = g.ReadAsArray()
g = gdal.Open("boulderBand3.tif") # red band
band3 = g.ReadAsArray()
g = gdal.Open("boulderBand4.tif") # near infrared band
band4 = g.ReadAsArray()

# ****************************************************************

# now we could do math on these arrays using numpy if we wanted to - yay!

# ****************************************************************

# these variables will get information about the input Tiff so we can
# write out our new Tiff into the correct geographic space and with correct row/column dimensions

geo = g.GetGeoTransform()  # get the datum
proj = g.GetProjection()   # get the projection
shape = band1.shape        # get the image dimensions - format (row, col)

# ****************************************************************

# here we write out the new image, as a multiband Tiff

driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.Create( "boulderMulti.tif", shape[1], shape[0], 4, gdal.GDT_Byte)
                                                         # here we set the variable dst_ds with 
                                                         # destination filename, number of columns and rows
                                                         # 4 is the number of bands we will write out
                                                         # gdal.GDT_Byte is the data type - integers
dst_ds.SetGeoTransform( geo ) # set the datum
dst_ds.SetProjection( proj )  # set the projection

# probably could loop this so we wouldn't have to have so many lines...can you imagine how??

dst_ds.GetRasterBand(1).WriteArray( band1 )  # write numpy array band1 as the first band of the multiTiff - this is the blue band
stat = dst_ds.GetRasterBand(1).GetStatistics(1,1)  # get the band statistics (min, max, mean, standard deviation)
dst_ds.GetRasterBand(1).SetStatistics(stat[0], stat[1], stat[2], stat[3]) # set the stats we just got to the band

dst_ds.GetRasterBand(2).WriteArray( band2 ) # write numpy array band1 as the first band of the multiTiff - this is the green band
stat = dst_ds.GetRasterBand(2).GetStatistics(1,1)  # get the band statistics (min, max, mean, standard deviation)
dst_ds.GetRasterBand(2).SetStatistics(stat[0], stat[1], stat[2], stat[3]) # set the stats we just got to the band

dst_ds.GetRasterBand(3).WriteArray( band3 ) # write numpy array band1 as the first band of the multiTiff - this is the red band
stat = dst_ds.GetRasterBand(3).GetStatistics(1,1)  # get the band statistics (min, max, mean, standard deviation)
dst_ds.GetRasterBand(3).SetStatistics(stat[0], stat[1], stat[2], stat[3]) # set the stats we just got to the band

dst_ds.GetRasterBand(4).WriteArray( band4 ) # write numpy array band1 as the first band of the multiTiff - this nearInfrared
stat = dst_ds.GetRasterBand(4).GetStatistics(1,1)  # get the band statistics (min, max, mean, standard deviation)
dst_ds.GetRasterBand(4).SetStatistics(stat[0], stat[1], stat[2], stat[3]) # set the stats we just got to the band

# that's it, close Python and load your new multiband image into QGIS !!