基于MOD09Q1数据批量计算NDVI

基于MOD09Q1数据批量计算NDVI

通过MRT处理好b01和b02波段后，分别存储至两个不同的文件夹（b01和b02）。

接下来打开在arcgis自带的python2.7中键入以下代码：（如果用pycharm和spyder的朋友需要把环境配置为python2.7并且安装arcpy包）

'''
# -*- coding: utf-8 -*-
import os,arcpy,datetime
from arcpy.sa import *

#存放文件夹的位置
Red_path = r'C:\Users\adminstor\Desktop\b01'#红色波段存放文件夹
NIR_path = r'C:\Users\adminstor\Desktop\b02'#近红外波段存放文件夹
path =r'C:\Users\adminstor\Desktop'#输出结果存放文件夹

# set the intermediate data folder
intermediateDataPath = path+"\\"+"IntermediateData"
# set result data folder
resultDataPath = path+"\\"+"Result"

# determine if the folder exists
if os.path.exists(intermediateDataPath):
    print "IntermediateData floder exists"
else:
    # create a intermediate data floder
    arcpy.CreateFolder_management(path, "IntermediateData")    
if os.path.exists(resultDataPath):
    print "Result floder exists"
else:   
    # create a result floder
    arcpy.CreateFolder_management(path, "Result")

print "Under calculation......"

arcpy.env.workspace = Red_path  #所在文件夹
rasters = arcpy.ListRasters("*", 'tif')  #栅格数据格式设置
for raster in rasters:
    print(raster)
    # set workspace
    arcpy.env.workspace = intermediateDataPath
    arcpy.env.overwriteOutput = True
    # Check out the ArcGIS 3D Analyst extension license
    arcpy.CheckOutExtension("3D")
    NIRraster = NIR_path +'\\'+raster[:-5]+"2.tif"
    NDVIoutPath = resultDataPath+'\\'+ raster[:-16]+ "_NDVI.tif"
    print(NIRraster)
    # Converted to floating-point data
    arcpy.Float_3d(Red_path+'\\'+ raster, "floatRedBand.tif")
    arcpy.Float_3d(NIRraster, "floatNIRBand.tif")
    #计算NDVI分子
    arcpy.Minus_3d("floatNIRBand.tif", "floatRedBand.tif", "outminus.tif")
    #计算NDVI分母
    arcpy.Plus_3d("floatNIRBand.tif", "floatRedBand.tif", "NDVIfenmu.tif")
    #计算NDVI
    arcpy.Divide_3d("outminus.tif", "NDVIfenmu.tif", NDVIoutPath)#计算NDVI
    print(raster+"  has done")  
print("All done")
#清理workspace中的缓存数据        
for i in os.listdir(intermediateDataPath):
    path_file = os.path.join(intermediateDataPath,i)
    if os.path.isfile(path_file):
        os.remove(path_file)

算好了以后就可以拿来批量裁剪后使用了，批量裁剪可以在arcgis中model building或arcpy中进行。