利用Arcpy模块进行干旱严重度指数DSI的批量计算（栅格批量计算）

        进行干旱研究的朋友可能会经常使用到干旱严重度指数（DSI）的计算，在arcmap中我们可以利用栅格计算器轻松实现，但当文件量较多时，栅格计算器就显得力不从心，这时候我们便可以使用arcpy代码进行批量计算,此外，对此代码根据相应公式修改便可以计算其他指数。以下是代码

# coding:utf-8

# 程序说明：使用ARCGIS自带的python2.7编译器。

# 本程序对文件夹中的文件名有要求，需要将文件名统一成0，1，2的递进数列形式，方便识别文件。

import os

import arcpy

from arcpy import env

from arcpy.sa import *

import string

arcpy.CheckOutExtension("spatial")

arcpy.gp.overwriteOutput = 1

arcpy.env.overwriteOutput = 1

#arcpy.CheckOutExtension("spatial")

#arcpy.gp.overwriteOutput = 1

fir_Path = "**********"  # NDVI数据

sec_Path = "**********"  # ET数据

thi_Path = "**********"# PET数据

OutPath = "**********"  # 输出结果1路径

n = 1  # 起始编号

for i in range(0, 12):  # 总文件数量。我这边示例计算12个文件

    OutName = OutPath + "DSI_" + str(i+1) + ".tif"
    NDVI = fir_Path + str(n) + ".tif"  # ndvi描述
    ET = sec_Path + str(n) + ".tif"  # et文件名描述
    PET = thi_Path + str(n) + ".tif"  #pet文件名描述
    ep = arcpy.sa.Float(Raster(ET)/Raster(PET))
    meanValueInfo_NDVI = arcpy.GetRasterProperties_management(NDVI,"MEAN")#获取NDVI的平均值
    mean_NDVI  = float(meanValueInfo_NDVI.getOutput(0))#将获取的mean转换成数值
    STD_NDVIInfo= arcpy.GetRasterProperties_management(NDVI,"STD")#获取NDVI的标准差
    STD_NDVI= float(STD_NDVIInfo.getOutput(0))#将获取的std转换成数值
    mean_EPInfo = arcpy.GetRasterProperties_management(ep,"MEAN")
    mean_EP = float(mean_EPInfo.getOutput(0))
    STD_EPInfo = arcpy.GetRasterProperties_management(ep,"STD")
    STD_EP = float(STD_EPInfo.getOutput(0))
    Z_NDVI = arcpy.sa.Float((Raster(NDVI)-float(mean_NDVI))/ float(STD_NDVI))
    Z_EP = arcpy.sa.Float(ep- float(STD_NDVI)) / float(STD_EP)
    Z = arcpy.sa.Float(Z_NDVI + Z_EP)
    Z_mean = float(Z_meanInfo.getOutput(0))
    Z_STDInfo = arcpy.GetRasterProperties_management(Z,"STD")
    Z_STD = float(Z_STDInfo.getOutput(0))
    DSI = arcpy.sa.Float((Z - Z_mean) / Z_STD)
    DSI_save(str(OutName))
    print(str(i+1) + "  has done")
    n = n + 1
    i = i + 1

print("all over")