Python批處理MODIS_NDVI數據

NDVI是歸一化植被指數，常用於檢測植被生長狀態、植被覆蓋度和消除部分輻射誤差等，范圍為-1到1，負值表示地面覆蓋為雲、水、雪等；0表示有岩石或裸土等；正值，表示有植被覆蓋，且隨覆蓋度增大而增大；

此次所用NDVI數據來源於NASA官網的MOD13A3產品，前期用MRT對所下載的數據進行了幾何校正、拼接、投影和格式轉換，所得NDVI數據為tif格式。GDAL也能用於處理柵格數據，但是自己正好看到arcpy所以使用該包進行嘗試，效果很不錯。

自己電腦正常使用的是python3，但是arcgis支持的是python2，因此本文使用安裝arcgis時自帶的python2進行計算，首先對C:\Python27\ArcGIS10.6中的python.exe和pythonw.exe重命名為python2.7.exe和pythonw2.7.exe，並配置環境變量，完成后在cmd中測試，輸入python2.7，顯示python版本，則證明python2可用。如下圖。

 

有了可用的python2環境后，編寫代碼：

# -*- coding: utf-8 -*-
import os
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.env.workspace = "e:/MODISNDVI/BYNDVI"
rasters = arcpy.ListRasters("*","tif")
out_path = "e:/MODISNDVI/BYNDVI_cl/"
for raster in rasters:
    
    (filepath, fullname) = os.path.split(raster)
    (prename, suffix) = os.path.splitext(fullname)
    print(prename)
    arcpy.CheckOutExtension("ImageAnalyst") #檢查許可
    arcpy.CheckOutExtension("spatial") #檢查許可
    whereClause = "VALUE = -3000" #無效值
    outSetNull = SetNull(raster, raster, whereClause) * 0.0001 #去除無效值並乘以0.0001
    #outname=r"E:\MODISNDVI\BYNDVI\try1.tif" #輸出路徑
    outSetNull.save(out_path + prename + '_ndvi_qcwxz.tif') #保存數據
    print('over')

對於處理長時間序列的遙感影像，批量處理是遙感數據處理家常便飯的事。需要注意的是，批量輸出文件名要重命名，因此引用了OS包。程序中的數據是整型，需要除以10000，轉化為大氣表觀反射率，再進行公式計算。

嘗試直接使用pycharm，但是在setting中配置了python2.7.exe后沒有顯示存在arcpy包，因此退而求其次，使用cmd運行腳本，在cmd中進入腳本所在文件夾下，輸入python2.7 xxxx.py即可。

 

 去除無效值后的數據。

后期會參考https://blog.csdn.net/summer_dew/article/details/78368184的文章，計算生長季NDVI均值。