简介针对MODIS遥感数据批量预处理需求这份资料面向遥感数据分析与GIS方向的学习者提供一套Python脚本调用NASA MRT工具自动完成重投影、镶嵌与裁剪的完整示例。压缩包共9个文件约22.59MB包含两个核心Python脚本runmrt.py、runmrtresample.py、MRT参数配置文件sample.prm、两份MODIS HDF原始样例数据、一张处理后的TIFF成果图及其ovr/xml辅助文件另有docx格式使用说明结构清晰便于上手。目前已有1379人学习下载。通过阅读代码与文档可以掌握子进程调用外部工具、批量处理HDF数据、按需重采样及格式转换等关键技能适合作为入门MODIS批处理流程的实战参考尤其有助于理解命令行工具与Python自动化脚本的结合方式。1. python调用MRT批处理MODIS数据先打破一个预期MRTMODIS Reprojection Tool是NASA处理MODIS产品的官方工具负责投影转换、重采样、裁剪和格式转换但图形界面一次只能处理一个文件。面对一年期的NDVI时间序列几十上百个HDF文件在GUI里逐个处理根本不现实。其实Python和MRT之间没有高级API中间只隔一个subprocess调用——把resample.exe、mrtmosaic.exe跑起来再在Python层补齐文件遍历、prm生成、日志和异常处理。下文从确认电脑里的MRT能调用开始到写出可复用的批处理脚本再到投影参数、镶嵌、并发和产物验证。适合遥感数据处理工程师、GIS开发者以及需要定期更新MODIS数据集的团队有Python基础但没碰过MRT的读者按章节往下走就能复现。2. Python调用MRT前先解决三件事路径、可执行文件、subprocess封装2.1 电脑里没有mrt多半是没加PATH不是没装「电脑里没有mrt」是这类任务里出现频率最高的求助但九成情况不是MRT没装上而是它的bin目录没进系统PATH。MRT安装完成后默认只写注册表和开始菜单不像常规软件那样自动配置命令行环境所以在cmd里敲resample会得到「不是内部或外部命令」的提示。先到安装目录确认Windows下默认路径类似C:\MRT或D:\Program Files\MRTLinux下常见/usr/local/MRT或/opt/MRTbin子目录里能找到resample、mrtmosaic、swath2grid这些执行文件。把bin路径加进PATH当然可以但更推荐在Python脚本里用绝对路径常量这样环境差异只集中在一个配置开关上。import os MRT_ROOT rD:\MRT # Windows 安装根目录 LINUX_MRT_ROOT /usr/local/MRT # Linux 安装根目录 def resolve_mrt(exe_name): 返回MRT可执行文件的绝对路径找不到直接抛错。 if os.name nt: candidates [os.path.join(MRT_ROOT, bin, exe_name .exe)] else: candidates [os.path.join(LINUX_MRT_ROOT, bin, exe_name)] for path in candidates: if os.path.exists(path): return path raise FileNotFoundError(f未找到 {exe_name}请检查MRT安装目录)代码说明os.name在Windows下是ntLinux下是posix利用它自动拼出正确的可执行文件后缀resolve_mrt把「找不到程序」的失败提前到脚本开头暴露而不是等循环跑到一半才报FileNotFoundError。Python环境本身不用追求最新版Anaconda或python.org的安装包都行3.8以上即可脚本不依赖任何第三方库。Linux下还有两个环境变量要留意MRT_HOME和MRT_DATA_DIR。MRT靠它们定位内部的参数模板和HDF-EOS运行时库安装说明里会让用户source一个setenv脚本或手动export。这个动作发生在shell里Python通过subprocess启动的子进程会继承这些变量所以不需要在Python内再设置。Windows安装包一般会写入注册表但如果是绿色版迁移过来的机器这两个变量丢了会报「找不到数据文件」之类的错误排查时先检查环境变量。2.2 批处理常碰到的三个MRT可执行文件各管一段MRT不是单一程序而是一组命令行工具。按数据形态选工具的规则很简单格网化产品直接resample条带数据先swath2grid多景合并先mrtmosaic再做投影转换。可执行文件职责适用产品resample重投影、重采样、裁剪、HDF转GeoTIFFMOD13Q1、MOD09GA、MOD11A1 等格网产品mrtmosaic多景HDF镶嵌成单景输出仍是HDF相邻分幅、逐日多轨数据合并swath2grid条带观测转网格MOD021KM、MOD03 等L1/L2条带数据多数批处理任务只用到前两个。比如NDVI时间序列常用的MOD13Q1或地表反射率MOD09GA都已经按分幅格网组织好直接对每个文件跑resample。如果是逐轨的L1条带数据要先swath2grid转成网格再走resample链路长一步参数也多一组。另外需要明确边界MODIS数据大气校正不是MRT的职责MRT处理层面对应的是几何校正和重投影需要大气校正后的产品就去订MOD09系列或者用6S、ATCOR、MAIAC这类专门工具别指望在prm里能找出答案。2.3 用subprocess代替os.system调用resample旧教程常见os.system(Resample.exe -p xxx.prm)的写法能跑但不适合批处理拿不到进程退出码stdout和stderr混在终端里无法程序化处理路径带空格时还要额外转义。subprocess.run返回CompletedProcess对象returncode、stdout、stderr都是结构化字段可以直接做分支判断和日志落盘。import subprocess def run_mrt(cmd_list, timeout900): proc subprocess.run( cmd_list, capture_outputTrue, textTrue, timeouttimeout, encodinggbk, # Windows下MRT输出常为本地编码 errorsignore, # 个别字符解不开就忽略保留主线索 ) if proc.returncode ! 0: raise RuntimeError( fMRT退出码 {proc.returncode}\n fstdout: {proc.stdout}\n fstderr: {proc.stderr} ) return proc参数说明capture_outputTrue把程序输出收进内存几百个文件的批处理日志才干净timeout防止某个异常文件卡死整个循环单景250m重投影900秒通常足够处理1km产品可以放宽到1800encoding按平台选Windows下MRT的报错文本可能是GBK强行用utf-8解码会先抛UnicodeDecodeErrorerrorsignore是防御性设置排错时至少能看到退出码和绝大部分输出。返回码非0就抛异常由上层循环统一决定是终止还是记日志继续。3. 写一个能直接跑的批处理脚本HDF遍历、prm模板、resample调用3.1 脚本骨架glob收集HDF按文件名生成输出批处理的结构固定glob找出所有输入HDF循环内每个文件做「生成prm、调用resample、检查结果」。真正的复杂度不在调用MRT本身而在管理每个文件的差异点——输入路径、输出路径、波段选择。把这三件事拆成独立函数后续做并发、断点续传时不用改主流程。脚本在命令行直接运行也可以在VS Code里配置好解释器后按F5调试断点时能看清当前处理到哪个文件、prm长什么样。import glob import os IN_DIR rD:\MODIS\MOD13Q1 OUT_DIR rD:\MODIS\MOD13Q1_tif def collect_inputs(in_dir): 返回排序后的HDF列表让处理顺序稳定可复现。 hdf_list glob.glob(os.path.join(in_dir, *.hdf)) if not hdf_list: raise FileNotFoundError(f目录下没有HDF文件: {in_dir}) return sorted(hdf_list)说明glob的模式.hdf 覆盖MODIS V6产品常见命名如果目录里混有 .HDF 大写后缀补一行glob.glob(os.path.join(in_dir, .HDF))合并去重。sorted排序不是形式主义后续日志对照、断点续传判断「处理到哪了」都依赖稳定顺序。3.2 prm模板固定参数写成字符串用format填充变量prm是MRT的控制文件纯文本GUI和命令行共用同一种格式。批处理的核心技巧是维护一份prm模板每次循环只替换输入、输出文件名。模板里最容易写错的是SPECTRAL_SUBSET和SPATIAL_SUBSET前者是0/1序列长度必须与产品的SDS数量一致后者是裁剪范围的左上角和右下角经纬度必须按左上右下成对填写。# 注意SPECTRAL_SUBSET 括号里 1/0 的个数必须等于该产品的SDS数量 PRM_TEMPLATE INPUT_FILENAME {input_hdf} OUTPUT_FILENAME {output_tif} SPATIAL_SUBSET_UL_LAT {ul_lat} SPATIAL_SUBSET_UL_LON {ul_lon} SPATIAL_SUBSET_LR_LAT {lr_lat} SPATIAL_SUBSET_LR_LON {lr_lon} SPECTRAL_SUBSET {spectral} RESAMPLING_TYPE {resample_type} OUTPUT_PROJECTION_TYPE GEOGRAPHIC OUTPUT_PROJECTION_PARAMETERS ( 0.0 0.0 0.0 0.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 250.0 OUTPUT_GEOTIFF TRUE 参数说明RESAMPLING_TYPE可选NEAREST_NEIGHBOR、BILINEAR、CUBIC_CONVOLUTION分类产品用最近邻避免插值产生新值连续变量产品用双线性或三次卷积OUTPUT_PIXEL_SIZE的单位跟随OUTPUT_PROJECTION_TYPEGEOGRAPHIC下是度250m大约对应0.00225度换成UTM或SINUSOIDAL后是米这里最容易因为单位错误产出像元大小离谱的文件。拿不准SDS顺序时用MRT GUI打开一个样本HDF在Spectral Subset面板勾选需要的层另存prm后当作模板GUI生成的参数行一定合法脚本只做字符串替换。提示prm路径和所有输入输出目录统一用英文MRT对非ASCII路径和文件名的支持不稳定中文路径会导致部分版本直接打不开输入文件。3.3 完整循环生成prm、调用resample、统计成败import subprocess import os RESAMPLE resolve_mrt(resample) def process_single(hdf, out_dir, prm_template): base os.path.basename(hdf) out_tif os.path.join(out_dir, base.replace(.hdf, .tif)) prm_path os.path.join(out_dir, base.replace(.hdf, .prm)) with open(prm_path, w, encodingutf-8) as f: f.write(prm_template.format( input_hdfhdf, output_tifout_tif, ul_lat40.5, ul_lon110.2, # 示例裁剪范围左上角 lr_lat30.1, lr_lon123.5, # 示例裁剪范围右下角 spectral( 1 0 0 0 0 0 0 0 0 0 0 0 0 ), # 按产品SDS数调整 resample_typeNEAREST_NEIGHBOR, )) proc subprocess.run( [RESAMPLE, -p, prm_path], capture_outputTrue, textTrue, timeout900, ) return proc.returncode, out_tif def batch_run(in_dir, out_dir, prm_template): os.makedirs(out_dir, exist_okTrue) hdf_list collect_inputs(in_dir) ok, fail 0, [] for hdf in hdf_list: try: rc, out_tif process_single(hdf, out_dir, prm_template) if rc 0 and os.path.exists(out_tif): ok 1 else: fail.append((hdf, rc)) except Exception as exc: fail.append((hdf, str(exc))) print(f成功 {ok} / {len(hdf_list)}失败 {len(fail)}) for item in fail: print(失败:, item) return fail if __name__ __main__: batch_run(IN_DIR, OUT_DIR, PRM_TEMPLATE)逻辑说明resample用-p参数读prmprm里的INPUT_FILENAME和OUTPUT_FILENAME决定输入输出因此循环里只改prm文件内容不需要拼一长串命令行参数成功判定用returncode加「tif确实存在」双条件因为MRT偶发返回0但未写文件多半与磁盘空间或输出目录权限有关每次生成的prm留在out_dir里一是留作处理记录二是排错时可以直接用GUI打开同一份prm复现问题。整套流程不打开MRT图形界面这也是批处理与GUI操作在可维护性上的本质差别。4. 投影参数、批量镶嵌和排错MRT批处理最容易翻车的三个地方4.1 投影参数怎么写GEOGRAPHIC最省事UTM和ALBERS各有边界投影选择直接影响OUTPUT_PIXEL_SIZE的单位和OUTPUT_PROJECTION_PARAMETERS的写法。批处理里三种投影最常见对比如下。投影类型OUTPUT_PIXEL_SIZE单位参数行写法典型场景GEOGRAPHIC度15个0占位与站点经纬度对齐、全球范围拼接SINUSOIDAL米15个0占位保持MODIS原始几何、定量反演UTM米参数行全0带号由MRT根据输入范围自动计算与矢量、DEM、地形图配准ALBERS等积圆锥米双标准纬线、中央经线、原点纬度、假东假北按序填写全国或省域范围制图拼接GEOGRAPHIC是最省心的批处理选型参数行全0像元单位是度输出直接带WGS84经纬度坐标。UTM的带号一般不用手填MRT按输入数据范围的几何中心自动确定但要注意多景拼接后东西跨度一旦过大仍按中心带输出会导致边缘变形这种情况换ALBERS或GEOGRAPHIC更合理。ALBERS在国内制图常用双标准纬线25N/47N、中央经线105E这组配置参数行依次填25.0、47.0、105.0、0.0后面补0到15个参数MRT严格按位置解析顺序不能错。4.2 批量镶嵌mrtmosaic和resample的顺序不能反研究区跨两个分幅比如h25v05和h26v05时每天一个产品对应两景HDF正确的做法是先mrtmosaic镶嵌成中间HDF再对整个镶嵌结果做一次resample。反过来的常见误用是每景先单独转GeoTIFF再用GDAL或桌面GIS去拼接这样重采样做了两遍接缝处的像素也被二次插值耗时和精度都吃亏。import re import subprocess MOSAIC resolve_mrt(mrtmosaic) def group_by_date(hdf_list): 按文件名中的DOY字段如A2020021把多分幅HDF分组。 groups {} for hdf in hdf_list: m re.search(rA(\d{7}), os.path.basename(hdf)) if m: groups.setdefault(m.group(1), []).append(hdf) return groups def mosaic_and_resample(hdf_group, out_dir, prm_template): list_file os.path.join(out_dir, mosaic_input.txt) with open(list_file, w, encodingutf-8) as f: for hdf in hdf_group: f.write(hdf \n) merged_hdf os.path.join(out_dir, merged.hdf) subprocess.run([MOSAIC, -i, list_file, -o, merged_hdf], checkTrue) rc, out_tif process_single(merged_hdf, out_dir, prm_template) os.remove(merged_hdf) # 中间镶嵌产物体积大转完即删 return out_tif参数说明mrtmosaic的-i参数指向一个纯文本清单文件每行一个HDF绝对路径-o指定镶嵌输出输出仍是HDF格式必须继续交给resample做投影和格式转换中间HDF保持原始分辨率MOD13Q1这种250m产品两景拼接后可能上GB处理完随手删除。group_by_date用正则A(\d{7})提取文件名里的儒略日字段把同一天的不同分幅归到一组如果产品命名不同换成对应的日期字段即可。这样主循环就变成了按日期分组对每组调用mosaic_and_resample产出该日期一景GeoTIFF。4.3 批处理常见报错对照表现象根因处理方式返回码非0提示SDS数量不匹配SPECTRAL_SUBSET长度与产品SDS数不符用GUI打开样本文件按实际SDS数改模板输出tif存在但全黑或范围错SPATIAL_SUBSET的经纬度写反确认UL是左上角、LR是右下角纬度值UL大于LR打不开含中文或空格的路径MRT对非ASCII路径支持不稳定输入输出全换英文路径文件名也改成ASCIILinux下报找不到HDF-EOS相关库MRT_HOME、MRT_DATA_DIR未设置source安装脚本或手动export后重跑不同产品复用同一个prm模板失败各产品的SDS数量与名称不同按产品各建一份模板不要共用处理速度异常慢且CPU占用低投影类型与像元单位不匹配输出了超大范围核对OUTPUT_PIXEL_SIZE单位与投影是否对应排错顺序建议固定先看returncode再看stderr最后一行最后回查prm。returncode为0却没产出基本锁定磁盘空间或输出路径权限returncode非0且报错指向prm某行直接用MRT GUI打开这份prm文件复现GUI的报错信息往往比命令行更具体。上述表格覆盖了批处理起步阶段绝大多数的失败场景。4.4 断点续传以日志为准跳过已完成批处理几百个文件时中断是常态断电、超时、人工干预都可能发生。断点续传的逻辑不需要复杂每处理完一个文件就往日志里追加一行重跑时先读日志已记录的文件直接跳过。判断依据用日志而不是文件是否存在因为MRT有时留下半成品tif靠文件存在判断会误认为已处理。def batch_run_resumable(in_dir, out_dir, prm_template, log_path): os.makedirs(out_dir, exist_okTrue) hdf_list collect_inputs(in_dir) log_keys set() if os.path.exists(log_path): with open(log_path, encodingutf-8) as f: log_keys {line.strip().split(,)[0] for line in f if line.strip()} for hdf in hdf_list: base os.path.basename(hdf) if base in log_keys: continue try: rc, out_tif process_single(hdf, out_dir, prm_template) done (rc 0 and os.path.exists(out_tif)) except Exception: done False with open(log_path, a, encodingutf-8) as f: f.write(f{base},{OK if done else FAIL}\n) f.flush()说明日志文件每处理完一条立即flush即使进程被强杀已完成的记录也已落盘重跑时把日志里的文件名读进一个set跳过已完成项只处理失败和未开始的。这样脚本就能放心交给任务计划或无人值守环境不需要担心跑到一半出问题就前功尽弃。5. MRT批处理提速与收尾多进程并发、GDAL验证和bat定时运行5.1 用multiprocessing按文件粒度并行每个resample进程独立处理数据、互不共享状态天然适合进程池并行。但MRT是CPU和IO双密集并发数不建议超过物理核数一半机械硬盘上并发过高会让磁头频繁寻道实测4个进程往往比8个更快SSD上可以适当放宽。imap_unordered按完成顺序返回结果适合观察整体进度。from multiprocessing import Pool, cpu_count def worker(hdf): try: rc, out_tif process_single(hdf, OUT_DIR, PRM_TEMPLATE) if rc 0 and os.path.exists(out_tif): return os.path.basename(hdf), True return os.path.basename(hdf), False except Exception as exc: return os.path.basename(hdf), ferror: {exc} if __name__ __main__: hdf_list collect_inputs(IN_DIR) workers min(4, cpu_count()) with Pool(processesworkers) as pool: for name, status in pool.imap_unordered(worker, hdf_list): print(name, status)参数说明workers取min(4, cpu_count())是保守起步值机器核心多可以逐步上调观察单批耗时Pool上下文管理器确保所有子进程结束再退出避免残留resample进程占用文件锁或造成后续任务冲突。5.2 用GDAL验证输出的空间参考处理完的产物不要直接拿去用先做一轮机器验证。GDAL读GeoTIFF元数据比看图判断可靠得多重点检查四件事能否正常打开、左上角坐标、像元尺寸、波段数。MRT偶尔会输出缺少空间参考的文件这类tif在桌面软件里能显示但和矢量叠加时位置全错。from osgeo import gdal def verify_tif(path): ds gdal.Open(path) if ds is None: return {ok: False, 原因: 无法打开} gt ds.GetGeoTransform() left, top gt[0], gt[3] right left gt[1] * ds.RasterXSize bottom top gt[5] * ds.RasterYSize return { ok: True, 左上角: (round(left, 6), round(top, 6)), 右下角: (round(right, 6), round(bottom, 6)), 像元尺寸: (gt[1], abs(gt[5])), 波段数: ds.RasterCount, }验证脚本放在批处理末尾遍历所有输出tif把okFalse的条目打印出来和日志里的FAIL记录交叉对比就能区分「MRT失败了」和「MRT说成功但产物不可用」两类问题。5.3 包一层bat壳交给任务计划程序最后把整个流程包成一个bat启动壳好处有两个定时任务直接指向bat不暴露Python解释器路径日志统一追加到一个文件排障只看单一日志。echo off chcp 65001 nul cd /d D:\MODIS C:\Python39\python.exe batch_mrt.py run.log 21说明chcp 65001把cmd代码页切到UTF-8避免Python打印的中文日志在重定向文件里乱码cd /d先切到脚本目录让脚本里的相对路径稳定表示追加写入21把stderr合并到stdout错误信息不会丢失。配合Windows任务计划程序设置每日定时执行新增的HDF文件会被自动处理已处理的按断点续传日志跳过整个增量流程可以长期挂着跑。本文还有配套的精品资源点击获取