# 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 !!