用Arcpy批量进行植被覆盖度计算

近期准备计算某地区的植被覆盖度，所用数据为MODIS13A3，前期经过了MRT实现批量投影、裁剪及格式转换，并利用arcpy进行了异常值处理及浮点转换，现在为标准tif格式数据。由于暂时没找到好的办法可以像ENVI一样进行置信度区间取值，所以打算采用经验值计算植被覆盖度，NDVIveg取值0.7，NDVIsoil取值0.05，当NDVI大于0.7时，植被覆盖度为1，当NDVI小于0.05时，植被覆盖度为0，中间时采用二分法公式计算：

vfc = [(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)]

 

事先在记事本中写好代码，arcgis分段计算利用Con函数实现，具体代码如下：

import sys
reload(sys)
sys.setdefaultencoding('utf-8')
import os
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.env.workspace = "E:/MODISNDVI/ceshi"
rasterlist = arcpy.ListRasters("*","tif")
output_path = "E:/MODISNDVI/BYVFC/"
for raster in rasterlist:
   print str(raster)
   (filepath, fullname) = os.path.split(raster)
   (prename, suffix) = os.path.splitext(fullname)
   #Raster(raster)读取已有的栅格数据来建立栅格对象
   inRaster = arcpy.Raster(raster)
   outCon = Con(inRaster < 0.05,0, Con((inRaster >= 0.05) & (inRaster  <=  0.7), (inRaster-0.05)/(0.7-0.05), 1))

   outCon.save(output_path + prename + '_vfc.tif')
print("OK!")

 

打开Arcgis，点击工具箱旁边的小图标，进入python界面，将写好的代码复制到里面，按enter键运行代码

 

PS：如果有会用arcgis实现置信去见取值的大神，请指教我一下哈，IDL对于我而言实在是太晕了