MODIS数据下载、HDF4定标与MRT重投影实战

简介：一份围绕MODIS数据入门与应用的文档，面向遥感、地理信息、环境监测等方向的学生和初级从业者，帮助梳理从传感器背景到数据获取与处理的关键知识。压缩包内仅含1个doc文件，约503KB，体量轻便，适合快速浏览和按主题查阅。内容涵盖NASA地球科学事业与EOS计划、Terra/Aqua卫星分工、MODIS 36个波段及空间分辨率等技术指标，并延伸到产品下载流程、HDF格式读取，以及ENVI、ArcGIS等工具的处理思路。文中还涉及预处理、校正、合成与分析等环节，可辅助理解植被指数、地表温度等常见应用。已有304人学习，适合用作MODIS入门参考和课堂/项目前的基础资料。

1. 600MB 的 HDF 拖进软件只显示一片灰：MODIS数据到底该从哪儿下手

第一次接触 MODIS数据的人，几乎都会卡在同一个画面上：文件下下来了，双击打不开，拖进 GIS 软件只弹出一个灰蒙蒙的空框，图层列表里全是EV_1KM_RefSB、sur_refl_b01、250m_16_days_NDVI这类没头没尾的名字。这不是文件坏了，而是 MODIS 的组织方式和普通影像完全不同——一个.hdf里塞着几十个科学数据集、属性表和定标参数，波段被拆成一个个“抽屉”，你不点名去取，它就什么都不给你看。

这份资料覆盖的范围比一个 HDF 文件大得多：从 Terra/Aqua 双星的轨道设计、36 个波段的物理含义，到 L0~L4 五级产品体系、44 种标准产品的命名规则，再到检索下单、批量下载、MRT 重投影与 Albers 参数配置，最后落到定标和有效值处理。适合第一次做植被指数、地表温度或土地覆被时间序列的人，也适合被 HDF4 折磨过一轮、想把下载与批处理固化成脚本的老手。往下读之前先记住一句话：MODIS 的坑，八成不在算法，而在“你下错产品”和“你没做定标”。

2. 36 个波段与 44 种产品：先定产品谱系再谈下载

绝大多数翻车都不是从代码开始的，是从检索页开始的。选了MOD021KM想做 NDVI，或者下了MOD13Q1想反演气溶胶，后面所有工作都是白费。这一章要做的事只有一件：在看到下载按钮之前，先能说清自己要的是哪一级、哪个短名、哪个分辨率桶。

2.1 三个分辨率桶与两类波段：为什么不能“全波段一把梭”

MODIS 的 36 个波段不是平均分布的，它被硬切成了三个分辨率桶：250m 只有波段 1、2；500m 是波段 3~7；剩下 8~36 全部是 1km。这个设计直接决定了很多产品为什么长成现在这样——MOD13Q1的 NDVI 是 250m，因为它只依赖波段 1、2；MOD11A2的地表温度是 1km，因为热红外波段本来就只有 1km。

波段还分两类：1~19 是反射光谱段，探测的是地表和云对太阳短波辐射的反射，物理量是反射率；20~36 是发射光谱段，探测的是地物自身的热辐射，物理量是亮温。这两类数据的定标方式、有效值范围、甚至可视化方式都不一样。把波段 31 当反射率去算归一化指数，结果会是一张全黑或者全白的图，很多人到这一步还以为是自己算错了。

常用波段的用途可以对着下面这张表看，不用背全 36 个，记住七个就够日常用：

表里能读出一个规律：凡是做植被、土地覆盖、积雪的，基本都落在 1~7；凡是做温度、火点、云的，基本都落在 20 以后。这也解释了为什么陆地产品多是 250m/500m，大气和海洋产品多是 1km。

提示：波段 13 有 13h 和 13l 两个（高增益与低增益），波段 14 同理。做叶绿素荧光时选错增益，数值会差一个量级，且不会有任何报错。

2.2 L0 到 L4：44 种标准产品里真正会用到的十来个

MODIS 标准产品按处理深度分五级。L0 是地面站直接收到的原始比特流；L1A 是重建、时间配准、加了辅助数据之后的产物；L1B 是做过辐射定标的标准产品，也就是常说的MOD02；L2 以上是在 L1B 基础上做应用处理得到的。这里有个容易被忽略的关键差别：L1B 和 L2 是沿轨条带（swath），每一景跟着卫星走，是斜的；L3 和 L4 是投影到正弦格网上的瓦片（tile），是正的、可按瓦片号索引的。

