基于Python（Arcpy）批量计算MNDWI_mndwi批量计算水体边界-CSDN博客

本文链接：https://blog.csdn.net/guoshihui123/article/details/127790706

一、摘要：

本文使用Python实现MNDWI指数的计算，当然如果你需要j计算例如NDVI、NDBI等其他指数，也可以参考本文，只需要将代码中的公式以及用到的数据替换即可。

本文涉及的内容有：

（1）使用矢量文件，批量裁剪，得到研究区范围的栅格数据；

（2）批量去除多个栅格数据中的无效值；

（3）批量计算MNDWI指数（计算例如NDVI、NDBI等其他指数，只需要把数据和公式替换即可，后面会说到）

二、数据说明

我用到的数据是MOD09A1数据，因为我要计算的是MNDWI，所以需要用到的是绿波段和中红外波段，分别对应MOD09A1数据的第４波段和第６波段。

MOD09A1数据的第4波段和第6波段的有效值范围是（-100，16000），填充值是-28672；Scale Factor（比例因子）为0.0001，这些可以在官方说明文件里找到。

三、数据处理与计算

（1）批量裁剪

裁剪这一步比较容易理解，注意看代码中的注释就行，不多说，直接上代码！

# -*- coding: UTF-8 -*-
import arcpy
import glob
import os


arcpy.CheckOutExtension('Spatial')

# 输入数据所在的文件夹
inws = r"E:\MOD09A1\withdraw_result\2020_result"

# 结果保存的位置
outws = r"E:\MOD09A1\2020_Clip"

# 选择需要用来裁剪的矢量数据，即研究区范围的shp文件
mask = r"E:\MOD09A1\Area_shp\Study_Area.shp"

# 将所有tif文件读取到rasters中
rasters = glob.glob(os.path.join(inws, "*.tif"))

# 开始循环裁剪
for ras in rasters:
    # nameT：输出文件名
    nameT = os.path.basename(ras) + "_clp.tif"
    outname = os.path.join(outws, nameT)  # 合并输出文件名+输出路径
    out_extract = arcpy.sa.ExtractByMask(ras, mask)  # 执行按掩模提取操作
    out_extract.save(outname)  # 保存数据
    print os.path.basename(ras) + " ---- 裁剪完成"
print "         ---- 全部裁剪完成 ----         "

（2）去除无效值，并乘上比例因子

关于我使用的数据，无效值是-28672，所以在处理的时候，只需要将栅格值中的-28672置为空即可；至于比例因子，如果你得数据没有这个东东就不用管，我的数据则需要乘上0.0001。了解了这些，就看代码吧，注意看注释。

# -*- coding: UTF-8 -*-
#用来检验栅格值是否都在有效范围内（-100，16000），并乘上Scale Factor（0.0001）
# -*- coding: UTF-8 -*-
import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("ImageAnalyst")  # 检查许可
arcpy.CheckOutExtension("spatial")  # 检查许可

arcpy.gp.overwriteOutput = 1
arcpy.env.overwriteOutput = 1

arcpy.env.workspace = "E:/MOD09A1/Band_6/"   # 工作空间
rasters = arcpy.ListRasters("*", "tif")  # 遍历工作空间中的tif格式数据
outPath = 'E:/MOD09A1/Band_6_OUT/'  # 输出文件路径

whereClause = "VALUE = -28672"   # 设置异常值满足的条件

# 循环rasters中的所有影像，进行去异常值操作
for ras in rasters:
    outrasters = outPath + str(ras)  # 更改输出栅格名字
    outSetNull = SetNull(ras, ras, whereClause)  # 去除异常值
    ra_DY = Times(outSetNull, 0.0001)  #乘上 Scale Factor
    ra_DY.save(outrasters)  #保存结果
    print(str(ras))

print("La La La, All is OK!")

(3)计算MNDWI

计算MNDWI需要用到第4波段和第6波段，我将五年的所有第4波段放在Band_4_OUT文件夹内，将五年的所有第6波段放在Band_6_OUT文件夹内（这一步如果不想手动，Python也能实现）。

接下来就可以开始批量计算MNDWI了！

这里有个需要注意的问题是，如果需要计算例如NDVI、NDBI等其他类型的指数，需要将代码中的公式替换，数据也要换。具体在哪里换，注意看注释！

# -*- coding: UTF-8 -*-
#用来批量计算MNDWI，用到的绿波段和中红外波段分别存在两个不同的文件夹内
import arcpy
import os

arcpy.CheckOutExtension("spatial")  # 检查许可
arcpy.gp.overwriteOutput = 1
arcpy.env.overwriteOutput = 1

arcpy.env.workspace = "E:/MOD09A1/"   # 工作空间

path_1="E:/MOD09A1/Band_4_OUT"  #所有绿波段存放的位置
path_2="E:/MOD09A1/Band_6_OUT"  #所有中红外波段存放的位置

files_1=os.listdir(path_1)  #得到所有的绿色波段文件
files_2=os.listdir(path_2)  #得到所有的中红外波段文件

outpath="E:/MOD09A1/MNDWI/"  #用来保存结果的文件夹

for file1 in files_1:
    raster_1=arcpy.Raster(path_1+'/'+file1) #读取到当前参与循环的绿色波段

    for file2 in files_2:  #读取参与循环的中红外波段
        if file2.split('.',-1)[1]==file1.split('.',-1)[1]:  #判断当前中红外波段与绿色波段是否为同一时间的数据，如果是，则计算
            print(file1.split('.',-1)[1],file2.split('.',-1)[1])
            raster_2=arcpy.Raster(path_2+'/'+file2)
            outrasters = outpath + str(file2.split('.',-1)[1]) +'_'+"MNDWI.tif" # 更改输出栅格名字
            MNDWI=(raster_1-raster_2)/(raster_1+raster_2)  #这里就是计算MNDWI的公式，如果需要计算其他指数，把这里的公式替换即可
            MNDWI.save(outrasters)　＃保存结果

print ("Finish!")

本文到这里就结束了，希望对你有所帮助！