前段时间接手一个活儿，对方扔过来一块硬盘，里面躺着四百多景 MODIS 的 HDF 文件，要求统一转成 WGS84 下的 UTM 投影，再按研究区边界裁好。随手打开其中一个，投影信息里写着 Sinusoidal，坐标是 X/Y 的米制值，扔进常规 GIS 桌面软件能显示，可一旦叠上自己的矢量边界，位置偏得让人怀疑人生。这不是软件坏了，正弦投影本来就不是给人直接用的——它把全球切成了一格一格的瓦片，每块瓦片自带一套局部坐标，跨瓦片拼接时要先还原到一个统一的平面坐标系里。这时候绕不开一个老工具：MRT（MODIS Reprojection Tool）。它年纪不小，界面朴素，命令行参数也不友好，但在“重投影 + 批量”这个组合上，很长时间里没有更省事的替代品。这篇就聊聊怎么用 MRT 批处理 MODIS 数据，从参数文件的字段含义，到批量脚本怎么组织，再到那些跑起来才发现的坑。

1. 正弦投影瓦片为什么不能直接用：MRT的存在理由与当下的替代关系

1.1 从镜头到瓦片：MODIS的正弦网格是怎么切出来的

要理解 MRT 存在的意义，得先明白手上的 HDF 文件到底是什么形态。MODIS 的观测数据并不是按经纬度规则排列的，它先被投影到一个以地球为中心的**正弦投影（Sinusoidal）**平面上，然后切成 36 列 × 18 行的网格，每一格叫一个瓦片（tile），编号形如 h26v05。h 表示水平方向第 26 列，v 表示垂直方向第 5 行。每块瓦片覆盖大约 10° × 10° 的经纬范围，内部是 1200 × 1200 或 2400 × 2400 的规则栅格。

这么设计的初衷很实在：正弦投影在面积上等积，全球任何位置的像元代表的地表面积基本一致，做区域统计、算植被指数、反演地表温度时不会因为高纬度拉伸而失真。代价是坐标不是经纬度，而是一套以米为单位的平面坐标，原点在全球网格的中心附近。你拿到的文件里，每个像元对应的 X/Y 值只是相对某块瓦片的局部值，丢进需要经纬度的工具链里就没法用。

更麻烦的是跨瓦片的问题。研究区稍微大一点，就要用好几块瓦片拼起来，但不同瓦片的局部坐标原点不一样，直接叠加会出现接缝错位。所以流程上必须有一步把数据拉回到一个统一的地理或投影坐标系里，这一步就是重投影，MRT 干的头号活儿就是这个。

1.2 MRT能做的事和它做不到的事

MRT 的功能其实很聚焦，就是围绕 MODIS 的 HDF-EOS 数据做三件事：重投影、重采样、子集提取。它能把正弦投影的数据转到 UTM、Albers、Lambert、Geographic 等二十来种投影上，能按经纬度或投影坐标裁出一个矩形范围，能只保留需要的波段，输出格式可以是 GeoTIFF、HDF-EOS 或者原始二进制。配套的还有一个mrtmosaic工具，专门把多块瓦片拼成一块。

但它的边界也很清楚。MRT 只认HDF4-EOS 封装的文件，如果你的数据已经换成了 HDF5 封装，它会直接报错退出，连个像样的提示都不给。它不处理质量波段与主波段的联动，子集裁剪时你要自己保证 SCI（state/quality）波段和数值波段用的是同一套范围。它也不做投影后的重采样精度优化，双线性就是双线性，不会帮你考虑像元中心偏移。

还有一个现实问题：这个工具停留在 HDF4 时代，官方早已停止更新。所以现在做 MODIS 处理，实际路线往往是分层的——如果手上是传统的 MODIS 陆地产品（大部分仍然是 HDF4-EOS），MRT 依然能用；如果数据格式已经变了，或者要做更复杂的重投影组合，就得换用 GDAL 的gdalwarp，它对 HDF4 的支持依赖子数据集（subdataset）路径，写法是HDF4_EOS:EOS_GRID:文件名:网格名:波段名，灵活度高但参数更多。

1.3 什么情况下不该用MRT

不是所有场景都值得上 MRT。如果只是单景、单波段、一次性的处理，直接在 GDAL 里一行gdalwarp就能解决，没必要为了生成一个 prm 文件绕一圈。如果研究区只落在一块瓦片内部、而且后续分析对投影精度要求不高，甚至可以用 HDF 自带的经纬度子数据集做查找表，跳过重投影直接采样，速度快好几倍。