这个差别决定了要不要做几何校正。拿 L1B 直接出图，你会得到一条歪着的、带弯曲边缘的影像，必须先配合MOD03地理定位文件做几何纠正；而拿MOD13Q1这类 L3 瓦片，它已经在正弦投影下规整好了，直接重投影就行。新手最常走错的路径，就是想做时间序列却从 L1B 开始，结果光是几何校正和拼接就耗掉一周。

下面这张表是实际项目里出现频率最高的十来个短名，覆盖了陆地、大气、海洋三个方向：

前缀是本篇最该刻进脑子的一条规则：MOD是 Terra（上午星），MYD是 Aqua（下午星），MCD是两颗星融合。同一个产品名换成MYD前缀，过境时间就从上午 10:30 变成下午 13:30。做物候或者能量平衡时，上午星和下午星的地表温度能差好几度，混着用会得到一条根本解释不通的曲线。

选型逻辑可以归纳成一句话：做算法反演、需要原始辐射亮度或自定义指数，走 L1B + MOD03；做区域制图、时间序列、统计分析，走 L3/L4 瓦片产品。前者自由度高但工作量大，后者开箱即用但只有官方定好的那几十个变量。

2.3 文件名逐段拆解：从 MOD13Q1.A2020032.h25v04.006 读出全部元数据

MODIS 的文件名本身就是一份元数据，只是被压缩成了一串点分字段。以MOD13Q1.A2020032.h25v04.006.2020052184512.hdf为例：

• MOD13Q1：产品短名，MOD 是 Terra，13 表示植被指数，Q 表示 250m 分辨率档，1 表示产品系列内的编号。
• A2020032：A加四位年份加三位年积日。2020 年第 32 天，也就是 2 月 1 日。
• h25v04：正弦投影下的瓦片行列号，水平方向 h00~h35，垂直方向 v00~v17，每块在赤道处约 10°×10°。
• 006：Collection 版本号。不同版本之间的定标系数和算法可能改过，跨版本拼时间序列会引入系统性偏差。
• 2020052184512：生产时间，格式是yyyydddhhmmss，即 2020 年第 052 天 18 时 45 分 12 秒。
• .hdf：文件格式。

L1B 的文件名结构略有不同，会多出一段采集时刻，比如MOD021KM.A2020032.0345.061.2020032213456.hdf，中间的0345就是 UTC 03:45 的过境时刻，这个字段在做多星匹配时非常有用。

手工数位太容易出错，我一般直接写个函数把文件名解析成结构化字典，后面做批量索引、按日期筛选、按版本分组都靠它：

import re from datetime import datetime, timedelta # 兼容 L1B（带时刻、无瓦片）和 L3（带瓦片、无时刻）两种命名 PATTERN = re.compile( r'^(?P<product>[A-Z]{3}\d{2}[A-Z0-9]{0,3})' r'\.A(?P<year>\d{4})(?P<doy>\d{3})' r'(?:\.(?P<utc_hhmm>\d{4}))?' r'(?:\.(?P<tile>h\d{2}v\d{2}))?' r'\.(?P<version>\d{3})' r'\.(?P<prod_time>\d{13})' r'\.(?P<ext>hdf)$' ) def parse_modis_name(name): m = PATTERN.match(name) if not m: return None d = m.groupdict() # 年积日转真实日期：1 月 1 日是第 1 天，所以要减 1 d['acq_date'] = datetime(int(d['year']), 1, 1) + timedelta(days=int(d['doy']) - 1) # 生产时间 yyyydddhhmmss p = d['prod_time'] d['prod_datetime'] = ( datetime(int(p[0:4]), 1, 1) + timedelta(days=int(p[4:7]) - 1) ).replace(hour=int(p[7:9]), minute=int(p[9:11]), second=int(p[11:13])) return d print(parse_modis_name('MOD13Q1.A2020032.h25v04.006.2020052184512.hdf')) print(parse_modis_name('MOD021KM.A2020032.0345.061.2020032213456.hdf'))

这段代码的关键在两个地方。一是正则里的两组可选字段都用(?:...)?包起来，因为 L1B 有时间无瓦片、L3 有瓦片无时间，如果不加可选，两类文件必然有一类匹配不上。二是年积日的换算要减 1，datetime(year, 1, 1) + timedelta(days=n-1)才等于第 n 天，写成days=n会整体偏移一天，这种错在按日聚合时特别隐蔽——数据一天不缺，但每一天都对不齐。

