MODIS HDF批量重投影：MRT配置与批处理实战指南

简介Mrt_bat.zip 是一份面向遥感、GIS 初学者的 MODIS 数据批处理脚本核心价值在于借助 NASA 的 MRTMODIS Reprojection Tool工具将原始 HDF-EOS 格式的 MODIS 产品批量执行重投影、格式转换与数据集提取解决手动逐景处理效率低下的问题。压缩包内仅含 1 个 Mrt_bat.m 文件大小约 1KB由 MATLAB 编写用户可直接修改输入路径、输出目录、投影参数及感兴趣区等关键选项即可驱动 MRT 完成自动批处理。对于需要处理多时相、多区域地表温度、植被指数等产品的科研人员来说该脚本能节省大量时间并便于后续 GIS 软件集成分析。资源已有 417 人学习下载适合具有一定 MATLAB 和遥感基础、希望搭建可复用处理流程的研究者使用。1. 下载到的 MODIS HDF 不能直接用Mrt_bat.zip 到底帮你解决了什么从 NASA 或 USGS 下载的 MOD11A1 这类 MODIS 产品打开 HDF 文件的那一刻你就会发现一个问题影像的行列号跟经纬度对不上放进 ArcGIS/QGIS 里和矢量边界完全错位。原因是 MODIS 原始数据用的是正弦投影Sinusoidal不是 WGS84 地理坐标而且一个 HDF 里塞了十几个子数据集SDS真正要的往往只有一个波段。Mrt_bat.zip 这个东西核心就是把 MRTMODIS Reprojection Tool的重投影、数据集提取和 Windows 批处理封装成一套固定套路输入一堆 HDF输出一组投影正确、只含目标波段的 GeoTIFF。适合手里攒了好几年 MODIS 数据、需要批量转投影做时间序列的从业者。新手按这套走能直接产出可用数据老手可以拿它当脚本模板改。2. 先跑通 MRT 再做批处理为什么选 MRT 而不是 GDAL以及 mrt.bat 环境配置2.1 为什么不用 GDAL 或 ArcGIS 批量重投影做 MODIS 重投影可选工具不少GDAL 的 gdalwarp、ArcGIS 的 Project Raster、ENVI 的 Georeference MODIS以及 NASA 官方出的 MRT。我大多数情况下会让新人直接用 MRT原因有三。第一MRT 是为 MODIS 量身定做的它认识 HDF 内部的 SDS 结构。GDAL 也能读 MODIS HDF但用 gdalwarp 之前你得先搞清楚要转哪个子数据集还要处理 HDF 里特有的缩放因子和填充值命令行写起来很长新手容易在子数据集选择上翻车。MRT 的 prm 参数文件里写一个 SDS_NAME 就能解决。第二MRT 一条命令同时完成重投影、数据集提取、范围裁剪、重采样和格式转换。ArcGIS 的 Project Raster 虽然也能转投影但提取波段和批量处理要靠 Model Builder 搭流程搭建和调试成本比写 prm 文件高得多。第三MRT 的参数文件是可复用的文本文件。换一天的数据只要改输入输出文件名投影、分辨率、重采样方法这些设置原样不动。做长时间序列比如 2000 年至今的 LST时这种可复制性极其重要。MRT 的缺点是它是个停止维护的 Java 程序界面老旧在高版本 Windows 和新的 JRE 上偶尔有兼容问题还有一个限制输入必须是 MODIS HDF不能拿普通 GeoTIFF 喂给 MRT 转投影。所以我的选型建议是手头是纯 MODIS HDF用 MRT如果是处理别人已经转好的 GeoTIFF或者要在 Linux 服务器上跑选 GDAL 更合适。两者的对比可以列成一张表。对比项MRTGDAL gdalwarp输入格式MODIS HDF-EOS几乎全部栅格格式数据集提取SDS_NAME 直接指定需要 /subdataset 语法批处理方式prm 模板 bat 循环shell/python 脚本维护状态已停止更新持续维护适合场景MODIS 时间序列批量处理通用栅格处理与服务器部署2.2 先把 mrt.bat 跑通Java 环境与 MRT_HOME 配置很多人拿到 MRT 压缩包解压之后双击 bin 目录下的 mrt.bat窗口一闪而过或者在命令行敲 mrt提示“不是内部或外部命令”。这通常不是 MRT 软件本身坏了而是两个前置条件没满足。第一个前置条件是 Java 环境。MRT 是 Java 写的运行 mrt.bat 前必须先装 JRE。这里要注意版本匹配MRT 4.1 是 2011 年左右发布的官方推荐 JRE 1.6 或 1.7。实测在 JRE 1.8 上部分功能也能跑但你要是装了很新的 JDK启动时容易报 UnsupportedClassVersionError 或者直接闪退。我建议装一个 JRE 1.7 专门给 MRT 用装完在命令行验证一下java -version看到版本号是 1.7.x 就说明 Java 基本没问题。如果系统里有多个 Java 版本最好在调用 MRT 的脚本里显式指定 Java 路径比如set JAVA_HOMEC:\Program Files\Java\jre7 set PATH%JAVA_HOME%\bin;%PATH%第二个前置条件是调用方式。mrt.bat 不会自动进 PATH直接敲 mrt 当然找不到。两个解决办法要么把C:\MRT\bin加进系统环境变量 PATH要么每次用完整路径调用。考虑到批处理脚本要反复调用我一般把 MRT_HOME 写成一个变量脚本里统一引用这样换机器只改一行。set MRT_HOMEC:\MRT %MRT_HOME%\bin\mrt.bat执行上面这行如果没有报错、出现 MRT 的版本信息就说明环境通了可以进入下一步写 prm。注意MRT 的解压路径不要带中文和空格。放在D:\Program Files\MRT这种带空格的路径下批处理脚本里经常出一些莫名其妙的问题直接把目录放在C:\MRT或D:\MRT最省心。3. 把重投影和波段提取写进一个 .prm六个必须手动改的字段3.1 prm 是所有重投影参数的唯一入口MRT 的 prm 文件是一个纯文本参数文件里面每行是一个配置项格式是配置项 值。GUI 界面你点选的所有内容最终都会生成一个 prm 文件反过来你手工写一个 prm 文件用 mrt.bat 加载它就能执行同样的操作。这也正是 Mrt_bat.zip 这类批处理脚本能成立的前提重投影的所有设置都集中在一个可文本替换的文件里批处理脚本只需要负责把文件名替换进去。初学的人容易把 prm 想复杂其实它就是个 key-value 配置文件。我建议不要用 GUI 去生成 prm而是自己新建一个模板文件把常用的配置写死每次只改输入输出文件名和关键参数这样既可控又可追溯。3.2 六个必须手动改的字段一份面向 MODIS 地表温度产品MOD11A1的 prm核心字段就六个。每个字段的作用和设置逻辑如下。第一个是INPUT_FILENAMES指定输入 HDF 的完整路径。多个输入文件用英文逗号分隔但批处理场景下我建议一次只处理一个文件这样出错时容易定位。路径里不要出现中文这是 MRT 的老毛病。第二个是SDS_NAME这就是“数据集提取”的关键字段。MODIS HDF 里通常有十几个子数据集比如 MOD11A1 包含 LST_Day_1km、QC_Day、Day_view_time、Emis_31、Emis_32 等。全转出来既慢又占磁盘正确做法是只写你需要的那个。比如提取白天地表温度就写LST_Day_1km需要质量控制波段就加上QC_Day多个 SDS 用逗号分隔。这一步能把输出数据量缩小到原来的几十分之一是整个流程里收益最大的一处。第三个是OUTPUT_FILENAME输出文件路径。后缀决定格式MRT 会按后缀自动判断写 .tif 输出 GeoTIFF写 .hdf 输出 HDF-EOS。批处理里这个字段每次都要替换成新的。第四个是OUTPUT_PROJECTION_TYPE投影类型。最常用的是GEO等经纬度适合和矢量数据叠图和UTM适合面积量算和区域统计。如果要和已有数据一致用 ArcGIS 里那个数据的投影参数填到这里。MODIS 原生的正弦投影用SIN但一般不会作为输出选项因为转出来还是对不上经纬度。第五个是RESAMPLING_TYPE。MRT 提供三种重采样NEAREST_NEIGHBOR、BILINEAR、CUBIC_CONVOLUTION。选择逻辑是连续变量地表温度、反射率用双线性离散变量质量控制波段、土地利用类型用最邻近因为双线性插值会插出不存在的类别值追求平滑视觉效果可以选三次卷积但计算量大且可能产生过冲。我做 LST 时间序列的习惯是温度波段用 BILINEARQC 波段用 NEAREST_NEIGHBOR。第六个是OUTPUT_PIXEL_SIZE输出像素大小。这个值受投影类型影响投影选 GEO 时单位是度0.01 度约等于 1 公里和 MOD11A1 原始分辨率一致投影选 UTM 时单位是米填 1000 表示 1 公里。很多人在这里填错填了 0.01 又选 UTM结果输出影像只有巴掌大。其他字段如SPATIAL_SUBSET_UL_CORNER和SPATIAL_SUBSET_LR_CORNER用于裁剪范围不需要裁剪就先不写OUTPUT_PRODUCT_FORMAT保持默认的 GeoTIFF 即可。3.3 一份可直接抄的 prm 模板下面是一份我常用的 MOD11A1 批处理 prm 模板占位符用 INPUT_FILE 和 OUTPUT_FILE 代替方便批处理脚本做文本替换。INPUT_FILENAMES INPUT_FILE SDS_NAME LST_Day_1km,QC_Day OUTPUT_FILENAME OUTPUT_FILE OUTPUT_PROJECTION_TYPE GEO OUTPUT_PROJECTION_PARAMETERS RESAMPLING_TYPE BILINEAR OUTPUT_PIXEL_SIZE 0.01 OUTPUT_PRODUCT_FORMAT GeoTIFF这份模板实现的效果是输入一个 MOD11A1 HDF提取白天地表温度和质量控制两个子数据集输出 WGS84 等经纬度投影、像素 0.01 度的 GeoTIFF。分辨率 0.01 度对应的地面距离约 1.1 公里比原始 1 公里略粗一点但好处是每个像素对应的经纬度是规整的后续做时序分析时按像素坐标直接读取即可。OUTPUT_PROJECTION_PARAMETERS这里留空即可GEO 投影不需要额外参数。如果你的数据要转到 UTM则需要在这一行填上中央经线等参数格式是0, 0, 0, 0, 0, 0这样的六个数值具体数值能通过 ArcGIS 里创建自定义投影的方式查出来。SDS_Name 里同时写了 LST_Day_1km 和 QC_Day意味着输出 GeoTIFF 会包含两个波段。后续用 Python 读文件时波段 1 是温度波段 2 是质量控制标志。我个人习惯把 QC 波段一起转出来因为后续过滤云污染和异常像元时要用它省得再回去读原始 HDF。4. 从单文件到批量重投影Windows 批处理里 for 循环、计数器与临时 prm 生成4.1 单文件先跑通再谈批量批处理脚本之前先用单文件验证 prm 模板正确。把模板里的 INPUT_FILE 和 OUTPUT_FILE 替换成真实路径保存为 test.prm然后执行C:\MRT\bin\mrt.bat -d -p test.prm参数-d表示执行数据处理-p指定 prm 文件。如果这一步跑通输出目录里出现预期的 GeoTIFF再用 gdalinfo 或 QGIS 打开检查投影正确就可以进入批量环节。如果这一步就报错先不要急着写循环把 prm 字段逐个排查否则批量脚本会跟着错一片。4.2 写一个能用的 for 循环文件排序、变量自增与临时 prm批量处理的核心难点不是 mrt.bat而是 Windows 批处理本身的几个坑。常见需求是把一个目录下的所有 HDF 逐个处理并给每个文件生成对应的 prm批处理脚本如下echo off setlocal enabledelayedexpansion set MRT_HOMEC:\MRT set SRC_DIRE:\modis\hdf set OUT_DIRE:\modis\out set PRM_TMPLE:\modis\prm_template.txt if not exist %OUT_DIR% mkdir %OUT_DIR% set /a COUNT0 for /f delims %%f in (dir /b /on %SRC_DIR%\*.hdf) do ( set /a COUNT1 set IN_FILE%SRC_DIR%\%%f set OUT_FILE%OUT_DIR%\%%~nf.tif set TMP_PRME:\modis\tmp_!COUNT!.prm powershell -Command (Get-Content %PRM_TMPL%) -replace INPUT_FILE, !IN_FILE! -replace OUTPUT_FILE, !OUT_FILE! | Out-File -Encoding ascii !TMP_PRM! echo [!COUNT!] Processing %%f !MRT_HOME!\bin\mrt.bat -d -p !TMP_PRM! del !TMP_PRM! ) echo Done. Total %COUNT% files. pause这段脚本有几个关键点需要说明。setlocal enabledelayedexpansion是 Windows 批处理里最容易踩坑的一行。for 循环内部对COUNT这类变量做自增时如果不用延迟变量展开%COUNT%会在循环开始前就被解析成初始值导致计数器永远输出 0。开启延迟展开后用!COUNT!引用变量才能取到循环内的实时值。记住一条规矩凡是 for 循环内部需要动态读取的变量一律用!变量名!而不是%变量名%。dir /b /on %SRC_DIR%\*.hdf的作用是按名称排序输出 HDF 文件列表。/b只显示文件名不带其他信息/on按名称升序排序。这个细节能保证处理顺序和文件名排序一致对后续按日期拼接影像、检查缺哪些日期的数据很有帮助。不用for %%f in (*.hdf)而用for /f接命令输出的原因是为了得到稳定的排序结果。文件名的处理用了两个批处理变量修饰符%%~nf表示取文件名主体去掉 .hdf 扩展名这样输出文件名能和输入文件名保持对应比如MOD11A1.A2010001.h25v06.006.2010002.tif。%%f则是完整的文件名。临时 prm 的生成用的是一行 PowerShell 命令。Get-Content读取模板文件-replace做文本替换Out-File -Encoding ascii写回临时文件。这里强制指定 ascii 编码非常重要原因在避坑章节详细讲。整个循环的核心思路没有变每个 HDF 对应一个临时 prmmrt 执行完立刻删除不留垃圾文件。4.3 多 tile 数据的拼接mrtmosaic 与重投影的前后顺序上面脚本处理的是单景 HDF。如果研究区覆盖多个 MODIS tile比如 h25v05、h25v06、h26v05那每个 tile 转出来的 GeoTIFF 只是研究区的一块下一步要拼接。MRT 自带的 mrtmosaic 工具负责这个活。常见做法是先把多个 HDF 拼成一个大的 HDF再做重投影顺序不能反。dir /b /on %SRC_DIR%\*.hdf %SRC_DIR%\list.txt C:\MRT\bin\mrtmosaic.bat -i %SRC_DIR%\list.txt -o %SRC_DIR%\mosaic.hdf-i指定输入文件列表-o指定输出拼接文件。拼接完成后把 prm 的 INPUT_FILENAMES 改成 mosaic.hdf 的路径SDS_NAME 依旧保留重投影一次就能得到覆盖全研究区的完整 GeoTIFF。这里有一个选择先拼接再重投影还是先重投影再拼接。我的经验是先拼接后重投影。因为 HDF 里的数据在正弦投影网格上是规则排列的MRT 对原始网格做拼接更准确如果在重投影后的经纬度网格上再做一次重采样等于两次插值精度损失叠加LST 这类连续变量的数值会发生肉眼可见的偏差。另外要注意 mosaic 之前必须确保所有 HDF 的 SDS 结构完全一致比如都是 MOD11A1 且版本相同混入一个不同产品会导致 mosaic 直接失败。5. 批处理重投影的高频报错与排查内存溢出、编码、SIN 拼接错位五个实例5.1 双击 mrt.bat 闪退提示“mrt 不是内部或外部命令”现象在 cmd 里执行 mrt.bat系统提示找不到命令双击 mrt.bat 则窗口一闪而过什么都看不清。原因分两层。第一层是 mrt.bat 不在系统 PATH 里这是“不是内部或外部命令”的直接原因。第二层是即使加了 PATH双击 bat 文件也会因为执行完自动关闭窗口而看不到错误信息常见的背后原因是 Java 没装或者 JRE 版本太高mrt.bat 里的 java.exe 启动失败。解决不要双击改成在 cmd 窗口里手动执行这样错误信息能留在屏幕上。先确认 Java 可用java -version如果确认 Java 版本没问题就修改 mrt.bat在文件开头加一行pause双击后窗口会停住显示报错信息。看到UnsupportedClassVersionError就是 JRE 版本太新换 JRE 1.7看到ClassNotFoundException一般是 MRT 解压不完整重新解压。5.2 prm 第一行就报错parameter 无法解析现象mrt.bat 执行后终端报错说 prm 文件里的某个参数无法识别特别是用记事本编辑过 prm 之后必现。原因Windows 记事本保存文本文件时默认带 UTF-8 BOM字节顺序标记MRT 是 2011 年的 Java 程序它的配置解析器不认识 BOM 这三个字节把INPUT_FILENAMES读成了\ufeffINPUT_FILENAMES自然报错。这正是批处理脚本里Out-File -Encoding ascii存在的意义强制写没有 BOM 的纯 ASCII 文件从源头避开这个坑。解决如果已经手改过 prm用记事本打开后选择“另存为”编码选“ANSI”而不是“UTF-8”如果用 VS Code 编辑右下角把编码切到“GBK 或 UTF-8 with BOM”之外的其他选项。更稳妥的是别用记事本手工改 prm交给批处理脚本统一生成。5.3 批量跑到一半 OutOfMemoryError现象处理单景 HDF 正常但批量处理大范围几十个 tile时跑到第几个文件就报java.lang.OutOfMemoryError: Java heap space进程中断。原因MRT 默认的 JVM 堆内存上限是 1GB参数-Xmx1024m处理高分辨率或大范围数据时需要把整个 HDF 的对应 SDS 读入内存多波段同时处理时内存就崩了。解决在 bin 目录下的 mrt.bat 中找到启动 Java 的那一行把堆内存参数调大比如-Xmx2048m或-Xmx4096m前提是你的机器内存够用。另外一个更有效的策略是减少同时处理的数据量prm 里 SDS_NAME 只保留真正需要的波段不要把所有子数据集都转出来再者按 tile 分批处理不要一次把几百个文件全扔进循环。5.4 拼接出来的影像经度错位、出现很大的负坐标或 0 到 360 度异常值现象mrtmosaic 拼接后重投影得到的 GeoTIFF 在 GIS 软件里显示位置完全不对坐标范围出现类似-20037508这样的数值或者经度跑到 300 多度。原因MODIS 原始正弦投影的范围在 X 方向是 -11119505 到 11119505 米左右全球范围Y 方向类似。如果你在 mosaic 阶段没有给目标投影mosaic 输出的 HDF 仍然保持正弦投影的坐标此时再交给重投影步骤投影转换参数一旦没配对比如 prm 里选了 GEO 但没指定基准面坐标就会错乱。还有一个常见手误是OUTPUT_PIXEL_SIZE填了 0.01 但投影类型还是 SIN这时输出的“0.01”被当成米而不是度结果影像只有几米宽完全显示不出来。解决mosaic 之后不要手动拼凑投影信息直接在 prm 里设置目标投影为 GEO并确保OUTPUT_PROJECTION_PARAMETERS不填任何正弦投影参数。处理完用 gdalinfo 查看坐标范围正常的全球 1 公里 LST 数据应为经度 -180 到 180、纬度 -90 到 90 附近看到异常范围立刻检查参数。5.5 提取出的 LST 数值像“假温度”直接当摄氏度用差得离谱现象重投影和数据集提取都成功LST 波段在 GIS 里显示的数值范围在 15000 到 16000 左右直接当温度读完全不像地表温度。原因MODIS 的 LST 产品存储的是缩放后的 DN 值MOD11A1 的缩放因子是 0.02真实温度开尔文 DN 值乘以 0.02。所以 16000 乘以 0.02 等于 320 开尔文减去 273.15 之后约 47 摄氏度这才是真实地表温度。MRT 做的是几何重投影和格式转换它不做定标不会帮你乘缩放因子。这也是“MODIS 下载的地表温度数据能不能直接用”这个问题的答案几何上转了投影可以直接叠图但数值必须先乘 scale factor 再做开尔文转摄氏度。解决在 Python 或 GDAL 里处理时加上定标逻辑。这步也可以写进批处理后续环节用 gdal_calc 或 Python 脚本统一处理from osgeo import gdal import numpy as np ds gdal.Open(rE:\modis\out\MOD11A1.A2010001.tif) band ds.GetRasterBand(1) dn band.ReadAsArray().astype(np.float32) lst_kelvin dn * 0.02 lst_celsius lst_kelvin - 273.15注意质量控制波段不参与这个换算QC 是位标志直接读取判断即可。6. 重投影做完不算完用 gdalinfo 和像素统计验证成品数据6.1 用 gdalinfo 检查投影是否真的转对了批处理跑完几十个文件不要急着拿去用先抽一个文件做投影验证。GDAL 自带的 gdalinfo 是最快的检查工具它能输出投影信息、坐标范围和像素大小。gdalinfo E:\modis\out\MOD11A1.A2010001.tif重点看两段输出。第一段是Coordinate System is如果显示GEOGCRS[WGS 84]说明投影转对了如果还显示Sinusoidal或者一大串PROJCS[...]里带Sinusoidal字样说明 prm 里的投影类型没生效回查第 3 章字段。第二段是Pixel SizeGEO 投影下应显示(0.010000000000000, 0.010000000000000)左右如果显示(926.625433, 926.625433)这类接近千米的米制数值说明输出还是正弦投影。6.2 批量对比像素统计快速筛出坏文件投影验证只解决“转没转对”的问题“转得好不好”要靠统计检查。我处理时间序列数据时有个习惯处理完一批立即用一个 Python 脚本遍历所有输出 GeoTIFF输出每个文件的均值、标准差、最小值和最大值然后目测序列是否平稳。from osgeo import gdal import glob import numpy as np files sorted(glob.glob(rE:\modis\out\*.tif)) for f in files: ds gdal.Open(f) band ds.GetRasterBand(1) data band.ReadAsArray() valid data[data 0] # 过滤填充值 if valid.size 0: print(f, 全空) continue lst valid * 0.02 - 273.15 print(f{f}: mean{lst.mean():.2f}, std{lst.std():.2f}, fmin{lst.min():.2f}, max{lst.max():.2f})把每个文件的统计结果按日期顺序排开看正常情况下地表温度均值应该随季节平滑变化不会出现相邻两天从 20 度跳到 50 度再跳回来的剧烈震荡。如果有文件均值明显偏离整条序列优先怀疑三个原因HDF 源数据有云污染没过滤干净、QC 波段显示该像元为无效值、重投影参数在该文件上异常。我用这套方法筛出过好几次下载源文件损坏的问题——MRT 处理损坏的 HDF 不一定会报错但输出统计会异常这也是我坚持批量处理完必须跑一遍统计的原因。最后一个建议prm 模板和批处理脚本里保留处理日期和参数说明的注释每批数据的 prm 不要覆盖旧文件。这样哪一批数据用了什么投影类型、什么重采样方法三个月后翻脚本还能对上。重投影这个活不复杂但数据几十个文件、处理跨好几天时参数可追溯就是后悔药。希望帮到你。本文还有配套的精品资源点击获取