Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
A script to process TM data to calculate NDVI
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 3 11:56:48 2011
The main program. The program takes 3 files as input:
#. The red band (any GDAL-compatible format will do)
#. the near-infrarred band (again as above)
#. An output file name
dditionally, one can specify the output format. By default, we use GeoTIFF.
@author: Jose Gómez-Dans (UCL/NCEO) - j.gomez-dans@ucl.ac.uk
"""
import os
import sys
import optparse
import numpy as np
from osgeo import gdal
def calculate_ndvi ( red_filename, nir_filename ):
"""
A function to calculate the Normalised Difference Vegetation Index
from red and near infrarred reflectances. The reflectance data ought to
be present on two different files, specified by the varaibles
`red_filename` and `nir_filename`. The file format ought to be
recognised by GDAL
"""
g_red = gdal.Open ( red_filename )
red = g_red.ReadAsArray()
g_nir = gdal.Open ( nir_filename )
nir = g_nir.ReadAsArray()
if ( g_red.RasterXSize != g_nir.RasterXSize ) or \
( g_red.RasterYSize != g_nir.RasterYSize ):
print "ERROR: Input datasets do't match!"
print "\t Red data shape is %dx%d" % ( red.shape )
print "\t NIR data shape is %dx%d" % ( nir.shape )
sys.exit ( -1 )
passer = np.logical_and ( red > 1, nir > 1 )
ndvi = np.where ( passer, (1.*nir - 1.*red ) / ( 1.*nir + 1.*red ), -999 )
return ndvi
def save_raster ( output_name, raster_data, dataset, driver="GTiff" ):
"""
A function to save a 1-band raster using GDAL to the file indicated
by ``output_name``. It requires a GDAL-accesible dataset to collect
the projection and geotransform.
"""
# Open the reference dataset
g_input = gdal.Open ( dataset )
# Get the Geotransform vector
geo_transform = g_input.GetGeoTransform ()
x_size = g_input.RasterXSize # Raster xsize
y_size = g_input.RasterYSize # Raster ysize
srs = g_input.GetProjectionRef () # Projection
# Need a driver object. By default, we use GeoTIFF
if driver == "GTiff":
driver = gdal.GetDriverByName ( driver )
dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
gdal.GDT_Float32, ['TFW=YES', \
'COMPRESS=LZW', 'TILED=YES'] )
else:
driver = gdal.GetDriverByName ( driver )
dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
gdal.GDT_Float32 )
dataset_out.SetGeoTransform ( geo_transform )
dataset_out.SetProjection ( srs )
dataset_out.GetRasterBand ( 1 ).WriteArray ( \
raster_data.astype(np.float32) )
dataset_out.GetRasterBand ( 1 ).SetNoDataValue ( float(-999) )
dataset_out = None
if __name__ == "__main__":
arg_parser = optparse.OptionParser()
arg_parser.add_option( '-r', '--red', dest="red_fname", \
help="The RED data" )
arg_parser.add_option( '-n', '--nir', dest="nir_fname", \
help="The NIR data" )
arg_parser.add_option( '-o', '--output', dest="out_fname", \
help="The output dataset" )
arg_parser.add_option( '-f', '--format', dest="out_format", \
default="GTiff", help="Output format" )
options, extra_junk = arg_parser.parse_args ()
if not os.path.exists ( options.red_fname ):
print "ERROR: The red filename %s does not exist" % options.red_fname
sys.exit ( -1 )
if not os.path.exists ( options.nir_fname ):
print "ERROR: The nir filename %s does not exist" % options.nir_fname
sys.exit ( -1 )
if os.path.exists ( options.out_fname):
print "ERROR: The output filename %s does already exist" % \
options.out_fname
print "\t Select a different one, or delete the file."
sys.exit ( -1 )
c_ndvi = calculate_ndvi ( options.red_fname, options.nir_fname )
save_raster ( options.out_fname, c_ndvi, options.red_fname, \
driver=options.out_format )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.