用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對於我而言實在是太暈了