有了这个解析结果，按年份、按瓦片、按版本分组或者筛出某段日期内的文件，都只是普通的列表过滤，不用再靠文件名切片数位置。

2.4 用 gdalinfo 验货：下载完先看 SDS 清单再写脚本

文件到手之后不要急着写批处理，先花两分钟验一下货。GDAL 对 HDF4 的支持是开箱即用的，一条gdalinfo就能把文件里所有子数据集的完整路径打印出来：

# 列出 HDF4 文件内所有子数据集（SUBDATASET） gdalinfo MOD13Q1.A2020032.h25v04.006.2020052184512.hdf | grep -A2 SUBDATASET | head -40 # 只看某个具体子数据集的元信息：尺寸、投影、类型、NODATA gdalinfo 'HDF4_EOS:EOS_GRID:"MOD13Q1.A2020032.h25v04.006.2020052184512.hdf":MODIS_Grid_16DAY_250m_500m:250m_16_days_NDVI'

输出里最该关注四项：Size is给出行列数，正弦瓦片通常是 4800×4800 这个量级；Coordinate System is给的是正弦投影，带一个中央经线参数；Type=一般是Int16，说明它存的是缩放后的整数而不是浮点物理量；NoData Value=决定后续掩膜怎么写。

子数据集路径是有固定语法的：HDF4_EOS:EOS_GRID:"文件路径":网格名:数据集名。网格名和数据集名在不同产品、不同版本里会变，所以脚本里不要硬编码，正确的做法是先用GetSubDatasets()把清单取出来，再用关键字过滤：

from osgeo import gdal fp = 'MOD13Q1.A2020032.h25v04.006.2020052184512.hdf' ds = gdal.Open(fp) for path, desc in ds.GetSubDatasets(): print(desc)

desc一般是[4800x4800] 250m_16_days_NDVI MODIS_Grid_16DAY_250m_500m这种格式，用if 'NDVI' in desc就能精准定位，比手写完整路径稳得多。

注意：同一个产品在不同 Collection 版本下，网格名可能不同。脚本里写死MODIS_Grid_16DAY_250m_500m在换版本后会静默失败（返回 None），而不是报错。所有硬编码路径前，都建议先跑一遍清单打印。

3. 从检索页到硬盘：1B 与 L3 产品的两条下单路线和批量脚本

搞清楚要什么之后，下一个问题是去哪里拿。MODIS 的下载入口不是统一的，而是按产品类型分在两个体系里。走错门不是拿不到数据，而是根本搜不出结果，很多人到这一步会怀疑是自己网络或者账号有问题。

3.1 两类门户的分工：swath 产品走订单，tile 产品走检索

L1B 和 L2 级产品（MOD02、MOD03、MOD04、MOD06、MOD35等）沿轨条带组织，检索维度是时间加空间范围，走的是大气产品那条线。这类门户的检索界面通常要求你填起止时间、经纬度框，然后列出覆盖这个范围的若干条轨道段，勾选后提交订单。

L3/L4 产品（MOD09、MOD11、MOD13、MOD15、MOD17、MCD12Q1等）按正弦瓦片组织，检索维度是时间加瓦片号。这类门户一般提供地图点选或者经纬度反查瓦片号的功能，选中瓦片和时间范围后直接给出一批文件链接，不需要走订单流程。

两条路线的差别不只是界面。swath 产品的单景覆盖范围是一条窄带，一天之内同一地区可能被多条轨道部分覆盖，要做完整覆盖就得把相邻轨道段都下下来再拼接；瓦片产品是规整的 10°×10° 网格，同一瓦片同一天只有一景，取数逻辑简单得多。这也是为什么做区域时间序列的人都优先选瓦片产品。

实际检索时还有一个绕不开的问题：时区。门户上显示和填写的时间全是 UTC，北京时间比 UTC 快 8 小时。要拿到北京时间 1 月 16 日全天的数据，UTC 检索范围应该是 1 月 15 日 16:00 到 1 月 16 日 16:00，而不是 1 月 16 日 00:00 到 24:00。按后者填，会丢掉北京时间当天前 8 小时的数据，而且丢得毫无痕迹——文件数看着正常，只是早上那几景没了。