真正适合 MRT 的，是大批量、参数统一、需要可重复的场景。比如几十上百景同一产品的时序数据，全部要转到同一个投影、同一个范围、同一个像元大小，这时候把参数固化成模板文件，批处理跑一遍，比每次手点 GUI 靠谱得多。反过来说，如果你的每景数据参数都不一样，或者需要做自适应裁剪，MRT 的模板机制反而会成为负担。

提示：判断标准很简单——如果需要处理的文件数超过 10 个，且投影参数完全相同，就值得写脚本；低于这个数量，手工操作反而更快。

2. 把GUI里点出来的东西固化成prm文件：字段逐项拆解

2.1 一个能跑通的prm长什么样

MRT 的核心是参数文件（.prm）。GUI 界面上点的每一个选项，最终都会落成这个纯文本文件里的一行键值对。批处理的关键，就是把这一份文件当成模板，反复替换其中随文件变化的部分。

一份最小可用的 prm 大致长这样：

INPUT_FILENAME = D:\modis\MOD09A1.A2020001.h26v05.006.hdf SPECTRAL_SUBSET = ( 1 0 1 0 0 0 0 0 0 0 0 0 0 ) SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG SPATIAL_SUBSET_UL_CORNER = ( 40.0 100.0 ) SPATIAL_SUBSET_LR_CORNER = ( 20.0 120.0 ) OUTPUT_FILENAME = D:\out\MOD09A1.A2020001.h26v05.tif OUTPUT_FILETYPE = GEOTIFF RESAMPLING_TYPE = NEAREST_NEIGHBOR OUTPUT_PROJECTION_TYPE = UTM OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ) DATUM = WGS84 UTM_ZONE = 49 OUTPUT_PIXEL_SIZE = 500

逐项说一下。INPUT_FILENAME和OUTPUT_FILENAME是唯一需要随文件变化的两行，其余在整批任务里保持不变。SPECTRAL_SUBSET是一串 0/1，位置对应波段序号，1 表示保留。MODIS 产品波段数不一样，反射率产品有 7 个主波段加若干质量波段，最多能到十几个，所以要按实际产品的波段列表来填，长度对不上会报错。

SPATIAL_SUBSET_TYPE有两种取值：INPUT_LAT_LONG表示用经纬度给范围，OUTPUT_PROJ_COORDS表示用输出投影的平面坐标给范围。前者更通用，因为经纬度范围跟投影无关；后者在做精确对齐时更准，但要先算好四个角点的平面坐标。

OUTPUT_FILETYPE一般填GEOTIFF，需要保留 HDF 结构就填HDF-EOS，做自定义解析可以填RAW_BINARY。OUTPUT_PIXEL_SIZE是输出像元大小，单位和输出投影一致，UTM 下是米。

2.2 投影参数那一串0到底该填什么：UTM与Albers的取值推演

OUTPUT_PROJECTION_PARAMETERS那一长串 0 是最容易让人懵的地方。它其实是一个按投影类型解释的位置参数数组，不同投影用到的位置不一样，用不到的填 0。

选 UTM 时，这串参数基本全是 0，因为 UTM 的关键信息在UTM_ZONE和DATUM两个独立字段里。中央经线由分带号自动推出，公式是中央经线 = 带号 × 6 - 183。带号 49 对应的中央经线就是 49 × 6 - 183 = 111°。这个公式值得记一下，后面判断研究区该用哪个带的时候天天用得到。

选 Albers 等积圆锥投影时，参数数组就派上用场了。一般需要填：中央经线、原点纬度、第一标准纬线、第二标准纬线。假设研究区在东经 95° 到 115° 之间、北纬 20° 到 45° 之间，一套常见取法是中央经线取 105°（范围中点），原点纬度取 0°，两条标准纬线分别取 25° 和 47°（大约覆盖范围上下各外扩一点）。这样做的理由是：两条标准纬线之间的区域变形最小，把它们放在研究区南北边界略往外，能让中间地带的面积误差压到千分之几。

注意：prm 里参数数组的具体顺序，不同版本的模板可能略有差异，最稳的做法是在 GUI 里配好一次，让它生成一份 prm，然后照着那份文件的顺序抄。我见过不止一次因为顺序错位导致输出投影完全跑偏的情况，而且工具不会报错，只是默默给你一个错误的坐标系。

