做遥感的人,谁没跟MODIS打过交道呢。这颗中分辨率卫星数据免费开放、时间序列长、波段覆盖从250米到1公里,地表温度、植被指数、反照率这些常用产品全都有。可一到预处理阶段就犯了难:下载下来的HDF文件默认是Sinusoidal正弦投影,不是我们习惯的经纬度或UTM坐标;研究区稍微大一点,一景影像往往盖不全,得把同一时间两个甚至四个相邻条带拼起来。这套工作要是靠手工在专业软件里一景一景地导入、投影、拼接,一次两次还行,连做十天八天的逐日数据,人能累到怀疑人生。这时候MODIS Reprojection Tool(MRT)就该出场了。
MRT是专门给MODIS数据做重投影、拼接、重采样和格式转换的小工具。它吃进去的是MODIS的HDF原始文件,输出的是带地理坐标的GeoTIFF,支持一个命令把拼接和转投影一起完成,而且自带批处理机制,能让你把几十上百个参数文件丢给电脑自己跑。这篇文章就围绕"MRT批处理拼接+转投影"这条主线,把这些年的实操经验梳理一遍,适合刚接触MODIS数据、还没找到顺手的预处理工具的人,也适合已经装了MRT但只会点图形界面、想转命令行做批量生产的人。全文不讨论复杂的算法,只讲怎么真正把活干完。
1. 为什么MODIS数据要过一遍MRT
1.1 MODIS原始数据的"麻烦"
先花一点点时间把MODIS数据究竟麻烦在哪讲清楚。MODIS原始产品是HDF格式文件,文件名的形如:MOD13Q1.A2014185.h25v05.006.2014193022529.hdf。
拆开来看:
- MOD13Q1是产品名,这里是以植被指数产品为例;
- A2014185是日期,意思是2014年第185天;
- h25v05是条带号,h表示水平方向第二十五块,v表示垂直方向第五块;
- 006是数据版本。
MODIS全球数据按一种叫Sinusoidal正弦投影的坐标系来分块存储,全球被切割成许多个约十度乘十度的瓦片(tile),每条瓦片有一个像h25v05这样的编号。这种设计对全球产品归档非常友好,可对我们做区域研究的人就不太友好了:首先,研究区几乎不可能刚好落在一个瓦片里,通常横跨两个甚至四个瓦片,这就有了拼接需求;其次,Sinusoidal投影并非我们常用的地理坐标或UTM坐标,直接拿进GIS里跟矢量边界叠加,坐标对不上,量算面积也不直观;再者,MODIS每个产品文件里不止一层数据,比如地表温度产品里有白天温度、夜晚温度、质量控制标志等,需要按需提取。
此外,原始HDF文件在很多通用软件里没法直接识别。虽然现在也有一些商业软件能打开,但打开之后默认还是Sinusoidal,你还是得转投影。与其在项目后期反复调整,不如在预处理阶段就用MRT把拼接、投影、格式一起搞定。MRT存在的意义,就是把上述这些麻烦一次性解决。
1.2 MRT一次性解决的三个问题
用MRT做MODIS批处理,核心解决三个问题。
第一是拼接(Mosaic)。MRT允许在同一个参数文件里列出多个输入HDF文件,处理时自动把它们放到同一输出影像里。比如研究区跨h25v05和h25v06两个瓦片,就在输入文件列表里把这两个文件都写上,输出就是一景完整覆盖研究区的影像。要注意的是,MRT拼接时重叠区域的值以文件列表中靠后的文件为准,所以尽量不要在列表里放同一时间、同一位置相互重叠的文件,否则后面的会覆盖前面。
第二是重投影。MRT支持的输出投影类型很多,包括常用的地理经纬度坐标(GEO)、通用横轴墨卡托坐标(UTM)、阿尔伯斯等积圆锥(AEA)、兰伯特等角圆锥(LAM),以及保留原始正弦投影(SNS)。绝大多数区域研究用GEO或UTM就够,遇到大范围不规则区域时用AEA做面积统计更合适。重投影不是简单地给影像贴个新坐标,它涉及像元重采样,所以MRT同时提供了最邻近、双线性、三次卷积三种重采样方法。
第三是格式和波段的重组。MRT可以把HDF输出成带地理参考的GeoTIFF、TIFF等格式,还能在参数文件里指定保留哪些波段/科学数据集。很多人忽略最后这一点,其实它能让后面的处理量小不少。比如1公里地表温度产品原始文件里有很多质量控制波段,你如果只关心白天地表温度,完全可以在MRT阶段就把其他层剔除,输出一个单波段GeoTIFF,处理速度会快很多。
一句话总结:MRT像一个预处理流水线,输入是"原始HDF列表",输出是"拼接好、投影正确、按需裁剪的GeoTIFF",批处理脚本负责把这条流水线成百上千次地跑起来。
2. 环境准备与数据整理
2.1 MRT的安装与运行环境
MRT本身是一个基于Java的桌面/命令行工具,安装包里通常会自带一个Java运行环境,所以一般情况下装完就能用。安装时有几点值得注意:
- 安装目录不要带中文、不要有空格,尽量用如
D:\MRT这样简单的路径。这一点非常重要,因为MRT在处理文件时对路径里的特殊字符非常脆弱。 - 装完后核心有三个东西:图形界面启动程序、命令行程序resample、批处理程序mrtbatch。图形界面适合用来生成参数文件模板,命令行和批处理程序才是批量生产的主力。
- 为了在任何目录下都能直接调用命令,可以把MRT的bin目录加到系统环境变量PATH里。否则每次都写完整路径也行,脚本里写清楚就行。
我的习惯是,把MRT装到一个固定位置,然后在工作目录下放一个小工具目录,把所有脚本和prm参数文件放进去。这样即使换电脑,只要把工作目录拷过去,环境变量重新配一次就能接着跑。
另外多说一句,MRT的图形界面和命令行是共用同一套参数体系的。你在图形界面里配置好一次,保存出来的prm文件拿到命令行跑完全没问题。所以不要觉得"只用命令行就不需要GUI",第一次配参数时用GUI反而最不容易出错。
2.2 数据准备与目录规划
处理之前的数据整理往往比操作本身更影响效率。我的建议是按照这样的目录结构来放数据:
D:\modis\ ├── raw\2014185\ 原始HDF ├── raw\2014186\ ├── prm\ prm参数文件 ├── out\2014185\ 输出GeoTIFF ├── logs\ 批处理日志两件事务必提前做。
第一件是核对数据完整性。MODIS产品在网络上分块下载时容易碰到下载中断导致的文件损坏。最简单的方法是看文件大小是否和产品说明一致,但更靠谱的是在批量处理前先写个小脚本把所有文件列出来,检查文件名里日期、条带号是否齐整套。别嫌这一步多余,我最开始跑批量时,连续输出了几十景"缺东少西"的影像,排查了很久才发现是源文件里混进了几个下载损坏的文件。
第二件是明确你要的输出规格。输出投影是什么坐标系?UTM的话涉及带号,得根据研究区经度确定;如果研究区跨多个带,不如直接用Albers或经纬度投影,省得后面再去拼接两个不同投影的文件。输出分辨率是多少?MODIS产品一般有250m、500m、1km几种,转投影后这个分辨率最好与原始分辨率一致,否则就人为改变了数据精度。重采样选什么方法?这取决于后续分析,分类、阈值提取这类要保持原始像元值的用最邻近,做连续地表参数反演、插值这类可以选择双线性或三次卷积。
这些问题在开始跑批处理之前就想清楚,比跑完几百张图再回头返工要划算得多。
3. 写出一个合格的prm参数文件
3.1 prm文件的组成
MRT批处理的核心是一份以.prm结尾的参数文件。你可以用图形界面手动生成,也可以手写,但我的建议是第一次务必用图形界面生成一次,再拿生成的prm当模板去改。原因很简单:prm里有很多编号类参数,比如投影参数那十几个数字,不同投影方式意义完全不同,手工记忆成本太高,而图形界面会按你选的投影自动填好。
一份典型的prm文件长这样:
INPUT_FILENAMES = MOD13Q1.A2014185.h25v05.006.2014193022529.hdf,MOD13Q1.A2014185.h25v06.006.2014193022530.hdf INPUT_PATH = D:\modis\raw\2014185 SPECTRAL_SUBSET = ( 1 1 1 1 1 1 ) SPATIAL_SUBSET_TYPE = 0 SPATIAL_SUBSET_UL_CORNER = ( 0 0 ) SPATIAL_SUBSET_LR_CORNER = ( 0 0 ) OUTPUT_PROJECTION_TYPE = UTM OUTPUT_PROJECTION_PARAMETERS = ( 50 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ) RESAMPLING_TYPE = NEAREST_NEIGHBOR OUTPUT_PIXEL_SIZE = 250 OUTPUT_FORMAT = GeoTIFF OUTPUT_FILENAME = MOD13Q1.A2014185.250m.UTM.tif OUTPUT_PATH = D:\modis\out\2014185提醒一下,上面这个prm示例里SPECTRAL_SUBSET的波段数和OUTPUT_PROJECTION_PARAMETERS的具体数字只是示意,不同产品的波段数完全不同,投影参数也要以GUI导出为准,千万不要拿着别人模板里的数字硬套,否则会出现波段数不匹配、投影参数错误这类问题。
逐个字段说一下:
- INPUT_FILENAMES:输入文件列表,多个文件用英文逗号分隔,这就是实现拼接的地方。
- INPUT_PATH:输入文件所在目录。写prm时要注意MRT对路径分隔符的兼容性,统一用反斜杠或正斜杠,别混用。
- SPECTRAL_SUBSET:输出哪些波段/科学数据集。不同产品的波段数不同,建议用GUI自动生成。
- SPATIAL_SUBSET_TYPE:空间子集类型,0表示使用输入文件的覆盖范围,1表示指定左上角和右下角经纬度,2表示指定条带号。批处理时如果希望所有输出影像范围严格一致,可以设成1并固定两个角点坐标。
- OUTPUT_PROJECTION_TYPE:输出投影类型。
- OUTPUT_PROJECTION_PARAMETERS:投影参数集合,UTM、AEA等投影的核心参数都在这一串数字里。区域研究最常用的UTM投影,第一个参数通常就是投影带号,后面那些默认值在不同版本里表现不太一样,所以最稳妥的办法是用GUI生成。
- RESAMPLING_TYPE:重采样方法。
- OUTPUT_PIXEL_SIZE:输出像元大小,单位与投影单位一致。GEO下这个数字的单位是度,UTM下是米。
- OUTPUT_FORMAT:输出格式,常用GeoTIFF。
- OUTPUT_FILENAME和OUTPUT_PATH:输出文件名与路径。
3.2 两个容易踩坑的细节
第一个坑是拼接到最终输出的范围控制。MRT做拼接时,如果不设置空间子集,默认输出范围是所有输入文件的并集,也就是说如果两个瓦片一个在纬度30度、一个在纬度40度,输出结果会包含两者之间的"空白带",这个空白带可能填满NoData。这倒不算错误,但会让你的输出影像四周出现大片无效值,文件也会膨胀。我的做法是先用数据预览弄清研究区的经纬度范围,然后把SPATIAL_SUBSET_TYPE设为1,把左上角、右下角坐标写死,这样不管输入几个瓦片,输出范围都是固定的,后续做时间序列的时候影像才能一景对一景地对齐。
第二个坑是重采样方法不要乱选。做批处理时很多人习惯全部用最邻近,因为速度快;也有人追求平滑全部用三次卷积。但我见过不少失败的案例:明明是地表温度产品,却用了三次卷积重采样,结果在云边缘产生了非常可疑的异常像元。我的经验是,如果后续要做分类、阈值提取、像元级统计,默认最邻近就好;如果是做连续变量插值、制图展示,可以用双线性或三次卷积,但也一定要检查重采样后影像的数值范围有没有超出原产品的有效范围。MRT不会替你判断这个,输出影像里出现异常大值,往往不是算法问题,而是重采样方法选错了。
3.3 用GUI快速生成模板prm
这里再说说具体操作。启动MRT图形界面后,先打开一个HDF文件,此时右侧会列出所有可用的科学数据集。勾选你需要的层级,然后在Projection参数区选择投影类型,设置好重采样方法、输出像元大小、输出格式和路径,最后在File菜单里保存参数文件。这个界面里所有参数都会同步到保存的prm文件里。
我实际操作时,会先保存一个模板prm,然后用文本编辑器或者Python脚本去替换文件名、日期这些变量。比如模板里INPUT_FILENAMES只有一行,我就在脚本里把它替换成当天需要拼接的两个瓦片。这样,图形界面只用来做一次参数设计,剩下的批量工作完全是文本级的替换,速度快且不容易出错。
4. 批量处理脚本怎么写
4.1 思路:模板prm加脚本替换
真正进入批量生产阶段,我的核心工作流是"模板prm + 脚本批量生成 + mrtbatch排队执行"三步走。
第一步,在前面提到的GUI里生成一个模板prm,把投影、重采样、格式这些固定参数都设置好。第二步,写一个脚本,按日期遍历原始数据目录,根据每天的瓦片组合生成对应的prm文件。第三步,调用mrtbatch命令,让程序依次处理所有prm文件。整个过程里我只需要维护模板和脚本,新增一个日期的数据不过是一句话的事。
这里有一个很重要的取舍:prm文件里是否写死OUTPUT_FILENAME。如果写死了,每次运行都会生成同一个文件名;如果模板里删掉这个字段,命令行就可以用-o和-p参数统一控制输出目录和前缀,更灵活。我的模板里默认不写OUTPUT_FILENAME,统一交给命令行管理。
4.2 Windows下的批处理示例
如果你在Windows下工作,最简单的方案是写一个bat脚本遍历目录下的所有prm文件,逐个交给mrtbatch执行:
@echo off set MRT_BIN=D:\MRT\bin set PRM_DIR=D:\modis\prm set OUT_DIR=D:\modis\out cd /d %MRT_BIN% for %%f in (%PRM_DIR%\*.prm) do ( echo Processing %%f mrtbatch "%%f" -o %OUT_DIR% -p out_ if errorlevel 1 echo %%f failed >> D:\modis\logs\error_log.txt ) echo All done! pause注意这里有两个细节。第一,mrtbatch后面的-o指定输出目录,-p指定输出文件前缀。如果prm文件里已经写了OUTPUT_FILENAME,不同prm就有可能出现同名输出互相覆盖的情况,所以批处理模式下我建议模板prm里不要写死文件名。第二,if errorlevel 1用来把失败的任务记录到日志里,这个习惯能让你在跑完几百个任务后迅速定位失败的那几个,而不是靠肉眼翻控制台输出。
4.3 Linux下的Bash与Python示例
Linux服务器上跑批处理更常见,尤其是数据量大、需要长期跑的时候。先用Bash版本:
#!/bin/bash MRT_BIN=/opt/MRT/bin PRM_DIR=/data/modis/prm OUT_DIR=/data/modis/out LOG_DIR=/data/modis/logs mkdir -p "$OUT_DIR" "$LOG_DIR" for prmfile in "$PRM_DIR"/*.prm; do echo "Processing $prmfile" "$MRT_BIN/mrtbatch" "$prmfile" -o "$OUT_DIR" -p out_ \ >> "$LOG_DIR/mrt.log" 2>>&1 status=$? if [ $status -ne 0 ]; then echo "$(date) :: $prmfile failed" >> "$LOG_DIR/error.log" fi done如果要用Python做更灵活的控制,比如根据文件名里的日期自动生成prm,可以这样:
import os import re import subprocess from pathlib import Path MRT_BIN = r"D:\MRT\bin" RAW_DIR = Path(r"D:\modis\raw") PRM_DIR = Path(r"D:\modis\prm") OUT_DIR = Path(r"D:\modis\out") PRM_DIR.mkdir(exist_ok=True) template = Path("template.prm").read_text() for raw_dir in sorted(RAW_DIR.iterdir()): if not raw_dir.is_dir(): continue date_id = raw_dir.name hdf_files = sorted(raw_dir.glob("*.hdf")) if len(hdf_files) == 0: print(f"no hdf in {raw_dir}") continue input_list = ",".join(f.name for f in hdf_files) prm_text = template.replace("__INPUT_FILENAMES__", input_list) prm_path = PRM_DIR / f"{date_id}.prm" prm_path.write_text(prm_text) print(f"generated {prm_path}")然后外面再套一个遍历prm文件调用mrtbatch的循环就行了。这样做的最大好处是,日期、条带组合全都由文件名自动推断,不需要每天手写一整份prm文件。我实际处理某区域生长季逐日MODIS数据时,一百多天的任务就是用这种方式一次性生成的,中间只有一天因为源文件下载不完全报错,其他全部顺利跑完。
4.4 关于并行处理的体会
mrtbatch本身是单任务执行的,但你可以同时开多个命令行进程来加速。做法很简单:先按日期把prm文件分成几组,每组开一个终端窗口跑;或者用xargs -P 4这类并行工具来控制同时运行的任务数。不过有一个重要前提,就是磁盘IO和内存足够。MRT处理一个瓦片需要先把原始HDF读进内存,输出GeoTIFF又要写盘,如果同时跑四五个任务,普通机械硬盘会先扛不住。我建议先拿一个prm测一下单任务耗时和内存占用,再决定并行度;SSD上并行4个任务通常比较稳,机械硬盘老老实实两个就行。
5. 常见问题与排查经验
5.1 程序报错与内存问题
MRT最常见的报错是Java虚拟机内存耗尽,表现为控制台直接抛出内存异常然后终止。这种问题多发于一次拼接多景250m分辨率影像的场景。解决办法有几个方向:一是减少单次输入文件数量,别贪心,一个prm最多放三四个瓦片;二是如果不需要全分辨率输出,把输出像元尺寸调大一些;三是必要时给Java虚拟机分配更多内存,但MRT自带的启动脚本不一定支持直接调整内存参数,所以调整优先级不如前两者。
还有一种隐蔽的崩溃是prm文件编码问题。MRT对prm文件的编码要求比较严格,如果文件不是纯英文ASCII码,处理时会出各种奇怪的错误。所以写prm或者生成prm的脚本里不要出现中文注释,路径和文件名也一律用英文数字。
5.2 中文路径与文件名问题
因为MODIS文件名是小写字母加数字,路径多半能保持英文,问题常常出在用户自己的工作目录上。比如实验数据存在"桌面/实验二/影像"这样的目录里,MRT一碰到中文目录就罢工。这个问题的根源是MRT底层对本地字符集支持不完善,换电脑、换系统版本都可能改变表现。排查方法也很简单:如果某个prm单独处理就报错,先把它移到磁盘根下一个纯英文目录再试。
Windows下另一个隐形坑是路径分隔符。在prm文件里写INPUT_PATH时,如果用D:\modis\没问题,可以保留反斜杠;但同样一份prm拿到Linux下跑就完全不认。所以我在跨平台跑批时,统一用正斜杠,比如D:/modis/raw/2014185,这样Windows和Linux都能用。这一点可能是我实际踩过最多次的坑,值得专门记录。
5.3 输出影像"不对"的检查清单
花了几个小时跑完批处理,最怕的是打开输出影像发现坐标不对、范围不对、或者一堆NoData。我给自己列了一个检查清单,每次批量跑完先抽查几景:
- 输出影像有没有地理坐标信息?可以用GIS软件看文件属性,或者用命令行工具读一下GeoTIFF的坐标标签。如果连投影信息都没有,说明prm里投影设置有遗漏。
- 影像范围和研究区对不对?重点看输出范围是否包含研究区,尤其注意MRT在指定经纬度子集时经度跨度的写法,别把输入文件的原始Sinusoidal坐标和经纬度坐标搞混。子集的范围单位是度,不是米。
- 投影/重采样是否导致影像变形?转UTM后视场应该基本方正,如果影像被拉伸得厉害,多半是输出像元尺寸写错,或者投影参数没设置对。
- 无效值范围是否合理?打开一张影像的属性表,看NoData和最大值最小值是否在合理范围内。如果出现类似-28672这种MODIS原始填充值,说明MRT没有把它识别成NoData,这部分区域在后续统计时如果不做掩膜会产生假数据。
5.4 源文件损坏的快速识别
批量处理前花十分钟做一个文件完整性检查,能省下后面几小时的麻烦。我一般用Python脚本统计每个目录下HDF文件的数量和大小,然后跟官方产品说明对比。若有某个文件大小明显小于同类,就单独挑出来做连通性测试。这一点特别适合处理逐日数据,因为一天一个文件,出问题往往只影响一天,但要靠人眼在几百个文件里找出缺漏,成本太高。
6. 一点实操体会
最后说几点个人心得。MRT虽然是个很多年前的工具,看起来界面也不现代,但它在批处理MODIS数据这件事上非常可靠。我到现在仍然会在新环境里装一个,因为它的参数体系、prm文件规划方式是很经典的"把规则写进配置、把重复交给脚本"的思路,这套思路即使以后换到其他工具也一样适用。如果你已经掌握了prm和脚本这套做法,将来哪怕不用MRT而改用其他更现代的库,也只是换一个执行命令的问题,规划和排查的框架完全可以复用。
另外提醒一句,MRT之后官方提出了功能更全面的替代工具HEG,但在做"拼接+转投影"这个基础操作上,MRT的prm机制依然简单直接。如果你只是想快速拿到一套预处理好的MODIS影像,完全不必纠结工具是不是最新,能稳定产出就是好工具。
批处理这件事,做得多了以后你就会明白:真正耗时间的不是执行本身,而是前期的参数设计、数据整理、错误排查。把这些前期规范做扎实,几百景影像也就是脚本里跑一圈的事。