ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

MODIS HDF批处理重投影:MRT参数配置与GDAL裁剪实战

MODIS HDF批处理重投影:MRT参数配置与GDAL裁剪实战 简介这份资源面向从事遥感与地球科学分析的科研人员及学生提供基于MATLAB驱动MRTMODIS Reprojection Tool对MODIS数据进行批处理重投影与数据集提取的脚本工具可解决多产品批量转换投影、裁剪镶嵌及按感兴趣区域筛选波段的效率问题。压缩包内共1个文件为m格式的MATLAB脚本体积约1KB主要用于指定输入路径、输出格式与投影参数实现自动化处理流程。目前已有418人学习下载适合具备一定遥感、GIS背景与MATLAB编程经验的用户参考。借助该脚本读者可快速搭建MODIS批处理工作流将HDF-EOS数据转换为GeoTIFF或ENVI等常用格式并适配UTM、Lambert等投影系统从而减少重复手工操作提升大规模遥感数据分析与集成到GIS软件中的效率。1. 从一堆 HDF 到一张能用的投影栅格MODIS 批处理重投影到底在解决什么如果你手头有一批 MODIS 的 HDF 文件比如 MOD13Q1、MOD11A2、MCD43A4 这类产品想把它们统一重投影、裁出研究区、再提取成一张张能直接进 GIS 或喂给模型的栅格那你大概率会撞上同一个问题单文件用 MRT 点点鼠标还行几十上百个文件就彻底变成体力活。标题里的Mrt_bat.zip说的就是这件事——用 MRT 的批处理能力把 MODIS 数据从原始正弦投影Sinusoidal批量重投影成地理坐标或你指定的投影顺带做子集提取。它解决的不是“能不能转”而是“怎么一次转一批、参数怎么统一、转完怎么验证没翻车”。适合做遥感时序分析、地表参数反演、生态监测的从业者尤其是那些被 HDF 的投影和波段层数折磨过的人。2. 先把 MRT 的批处理逻辑吃透为什么不是点开 GUI 就完事2.1 MRT 到底做了什么批处理模式又省掉了哪一步MRTMODIS Reprojection Tool是专门处理 MODIS 层级数据的工具核心能力有三块重投影、镶嵌、子集提取。GUI 模式下你选输入文件、选输出投影、选波段、点 Run它背后其实是在拼一条命令行再调用 Java 引擎执行。批处理模式-p参数指定参数文件就是把这条命令行固化成一个文本文件让工具循环读文件列表、套同一套参数、批量输出。这里有个容易被忽略的点MRT 的批处理不是“多线程并发”它是顺序处理。你给它一个文件列表它一个一个转。所以文件多的时候瓶颈往往在磁盘 I/O 和 Java 堆内存而不是 CPU。我一般会先把文件按产品类型和年份分目录避免一次塞几百个文件进去既好排查也方便中途断了续跑。另一个关键认知MRT 对 MODIS 正弦投影的识别是内置的你不需要手动指定源投影参数。但输出投影必须显式写清楚否则默认输出还是正弦投影等于白转。很多人第一次跑完发现“怎么还是斜的”就是漏了输出投影配置。2.2 参数文件.prm的字段拆解哪几行决定成败MRT 的批处理参数文件是一个纯文本结构固定。下面是一个我常用的模板输出 Albers 等积投影做中国区域研究时比较稳# 输入文件列表 INPUT_FILES C:\modis\list.txt # 输出文件目录 OUTPUT_FILE C:\modis\out # 输出投影类型 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 ) # 重采样方法 RESAMPLING_TYPE NEAREST_NEIGHBOR # 输出像素大小 OUTPUT_PIXEL_SIZE 250 # 波段子集按 SDS 名称 SPECTRAL_SUBSET ( 1 ) SPATIAL_SUBSET_TYPE OUTPUT_PROJ_COORDS SPATIAL_SUBSET_UL_CORNER ( 70.0 10.0 ) SPATIAL_SUBSET_LR_CORNER ( 140.0 55.0 ) # 输出格式 OUTPUT_TYPE GEOTIFF逐段说。INPUT_FILES指向一个文本文件里面每行一个 HDF 绝对路径这是批处理的核心入口。OUTPUT_PROJECTION_TYPE选 ALBERS 或 GEOGRAPHIC做区域分析我倾向 Albers面积不变形。OUTPUT_PROJECTION_PARAMETERS那串数字对应椭球参数、中央经线、标准纬线等顺序不能乱写错一位投影就偏到姥姥家。RESAMPLING_TYPE对分类产品如土地覆盖必须用 NEAREST_NEIGHBOR对连续值如 NDVI、LST可以用 BILINEAR但时序分析里我仍然常用最近邻避免引入虚假值。SPECTRAL_SUBSET里的数字对应 SDS 索引不是波段号得先用gdalinfo或 MRT 的-s参数查清楚。SPATIAL_SUBSET_TYPE设成 OUTPUT_PROJ_COORDS 表示你给的角点坐标是输出投影下的不是原始经纬度这个坑后面会细说。2.3 用命令行把参数文件跑起来最小可复现流程假设 MRT 装在C:\MRT\bin参数文件叫batch.prm文件列表叫list.txt。打开命令行cd C:\MRT\bin java -jar MRTBatch.jar -d C:\modis\input -p C:\modis\batch.prm -o C:\modis\out如果你用的是mrtmosaic或旧版resample命令会略有不同但核心就这三样输入目录、参数文件、输出目录。跑之前先确认list.txt里路径没有中文和空格MRT 对这两样东西的容忍度极低路径带空格会直接报“file not found”带中文可能静默失败。我一般会把数据放在纯英文短路径下比如D:\work\modis。跑完后检查输出目录正常情况每个输入 HDF 对应一个 GeoTIFF文件名带_reproject或类似后缀。如果某个文件没输出先看控制台有没有Error字样再单独拿那个文件跑一次 GUI通常能定位到是 SDS 索引不对还是角点超范围。3. 从 HDF 到研究区栅格子集提取与批量裁剪的完整链路3.1 空间子集 vs 波段子集别把两种提取混成一件事MRT 的子集提取分两个维度空间子集SPATIAL_SUBSET和波段子集SPECTRAL_SUBSET。空间子集是按经纬度或投影坐标裁一个矩形范围波段子集是按 SDS 名称选层。很多人以为选了空间子集就自动只输出研究区结果发现输出还是全球幅面只是被裁了一刀——这没错但如果你没设SPATIAL_SUBSET_TYPE它可能按输入投影坐标裁裁出来的位置和你想要的经纬度范围对不上。我的习惯是先在 MRT 里做重投影 波段子集输出整幅 GeoTIFF再用 GDAL 做精确裁剪。这样分工明确MRT 负责投影和格式转换GDAL 负责按矢量边界裁。因为 MRT 的空间子集只支持矩形遇到不规则研究区还是得靠 GDAL。3.2 用 GDAL 做二次裁剪一条命令裁一个目录假设你已经用 MRT 转出了一批 GeoTIFF现在要用一个 Shapefile 边界裁剪。写个简单循环# 遍历目录下所有 tif按边界裁剪并输出到 clipped 目录 for %f in (C:\modis\out\*.tif) do ( gdalwarp -cutline C:\boundary\study_area.shp ^ -crop_to_cutline ^ -dstnodata -9999 ^ -co COMPRESSLZW ^ C:\modis\out\%~nf.tif ^ C:\modis\clipped\%~nf_clip.tif )-cutline指定矢量边界-crop_to_cutline让输出范围紧贴边界外接矩形-dstnodata设无效值-co COMPRESSLZW做无损压缩省空间。注意%~nf是取文件名不带扩展名Windows 批处理里的写法Linux 下换成basename即可。这一步跑完你拿到的就是投影统一、范围一致、带无效值标记的栅格可以直接进时序分析。3.3 批量提取像元值把栅格变成表格的最后一步如果你要做统计分析比如提取每个像元的 NDVI 时间序列可以用gdallocationinfo或 Python 的 rasterio。下面是一个 Python 片段遍历裁剪后的栅格按点坐标提取值import rasterio import glob import pandas as pd points [(116.3, 39.9), (121.5, 31.2)] # 示例坐标经纬度 records [] for tif in glob.glob(rC:\modis\clipped\*.tif): with rasterio.open(tif) as src: for x, y in points: # 注意这里坐标需与栅格投影一致若栅格是 Albers 需先转换 row, col src.index(x, y) val src.read(1)[row, col] records.append({file: tif, x: x, y: y, value: val}) df pd.DataFrame(records) df.to_csv(extracted_values.csv, indexFalse)关键在src.index(x, y)那行它要求你给的坐标和栅格投影一致。如果栅格是 Albers而你手里是经纬度得先用 pyproj 转一下否则提取的位置会偏出几百米甚至几公里。这个坑我在做跨区域对比时踩过明明点落在农田里提取出来却是水体值查了半天才发现是投影没对齐。4. 避坑与排查批处理重投影里最容易翻车的 5 个地方4.1 输出投影参数写错结果整体偏移现象转出来的栅格能打开但和底图对不上整体平移了几百公里。原因OUTPUT_PROJECTION_PARAMETERS里的中央经线或标准纬线填错或者参数顺序搞混。解决拿一个已知控制点在原始 HDF 和输出栅格上分别读坐标反算偏移量或者直接用 MRT GUI 配一次把生成的参数文件抄下来对比。4.2 文件列表里有隐藏字符或空行现象批处理跑到一半报“cannot open file”但手动打开那个 HDF 没问题。原因list.txt从 Excel 或网页复制时带了不可见字符或者末尾有空行。解决用notepad打开显示所有字符删掉空行和 BOM或者在生成列表时用脚本写别手敲。4.3 重采样方法选错分类产品出现新类别现象土地覆盖产品重投影后类别数变多了出现了原本不存在的过渡类。原因用了 BILINEAR 或 CUBIC 对分类值插值。解决分类产品一律 NEAREST_NEIGHBOR连续值才考虑双线性。如果不确定先拿一个文件两种方法各跑一次对比直方图。4.4 内存不足导致 Java 崩溃现象跑到第几十个文件时控制台报OutOfMemoryError进程直接退出。原因MRT 默认 Java 堆内存较小处理大文件或长列表时不够。解决在启动命令前加-Xmx4g比如java -Xmx4g -jar MRTBatch.jar ...根据机器内存调整一般 4G 够用。4.5 输出文件名冲突后一个覆盖前一个现象输出目录里文件数比输入少某些文件被覆盖。原因MRT 默认按输入文件名生成输出名如果不同子目录下有同名 HDF就会冲突。解决在list.txt里保证文件名唯一或者输出目录按产品/年份分文件夹别全堆在一起。5. 进阶技巧用 Python 驱动 MRT 并自动校验输出5.1 自动生成参数文件和文件列表手动维护list.txt和.prm在文件多的时候很烦。我一般用 Python 扫目录、生成列表、再拼参数文件import os input_dir rD:\work\modis\raw output_dir rD:\work\modis\out list_file rD:\work\modis\list.txt prm_file rD:\work\modis\batch.prm hdf_files [os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.endswith(.hdf)] with open(list_file, w) as f: f.write(\n.join(hdf_files)) prm_content fINPUT_FILES {list_file} OUTPUT_FILE {output_dir} 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 ) RESAMPLING_TYPE NEAREST_NEIGHBOR OUTPUT_PIXEL_SIZE 250 SPECTRAL_SUBSET ( 1 ) OUTPUT_TYPE GEOTIFF with open(prm_file, w) as f: f.write(prm_content)这样每次换数据目录改两个变量就行不用碰参数文件内容。注意SPECTRAL_SUBSET里的索引要根据产品调整MOD13Q1 的 NDVI 通常是第 1 个 SDS但不同版本可能变跑之前用gdalinfo确认。5.2 跑完后自动校验文件数、投影、范围三查批处理最怕“跑完了但结果是错的”。我习惯加一段校验脚本import glob import rasterio out_files glob.glob(os.path.join(output_dir, *.tif)) print(f输出文件数: {len(out_files)}) for tif in out_files[:3]: # 抽查前三个 with rasterio.open(tif) as src: print(tif, src.crs, src.bounds)看src.crs是不是你指定的投影src.bounds是不是落在研究区范围内。如果 CRS 显示None或还是 Sinusoidal说明参数没生效如果 bounds 跑到南极或太平洋说明角点坐标写反了。这一步花不了一分钟但能省掉后面几小时的返工。5.3 一个我常犯的错忘了检查 NoData 值MRT 输出的 GeoTIFF 默认 NoData 可能是 0 或 -9999取决于产品。如果你后续做统计没屏蔽 NoData均值会被拉低。我现在的习惯是裁剪时统一用-dstnodata -9999然后在 Python 里读的时候先src.read(1, maskedTrue)让无效值自动被 mask 掉。这个习惯帮我省过好几次“数据看着正常但统计结果离谱”的排查时间。希望帮到你。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进