提示：做长时间序列时，我习惯把检索范围前后各多放一天，下载完再按解析出的采集时刻裁掉多余的。宁可多下几景，也不要在边界上少一天。

3.2 空间范围、瓦片与覆盖完整性

经纬度框看着简单，坑却不少。第一个坑是跨 180° 经线：如果研究区横跨日界线，用minLon=170, maxLon=-170这种写法，大部分检索接口会直接返回空结果。正确做法是拆成两个框分别检索，或者干脆改用瓦片号检索。

第二个坑是研究区跨瓦片。一个 10°×10° 的瓦片对中等尺度的研究区来说往往不够，比如经度 111°–120°、纬度 35°–43° 这样的范围，至少要覆盖 2×2 共四块瓦片，边缘还可能再搭上一块。少下一块的表现是拼接结果边缘缺一条，而缺的那条往往正好是研究区边界。

第三个坑是瓦片号手算。正弦投影下瓦片的经度跨度随纬度变化，赤道处是 10°，高纬度会变宽，所以用简单的等经纬度除法算出来的瓦片号经常差一格。可靠的做法是在检索页的地图上直接框选，让系统把覆盖的瓦片列出来，把这几个瓦片号记下来，后面所有检索复用同一组。

3.3 批量下载脚本：令牌认证与断点续传

门户上的“一键下载”适合拿几景试试水，真正做项目必须走脚本。现在的数据门户普遍用 Earthdata 账号登录，脚本里用账号密码硬编码既不安全也不稳定，标准做法是在账号设置里生成一个令牌（token），通过 HTTP 头或者 cookie 传给下载工具。

#!/usr/bin/env bash set -euo pipefail # 数据门户的归档根路径，具体值以检索页给出的下载链接前缀为准 BASE="https://<data-portal-host>/archive/allData" PRODUCT="MOD13Q1" VERSION="006" YEAR="2020" DOY="032" # 令牌放在环境变量里，不要写进脚本提交到仓库 TOKEN="${EARTHDATA_TOKEN:?请先导出 EARTHDATA_TOKEN}" # 断点续传 + 失败重试 + 限速，避免被服务端掐断 wget --header="Authorization: Bearer ${TOKEN}" \ --continue \ --tries=5 \ --waitretry=10 \ --timeout=60 \ --limit-rate=8m \ --no-check-certificate \ -P ./raw \ "${BASE}/${PRODUCT}/${VERSION}/${YEAR}/${DOY}/" # 只保留 hdf 主文件，清掉同目录下的 .xml 元数据（如果不需要溯源） find ./raw -name '*.hdf.xml' -delete echo "已下载 $(ls -1 ./raw/*.hdf | wc -l) 个 HDF 文件"

脚本里有几个参数值得单独说。--continue是断点续传的开关，MODIS 单景动辄几百兆，一次下完的概率并不高；--tries=5 --waitretry=10让它在遇到 5xx 时退避重试，而不是直接放弃；--limit-rate=8m是血泪经验，并发拉满容易被服务端限流，反而更慢，压到 8MB/s 左右反而稳定。

如果要下的是按订单生成的 swath 产品，脚本思路一样，只是 URL 列表要先从订单页导出成文本文件，再用wget -i list.txt批量拉。

# 从订单页导出的链接清单逐行下载，每行一个 URL while IFS= read -r url; do [ -z "$url" ] && continue wget --header="Authorization: Bearer ${EARTHDATA_TOKEN}" \ --continue --tries=5 -P ./raw "$url" done < order_links.txt

while read这种写法看起来笨，但比xargs -P的好处是失败可感知：某一行断了，日志里能定位到具体是哪个文件，重跑时--continue会接着下。并行下载在文件数上千时确实快，但一旦出现半截文件，排查成本远高于省下的那点时间。

3.4 下载完先做完整性检查

脚本跑完不算结束，还要验一遍。三个检查动作：文件数量对不对、文件大小有没有异常的、GDAL 能不能打开。

# 1) 逐个文件做 HDF4 结构校验，打不开的会报错 for f in ./raw/*.hdf; do gdalinfo "$f" > /dev/null 2>&1 || echo "损坏: $f" done # 2) 找出明显偏小的文件（半截下载） find ./raw -name '*.hdf' -size -1M -print # 3) 统计文件数，和检索页给出的数量对一下 ls -1 ./raw/*.hdf | wc -l

