一、摘要:
本文使用Python实现MNDWI指数的计算,当然如果你需要j计算例如NDVI、NDBI等其他指数,也可以参考本文,只需要将代码中的公式以及用到的数据替换即可。
本文涉及的内容有:
(1)使用矢量文件,批量裁剪,得到研究区范围的栅格数据;
(2)批量去除多个栅格数据中的无效值;
(3)批量计算MNDWI指数(计算例如NDVI、NDBI等其他指数,只需要把数据和公式替换即可,后面会说到)
二、数据说明
我用到的数据是MOD09A1数据,因为我要计算的是MNDWI,所以需要用到的是绿波段和中红外波段,分别对应MOD09A1数据的第4波段和第6波段。
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!")
本文到这里就结束了,希望对你有所帮助!