Python遥感图像处理应用篇：Arcpy遥感图像EVI指数计算批量处理

Posted 2021-05-17 空中旋转篮球

1.MOD09A1数据

采用MRT工具对MOD09A1数据进行初步处理，得到8day7波段数据，数据处理如下：

2.EVI指数计算公式

evi指数计算公式采用近红外、红光、蓝光三个波段进行计算，公式如下：

MODIS数据波段介绍如下：

MODIS sensor has 36 spectral bands, seven of which are designed for the study of vegetation and landsurfaces: blue (459–479 nm), green (545–565 nm), red(620–670 nm), near infrared (NIR1: 841–875 nm; NIR2:1230–1250 nm), and shortwave infrared (SWIR1: 1628–1652 nm, SWIR2: 2105–2155 nm). Daily global imageryis provided at spatial resolutions of 250-m (red andNIR1) and 500-m (blue, green, NIR2, SWIR1, SWIR2)
这里选择第1、2波段和第3波段进行计算。

3.Modis数据整理

把MOD09A1处理出来的各波段数据分别按照波段对应放在下面3个文件夹里面：

4.Arcpy进行EVI指数计算代码

代码如下：

# -*- coding: utf-8 -*-
import os,arcpy,datetime
from arcpy.sa import *
#存放文件夹的位置
Blue_path = r'H:\\NDVITest03'#蓝色波段存放文件夹
Red_path = r'H:\\NDVITest01'#红色波段存放文件夹
NIR_path = r'H:\\NDVITest02'#近红外波段存放文件夹
path =r'H:\\NDVITestOut'#输出结果存放文件夹

# 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 "-----------------------------------------------------------"
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
    # Caculate EVI
    # Check out the ArcGIS 3D Analyst extension license
    arcpy.CheckOutExtension("3D")
    NIRraster = NIR_path +'\\\\'+raster[:-5]+"2.tif"
    Blueraster = Blue_path +'\\\\'+raster[:-5]+"3.tif"
    EVIoutPath = resultDataPath+'\\\\'+ raster[:-24]+ "_EVI.tif"
    print(NIRraster)
    # Converted to floating-point data
    arcpy.Float_3d(Red_path+'\\\\'+ raster, "floatRedBand.tif")
    arcpy.Float_3d(NIRraster, "floatNIRBand.tif")
    arcpy.Float_3d(Blueraster, "BlueBand.tif")
    #后续操作在workspace中进行不需要设置具体路径
    #计算EVI分子
    arcpy.Minus_3d("floatNIRBand.tif", "floatRedBand.tif", "outminus.tif")
    arcpy.Times_3d("outminus.tif", 2.5, "Fenzi.tif")
    #计算EVI分母
    arcpy.Times_3d("floatRedBand.tif", 6, "RedTimes6.tif")
    arcpy.Times_3d("BlueBand.tif", 7.5, "BlueTimes75.tif")
    arcpy.Plus_3d("floatNIRBand.tif", "RedTimes6.tif", "NIRPRed.tif")
    arcpy.Minus_3d("NIRPRed.tif", "BlueTimes75.tif", "NPRPB.tif")
    arcpy.Plus_3d("NPRPB.tif", 10000, "Fenmu.tif")
    #计算EVI
    arcpy.Divide_3d("Fenzi.tif", "Fenmu.tif", EVIoutPath)#计算EVI
    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)

print "Finish!"