半截文件是批量下载里最阴的一种：它存在、它有大小、文件名也对，但gdalinfo打不开或者子数据集缺一半。如果不在入库前筛出来，等跑到批处理第 300 个文件时再崩，重跑成本会翻好几倍。我的习惯是下载完立刻跑一遍上面的循环，把坏文件挪到./bad目录，重新下一遍。

4. HDF4 拆箱：子数据集、定标系数与 MRT 重投影

文件完好、路径也对，接下来才是真正的处理。这一步要解决三件事：怎么把数据从 HDF 里取出来、怎么把整数 DN 变成有物理意义的反射率或亮温、怎么把正弦投影转成能和其他数据叠加的投影。

4.1 HDF4 与 HDF-EOS：一个文件里的三层结构

MODIS 用的是 HDF4 格式（部分新产品线已转向 HDF5 或 NetCDF，但经典产品仍是 HDF4）。HDF4 本身是个通用容器，NASA 在它上面加了一层约定，叫 HDF-EOS，用来规范地理数据的组织方式。所以一个 MODIS 文件里其实有三层东西：

第一层是科学数据集（SDS），就是真正的栅格数据，每个波段、每个指数都是独立的 SDS，有自己的尺寸、类型、填充值。第二层是属性（attributes），分全局属性和 SDS 属性，定标系数scale_factor、add_offset、有效范围valid_range、填充值_FillValue都藏在这里，不读属性直接拿数组算，结果一定是错的。第三层是 Vdata 和 Vgroup，存储地理定位查找表、质量标志的位定义等辅助信息，元数据浏览工具能看到，普通处理一般碰不到。

这解释了为什么不能把.hdf当.tif用。TIFF 一个文件基本对应一张图，HDF 一个文件是一个小型数据库，你必须先“查表”再“取数”。

4.2 三条读取路径：GDAL、pyhdf、h5py 的适用边界

读取 MODIS 的 HDF4，常见的有三条路，各有各的适用场景，选错了会浪费很多时间。

第一条是 GDAL，也是我最推荐的默认路径。它把 HDF4 的子数据集抽象成了虚拟文件路径，读出来直接就是带投影和地理变换的栅格，后面的重投影、裁剪、统计全都能用 GDAL 生态一条龙做掉。缺点是必须先知道子数据集全名，靠GetSubDatasets()枚举。

from osgeo import gdal import numpy as np fp = 'MOD13Q1.A2020032.h25v04.006.2020052184512.hdf' ds = gdal.Open(fp) # 从清单里按关键字找目标子数据集，避免硬编码网格名 target = [p for p, d in ds.GetSubDatasets() if '250m_16_days_NDVI' in d][0] band = gdal.Open(target).GetRasterBand(1) arr = band.ReadAsArray().astype('float32') scale = band.GetScale() or 1.0 # 常见为 1e-4 offset = band.GetOffset() or 0.0 nodata = band.GetNoDataValue() # 先掩膜再定标，顺序反了会把填充值也乘上缩放系数 arr[arr == nodata] = np.nan ndvi = arr * scale + offset print('范围:', np.nanmin(ndvi), np.nanmax(ndvi)) # 正常应落在 -0.2 ~ 1.0

这段代码有两个细节。一是band.GetScale()在某些版本或某些子数据集上返回None，所以要给默认值 1.0，否则整个数组会变成 NaN。二是掩膜必须在定标之前做，如果先乘 1e-4 再比较填充值，-3000已经变成了-0.3，和nodata对不上，等于没掩。

第二条是 pyhdf，HDF4 的原生 Python 绑定。它的优势是能拿到 GDAL 抽象掉的东西，比如 Vdata 里的质量标志位定义、完整的属性字典。做质量控制、精度验证、或者要读一些冷门产品时，绕不开它。

from pyhdf.SD import SD, SDC hdf = SD('MOD11A2.A2020032.h25v04.006.2020052184512.hdf', SDC.READ) print(hdf.datasets().keys()) # 所有 SDS 名称 sds = hdf.select('LST_Day_1km') data = sds[:] # 读出为 numpy 数组 attrs = sds.attributes() # 定标系数、有效范围都在这里 print(attrs.get('scale_factor'), attrs.get('valid_range'))

第三条是 h5py，这条路走不通，必须说清楚。h5py 只能读 HDF5，遇到 HDF4 会直接抛OSError: Unable to open file (file signature not found)。网上不少教程把两者混为一谈，新手照着敲完报错，还以为是文件损坏。判断方法很简单：gdalinfo能打开且显示HDF4_EOS就是 HDF4，用 GDAL 或 pyhdf；显示HDF5才是 h5py 的战场。

