l8scale.py

# dirty l8 scaling in python
# Charlie Loyd, fall 2013
# note use of arbitrary 160000 scale factor in L48

from math import sin
from sys import argv, exit
from osgeo import gdal, gdalconst
import os.path
from numpy import *

infile_path = argv[1]
outfile_path = argv[2]

basename = os.path.basename(infile_path)

if os.path.exists(outfile_path):
	exit('Output file already exists!')

#infile = open(infile_path, 'r')
outfile = open(outfile_path, 'w')

mtl_path = os.path.dirname(infile_path) + '/' + basename[:21] + '_MTL.txt'

mtl = open(mtl_path).readlines()
metadata = {}
for line in mtl:
  line = line.split()
  try:
    metadata[line[0]] = line[2]
  except: pass

band = basename[basename.find('B')+1:basename.find('.')]

Mr = float(metadata['REFLECTANCE_MULT_BAND_' + band])
Ar = float(metadata['REFLECTANCE_ADD_BAND_' + band])
SE = radians(float(metadata['SUN_ELEVATION']))

print Mr, Ar, SE

input_image = gdal.Open(infile_path)
img_x = input_image.RasterXSize
img_y = input_image.RasterYSize
img_proj = input_image.GetProjection()
img_trans = input_image.GetGeoTransform()

DN = input_image.GetRasterBand(1).ReadAsArray().astype(float32)

toa = (((DN * Mr) + Ar) / sin(SE)) * 160000

output_driver = gdal.GetDriverByName("GTiff")
output_image = output_driver.Create(outfile_path, img_x, img_y, 1, gdalconst.GDT_UInt16, options = ['COMPRESS=LZW'])
output_image.SetGeoTransform(img_trans)
output_image.SetProjection(img_proj)
output_image.GetRasterBand(1).WriteArray(toa)