Python批处理MODIS数据并计算NDVI

2022-07-25 13:01:39

1.登录NASA官网下载MOD13Q1数据，红框标出来的是筛选条件，我的筛选条件列出符合要求的文件如下：

NASA官网：https://ladsweb.modaps.eosdis.nasa.gov/

MODIS数据的介绍：https://www.cnblogs.com/cuteshongshong/articles/3622855.html

2.利用指定的MRT工具对MODIS数据进行批处理

MRT下载地址：链接：https://pan.baidu.com/s/1aqD4UAhPQAWq83zqsR3_2w
提取码：uv43
复制这段内容后打开百度网盘手机App，操作更方便哦

MRT安装见：https://blog.csdn.net/gisboygogogo/article/details/75784080

（1）数据准备

安装在 E:\r_MRT\bin目录下，运行E:\r_MRT\bin\Modistool.bat 进入MRT GUI界面，选择一副影像制作批量处理需要的*.prm文件。

影像中参数设置如图，需要注意的是，在设置输出影像时需要确定输出影像的格式如D:\ProfessionalProfile\MODISrelevant\pro\A2020161.tif ，最重要的是要点saveparameter file 保存A2020161.prm文件，保存后不需要run，直接退出MRT GUI即可。

将保存的A2000049.prm文件放到需要处理的MODIS *.hdf格式的影像数据的文件目录中，如D:\ProfessionalProfile\MODISrelevant\pro中。

我选择的是两景影像，MRT可以将其进行图像拼接、图像投影格式转换等多种操作。

上图中的各项参数都可以自己选择，包括输出类型、重采样类型、输出投影类型；输出的像素可写可不写，都是输出默认的。

有一个问题：每次选取某一个或者某几个波段处理的时候，总是会处理全部。（希望大神可以不吝赐教~）

（2）CMD命令实现对数据的批处理

运行cmd 命令，将工作目录设置到 E:\r_MRT\bin 中，即MRT安装目录中的 bin 文件夹中。

输入 java -jar MRTBatch.jar -d D:\ProfessionalProfile\MODISrelevant\pro -p D:\ProfessionalProfile\MODISrelevant\pro\A2020161.prm -o D:\ProfessionalProfile\MODISrelevant\pro 其中，-d 表示的是影像数据存储的目录，-p 表示经过MRT GUI处理的prm文件路径，-o 表示输出路径。这串命令表示的是对所有的影像数据批处理得到每个影像的拼接和重采样的 prm 文件。运行成功并得到所有影像的 prm 文件后，继续在输入 MRTBatch.bat（进行批处理）命令，执行这个bat文件，即进行影像的批处理。

在ENVI中打开，显示如下：

这是两景影像拼接在一起的，这样就和在NASA网站看到的“瓦片TILE”似的（下下图）呈现方式差不多，刚才处理的时候选择的投影方式是UTM-47，原始数据使用“sinusoidal 正弦曲线投影”。

3.Python计算NDVI

我们刚才处理的影像输出很多文件，有NDVI也是NIR和RED，正好可以对比一下Python计算的NDVI是否和MRT软件输出的有差异。

得到NIR和RED文件：

代码：


import numpy as np
from osgeo import gdal
import glob



# 出错：SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 36-37: malformed \N character escape
# 出错原因：没有在目录位置之前加r
# glob函数参考：https://blog.csdn.net/wc781708249/article/details/78185839
list_tif = glob.glob(r'D:\ProfessionalProfile\MODISrelevant\NIRandRED\*.tif')
out_path = r'D:\ProfessionalProfile\MODISrelevant\file'

for i in range(0,len(list_tif),2):

    NIR = gdal.Open(list_tif[i+1])
    RED = gdal.Open(list_tif[i])

    # (filepath, fullname) = os.path.split(NIR)
    # (prename, suffix) = os.path.splitext(fullname)

    red = NIR.ReadAsArray() * 0.0001
    nir = RED.ReadAsArray() * 0.0001

    # https://www.cnblogs.com/beile/p/10789333.html
    try:
    #   ndvi = (nir - red) / (nir + red+0.000000000001)加一个小的数能避免被0除的问题
        ndvi = (nir - red) / (nir + red)
    # ImportError: No module named XXX，捕获该类错误。
    except ZeroDivisionError:
        print('ZeroDivisionError has be found!')

    # 将NAN转化为0值
    nan_index = np.isnan(ndvi)
    ndvi[nan_index] = 0
    ndvi = ndvi.astype(np.float32)
    # 将计算好的NDVI保存为GeoTiff文件
    gtiff_driver = gdal.GetDriverByName('GTiff')
    # 批量处理需要注意文件名是变量，这里截取对应原始文件的不带后缀的文件名
    out_ds = gtiff_driver.Create(out_path + 'MOD13Q1A2020'+ str(i) + '.tif',
                                 ndvi.shape[1], ndvi.shape[0], 1, gdal.GDT_Float32)
    #  将NDVI数据坐标投影设置为原始坐标投影
    out_ds.SetProjection(NIR.GetProjection())
    out_ds.SetGeoTransform(NIR.GetGeoTransform())
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(ndvi)
    out_band.FlushCache()

运行结果：

下图是MRT处理的NDVI影像：

4.遇到的问题

（1）RuntimeWarning: divide by zero encountered in true_divide

将 ndvi = (nir - red) / (nir + red) 改为：ndvi = (nir - red) / (nir + red+0.000000001)，加一个很小的数就能避免，但是这样容易造成误差。

（2）MOD13Q1.A2020241.250m_16_days_red_reflectance，第241天的NDVI图出不来。

单独读取试了一下，同样的代码，第一次运行就出来下图结果，第二次就还是无法写入ndvi值。搞不懂，有知道的大神，希望多指点。

5.参考

（1）https://blog.csdn.net/qq_43177210/article/details/107283776

（2）https://blog.csdn.net/XBR_2014/article/details/85041757?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.control&dist_request_id=&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.control

码农公寓