4.3 MRT 命令行：先拼接后重投影的固定动作

MODIS 的正弦投影和常见的等经纬度、Albers、UTM 都对不上，必须重投影才能和矢量边界、其他遥感数据叠加。官方工具 MRT（MODIS Reprojection Tool）是这条链路上最稳的一个，虽然有图形界面，但批量处理必须用命令行。

MRT 的核心逻辑是“先拼接、后投影”。多幅同一天、相邻瓦片的影像，mrtmosaic会先按原正弦投影拼成一整块，resample再做投影转换。这个顺序不能反，反过来先投影再拼接，接缝处会因为重采样误差出现明显的条带。

# 第一步：把同一天的四块瓦片列进清单，做拼接 cat > tile_list.txt <<'EOF' MOD13Q1.A2020032.h25v04.006.2020052184512.hdf MOD13Q1.A2020032.h26v04.006.2020052184512.hdf MOD13Q1.A2020032.h25v05.006.2020052184512.hdf MOD13Q1.A2020032.h26v05.006.2020052184512.hdf EOF # 光谱子集按需选，这里保留 3 个 SDS；顺序与文件内 SDS 顺序一致 mrtmosaic -i tile_list.txt -o mosaic_2020032.hdf -s "1 1 1" # 第二步：用参数文件做重投影 resample -p albers.prm -i mosaic_2020032.hdf -o mosaic_2020032_albers.tif

-s "1 1 1"这个参数容易被忽略，它的意思是“每个输入文件的第 1 个、第 1 个、第 1 个 SDS 参与拼接”，位置对应输出波段。如果拼的是MOD13Q1，它的 SDS 顺序里 NDVI 是第 1 个、EVI 是第 2 个，想要两个都保留就得写成-s "1 1 2 2"（四个文件各取前两个）。写错的后果不是报错，而是输出波段数和你预期的不一样。

albers.prm是纯文本的投影参数文件，内容如下：

INPUT_FILENAME = mosaic_2020032.hdf SPECTRAL_SUBSET = ( 1 1 1 1 ) SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG SPATIAL_SUBSET_UL_CORNER = ( 43.0 111.0 ) SPATIAL_SUBSET_LR_CORNER = ( 35.0 121.0 ) OUTPUT_FILENAME = mosaic_2020032_albers.tif RESAMPLING_TYPE = BILINEAR OUTPUT_PROJECTION_TYPE = ALBERS OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 25.0 47.0 105.0 0.0 0.0 0.0 ) DATUM = WGS84 OUTPUT_PIXEL_SIZE = 250

参数串的顺序是(SMajor SMinor STDPR1 STDPR2 CenMer OriginLat FE FN)。前两个填 0.0 时，MRT 会自动套用 WGS84 的长短半轴（6378137.0 和 6356752.314245 米），不用手打那一长串。STDPR1=25、STDPR2=47、CenMer=105是国内中低纬度区域制图最常用的一组，双标准纬线选在 25° 和 47°，中央经线 105°，这样全国范围内的面积变形能控制得比较均匀。如果是单个省份，标准纬线可以按该省南北跨度收窄，中央经线取省域中心经度，变形会更小。

RESAMPLING_TYPE选BILINEAR是连续量（NDVI、反射率、温度）的常规做法；如果是土地覆盖这类分类产品，必须换成NEAREST_NEIGHBOR，双线性会把类别码插值成中间值，MCD12Q1的类别 1 和 2 插出来变成 1.5，等于凭空造了个不存在的类型。

注意：SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG时，UL 是左上角、LR 是右下角，顺序写反（先 LR 后 UL）不会报错，但输出可能是一张空图或者尺寸为 0 的文件。

5. 五个翻车现场：MODIS 下载与处理里的避坑清单

这一章记录的都是真实踩过、且反复有人踩的坑。每一条按“现象—原因—解决”写，遇到对应症状可以直接对号入座。

5.1 现象：L1B 打开是一条歪的带子，边缘发黑；原因：漏下 MOD03；解决：配对下载并做几何校正

MOD021KM是沿轨条带数据，它的每个像元在文件里是按扫描顺序排列的，没有地理坐标，只有仪器坐标。想让它变成正射影像，必须配合MOD03地理定位文件，里面存着每个像元的经纬度查找表。漏下MOD03的表现就是：能打开、能出图、但影像斜着弯着，边缘一圈黑。

