ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

MRT批处理MODIS重投影:从HDF到GeoTIFF的完整实操与避坑指南

MRT批处理MODIS重投影:从HDF到GeoTIFF的完整实操与避坑指南 简介这份zip压缩包面向遥感、地理信息领域的科研与工程人员解决MODIS数据批量重投影与定制化区域提取的繁琐问题。包内依托NASA官方MRT工具通过一个MATLAB脚本Mrt_bat.m串联输入路径、输出格式、投影参数等指令实现HDF-EOS格式到GeoTIFF/ENVI等常见格式的自动转换并支持按感兴趣区域或波段筛选数据适合有一定遥感和MATLAB基础的学习者参考。整个资源仅1个文件大小1KB脚本短小精悍便于直接阅读与二次修改。目前已有417人学习浏览验证了其在MODIS批处理场景下的实用价值。读者可借助该脚本快速理解MRT批处理流程掌握重投影、数据集提取的自动化写法从而迁移到大规模地表温度、植被指数等产品的处理中显著提升遥感数据预处理效率。1. MRT批处理MODIS重投影不是玄学先把这套工具链看明白干遥感的人迟早会撞上同一个坎从LP DAAC或NASA Earthdata下载的MODIS标准产品清一色是HDF格式、Sinusoidal投影文件名带着h21v04这种行列号直接丢进GIS里根本没法跟当地矢量数据叠到一块儿。Mrt_bat.zip这类批处理脚本解决的正是这个环节——用MRT把一堆MODIS HDF逐个重投影成GeoTIFF再挑出自己真正要用的子数据集比如地表温度LST或者NDVI一次性批量跑完省掉手动在MRT GUI里点几百次鼠标的时间。这套做法适合做长时间序列LST、植被指数或者地表反射率研究的人也适合接了区域制图项目的工程师。它不解决下载问题也不解决后续建模问题只管把“原始的、投影奇怪的、波段混杂的”MODIS数据变成“能直接用的、和你的研究区坐标系一致的”栅格文件。2. MODIS数据的“黑匣子”先拆开HDF结构、Sinusoidal投影与MRT三个组件2.1 MODIS标准产品为什么是HDFSinusoidal不重投影就没法用MODIS标准产品MOD11A1、MOD13Q1、MOD09GA这些出厂时使用 Sine 投影也就是 Sinusoidal 投影空间网格按全球等面积划分成水平h和垂直v的瓦片每个瓦片大约1200×1200公里像元分辨率按产品不同从250米到1公里不等。NASA这么做是为了全球拼接方便但到了区域应用场景就非常难受你研究区在内蒙古中东部可能横跨h25v04、h26v04两块瓦片而且正弦投影下瓦片接缝是斜的没法直接和UTM坐标系的土地利用数据、气象站点数据做空间运算。另一个更隐蔽的问题是HDF内部结构。一个MOD11A1文件里面不是只有一层栅格而是塞了一大堆科学数据集SDSScientific Data Set比如LST_Day_1km、QC_Day质量控制、Day_view_time、Day_view_angle等等。有的SDS是16位整型有的是8位整型还有的是浮点。你真正想用的可能就一个子数据集但用ArcGIS直接拖进去只能看到一个多波段栅格波段名也不是SDS原名处理起来一头雾水。MRT的核心价值就是把“HDF里的某个SDS”提出来、重投影、转成GeoTIFF一步到位。有关“modis下载地表温度数据可以直接用吗”这个问题答案很明确不能直接用。MOD11A1的LST_Day_1km是16位无符号整型存的是开尔文温度乘以0.02之后的值无效值比如云遮挡区域用0或特殊值表示而且投影是Sinusoidal。直接用意味着坐标不对、数值不对、单位不对三错叠加。MRT批处理解决的就是投影和子集提取这两个问题数值定标通常还要在后续Python或ENVI里补一步乘以0.02、把无效值过滤掉。2.2 MRT的三个落地文件resample.exe、prm参数文件、Java环境MRTMODIS Reprojection Tool是NASA官方发布的MODIS重投影工具虽然官方早就停止更新、出了替代品但它在批处理和MODIS专项支持上依然有大量存量用户。MRT安装完后实际干活的是三个东西。第一个是bin目录下的resample.exe这是命令行核心程序功能是读取一个HDF输入文件、套用参数文件、输出一个GeoTIFF。第二个是参数文件.prm文本格式里面写明了输入文件、输出文件、投影类型、重采样方法、要输出的SDS列表、空间范围等。第三个是Java运行环境因为MRT的图形界面和部分库依赖Java而且是很老的32位Java这个坑后面会细说。用GDAL做MODIS重投影的对比值得先摆出来因为这决定了你是不是非要折腾MRT对比项MRT (resample.exe)GDAL (gdalwarp)MODIS HDF子数据集识别直接支持SDS名称提取需要手写HDF4::MODIS_...的子数据集路径正弦投影参数内置无需自己配需要指定Sinusoidal投影参数或依赖prm批处理bat脚本for循环即可同样可以bat或Python调用安装难度老软件Java环境配置麻烦用conda/OSGeo4W安装简单维护状态已停止更新持续维护我自己的经验是如果机器上已经装了GDAL、又是临时转几个文件直接用gdalwarp更省事但如果要批量处理上百个HDF、还要提取指定SDS并按原文件名输出MRT这套批处理逻辑更顺因为它一个prm文件就能把“提哪个波段、用什么投影、怎么重采样”全部固定下来循环里只换输入输出文件名就行。3. 单条命令先跑通再批处理MRT resample的输入输出与prm文件的最小用法3.1 最小可用示例一条resample命令把MOD11A1重投影成GeoTIFF先别急着写循环第一步是手动跑通一条命令。假设你手头有一个MOD11A1文件放在E:\modis\hdf\MOD11A1.A2019001.h25v04.006.2019031154523.hdfMRT安装在E:\MRT先创建一个最小prm文件我一般叫它default.prm内容如下INPUT_FILENAME E:\modis\hdf\MOD11A1.A2019001.h25v04.006.2019031154523.hdf OUTPUT_FILENAME E:\modis\out\MOD11A1.A2019001.h25v04.006.2019031154523_LST.tif RESAMPLING_TYPE NEAREST_NEIGHBOR OUTPUT_PROJECTION_TYPE UTM OUTPUT_PROJECTION_PARAMETERS 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 OUTPUT_PIXEL_SIZE 1000.00000 OUTPUT_SDS_NAME LST_Day_1km然后打开cmd切到工作目录执行E:\MRT\bin\resample.exe -i E:\modis\hdf\MOD11A1.A2019001.h25v04.006.2019031154523.hdf -o E:\modis\out\MOD11A1.A2019001.h25v04.006.2019031154523_LST.tif -p E:\modis\default.prm这段命令的逻辑是-i指定输入HDF-o指定输出GeoTIFF-p指定参数文件。resample.exe会先用-p读入投影和子集设置再用-i和-o覆盖prm文件里写的输入输出路径。这样做的意义是同一个prm可以复用于任意多个HDF不用每个文件都去改prm里的文件名批处理时只需要循环里改-i和-o的参数值。参数文件里每一行都有讲究RESAMPLING_TYPE是重采样算法处理LST这类连续场我用NEAREST_NEIGHBOR因为双线性会对无效值边缘做插值、把云污染区域的值抹开OUTPUT_PROJECTION_TYPE UTM表示输出为UTM投影UTM带号在OUTPUT_PROJECTION_PARAMETERS里通过字符串指定OUTPUT_PIXEL_SIZE 1000是输出像元大小单位米对应MOD11A1的1公里分辨率。OUTPUT_SDS_NAME指明我们要的SDS这里只提取LST_Day_1km一个子数据集输出就是一个单波段GeoTIFF而不是整个HDF的多波段堆叠。跑完之后用GIS软件打开输出文件确认坐标系是WGS84/UTM、范围和研究区对得上、像元大小是1000米。这一步过了批处理才有意义否则脚本跑一百个文件也是白跑。3.2 批量处理才是日常for循环遍历HDF并按原文件名生成输出单条命令通了之后批处理脚本的核心就是一个for循环。下面这个bat脚本是最常见的写法我实际项目中改个路径就能用echo off setlocal enabledelayedexpansion set MRT_BINE:\MRT\bin\resample.exe set PRM_FILEE:\modis\default_utm.prm set IN_DIRE:\modis\hdf set OUT_DIRE:\modis\out cd /d %IN_DIR% if not exist %OUT_DIR% mkdir %OUT_DIR% for %%f in (*.hdf) do ( echo [Processing] %%f %MRT_BIN% -i %%f -o %OUT_DIR%\%%~nf_LST.tif -p %PRM_FILE% if !errorlevel! 0 ( echo [OK] %%~nf_LST.tif ) else ( echo [FAILED] %%f ) )这个脚本的逻辑分四段前三行定义MRT路径、prm路径、输入输出目录cd /d %IN_DIR%保证在HDF所在目录执行因为for循环里%%f只展开文件名、不带路径循环体对每个.hdf文件调用resample.exe输出文件名用%%~nf取原始文件名前缀再拼“_LST.tif”后缀比如原文件是MOD11A1.A2019001.h25v04.006.2019031154523.hdf输出就是MOD11A1.A2019001.h25v04.006.2019031154523_LST.tif保证和输入一一对应。这里有两个关键点容易踩坑。第一bat文件里的for变量必须写成%%f如果你直接在cmd命令行里测试循环则写成%f两者在批处理文件里混用会导致变量不展开、循环只跑一次或者报语法错。第二setlocal enabledelayedexpansion这行不是摆设后面!errorlevel!必须在延迟变量展开环境下才能取到每条命令执行后的返回值如果不加这行errorlevel会被当作文本原样输出批处理结果要么全显示OK、要么全显示FAILED根本没法判断单个文件是否成功。另外要提醒的是prm文件里的INPUT_FILENAME和OUTPUT_FILENAME在命令行指定了-i/-o时会被覆盖所以一套prm可以应对所有文件。但如果prm里有SPATIAL_SUBSET或SPECTRAL_SUBSET这类按文件内容变化的设置批处理前一定要确认所有输入HDF的产品类型和波段结构一致。混入一个MOD13Q1到全是MOD11A1的文件夹里脚本不会报错但输出结果会缺波段或者全黑。4. 数据集提取与投影参数设置prm文件里那几行怎么填4.1 SPECTRAL_SUBSET与OUTPUT_SDS_NAME只提取你要的波段别把QC带出去很多人拿到MRT后只看GUI界面忽略prm文件里一个非常实用的字段SPECTRAL_SUBSET。它控制的是“HDF里的哪些SDS参与输出”。以MOD11A1为例它的SDS顺序大致是LST_Day_1km、QC_Day、Day_view_time、Day_view_angle、LST_Night_1km、QC_Night、Night_view_time、Night_view_angle共8个。如果你在GUI里直接勾选输出GeoTIFF默认全选结果是一个8波段文件每个波段还带不同单位的数值后续处理时还得自己记波段顺序非常被动。更稳的做法是在prm里显式指定OUTPUT_SDS_NAME或者用SPECTRAL_SUBSET精确控制。实践里我对地表温度产品会这样写prmINPUT_FILENAME OUTPUT_FILENAME RESAMPLING_TYPE NEAREST_NEIGHBOR OUTPUT_PROJECTION_TYPE UTM OUTPUT_PROJECTION_PARAMETERS 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 OUTPUT_PIXEL_SIZE 1000.00000 SPECTRAL_SUBSET ( 1 0 0 0 1 0 0 0 )SPECTRAL_SUBSET括号里的数字按SDS顺序排列1表示输出、0表示丢弃。上面的配置表示只保留第1个LST_Day_1km和第5个LST_Night_1kmSDS输出是一个双波段GeoTIFF波段顺序固定后续做白天地温序列和夜间地温序列就能直接按波段索引取数。如果你是做NDVI时序MOD13Q1里面SDS顺序大约是NDVI、EVI、VI_Quality、pixel_reliability等想同时保留NDVI和EVI就写成( 1 1 0 0 ... )。这里必须强调不同MODIS产品的SDS顺序并不是完全一致的同一个产品不同版本Collection 5 和 Collection 6也可能有差异。所以我每拿到一个新产品第一步就是用MRT GUI里的“打开HDF”看一下SDS列表或者用HDFView确认顺序然后才写SPECTRAL_SUBSET。凭记忆写括号里的01串翻车概率极高输出一个波段全对不上、数值范围也不对。4.2 投影参数选择UTM带号、Albers还是WGS84经纬度重采样算法怎么定OUTPUT_PROJECTION_TYPE是prm里最值得花时间的一项。MRT支持的类型包括UTM、Albers Conical Equal Area、Lambert Azimuthal Equal Area、Geographic经纬度等。选哪个取决于你的研究区形状和后续用途。如果是做省级或县级的地表温度产品我一般选UTM。中国从东到西跨了UTM 43到52区东部省份用WGS84/UTM 50N内蒙古中西部用UTM 49N或48N新疆要用UTM 45N-44N。UTM带号信息写在哪里写在OUTPUT_PROJECTION_PARAMETERS这一行里。MRT的规则是UTM类型的参数串中有一个数字表示带号通常放在靠前的位置。比如OUTPUT_PROJECTION_TYPE UTM OUTPUT_PROJECTION_PARAMETERS 49.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000第一个49代表UTM 49N带号后面参数全部置0。注意UTM带号和中央经线的对应关系是“带号×6-183”49N带范围是东经108到114度正好覆盖内蒙古中东部和京津冀以西。如果你的研究区跨了两个带不想分带处理可以不用UTM改用Albers等面积投影这样整个研究区在一个坐标系下、面积也不会变形。写Albers投影时参数行里要填中央经线、中央纬线、双标准纬线一般中央经线取研究区中心经度双标准纬线取研究区南北边界的1/6和5/6位置。这个可以直接用MRT GUI里的下拉框选Albers然后拖动参数GUI会帮你生成参数串再复制到prm里用。RESAMPLING_TYPE这一项我多说一句。默认的NEAREST_NEIGHBOR适合大多数MODIS产品因为它不会改变原始像元值分布对后续统计分析和数值定标友好。双线性插值适用于连续地表变量比如反射率但代价是会平滑边界、把无效值扩散到有效像元周边。三次卷积的视觉效果最好但计算最慢且对LST这种噪声较大的产品没有实际收益。还有一点如果你打算把重投影后的数据用来做像元尺度的时序分析不要用双线性因为每次处理同一个原始像元的插值权重不同会导致时间序列里出现和真实地表变化无关的抖动。5. MRT批处理避坑指南电脑运行不了mrt指令的五个常见原因5.1 现象双击bat一闪而过命令行却正常原因bat文件所在目录或输出路径里含有中文或空格for循环里的路径带引号后resample.exe读prm文件里的路径时出现解析错位。另一个常见原因是bat用了 无延迟变量 方式读errorlevel导致脚本逻辑混乱提前退出。解决所有目录统一改成英文路径bat脚本开头加一行cd /d %~dp0让脚本先切到自己所在目录调用resample.exe时路径%MRT_BIN%加双引号防止带空格路径被拆开。5.2 现象循环只处理了第一个文件或者完全不循环原因bat脚本里用了%f而不是%%f。在cmd命令行里执行for %f in (*.hdf) do ...是合法的但同样的语句写进.bat文件里必须改成%%f否则bat会把%f当环境变量展开成空值循环体执行时变量是空的。解决检查脚本里的for变量是不是双百分号。直接说吧我见过最多的翻车就是这一个。5.3 现象提示“计算机中丢失javaw.exe”或“指令引用的内存不能为read”原因MRT的运行依赖32位Java运行时环境新版64位JDK直接跑不了老MRT此外MRT的resample.exe本身是32位程序在64位Windows上调用32位Java时路径不一致就会出现这类提示。解决安装32位JDK 1.8x86版本装完把它的bin目录加到PATH最前面并在系统环境变量里新建JAVA_HOME指向32位JDK安装根目录。如果机器上同时有64位和32位Java确保resample.exe运行时PATH里先找到的是32位版本。建议用命令行执行java -version确认运行的是32位版本。5.4 现象输出GeoTIFF能打开但全黑或者数值范围完全不对原因OUTPUT_SDS_NAME写错或者SPECTRAL_SUBSET的01串顺序和实际HDF里的SDS顺序不一致。比如你想提取地表温度结果0/1配置把QC数据输出成了主波段QC是8位整型数值范围0-255全图偏黑看起来像全黑。另外如果直接查看15位整型的原始DN值而没有做尺度因子换算LST值范围是7500-13100对应开尔文×0.02在GIS里自动拉伸显示也会接近全黑或全白。解决先用MRT GUI把HDF打开一次看SDS列表的真实顺序和名字输出tif后在GIS里查看波段属性确认位深和有效值范围。别忘了LST数据要乘以0.02才是开尔文减去273.15才是摄氏度这是定标步骤MRT不会替你做。5.5 现象批处理跑完输出文件的时间戳不对或者有些文件根本没生成原因for循环里echo [Processing]之后立刻调resample.exeMRT处理单个文件耗时几十秒甚至几分钟bat脚本本身没有等进程结束就继续循环这种情况其实不会发生因为普通调用是阻塞式的。更常见的原因是bat脚本里用了start命令调用resample.exe导致脚本不等处理完成就进入下一轮循环多个MRT进程同时写同一个输出文件互相覆盖。解决不要在bat里用start调用resample.exe直接写完整路径加参数调用等它执行完再循环下一个。如果确实想并行处理建议分批开多个cmd窗口各跑各的文件夹而不是在同一个循环里用start。这些避坑经验总结成一句MRT批处理95%的问题不出在算法上而是出在Windows环境、bat语法和路径字符上。处理数据前先拿一个文件把整条链路跑通再扔给for循环能省下大量反复试错的时间。6. 结果验证这一关不能省从“能跑”到“结果可靠”的检查手段脚本跑完一堆GeoTIFF以后我会做的第一件事不是直接进建模而是拿两三个文件做交叉验证。做法分两步第一步是打开QGIS或ArcGIS Pro叠加同一地区的矢量边界检查影像范围和边界是不是对齐、有没有横向错位、像元大小是不是和预期一致。第二步是用Python快速抽检数值属性这里有个简单的脚本逻辑读一个输出tif统计无效值比例和有效值范围from osgeo import gdal import numpy as np ds gdal.Open(rE:\modis\out\MOD11A1.A2019001.h25v04_LST.tif) band ds.GetRasterBand(1) arr band.ReadAsArray() valid arr[arr 0] # MOD11A1中0是无效值 print(Shape:, arr.shape) print(Valid pixel count:, len(valid)) print(LST DN min/max:, valid.min(), valid.max()) print(LST Celsius range:, valid.min() * 0.02 - 273.15, valid.max() * 0.02 - 273.15)这段代码的逻辑是用GDAL打开输出tif读取第一波段为numpy数组把等于0的像元视为无效值剔除输出有效像元的最小/最大DN值再乘以0.02并转成摄氏度。做完这个抽检你就能判断MRT输出的数据能不能进时序分析——比如冬季地表温度的DN值对应摄氏度在-25到-5之间是正常的如果出来个80°C那一定是SDS提错了或者尺度因子填错了。我习惯在正式批量处理前拿一个文件跑通并验证确认无误后再提交全部任务处理完再从头尾各抽一个文件复查一遍数值和坐标范围。这是从早期被“全黑tif”坑过的血泪经验换来的习惯流程越固定翻车概率越小。MRT批处理这条链路跑通后收益是很稳定的——以后每次拿到新的MODIS批量数据改个路径、确认一下SDS顺序剩下的交给脚本就行。希望帮到你。本文还有配套的精品资源点击获取
返回列表