Python批处理MODIS_NDVI数据

2022-11-05 12:55:54

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均值。

码农公寓

相关文章