简介这份资源面向从事遥感与地球科学分析的科研人员和学生提供基于MATLAB驱动MRTMODIS Reprojection Tool对MODIS数据进行批处理重投影与数据集提取的脚本工具帮助解决大批量HDF-EOS格式数据在投影转换、区域裁剪和变量筛选中的效率问题。压缩包内共1个文件为m脚本类型体积约1KB核心是用于自动化调用MRT命令行、指定输入输出路径与投影参数的批处理脚本可一次性完成多产品重投影、镶嵌与感兴趣区域提取减少重复手工操作。目前已有418人学习下载适合具备一定遥感、GIS背景及MATLAB编程经验的用户参考能快速搭建可复用的MODIS预处理工作流将原始数据转换为GeoTIFF或ENVI等通用格式便于后续在GIS软件中集成分析也为气候研究、环境监测等场景下的数据准备提供实用起点。1. 从一堆 HDF 到一张能用的投影栅格MODIS 批处理到底卡在哪如果你手头有一批 MODIS 的 HDF 文件想把它们统一重投影、裁出研究区、再提取成 GeoTIFF 或表格那你大概率已经听过 MRT 这个名字。MODIS Reprojection Tool 就是干这件事的老牌工具配合一个批处理脚本能把几十上百景影像从正弦投影转到常见的等经纬或阿尔伯斯投影顺带做镶嵌和子集提取。听起来简单但真正上手的人都知道卡人的从来不是“会不会点按钮”而是批处理参数文件写错一个字段、瓦片编号对不上、投影参数和分辨率不匹配跑一晚上出来一堆空文件。这篇笔记就按我实际踩过的路把 MRT 批处理重投影和数据集提取的完整链路拆开讲从环境准备、参数文件写法、批处理脚本到提取像元值时怎么避开坐标偏移的坑尽量让第一次接触的人也能照着跑通熟手也能对照检查自己的参数边界。2. 先搞清楚 MRT 能做什么、不能做什么2.1 MRT 的定位与替代方案对比MRT 是专门为 MODIS 产品设计的重投影和子集工具支持正弦投影到地理经纬度、阿尔伯斯等距、兰勃特等角、墨卡托等常见投影的转换也能按经纬度范围或瓦片编号做空间子集还能顺带做波段选择和格式转换。它的优势在于对 MODIS 的 HDF 结构理解得比较透能自动识别波段、缩放因子和填充值不用自己从头解析元数据。但它的局限也很明显只认 MODIS 的 HDF 格式其他遥感数据它不管批处理靠参数文件驱动没有图形界面那么直观而且它依赖 Java 运行环境不同版本对 Java 版本还有要求。常见做法是如果只是少量文件用 MRT 的 GUI 跑一遍把生成的参数文件存下来再改成批处理模板。如果文件量大或者需要集成到自动化流程里就直接写参数文件加脚本循环。也有人用 GDAL 的gdalwarp替代但 GDAL 对 MODIS 正弦投影的处理需要手动指定projsinu lon_00之类的参数而且 HDF 里的子数据集要先用gdalinfo确认路径对新手来说反而更绕。所以如果你的数据就是 MODIS 标准产品MRT 仍然是省事的选择。2.2 环境准备与文件目录组织MRT 的安装包通常是一个压缩包解压后里面有bin、data、doc等目录。Windows 下直接运行bin里的可执行文件Linux 下需要给bin和lib里的文件加执行权限。我一般会把 MRT 放在一个不带空格和中文的路径下比如D:\tools\MRT或/opt/MRT因为参数文件里写路径时空格和中文经常导致解析失败这是血泪经验。数据目录建议按产品类型和年份分开比如MOD09A1_2020、MOD09A1_2021每个目录下放原始的 HDF 文件。输出目录单独建一个不要和输入混在一起否则批处理循环时容易把输出文件又当成输入读进去。另外MRT 在运行时会生成一些临时文件确保输出目录有写入权限。# Linux 下给 MRT 加执行权限 chmod x /opt/MRT/bin/* chmod x /opt/MRT/lib/* # 检查 Java 环境MRT 一般需要 Java 8 或更高 java -version这段命令做两件事一是给 MRT 的可执行文件和库文件加上执行权限否则运行时会报权限拒绝二是确认 Java 版本如果 Java 版本太低MRT 启动时会直接闪退或者报UnsupportedClassVersionError。参数说明chmod x后面的路径要换成你实际的 MRT 安装路径java -version输出的版本号如果低于 1.8建议先升级 Java。2.3 参数文件的结构与关键字段MRT 的批处理参数文件是一个纯文本文件通常以.prm结尾。它由若干行组成每行一个键值对格式是KEY VALUE。核心字段包括字段含义常见取值INPUT_FILENAME输入 HDF 文件路径绝对路径或相对路径SPECTRAL_SUBSET波段选择( 1 0 0 0 0 0 0 )表示选第1波段SPATIAL_SUBSET_TYPE子集类型INPUT_LAT_LONG或INPUT_PIXEL_LINESPATIAL_SUBSET_UL_CORNER左上角坐标经纬度或像元行列SPATIAL_SUBSET_LR_CORNER右下角坐标经纬度或像元行列OUTPUT_FILENAME输出文件名不带扩展名MRT 自动加.tifRESAMPLING_TYPE重采样方法NEAREST_NEIGHBOR、BILINEAR、CUBIC_CONVOLUTIONOUTPUT_PROJECTION_TYPE输出投影类型GEOGRAPHIC、ALBERS_EQUAL_AREA等OUTPUT_PROJECTION_PARAMETERS投影参数根据投影类型填写DATUM基准面WGS84等OUTPUT_PIXEL_SIZE输出像元大小根据投影单位填写这些字段里最容易翻车的是SPECTRAL_SUBSET的括号和数字个数。MODIS 不同产品的波段数不一样比如 MOD09A1 有 7 个波段SPECTRAL_SUBSET就要写 7 个数字选中的波段写 1不选的写 0。如果数字个数不对MRT 会报错或者输出空文件。另一个坑是OUTPUT_PROJECTION_PARAMETERS不同投影需要的参数个数不同地理经纬度只需要一个参数通常是 0阿尔伯斯需要多个写错了投影结果会扭曲。3. 写一个能复用的批处理参数模板3.1 从单文件参数到批处理模板最稳妥的做法是先用 MRT 的 GUI 跑一个文件把生成的.prm文件保存下来然后用文本编辑器打开把里面和具体文件相关的字段改成变量占位符。比如INPUT_FILENAME和OUTPUT_FILENAME改成{input}和{output}其他字段保持不变。这样你就得到了一个模板后续用脚本替换占位符生成每个文件的参数文件。我一般会保留一份template.prm里面除了输入输出路径其他都是固定值。如果研究区范围固定SPATIAL_SUBSET_UL_CORNER和SPATIAL_SUBSET_LR_CORNER也写死如果每个文件的范围不一样就把这两个字段也做成变量。重采样方法根据数据类型选分类数据用NEAREST_NEIGHBOR连续数据用BILINEAR或CUBIC_CONVOLUTION。输出投影如果选GEOGRAPHICOUTPUT_PROJECTION_PARAMETERS通常写( 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 )这是 MRT 的固定格式不用深究每个 0 的含义照抄就行。3.2 用 Python 生成参数文件并调用 MRT下面是一个 Python 脚本它会遍历输入目录下的所有 HDF 文件为每个文件生成一个.prm参数文件然后调用 MRT 的resample命令执行重投影。import os import subprocess from string import Template # 配置路径 mrt_bin /opt/MRT/bin/resample # MRT 可执行文件路径 input_dir /data/MOD09A1_2020 # 输入 HDF 目录 output_dir /data/output # 输出目录 template_file template.prm # 参数模板文件 # 读取模板 with open(template_file, r) as f: template_content f.read() # 研究区范围经纬度 ul_lon, ul_lat 110.0, 40.0 lr_lon, lr_lat 120.0, 30.0 # 遍历 HDF 文件 for filename in os.listdir(input_dir): if not filename.endswith(.hdf): continue input_path os.path.join(input_dir, filename) # 输出文件名去掉 .hdf 后缀加 _reproj base_name os.path.splitext(filename)[0] output_path os.path.join(output_dir, base_name _reproj) # 替换模板中的占位符 prm_content Template(template_content).substitute( inputinput_path, outputoutput_path, ul_lonul_lon, ul_latul_lat, lr_lonlr_lon, lr_latlr_lat ) # 写入临时参数文件 prm_file os.path.join(output_dir, base_name .prm) with open(prm_file, w) as f: f.write(prm_content) # 调用 MRT cmd [mrt_bin, -p, prm_file] result subprocess.run(cmd, capture_outputTrue, textTrue) if result.returncode ! 0: print(f处理 {filename} 失败{result.stderr}) else: print(f处理 {filename} 完成)这段脚本的逻辑很直接读模板、替换占位符、写参数文件、调 MRT。关键点在于Template的占位符要和模板文件里的变量名一致比如模板里写INPUT_FILENAME $input脚本里就传inputinput_path。参数说明mrt_bin要换成你实际的resample路径Windows 下是resample.exeul_lon等四个变量是研究区范围根据你的需求改subprocess.run的capture_outputTrue会捕获 MRT 的输出方便排查错误。如果 MRT 报错先看result.stderr里的信息常见的是路径不存在、波段数不对、投影参数格式错误。3.3 模板文件示例与字段解释下面是一个针对 MOD09A1 的模板文件示例输出地理经纬度投影选第 1 到第 3 波段重采样用最近邻。INPUT_FILENAME $input SPECTRAL_SUBSET ( 1 1 1 0 0 0 0 ) SPATIAL_SUBSET_TYPE INPUT_LAT_LONG SPATIAL_SUBSET_UL_CORNER ( $ul_lat $ul_lon ) SPATIAL_SUBSET_LR_CORNER ( $lr_lat $lr_lon ) OUTPUT_FILENAME $output RESAMPLING_TYPE NEAREST_NEIGHBOR OUTPUT_PROJECTION_TYPE GEOGRAPHIC OUTPUT_PROJECTION_PARAMETERS ( 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 ) DATUM WGS84 OUTPUT_PIXEL_SIZE 0.0025注意SPATIAL_SUBSET_UL_CORNER的顺序是纬度在前、经度在后这是 MRT 的约定写反了会裁到地球另一边。OUTPUT_PIXEL_SIZE的单位是度0.0025 度大约是 250 米对应 MOD09A1 的原始分辨率。如果输出投影换成阿尔伯斯OUTPUT_PROJECTION_PARAMETERS需要填更多参数具体格式可以在 MRT 的文档里找到或者用 GUI 生成一次后照抄。4. 重投影后的数据提取与坐标对齐4.1 用 GDAL 读取重投影结果并提取像元值MRT 输出的是 GeoTIFF每个波段一个文件文件名通常带波段后缀。提取像元值时可以用 GDAL 的 Python 绑定也可以用rasterio。我一般用rasterio因为它的 API 更简洁。下面是一个提取指定经纬度处像元值的例子。import rasterio import numpy as np # 打开重投影后的 TIFF tif_path /data/output/MOD09A1_2020_001_reproj.tif with rasterio.open(tif_path) as src: # 目标经纬度 lon, lat 115.0, 35.0 # 将经纬度转换为行列号 row, col src.index(lon, lat) # 读取该位置的像元值 value src.read(1, window((row, row1), (col, col1))) print(f像元值{value[0, 0]}) # 检查坐标参考系 print(fCRS: {src.crs}) print(f分辨率: {src.res})这段代码的核心是src.index(lon, lat)它把地理坐标转成栅格的行列号。注意src.index返回的是(row, col)不是(col, row)顺序别搞反。src.read(1, window...)读取第 1 波段在指定窗口内的值窗口大小是 1x1所以返回一个二维数组取[0, 0]就是单个像元值。如果src.crs不是EPSG:4326说明重投影时投影类型选错了需要检查参数文件里的OUTPUT_PROJECTION_TYPE。4.2 批量提取与表格输出如果要在多个点位上批量提取可以把点位存成 CSV循环读取每个点的经纬度然后写入结果表格。import csv import rasterio points [(115.0, 35.0), (116.0, 36.0), (117.0, 34.0)] tif_path /data/output/MOD09A1_2020_001_reproj.tif with rasterio.open(tif_path) as src: with open(extracted.csv, w, newline) as f: writer csv.writer(f) writer.writerow([lon, lat, value]) for lon, lat in points: row, col src.index(lon, lat) value src.read(1, window((row, row1), (col, col1)))[0, 0] writer.writerow([lon, lat, value])这个脚本把每个点的经纬度和提取值写进 CSV。参数说明points列表里是经纬度元组tif_path是重投影后的文件路径。如果点位很多建议先检查每个点是否在栅格范围内src.index对超出范围的点会返回超出行列号的值读取时会报错。可以在循环里加一个判断if 0 row src.height and 0 col src.width。4.3 坐标偏移的排查方法重投影后提取的值对不上最常见的原因是坐标参考系不一致。比如 MRT 输出的是地理经纬度但rasterio打开时如果没正确识别 CRSsrc.index就会按默认的像素坐标算导致偏移。排查方法是先打印src.crs和src.bounds确认 CRS 是EPSG:4326边界范围和研究区一致。如果 CRS 不对可以用rasterio.warp重新投影或者检查 MRT 参数文件里的DATUM和OUTPUT_PROJECTION_TYPE是否匹配。另一个常见问题是OUTPUT_PIXEL_SIZE设得太大或太小导致重采样后像元值被平滑或出现锯齿。对于分类数据OUTPUT_PIXEL_SIZE最好和原始分辨率一致重采样方法用NEAREST_NEIGHBOR对于连续数据可以适当放大像元但不要超过原始分辨率的 2 倍否则信息损失明显。5. 避坑与常见问题排查5.1 批处理中途报错但不知道哪个文件出问题现象脚本跑了一晚上早上发现输出目录里只有一部分文件日志里一堆错误但分不清是哪个文件。原因MRT 的错误信息有时只输出到标准错误而脚本没有捕获或者没有打印文件名。解决在循环里给每个文件加一个标识比如print(f正在处理 {filename})并且在subprocess.run之后检查returncode把stderr和文件名一起打印。更稳妥的做法是每个文件处理完后检查输出文件是否存在且大小大于 0否则记录到失败列表。5.2 输出文件为空或只有几 KB现象MRT 运行没有报错但输出的 TIFF 文件只有几 KB打开后全是 NoData。原因SPATIAL_SUBSET_UL_CORNER和SPATIAL_SUBSET_LR_CORNER的经纬度写反了或者研究区范围和影像实际覆盖范围没有交集。解决先用gdalinfo查看原始 HDF 的经纬度范围确认研究区在影像范围内。另外检查SPATIAL_SUBSET_TYPE是否设成了INPUT_LAT_LONG如果设成INPUT_PIXEL_LINE但填的是经纬度就会裁到错误位置。5.3 波段选择错误导致输出波段数不对现象明明选了 3 个波段输出却只有 1 个或者 7 个。原因SPECTRAL_SUBSET的括号里数字个数和产品波段数不匹配。比如 MOD09A1 有 7 个波段你只写了 3 个数字MRT 会按默认行为处理可能只输出第一个波段。解决查清楚产品的波段数SPECTRAL_SUBSET里写够数字选中的写 1不选的写 0。不确定的话先用 GUI 选一次看它生成的参数文件里SPECTRAL_SUBSET是怎么写的。5.4 投影参数格式错误导致重投影失败现象MRT 报错Invalid projection parameters或者输出结果扭曲。原因OUTPUT_PROJECTION_PARAMETERS的括号、数字个数或顺序不对。不同投影需要的参数不同地理经纬度通常 15 个 0阿尔伯斯需要更多。解决用 GUI 生成一次对应投影的参数文件直接复制OUTPUT_PROJECTION_PARAMETERS那一行。不要手动改数字除非你清楚每个参数的含义。5.5 路径中有空格或中文导致解析失败现象MRT 报错File not found但文件明明存在。原因参数文件里的路径包含空格或中文MRT 解析时把空格当成了分隔符。解决把所有路径改成不带空格和中文的英文路径或者用引号把路径括起来。我一般会在项目开始前就把数据目录改成纯英文无空格的路径省得后面折腾。6. 进阶把 MRT 批处理嵌入自动化流程的几个技巧如果你已经能跑通单次批处理下一步可以考虑把它嵌入更自动化的流程。我自己的习惯是用一个主脚本做调度MRT 只负责重投影和子集后续的提取、统计、入库用 Python 单独处理。这样分工明确出问题也容易定位。一个具体的技巧是不要每次都用 MRT 重新生成参数文件而是把参数文件模板和文件清单分开管理。文件清单可以是一个 CSV里面包含输入路径、输出路径、研究区范围等字段主脚本读 CSV 生成参数文件。这样修改研究区范围时只需要改 CSV不用动脚本代码。另一个技巧是并行处理。MRT 本身是单线程的但你可以用 Python 的multiprocessing或者concurrent.futures同时跑多个 MRT 进程。注意不要开太多一般设为 CPU 核心数的一半因为 MRT 运行时也会占内存。下面是一个简单的并行示例from concurrent.futures import ProcessPoolExecutor import subprocess def run_mrt(prm_file): cmd [/opt/MRT/bin/resample, -p, prm_file] result subprocess.run(cmd, capture_outputTrue, textTrue) return prm_file, result.returncode, result.stderr # 假设 prm_files 是已经生成好的参数文件列表 with ProcessPoolExecutor(max_workers4) as executor: for prm_file, code, err in executor.map(run_mrt, prm_files): if code ! 0: print(f{prm_file} 失败{err})这段代码用ProcessPoolExecutor并行调用 MRTmax_workers4表示同时跑 4 个进程。参数说明prm_files是参数文件路径列表run_mrt函数返回参数文件名、返回码和错误信息。并行处理能显著缩短大批量文件的处理时间但要注意磁盘 I/O 和内存占用如果输出目录在同一块磁盘上并行写入可能会成为瓶颈。最后说一个验证方法重投影完成后随机抽几个文件用gdalinfo检查 CRS、分辨率和范围再用rasterio提取几个已知点的值和原始 HDF 里的值对比。如果偏差在合理范围内比如最近邻重采样应该完全一致说明流程没问题。如果偏差很大回头检查投影参数和重采样方法。这个验证步骤我每次都会做虽然麻烦但能避免后面用错数据。希望帮到你。本文还有配套的精品资源点击获取