解决办法是在检索时把MOD03和MOD02按同一时刻成对下载，文件名里的时间字段（A2020032.0345）必须完全一致。做几何校正时，MRT 的 swath 模式或者 ENVI 的 GLT 校正都能用，前者更快，后者可控性更强。如果是做区域制图，其实更建议直接绕开 L1B，改用MOD09GA这类已经做过定位的 L2G 产品，省掉整个几何校正环节。

5.2 现象：时间序列中间莫名少几天；原因：按北京时间检索 UTC；解决：检索窗口整体前移 8 小时并加宽

按“2020-01-01 到 2020-12-31”去检索，拿到的实际是 UTC 的这一整年，对应北京时间 2020-01-01 08:00 到 2021-01-01 08:00。开头少了 8 小时，结尾多了 8 小时，中间还可能出现某天完全没有数据的情况——因为那颗星那天的过境时刻恰好跨了 UTC 日界，而你按 UTC 日聚合时把它归到了另一天。

解决办法是检索时把窗口前后各放宽一天，下载完用第 2 章那个解析函数把采集时刻提出来，统一换算成北京时间再做日聚合。另外，L3 产品的时间标签是合成周期的起始日，16 天合成的MOD13Q1标的是周期第一天，不是中间那天，做物候提取时这个偏移必须考虑进去，否则整条曲线会平移 8 天。

5.3 现象：脚本读文件报Unable to open file，或者返回None；原因：HDF4/HDF5 混淆、子数据集路径写错；解决：先枚举再取

三种情况会长得一模一样。第一种是拿 h5py 去读 HDF4 文件，报签名不匹配。第二种是gdal.Open()打开了主文件但忘了取子数据集，ReadAsArray()返回None，因为主文件本身不是栅格。第三种是子数据集路径里的网格名或数据集名拼错了一个字符，gdal.Open()静默返回None，不抛异常。

排查顺序建议固定下来：先用gdalinfo 文件名看它是 HDF4 还是 HDF5；再用GetSubDatasets()打印全部路径，从输出里复制而不是手打；最后确认路径里的引号是英文双引号，且文件路径部分用绝对路径。Windows 下的路径反斜杠还要转成正斜杠，否则 GDAL 的虚拟文件路径解析会失败。

5.4 现象：拼接后接缝处有条带，或者整块是 NoData；原因：先投影后拼接、或 UL/LR 顺序反了；解决：固定“先拼后投”顺序并检查角点

接缝条带的根因是重采样顺序。每块瓦片单独重投影，边缘像元在重采样时会因为邻域超出瓦片范围而被填成 NoData，拼起来就是一条缝。正确顺序是先mrtmosaic在正弦投影下拼成完整一块，再resample整体投影，接缝处的邻域信息是完整的。

整块 NoData 则多半是空间子集写错了。SPATIAL_SUBSET_UL_CORNER必须比SPATIAL_SUBSET_LR_CORNER的纬度大、经度小，写反了交集为空，MRT 会输出一个 0 字节或者全填充的文件，而且退出码是 0，看着像成功。我的习惯是在 prm 文件里把角点写成注释放在旁边，改参数时肉眼对一遍。

5.5 现象：算出来的 NDVI 全是几百上千的整数；原因：没做定标；解决：读 scale_factor 并先掩膜后换算

MODIS 的 L3 产品为了省存储，把浮点物理量压成了Int16。NDVI 存的是真实值 × 10000，地表温度存的是开尔文 × 50（MOD11A2的scale_factor是 0.02），反射率存的是真实值 × 10000。直接把数组拿去算或者画图，得到的是一张数值范围在 -3000 到 10000 之间的图，看着“有数据”，其实完全不对。

正确顺序是三步：先按valid_range和_FillValue做掩膜，再乘scale_factor加add_offset，最后才做指数计算或统计。MOD13Q1的 NDVI 有效范围是 -2000 到 10000，超出这个范围的都是填充值或者云污染像元。还有一点，MOD11A2的LST_Day_1km在海洋和常年冰雪区是填充值 0，如果直接参与均值统计会把结果拉低，务必先掩掉。

6. Albers 参数与批量流水线：让每一次重投影都对得上

投影参数这件事，写一次和写一百次是两种心态。单次处理时随便

本文还有配套的精品资源，点击获取