2.3 重采样方法不是随便选：分类数据与连续量的分界线

RESAMPLING_TYPE有三个可选值：NEAREST_NEIGHBOR、BILINEAR、CUBIC_CONVOLUTION。选错的后果比想象中严重。

最近邻是唯一能保持原始像元值不被改写的方案。它直接取最近的那个源像元的值，不做任何加权平均。所以只要你的数据是分类量——土地覆盖类型、质量标记、火灾掩膜、云掩膜——就必须用它。用双线性处理分类数据，会在类别边界上生成一堆原本不存在的中间值，比如类别 1 和类别 5 之间冒出个 3，而这个 3 在原始分类体系里可能代表完全不同的地物。

双线性适合连续量，比如反射率、地表温度、叶面积指数。它取周围四个像元按距离加权，结果更平滑，视觉上也更自然。代价是会轻微削弱极值，对于做异常检测的场景要慎用。

三次卷积更平滑，但会引入过冲，在数值边界附近可能产生超出原始值域的结果，比如 NDVI 本来最大 1.0，处理后冒出 1.05。除非明确需要高质量的重采样，一般不用。

还有一点容易被忽略：填充值在重采样时会参与计算。如果源数据边缘是填充值（比如 -28672），双线性会把这些值混进来，导致边缘出现一条诡异的过渡带。处理办法是在重投影前先把填充值掩掉，或者在输出后用阈值过滤。

2.4 空间子集与波段子集：少读一点数据就少等一会儿

空间子集看起来只是个裁剪功能，实际上对性能影响很大。MRT 在处理时是先按子集范围定位源数据块，再只读需要的部分，所以给一个紧凑的范围能显著减少 IO 和内存占用。我试过对同一景数据分别做全图和四分之一范围的裁剪，后者耗时大约只有前者的三分之一。

范围给的技巧是留一点余量。因为源数据是正弦投影、输出是 UTM，两个坐标系下的矩形不是简单对应的，边界处会有几个像元的错位。如果范围卡得死死的，输出边缘可能出现一条没有数据的空带。稳妥的做法是四个方向各外扩 0.1° 到 0.2°，输出后再裁一次。

波段子集的作用类似。MODIS 反射率产品动辄十几个波段，如果只算 NDVI，只需要红光和近红外两个波段，其余全填 0。少读一个波段就少一次解码，批量跑起来差别很明显。这里有个细节：SPECTRAL_SUBSET的括号里必须有空格分隔，写成( 1 0 1 0 )而不是(1,0,1,0)，后者会解析失败。

3. 批量执行的三种组织方式：模板替换、shell循环、Python调度

3.1 模板加占位符：最土也最不容易翻车

批处理最直接的做法，是准备一份带占位符的模板，然后对每个输入文件做一次文本替换，生成临时 prm，再调用命令行。模板里把两个文件名换成特殊标记：

INPUT_FILENAME = __IN__ OUTPUT_FILENAME = __OUT__ SPECTRAL_SUBSET = ( 1 0 1 0 0 0 0 0 0 0 0 0 0 ) ...

替换的时候要注意两件事。一是路径分隔符，Windows 下用反斜杠，Linux 下用正斜杠，混用会导致文件找不到。二是行尾符，prm 是纯文本，如果模板是在 Windows 下编辑的、拿到 Linux 上跑，行尾的\r可能被当成文件名的一部分，表现为“文件明明存在却提示打不开”。用dos2unix转一下，或者在编辑器里显式存成 Unix 格式。

这个办法的好处是完全可控，出了问题一眼能看出是哪一行不对。坏处是模板和脚本分离，改参数要动两个地方。

3.2 Windows批处理和Linux shell里几个必须注意的差异

Windows 下的循环大概是这样：

@echo off setlocal enabledelayedexpansion for %%f in (D:\modis\*.hdf) do ( set "IN=%%f" set "NAME=%%~nf" set "OUT=D:\out\!NAME!.tif" powershell -Command "(Get-Content template.prm) -replace '__IN__','!IN!' -replace '__OUT__','!OUT!' | Set-Content tmp.prm" resample -p tmp.prm if errorlevel 1 echo FAILED: !NAME! >> error.log )

几个坑点。%%f在批处理文件里是双百分号，在命令行里直接敲是单百分号，这个搞混了会一直提示语法错误。延迟变量展开要开enabledelayedexpansion，否则循环里!NAME!拿到的永远是第一次的值。路径里有空格的话，所有变量引用都要加引号。

