This might help get you started with gdal....Based on code I found here: http://www.gis.usu.edu/~chrisg/python/
'''
---------------------------------------------------------------------------
---------------------------------------------------------------------------
program: TOA to NDVI.py
Chris Bater
GIS and Remote Sensing Technologist
Wildfire Resource Information Unit
Wildfire Management Branch
Forestry Division
Alberta Sustainable Resource Development
9th floor Great West Life Building
9920 - 108 Street
Edmonton, AB
T5K 2M4
Ph: 780-422-0669
Cell: 780-446-1899
Email: chris.bater@gov.ab.ca
Date: September 2011
Purpose: Converts Landsat-5 TOA reflectance data to NDVI
---------------------------------------------------------------------------
----------------------------------------------------------------------------
'''
#----------------------------------------------------------------------
#Identify the directory containing the landsat data
workspace = r"C:\_LOCALdata\temp"
#----------------------------------------------------------------------
#----------------------------------------------------------------------
#Import and load modules
#----------------------------------------------------------------------
# Import system modules
import os, sys
from osgeo import gdal
import numpy as np
import numpy.ma as ma
#register the GDAL drivers
gdal.AllRegister()
driver = gdal.GetDriverByName('GTiff') #not sure if this is necessary
driver.Register()#not sure if this is necessary
print ""
print "All systems go"
print "Hit CTRL + c at any time to stop processing"
print ""
#----------------------------------------------------------------------
#Iterate through the landsat scenes and convert to TOA
#----------------------------------------------------------------------
try:
files = os.listdir(workspace)
for f in files:
if not f.endswith('tif'):
continue
name = f[:-4]
#open the landsat TOA file
TOA = gdal.Open(workspace + "/" + f)
#get the image size
rows = TOA.RasterYSize
cols = TOA.RasterXSize
bands = TOA.RasterCount
if bands != 6:
continue
print "processing " + f
#create an empty tif file
outData = driver.Create(os.path.join(workspace, name + "_NDVI.tif"), cols, rows, 1, gdal.GDT_Float32)
outBand = outData.GetRasterBand(1)
#assign the bands to variables
b3 = TOA.GetRasterBand(3)
b4 = TOA.GetRasterBand(4)
#set the initial blocksize for analysis
xBlockSize = 5000
yBlockSize = 5000
for i in range(0, rows, yBlockSize):
if i + yBlockSize < rows:
numRows = yBlockSize
else:
numRows = rows - i
for j in range(0, cols, xBlockSize):
if j + xBlockSize < cols:
numCols = xBlockSize
else:
numCols = cols - j
#read bands into arrays
d3 = b3.ReadAsArray(j, i, numCols, numRows).astype(np.float32)
d4 = b4.ReadAsArray(j, i, numCols, numRows).astype(np.float32)
#mask no data values (no data = 0)
m3 = ma.masked_equal(d3, 0)
m4 = ma.masked_equal(d4, 0)
del d3, d4
ndvi = (m4 - m3) /(m4 + m3)
del m3, m4
outBand.WriteArray(ndvi, j, i)
# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
outBand.SetNoDataValue(0)
stats = outBand.GetStatistics(0, 1)
# georeference the image and set the projection
outData.SetGeoTransform(TOA.GetGeoTransform())
outData.SetProjection(TOA.GetProjection())
#memory management
del TOA, outData, outBand
print 'script complete'
except Exception, e:
# If an error occurred, print line number and error message
tb = sys.exc_info()[2]
print "Line %i" % tb.tb_lineno
print e.message
