Statistical Processing of EVI and NDVI from MODIS in ArrcGIS 10 using Python

OK, ArcGIS got much better with version 10. I did not even mind having to re-write all of my scripts, since now, Python is really integrated. The only thing worrying me now is the performance and some file name issues (see below).

With this script, I can batch-process a number of Nasa HDF files. First, the data of the NDVI (Normalized Difference Vegetation Index) and the EVI (Enhanced Vegetation Index) band is extracted. Then we re-scale it (for storage, the data was multiplied by 1E4).
Finally, I create a zonal statistic for a 320 km buffer around Zurich.

	# This script processes the EVI/NDVI data for Francesco.
	#
	################################################################################
	# $Id: evi.py 6723 2011-08-08 12:58:26Z oderbolz $
	################################################################################
	# Import stuff #################################################################
	import os
	from glob import glob
	import sys
	import string
	import arcpy
	from arcpy.sa import *
	from arcpy import env
	# Check out the ArcGIS Spatial Analyst extension license
	arcpy.CheckOutExtension("Spatial")
	################################################################################
	inputPath='Y:\\oderbolz\\SatelliteData\\MODIS-EVI\\2011'
	inputPattern='*.hdf'
	outputPath='T:\\Projects\\OPS\\FC\\EVI\\output'
	# Buffer to use
	bufferFile='T:\\Projects\\OPS\\FC\\EVI\\buffer\\Zurich_Buffer.shp'
	################################################################################
	print("Storing output in " + outputPath)
	env.workspace = outputPath
	fileList=glob(os.path.join(inputPath,inputPattern))
	for fname in fileList:
	# Split file name
	currentInputFile=os.path.basename(fname)
	filePrefix, extension = os.path.splitext(currentInputFile)
	# Prefix must not contain all these dots (ArcGIS gets confused because it determined datatype by extension)
	filePrefix = string.replace(filePrefix,'.','_')
	print("Processing file " + currentInputFile)
	# Load data
	try:
	arcpy.ExtractSubDataset_management(fname,os.path.join(outputPath,filePrefix + '_ndvi_raw.tif'), "0")
	arcpy.ExtractSubDataset_management(fname,os.path.join(outputPath,filePrefix + '_evi_raw.tif'), "1")
	except:
	print "Extract Subdataset failed."
	print arcpy.GetMessages()
	sys.exit(1)
	print("Data loaded, converting back to NDVI and EVI...")
	# Convert back to EVI and NDVI
	# arcpy.env.mask = area #Extract area of interest
	ndvi = Raster(os.path.join(outputPath,filePrefix + '_ndvi_raw.tif')) * .0001 #Multiply by scale factor
	evi = Raster(os.path.join(outputPath,filePrefix + '_evi_raw.tif')) * .0001 #Multiply by scale factor
	ndvi.save(os.path.join(outputPath,filePrefix + '_ndvi.tif'))
	evi.save(os.path.join(outputPath,filePrefix + '_evi.tif'))
	try:
	print("Calculating mean in region of interest...")
	# Zonal statistics with buffer
	outZonalStatsNdvi = ZonalStatistics(bufferFile,"Name", os.path.join(outputPath,filePrefix + '_ndvi.tif'),"MEAN", "DATA")
	outZonalStatsEvi = ZonalStatistics(bufferFile, "Name", os.path.join(outputPath,filePrefix + '_evi.tif'), "MEAN", "DATA")
	except:
	print "Zonal statistics failed."
	print arcpy.GetMessages()
	sys.exit(1)
	# Save data
	outZonalStatsNdvi.save(os.path.join(outputPath,filePrefix + '_ndvi_stats.tif'))
	outZonalStatsEvi.save(os.path.join(outputPath,filePrefix + '_evi_stats.tif'))

view raw evi.py hosted with ❤ by GitHub

You notice that I replace all dots by underscores in line 48. This is needed because ArcGis otherwise has trouble identifying the datatype (which is done by extension).

Have fun!