Linux 下相对简单：

for f in /data/modis/*.hdf; do name=$(basename "$f" .hdf) sed -e "s|__IN__|$f|" -e "s|__OUT__|/data/out/$name.tif|" template.prm > tmp.prm resample -p tmp.prm || echo "FAILED: $name" >> error.log done

sed的替换分隔符建议用竖线而不是斜杠，因为路径里全是斜杠，用斜杠要转义到怀疑人生。文件名里有空格的话，$f一定要加双引号。

3.3 用Python把prm生成和进程调用串起来

文件数量上到几百之后，shell 的容错能力就不太够了，这时候用 Python 更省心。核心思路是：读模板、替换、写临时文件、调用外部进程、检查返回码。

import subprocess, pathlib, concurrent.futures tpl = pathlib.Path("template.prm").read_text(encoding="utf-8") src_dir = pathlib.Path(r"D:\modis") out_dir = pathlib.Path(r"D:\out") out_dir.mkdir(exist_ok=True) def run_one(hdf): out = out_dir / (hdf.stem + ".tif") if out.exists(): return ("skip", hdf.name) prm = tpl.replace("__IN__", str(hdf)).replace("__OUT__", str(out)) tmp = out_dir / (hdf.stem + ".prm") tmp.write_text(prm, encoding="utf-8") r = subprocess.run(["resample", "-p", str(tmp)], capture_output=True, text=True) if r.returncode != 0: return ("fail", hdf.name + " | " + r.stderr.strip()[:200]) return ("ok", hdf.name) files = sorted(src_dir.glob("*.hdf")) with concurrent.futures.ThreadPoolExecutor(max_workers=4) as ex: for status, msg in ex.map(run_one, files): print(status, msg)

这段有几处是刻意的设计。out.exists()的判断实现了断点续跑，四百景跑到一半崩了，重跑时已经完成的会直接跳过，不用从头来。capture_output=True把输出抓下来，失败时能记录错误摘要，而不是只看到一个返回码。用线程池而不是进程池，是因为真正的计算在外部进程里，Python 这一层只负责调度，线程足够。

提示：max_workers不要设太大。MRT 是单线程程序，每个实例都会把整景数据读进内存，四到六个并发在普通工作站上比较合适，再多就是拿内存换并不存在的速度。

3.4 拼接与重投影的先后顺序，决定了你要不要重算两遍

如果研究区跨了多块瓦片，就要面对顺序问题。MRT 配套的mrtmosaic只能处理同投影、同网格的输入，也就是说它拼接的是正弦投影下的原始瓦片，输出仍然是正弦投影。所以标准流程是：

1. 用mrtmosaic把同一天的多块瓦片拼成一块；
2. 再对拼好的文件跑resample做重投影。

mrtmosaic的命令行大概长这样：

mrtmosaic -i tilelist.txt -o mosaic.hdf -s "1 0 1 0 0 0 0 0 0 0 0 0 0"

-i指向一个文本文件，每行一个 HDF 路径；-o是输出；-s是波段子集，格式和 prm 里的SPECTRAL_SUBSET一致，但这里是用空格分隔的字符串，不带括号。

反过来先重投影再拼接也不是不行，那就得换 GDAL 的gdalbuildvrt加gdal_translate，好处是可以顺便做金字塔，坏处是每块瓦片都要单独重投影一次，边缘地带可能因为子集范围不一致而出现接缝。我的经验是：同一产品的时序数据，先拼后投；不同来源、不同分辨率的数据混在一起，先投后拼更可控。

4. 运行即失败、输出全黑、位置错位：三类故障的排查链路

4.1 双击没反应、命令行一闪而过：先看日志落在哪

MRT 的命令行版本在 Windows 下双击运行会瞬间关窗，什么提示都看不到。这不是程序崩溃，是它把信息输出到标准输出，窗口一关内容就没了。正确的排查方式是在命令行窗口里手动执行，或者把输出重定向到文件：

resample -p test.prm > run.log 2>&1

日志里最常见的三类信息：Unable to open input file是路径问题，Invalid projection parameter是参数数组长度或取值不对，Unsupported SDS是波段序号超出了文件实际波段数。

如果日志是空的、程序直接退出，大概率是运行库缺失。MRT 依赖 HDF4 的动态库，这些库要放在可执行文件同目录或者系统路径里。换了台机器就跑不起来，八成是这个原因。

4.2 输出影像一片黑或者全是填充值

这种情况通常是空间子集给错了坐标系。如果SPATIAL_SUBSET_TYPE写的是INPUT_LAT_LONG，但四个角点填的是投影坐标的数值（比如几十万米），工具会认为你要的范围在数据之外，输出就是全填充。反过来，写了OUTPUT_PROJ_COORDS却填了经纬度，结果也一样。

排查方法很直接：先用gdalinfo看一下源数据的经纬度范围，确认你给的角点确实落在里面。注意角点的顺序是左上和右下，纬度是北正南负，经度是东正西负。有人习惯性把纬度写成 20 到 40（南纬到北纬），在北半球研究区就会得到空结果。

另一种可能是波段子集全填了 0。SPECTRAL_SUBSET = ( 0 0 0 0 0 )意味着一个波段都不要，输出自然是空的，但工具不会拦你。

4.3 位置偏了几百米：分带、基准面和像元中心的问题

输出能打开、有数据，但叠上矢量边界发现整体偏移，这种问题最费时间。按可能性从高到低排：

分带号写错。UTM 带号差 1，中央经线就差 6°，位置偏移会到几百公里级别，很容易发现。但如果研究区正好在分带边界上，带号对了、中央经线选错，偏移就是几十公里。

基准面不匹配。DATUM填WGS84而参考矢量是另一套基准面，偏移从几十米到几百米不等。国内的业务数据经常遇到这个问题，需要做基准面转换而不是简单改个字段。

像元中心与角点约定。GeoTIFF 的坐标默认指向像元左上角，而有些软件按像元中心解释。差半个像元，500 米分辨率下就是 250 米的偏移。这个偏移是系统性的，整幅图往一个方向平移，很好辨认。

重采样后的范围外扩。前面提到子集范围要留余量，如果你留了余量又没在后续裁掉，输出的实际范围会比预期大一圈，做叠加分析时会被误认为错位。

排查建议：拿一个已知坐标的地面控制点，在输出影像上找对应像元，量一下偏差的方向和量级。方向和量级能直接指向上面四类原因中的哪一类。

4.4 路径里的空格、中文和超长文件名

这个问题看似低级，但踩的人非常多。prm 是纯文本解析的，路径里有空格时不会自动加引号，解析器会在空格处截断，于是它去找一个不存在的文件。中文路径的问题更隐蔽，取决于系统区域设置和文件编码，有时能读有时不能，换台机器就变了。

最省事的办法是把输入输出都放在一级短路径的无空格目录下，比如D:\mrt_in和D:\mrt_out，处理完了再挪走。文件名本身也别太长，MODIS 的原始文件名就有四十多个字符，加上临时后缀容易触及路径长度上限。

注意：不要在中文用户名下的桌面或文档目录里直接跑，那种路径里既有中文又有空格，还特别长，三个问题一次性凑齐。

5. 几百景数据无人值守跑完：并发、校验与后续衔接

5.1 并发怎么开才不把内存撑爆

MRT 单个实例处理一景 500 米分辨率的反射率产品，峰值内存大概在几百 MB 到 1 GB 之间，取决于子集范围大小。八核机器开八个并发，内存占用就可能到 6 到 8 GB，再加上系统本身的开销，16 GB 的机器会开始换页，速度反而下降。

我的经验值是并发数 = 物理核心数的一半，且不超过 6。跑之前先单独测一景，看耗时和内存峰值，再决定开几个。如果数据分布在机械硬盘上，并发数还要再降，因为多个进程同时随机读会让磁头来回跳，IO 成为瓶颈。

另一个技巧是按瓦片分组。同一个瓦片号的数据放在一起处理，文件的物理存储位置接近，读取效率更高。如果数据是按日期命名的，可以先按瓦片号建子目录做一次整理。

5.2 不开软件也能判断产物对不对

批量跑完几百个文件，不可能一个个打开看。可以写一个简单的校验脚本，用 GDAL 读元数据做批量检查：

from osgeo import gdal import pathlib for tif in pathlib.Path(r"D:\out").glob("*.tif"): ds = gdal.Open(str(tif)) if ds is None: print("UNREADABLE", tif.name); continue gt = ds.GetGeoTransform() band = ds.GetRasterBand(1) stats = band.GetStatistics(True, True) print(tif.name, ds.RasterXSize, ds.RasterYSize, round(gt[0], 1), round(gt[3], 1), round(stats[0], 2), round(stats[1], 2))

重点看三个指标。尺寸是否和预期一致，差太多说明子集范围写错了。左上角坐标是否在合理区间，全是一堆零或者极大值说明投影没生效。最小值和最大值是否落在物理合理范围内，比如反射率应该在 0 到 10000 之间（整数存储），如果出现 -28672 占了很大比例，说明填充值没处理好。

还可以补一条：文件体积。同一产品的不同日期，文件大小应该在一个窄范围内波动，如果某个文件只有几百字节或者比平均值大十倍，基本可以判定有问题。

5.3 输出之后交给GDAL做的那几件事

MRT 输出的 GeoTIFF 已经带了投影信息，但通常还缺几样东西，用 GDAL 补上很顺手。

建金字塔，让后续在 GIS 里缩放时不用实时重采样：

gdaladdo -r average out.tif 2 4 8 16

转压缩格式，四百景 GeoTIFF 占的空间相当可观，转成带压缩的格式能省一半以上：

gdal_translate -co COMPRESS=DEFLATE -co TILED=YES -co PREDICTOR=2 in.tif out.tif

PREDICTOR=2对整数型反射率数据效果明显，对浮点数据用PREDICTOR=3。

统一裁剪，如果前面为了保险留了余量，这时候用gdalwarp -te按精确范围裁掉：

gdalwarp -te 100.0 20.0 120.0 40.0 -te_srs EPSG:4326 in.tif out.tif

注意-te的坐标顺序是 xmin ymin xmax ymax，不要和 prm 里的左上/右下搞混。

6. 几个反复踩到的细节：值域、填充值和文件命名

6.1 填充值在重投影后会被改写

MODIS 产品用特定的整数表示无效值，反射率产品常见的是 -28672，地表温度产品是 0 或特定标记值。经过重投影和重采样之后，这些填充值不一定保持原样。最近邻会原样带过去，双线性会把填充值和有效值混在一起算出个中间数，看起来既不像填充值也不像正常值。

处理办法有两个。一是在 prm 里把空间子集卡得紧一些，尽量避开数据边缘的填充区；二是在输出后用gdal_calc或 numpy 做一次掩膜，把超出合理值域的像元统一设成 NoData。后者更彻底，但要多跑一遍。

值域本身也要注意。很多 MODIS 产品用整数存储并附带缩放因子，比如反射率要乘 0.0001 才是真实值。MRT 不会自动应用缩放因子，输出的还是原始整数。如果你在后续分析里忘了这一步，算出来的 NDVI 会是几千而不是 -1 到 1 之间。

6.2 命名里带上瓦片号和日期，后面能省很多事

批处理最容易出的低级错误是输出文件名冲突。不同日期的数据如果只用产品名命名，后一个会覆盖前一个。原始 HDF 文件名里已经包含了产品名、日期、瓦片号、处理批次，最省事的做法是直接把原文件名的主体部分搬过来，只换扩展名。

out_name = hdf.stem + ".tif"

这样输出就是MOD09A1.A2020001.h26v05.006.tif，日期和瓦片一目了然，排序后天然按时间排列。如果嫌长，至少保留日期和瓦片号两段，这两个是后续做时序分析时筛选文件的关键。

还有一个建议：把失败的记录单独写到一个日志文件里，格式包含文件名和错误摘要。跑完之后直接对着日志重跑失败的那几个，比全量重跑快得多。前面 Python 脚本里out.exists()的跳过逻辑配合这个日志，基本能实现无人值守。

6.3 工具链老化带来的取舍

MRT 是个老工具，用起来确实有这样那样的不便，但它在参数确定性上有优势——同样的 prm 文件，在任何机器上跑出来的结果逐像元一致，这对需要复现的研究场景很重要。GDAL 的gdalwarp更现代、更灵活，但默认参数下不同版本的重采样实现可能有细微差异。

我的实际取舍是：格式还是 HDF4-EOS、流程固定、需要可复现的批量任务，继续用 MRT；格式已经变了、需要复杂重投影组合、或者要和别的数据源对齐的任务，用 GDAL。两条路线不冲突，同一个项目里混着用也很正常——前半段用 MRT 批量重投影，后半段用 GDAL 做镶嵌、裁剪、统计。

真要说经验，最重要的一条是：别等跑完四百景才发现参数错了。先拿三景不同瓦片、不同季节的数据跑通全流程，把输出打开看一眼投影、范围、值域，确认没问题再全量开跑。这个前置验证花半小时，能省掉后面一整夜的返工。