简介针对MODIS数据重投影、镶嵌与重采样等批处理需求这套PythonMRT自动化示例兼具实用性和学习价值适合遥感数据分析、环境监测从业者及相关专业学生参考。压缩包共9个文件大小22.59MB包含两个Python处理脚本、两个HDF样例数据、一个TIFF结果及其ovr/xml辅助文件、一个prm模板配置和一份代码说明文档结构清晰便于直接运行或按需修改。目前已有1379人学习下载说明该方法已被不少同行验证过。借助这些代码与文档读者可掌握Python调用MRT工具链、理解HDF与TIFF栅格格式转换、配置投影与分辨率参数并批量完成MODIS数据预处理为后续遥感定量分析和环境监测提供可复用的技术路径。1. 为什么用Python调MRT而不是手动点GUINASA 的 MODIS 数据一景 HDF 动辄几百 MB一个年度合成产品下来是几十上百个 tileMRTMODIS Reprojection Tool的 GUI 虽然功能齐全但每次都要手动选文件、设投影、点 Resample遇到跨轨道、需要镶嵌的场景更是灾难。我在处理 GLASS17E011km 分辨率产品时深有体会一天的数据就要覆盖 h10v08、h10v09 等多个分幅手工处理费时不说参数一错就得全部重跑。正确的路子是让 MRT 的命令行接口跑批用 Python 做调度遍历 HDF 文件、动态生成 prm 参数文件、调用 MRT、检查返回码、收集输出 TIFF。这条流程在 MODIS 处理链里非常成熟本压缩包里的 runmrt.py、runmrtresample.py 和 sample.prm 就是一组可直接复用的骨架适合做遥感数据分析、环境监测和植被长势时序分析的人也适合被重复性预处理逼疯的 GIS 从业者。2. MRT 的工作方式与 sample.prm 配置解析2.1 MRT 到底在干什么MRT 不是简单的格式转换器。一个典型的 MODIS HDF 文件里包含多个科学数据集SDS比如反射率波段、质量标记层、云掩膜等。MRT 做的事是读取 SDS → 按指定投影参数重投影 → 按需重采样到目标分辨率 → 输出为 TIFF 等格式。如果输入是多个 tile 文件它还能把它们拼成一个镶嵌结果。理解这一点才能看懂 sample.prm 里为什么有SPECTRAL_SUBSET、SPATIAL_SUBSET_TYPE这些字段——做统计建模或深度学习不需要全部波段在 MRT 阶段直接把用不上的波段切掉能省大量磁盘 I/O 和后续处理时间。2.2 sample.prm 关键字段说明压缩包里的 sample.prm 是 MRT 的参数文件新版 MRT 虽然也有图形化生成 prm 的功能但批处理场景下你只需要维护这一份文本模板。它的典型结构如下INPUT_FILENAME E:\MODIS\GLASS17E01.V50.A2018001.h10v08.2020363.hdf SPECTRAL_SUBSET ( 1 1 1 1 1 1 1 1 ) SPATIAL_SUBSET_TYPE INPUT_LAT_LONG SPATIAL_SUBSET_UL_CORNER ( 50.0 70.0 ) SPATIAL_SUBSET_LR_CORNER ( 30.0 100.0 ) OUTPUT_FILENAME E:\MODIS\out\A2018001.h10v08.tif RESAMPLING_TYPE NEAREST_NEIGHBOR OUTPUT_PROJECTION_TYPE ALBERS OUTPUT_PROJECTION_PARAMETERS ( 0.0 0.0 25.0 47.0 105.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ) OUTPUT_PIXEL_SIZE 1000.0这段配置在 MRT 里叫 prm 文件旧版 MRT 只认这种固定格式的文本参数。我一般在处理区域尺度产品时把SPATIAL_SUBSET_TYPE设成INPUT_LAT_LONG配合两个角点控制裁切范围少切一点数据后面的重投影和镶嵌就快不少。2.3 重点参数怎么改参数含义常见改动INPUT_FILENAMEHDF 输入路径换成你自己的 MODIS 文件SPECTRAL_SUBSET要保留的 SDS 波段1 保留0 舍弃只用可见光波段就改成( 1 1 1 0 0 0 ... )SPATIAL_SUBSET_TYPE裁剪方式INPUT_LAT_LONG、NO_SUBSET等不裁剪可删除这两行RESAMPLING_TYPENEAREST_NEIGHBOR/BILINEAR/CUBIC_CONVOLUTION分类产品用最近邻地表温度等平滑产品用双线性OUTPUT_PROJECTION_TYPEALBERS/UTM/SIN/GEO中国区常用 ALBERS 或 UTMOUTPUT_PIXEL_SIZE输出像元大小米1km 数据写 1000想提高分辨率可以改 500注意NEAREST_NEIGHBOR是唯一能保住原始 DN 值含义的重采样方法质量标记层如 view zenith angle和反照率产品千万别用双线性否则会把无效值拉出像 1.2 这种物理上不存在的 DN 数值。2.4 MRT 命令行入口MRT 安装完有两个核心程序resample.exe和mrtmosaic.exe。前者做单景重投影/重采样后者做多景镶嵌。压缩包里的 runmrt.py 和 runmrtresample.py 对应的就是这两条链路的封装E:\MRT\bin\resample.exe -f -p E:\MODIS\sample.prm E:\MRT\bin\mrtmosaic.exe -i E:\MODIS\filelist.txt -o E:\MODIS\Mosaic.hdf-f代表强制覆盖输出文件-p指定 prm 参数文件mrtmosaic的-i接收一个纯文本列表每行一个 HDF 路径输出是还没重投影的镶嵌 HDF。你可以把 mrtmosaic 理解成只做拼接不做投影变换真正的重投影和重采样要再交给 resample 完成。3. runmrt.py 怎么封装 MRT 单景处理3.1 用 subprocess 调用外部程序Python 调用 MRT 最直接的方式是subprocess.run配合参数列表。runmrt.py 里的核心逻辑大致如下import subprocess import os def run_mrt(prm_path, exe_pathrE:\MRT\bin\resample.exe): cmd [exe_path, -f, -p, prm_path] ret subprocess.run(cmd, capture_outputTrue, textTrue) if ret.returncode ! 0: raise RuntimeError(fMRT failed: {ret.stderr}) return ret.stdout if __name__ __main__: prm rE:\MODIS\sample.prm out run_mrt(prm) print(out)这里有两个细节一是 MRT 某些版本会把正常日志也打到 stderr所以returncode 0不代表没有警告要配合后面的结果文件检查二是在 Windows 上如果 exe 路径含空格用列表形式传参而不是字符串拼接能避免 shell 转义问题。这是 Python 调用外部工具最常见的坑subprocess 里字符串命令的坑我在别的项目里踩过不止一次。3.2 动态生成 prm 而不是手写死单景执行的价值有限批处理的关键是每一景生成自己的 prm。我的做法是把 sample.prm 当模板用字符串替换动态生成临时 prm 文件import tempfile template open(rE:\MODIS\sample.prm, encodingutf-8).read() def make_prm(hdf_path, out_tif_path, out_dirE:/MODIS/out): prm_text template.replace( E:\\MODIS\\GLASS17E01.V50.A2018001.h10v08.2020363.hdf, hdf_path.replace(\\, \\\\) ) prm_text prm_text.replace( E:\\MODIS\\out\\A2018001.h10v08.tif, out_tif_path.replace(\\, \\\\) ) fd, tmp_path tempfile.mkstemp(suffix.prm, dirout_dir) with os.fdopen(fd, w, encodingutf-8) as f: f.write(prm_text) return tmp_path这里有个隐藏的坑MRT 对 Windows 路径里反斜杠的解析在不同版本里不一致。新版本直接接受E:\MODIS\xxx老版本可能把\M当成转义字符。稳妥做法是统一把单个反斜杠替换成\\\\或者干脆所有路径都用正斜杠E:/MODIS/xxxMRT 两种都能识别。压缩包里的代码使用说明.docx对这些路径细节也有对应说明跑之前值得翻一下。3.3 遍历 HDF 数据集并补上输出目录批量处理的重点在文件发现和输出目录管理。以下是最小可用的循环from pathlib import Path hdf_dir Path(rE:\MODIS\hdf) out_dir Path(rE:\MODIS\tiff) out_dir.mkdir(parentsTrue, exist_okTrue) for hdf in hdf_dir.glob(*.hdf): out_tif out_dir / (hdf.stem .tif) if out_tif.exists(): continue prm_file make_prm(str(hdf), str(out_tif)) try: run_mrt(prm_file) print(f[OK] {hdf.name} - {out_tif.name}) except RuntimeError as e: print(f[FAIL] {hdf.name}: {e}) finally: os.remove(prm_file)这段逻辑照顾了三件事输出已存在就跳过实现断点续跑单个 HDF 失败不影响整个批次临时 prm 用完删除避免堆积。真实生产环境里我会把失败的文件名追加到 error.log 里批处理结束后统一排查而不是盯着终端翻记录。3.4 关于电脑里没有 MRT怎么办有留言问电脑里没有 mrt——MRT 是 NASA 2012 年前后停止更新的老工具安装包里自带 Java 运行时但前提是系统里得有 Java JDK 1.6 及以上版本新一些的 JDK 11 也能跑。装完后建议手动把E:\MRT\bin加进系统 PATH否则subprocess.run在子进程里找不到resample.exe。我第一次跑这类脚本时习惯先打印环境和路径再进入批处理循环python -c import os; print(os.path.exists(rE:\MRT\bin\resample.exe))如果这条返回False先别急着查脚本逻辑把 MRT 路径和 Java 环境调通再继续。4. runmrtresample.py 与多景镶嵌重采样4.1 为什么要先镶嵌再重采样MRT 官方推荐的流程是用mrtmosaic.exe把多个 tile 拼成一个 HDF不做投影变换再用resample对镶嵌结果做重投影和重采样。压缩包里 runmrtresample.py 做的事情本质上就是这条链路的编排。先镶嵌再重投影的好处是所有 tile 在原始正弦投影SIN网格上天然对齐重投影只做一次避免两次插值引入的信息损失和拼接缝误差。def build_filelist(filelist_path, hdfs): with open(filelist_path, w, encodingutf-8) as f: f.write(\n.join(str(h) for h in hdfs)) filelist rE:\MODIS\filelist.txt build_filelist(filelist, sorted(hdf_dir.glob(*.hdf))) mosaic_cmd [rE:\MRT\bin\mrtmosaic.exe, -i, filelist, -o, rE:\MODIS\Mosaic_temp.hdf] subprocess.run(mosaic_cmd, checkTrue)这里强调sorted()的原因mrtmosaic 的输出顺序和输入顺序相关非确定性排序会导致镶嵌结果中同一像素可能取自不同输入尤其当相邻 tile 在填充值或云掩膜上存在差异时结果会对输入顺序敏感。用日期加 tile 编号排序至少保证同一批数据每次处理结果可复现。4.2 重采样脚本里该做什么重采样这一步和第三章的单景重投影几乎一样区别是输入变成Mosaic_temp.hdf而且 prm 里的SPATIAL_SUBSET不再裁剪改成覆盖整个研究区的范围mosaic_prm make_prm(rE:\MODIS\Mosaic_temp.hdf, rE:\MODIS\A2018001.NPP.tif) run_mrt(mosaic_prm)压缩包里出现的A2018001.NPP.tif、A2018001.NPP.tif.ovr和.aux.xml就是这一步之后的产物。.tif.ovr是金字塔文件GIS 软件加载大 tif 时自动读取.aux.xml描述统计信息和坐标参考方便 ArcGIS 或 QGIS 直接识别。如果你的下游是自己写代码读数组这两个文件直接删掉完全不影响像素值。4.3 批量跨时间序列处理处理一整年的 GLASS 数据常见做法是日期在外层循环globe 匹配在里层from datetime import date, timedelta start date(2018, 1, 1) for i in range(46): # GLASS 8天合成产品一年约46期 d start timedelta(days8 * i) doy d.timetuple().tm_yday files list(hdf_dir.glob(f*A{doy:03}*.hdf)) if len(files) 2: continue # 同一个日期的所有tile做镶嵌和重采样8 天合成产品就用 8 天步长扫日产品步长换成 1 即可。GLASS17E01 的文件名里有 V50版本号和 2020363处理日期用*A{doy:03}*只匹配产品日期部分避免把处理日期当产品日期这是这类时序批处理最容易看走眼的地方。4.4 内存与并发不要无脑 multiprocessingMRT 处理高分辨率镶嵌时内存占用不小尤其SPECTRAL_SUBSET保留了 7 个波段时多个 tile 镶嵌可能吃进 3 到 4GB 内存。我实际跑 2018 年全年数据时没有用multiprocessing.Pool(8)一次开 8 个 MRT 进程而是外层并行内层串行——同一天的不同 tile 用 mrtmosaic 处理不同日期之间用concurrent.futures最多开 2 到 3 个进程from concurrent.futures import ProcessPoolExecutor, as_completed def process_one_day(doy): return doy, process_glasses_data(doy) with ProcessPoolExecutor(max_workers3) as pool: futures [pool.submit(process_one_day, doy) for doy in doy_list] for fut in as_completed(futures): doy, ok fut.result() print(f{doy}: {done if ok else failed})并发数开到 3 就已经比串行快一倍左右且不触发 OOM开到 6 以上在 Windows 上容易因为磁盘 I/O 竞争和文件锁反而变慢。这类 IO 密集型批处理的通病就是并发开得越猛磁盘排队越严重最终耗时不一定下降。5. 验证结果与三个实用技巧5.1 用 gdalinfo 验证投影和分辨率重采样完别急着进模型先用 GDAL 检查输出gdalinfo E:\MODIS\A2018001.NPP.tif重点看Size、Origin、Pixel Size和你 prm 里写的是否一致。GLASS 数据经常出现 prm 里写了 ALBERS 但输出仍是 SIN 的情况原因多半是OUTPUT_PROJECTION_PARAMETERS里标准纬线填反了——Albers 等积投影要求第一条标准纬线小于第二条且不能相等填错时 MRT 不报错只是静默回退到 SIN 网格。5.2 检查 NoData 和有效值范围MRT 重投影后边缘会出现一圈 NoData 或 0 值。处理统计产品时我会用 Python 快速检查有效像元比例from osgeo import gdal import numpy as np ds gdal.Open(rE:\MODIS\A2018001.NPP.tif) arr ds.GetRasterBand(1).ReadAsArray() valid np.sum(arr 0) print(fvalid ratio: {valid / arr.size:.3f})如果有效比例低于 0.5多半是SPATIAL_SUBSET范围没盖全或 tile 缺景回到第 4 章的 filelist 和排序逻辑去查不要先怀疑算法。5.3 一个节省磁盘空间的小技巧Mosaic_temp.hdf中间产物体积通常是最终 tif 的两倍以上跑完直接删掉。如果确认不需要中间结果在 resample 成功后顺手清理临时 prm 和filelist.txtfor tmp in (rE:\MODIS\Mosaic_temp.hdf, rE:\MODIS\filelist.txt): if os.path.exists(tmp): os.remove(tmp)A2018001.NPP.tif.ovr也可以删掉让 ArcGIS 或 QGIS 在第一次打开时重新生成金字塔省去从 HDF 到最终产品这一路的临时文件占满磁盘的麻烦。这套批处理流程跑通一次之后后续换数据源、换投影、换分辨率都只需要改 sample.prm 里的五六个参数而已。本文还有配套的精品资源点击获取