MODIS数据选型、下载、HDF4读取、MRT重投影与QA筛云

简介：面向遥感与地理信息领域的初学者及科研应用人员，这份文档围绕MODIS数据的背景、下载与处理展开系统梳理，可帮助缺乏卫星遥感经验的读者快速建立从数据获取到分析应用的完整认知。包内仅含1个doc文档，压缩包约503KB，篇幅紧凑却覆盖了ESE与EOS计划脉络、Terra与Aqua双星互补观测机制、36个波段的光谱范围与空间分辨率、2330公里视场及705公里太阳同步极轨等关键技术指标，并给出波段设计目标与典型应用领域的对照。文档同时说明NASA官网的产品与时间范围选取流程、HDF格式文件的读取方式，以及借助ENVI、ArcGIS或开源工具完成预处理、辐射与几何校正、时序合成和植被指数、地表温度提取的思路。已有304人学习，适合作为课程作业、课题入门或环境变化监测项目的前期参考资料。

1. MODIS数据：1到2天覆盖全球的36波段，到底适合谁

如果你要把一个区域几年的植被长势、地表温度或者云量变化跑一遍，又不想在数据采购上耗预算，MODIS数据大概率是最先该试的那一套。Terra 和 Aqua 两颗太阳同步极轨卫星互相补位，上午星和下午星各自过境，1 到 2 天能把全球表面扫一遍；星上搭载的中分辨率成像光谱仪有 36 个波段，光谱范围从 0.405 微米一直铺到 14.385 微米，空间分辨率分 250 米、500 米、1000 米三档，幅宽 2330 公里。它做不了米级精细制图，也替代不了高光谱专题产品，但做区域趋势、长时序对比、大范围快速摸底，这套数据的性价比很难绕开。下面按产品选型、下载、HDF4 读取、MRT 投影、排错、合成验证的顺序，把这份资料里能直接落地的部分拆开讲。

2. 产品分级与文件名：MOD、MYD、MCD三个前缀背后的选型账

2.1 0级到4级：先分清你要的是原始观测还是算好的产品

MODIS 标准数据产品按处理深度分成 0 到 4 级，这个分级不是学术分类，它直接决定你下载的文件里装的是什么、后面要不要自己写算法。0 级是地面站直接收到的原始数据，一般用户碰不到；1A 级做了数据重建和时间配准，仍然是没有辐射校正的原始计数值；1B 级在 1A 基础上做了辐射校正和定标，输出的是辐亮度，MOD02 就是这一级。真正开始“能用”的是 2 级以后：2 级产品和 1B 具有相同空间分辨率、覆盖同一地理区域，比如 MOD09 表面反射率、MOD11 地表温度；3 级是把多天数据放进统一时空栅格做合成，比如 MOD13 的 16 天 NDVI；4 级则是跑完模型之后的结果，比如 MOD17 的初级生产力。

选型的判断依据很简单：如果你要自己做指数计算、自己控制云掩膜，就从 2 级甚至 1B 起步；如果你只想要一条能直接分析的时间序列，3 级产品省掉的预处理时间远超它损失的那点灵活性。我一般会先确认研究问题的时间尺度——日尺度看 2 级，旬月尺度看 3 级，年度碳收支看 4 级，别一上来就下 1B，几十景原始辐亮度堆在硬盘里，光定标和几何定位就能拖掉一周。

2.2 44种标准产品怎么倒着选

标准产品按内容分成大气、陆地、海洋三大类，合计 44 种。新手最容易犯的错是照着产品编号从 MOD01 往下翻，翻到哪算哪。更靠谱的做法是从应用倒推：先写清楚你要的输出是什么，再找对应的短名。下面这张表是日常用得最多的十几个，覆盖了土地覆被、植被、温度、火情、云和气溶胶几个方向。

表里没有列的，去产品手册里按“短名 + 版本号”检索，不要凭记忆猜。版本号很关键，同一产品不同 Collection 的处理基线不同，定标系数和 QA 位定义都可能变，混用会得到一条斜率对不上的时间序列。

2.3 文件名拆成七段，每一段都影响后面脚本

MODIS 的文件名不是随手拼的，拿MOD01.A2002032.1015.003.2002034184356.hdf举例，七段各有含义，写批处理脚本时经常要按段提取。

最容易翻车的是 c 段和 d 段。儒略日不是月日，032是 2 月 1 日不是 3 月 2 日；1015是 UTC 时间，换算到北京时间要加 8 小时，所以一景标着A2020032.1015的数据，实际对应北京时间 2 月 1 日 18 点 15 分。做日尺度统计时如果直接拿文件名里的日期当本地日期，跨年那几天必错。

2.4 正弦投影与h17v04：瓦片编号不是随便起的

大多数标准产品用正弦投影（Sinusoidal），赤道处按 10 度乘 10 度切网格，行号从左上角 (0,0) 起算，到右下角 (35,17)，全球共 36 列 18 行。文件里的h17v04就是水平第 17 列、垂直第 4 行，一景大致覆盖中纬度一个省级尺度的范围。瓦片编号的价值在于：同一h17v04的多年数据可以直接按像元对齐做时间序列，不需要先做空间重采样；跨瓦片研究才需要拼接，而拼接一定伴随投影和像元网格的变化。

这里有个常见误解：以为同一瓦片号在所有产品里范围完全一致。250 米、500 米、1 公里产品的瓦片覆盖范围是嵌套的，但像元数不同，MOD09GQ（250 米）和 MOD09GA（500 米/1 公里）的阵列尺寸差一倍，混着读之前先看清楚SDS的维度。投影参数、瓦片号、分辨率这三个信息凑齐，才能判断两景数据能不能直接逐像元运算。

3. 下载：从经纬度检索到wget批量拉取的四个步骤

3.1 检索条件写成“经纬度框+UTC时间段”

数据检索页通常提供按经纬度框选、按瓦片号选、按日期范围选三种方式。区域研究用经纬度框最直接，输入左上角和右下角坐标，再填时间范围。时间一定按 UTC 填：要取北京时间 1 月 16 日全天的数据，UTC 范围是 1 月 15 日 16:00 到 1 月 16 日 16:00，直接填 1 月 16 日 00:00 到 24:00 会漏掉前 8 小时。检索结果里会列出每景的文件名、大小、平台、版本，先别急着全选，用文件名里的 d 段确认过境时间是不是你要的上午星或下午星。

检索结果一般支持导出成文本清单，这个清单就是后面批量下载的输入。清单里通常带完整下载链接，但链接可能带会话参数，过几小时就失效，所以脚本里更稳的做法是只保留文件名，把数据目录作为变量拼前缀。L1B 数据记得把 MOD03 一起勾上，它是 1 公里、500 米、250 米产品做几何纠正的定位文件，漏下这一景，后面打开 MOD02 会发现没有经纬度，返工重下最费时间。

3.2 账号与token：把登录这一步提前做掉

数据平台现在基本都走统一登录，匿名 FTP 的老路子已经不好用。注册账号后在个人中心生成一个下载 token，之后所有请求带Authorization: Bearer头即可，不用每次输密码，也不会因为会话过期在半夜断掉。token 等同于密码，别写进公开仓库，脚本里用环境变量传入。下面这段是批量下载的骨架，把清单文件、数据目录、token 三个变量换掉就能跑。

#!/usr/bin/env bash # 批量下载 MODIS 数据：断点续传、失败重试、限速 set -u export ED_TOKEN="你的下载token" BASE_URL="https://<检索页复制下来的数据目录>" OUT_DIR="./modis_raw" mkdir -p "$OUT_DIR" while read -r name; do [ -z "$name" ] && continue # 跳过空行 [ -f "$OUT_DIR/$name" ] && continue # 已存在则跳过 wget -c -nc --tries=5 --timeout=60 --wait=2 \ --header="Authorization: Bearer ${ED_TOKEN}" \ -P "$OUT_DIR" "${BASE_URL}/${name}" || echo "FAILED $name" >> failed.txt done < urls.txt

-c是断点续传，网络中断后重跑不会从头拉；-nc防止覆盖已完成的文件；--tries=5和--timeout=60对付偶发的连接超时；--wait=2每次请求间隔两秒，批量几百景时能明显降低被限流的概率。失败的文件名单独记到failed.txt，第二天换时段重跑，比盯着终端一屏屏翻要省心。数据目录不要手拼，直接从检索结果里复制一层，不同产品、不同版本的目录结构不完全一样。

3.3 下载完先校验，别等读的时候才发现是HTML

下载中断最隐蔽的结果不是文件缺失，而是文件存在但只有几 MB，打开一看是登录失效后返回的 HTML 错误页。HDF4 文件有固定的文件头，正常文件开头是 HDF 标识，用几行 Python 就能批量筛一遍，顺便统计每个文件的 SDS 数量，数量为零的直接标红。

from pathlib import Path from pyhdf.SD import SD, SDC bad = [] for p in Path("modis_raw").glob("*.hdf"): try: h = SD(str(p), SDC.READ) n = len(h.datasets()) # SDS 数量 h.end() if n == 0: bad.append((p.name, "no SDS")) except Exception as e: bad.append((p.name, repr(e))) for name, reason in bad: print(name, reason) print("待重下:", len(bad))

SD打开失败基本就是文件截断或根本不是 HDF；能打开但datasets()为空，多半是下到了错误页。这个检查跑一遍几十秒，比后面在重投影环节报一个看不懂的错要划算。顺手把文件按“产品短名/年/儒略日”建目录，后面 MRT 批处理和质量筛选都按这个层级遍历，省得再写一套路径解析。

3.4 清单、日志、版本三样东西留好

批量下载最怕两周后回头找“当时到底下了哪几景”。我的习惯是每批数据在目录下留三个文件：原始清单urls.txt、下载日志failed.txt、以及一个manifest.csv，记录文件名、大小、下载时间、Collection 版本。看着繁琐，等审稿人或合作者问“你这批数据版本是不是统一”的时候，这份清单就是后悔药。尤其是 MOD 和 MYD 混用、或者跨 Collection 拼接时，版本不一致带来的反射率偏差不会报错，只会让趋势线的斜率悄悄偏掉。

4. HDF4结构与MRT重投影：SDS读取到Albers输出的完整链路

4.1 HDF4的层级：SDS、Vdata、Vgroup与属性

MODIS 的.hdf是 HDF4 格式，不是 HDF5，两者不兼容。HDF4 内部是层次结构：最常用的是 SDS（Scientific Dataset），也就是一维到多维的数组，反射率、温度、QA 都是 SDS；Vdata 存表格型数据；Vgroup 负责把相关数据集组织在一起；每个 SDS 自带属性，定标系数、填充值、有效范围都藏在属性里。用pyhdf.SD打开文件后先列一遍数据集名和维度，能避免“凭记忆 select 一个不存在的名字”这种低级错误。

from pyhdf.SD import SD, SDC import numpy as np f = "MOD09GA.A2020032.h17v04.006.2020040034512.hdf" hdf = SD(f, SDC.READ) for name, info in hdf.datasets().items(): print(name, info[0], info[1:]) # 名称、类型、维度

输出里会看到sur_refl_b01到sur_refl_b07、state_1km、solar_zenith、satellite_zenith等。info[0]是数据类型，info[1:]是各维长度。MOD09GA 的 500 米反射率维度大约是 1200×1200，1 公里 QA 是 600×600，两者不能直接逐像元运算，要么先重采样，要么用 500 米波段自己的 QA。读之前先看清楚维度，比读完之后发现形状对不上再回头查要快得多。

4.2 读反射率：scale_factor和_FillValue必须一起处理

MODIS 的反射率 SDS 存的是整型，不是物理值，直接读出来是几千上万的整数。每个 SDS 的属性里有scale_factor、add_offset、_FillValue、valid_range，正确顺序是先按填充值掩膜，再乘比例加偏移。顺序反了，填充值也会被定标成一个看似合理的负值，后面统计最小值时就会莫名其妙出现 -0.1 这种数。

sds = hdf.select("sur_refl_b01") raw = sds.get().astype("float32") attr = sds.attributes() scale = attr.get("scale_factor", 1.0) offset = attr.get("add_offset", 0.0) fill = attr.get("_FillValue", -28672) refl = raw * scale + offset # 定标 refl[raw == fill] = np.nan # 先掩膜再统计 print(refl.shape, np.nanmin(refl), np.nanmax(refl)) sds.end()

MOD09 的scale_factor通常是 0.0001，_FillValue是 -28672，valid_range在 -100 到 16000 之间。算完的反射率理论值域在 0 到 1 附近，出现大于 1 的值先别慌，看看是不是云边缘或观测几何极端的情况，再检查有没有把valid_range之外的像元算进来。对于 MOD11 的地表温度，scale_factor是 0.02，填充值是 0，单位是开尔文，要转摄氏度还得减 273.15，这些都在产品手册的属性表里，别靠猜。

4.3 geo_field与经纬度查找表：几何定位的两种来源

1B 数据本身不带完整的逐像元经纬度，MOD03 才是定位文件，里面有Latitude、Longitude、SolarZenith、SensorZenith等数据集，维度与 1 公里波段一致。做几何校正时可以把 MOD03 的经纬度读出来，按行列号对应到 MOD02 的像元上；也可以在 MRT 重投影时直接把 MOD03 作为定位输入。两种方式的区别在于：前者适合自己写校正算法、做精确的像元级配准；后者适合批量出图，MRT 内部会处理定位，速度快但可控参数少。

读 MOD03 的经纬度时注意维度顺序，HDF4 的数组是行优先，Python 读出来是[row, col]，和影像的上下左右方向一致；但有些产品会做翻转或偏移，拿一景有明显地物特征的数据（比如海岸线）叠一下就能确认。经纬度是浮点数组，文件体积不小，批量处理时不要一次性把几十景的经纬度都读进内存，按瓦片循环、读完即释放。

4.4 MRT：先拼接再投影，Albers参数一次配好

MRT（MODIS Reprojection Tool）负责重投影、拼接、重采样、裁剪和波段子集。它的处理顺序是先把同一投影下的多幅影像拼接，再做投影转换，这个顺序不能反：如果先投影再拼接，接缝处的像元网格对不齐，会出现一条明显的错位带。参数文件（.prm）里最关键的几项是输出投影类型、标准纬线、中央经线、基准面和重采样方式。国内区域制图常用 Albers 等积圆锥投影，双标准纬线 25 度和 47 度、中央经线 105 度是一套通用参数，省级范围可以按轮廓再调，比如东西跨度大的区域把中央经线往中间挪。

重采样方式要按数据类型选：反射率、温度这类连续量用双线性，土地覆被、QA 这类分类或位标志必须用最近邻，否则一个 3 类像元被插成 2.5 类，后面统计类别面积全是错的。输出像元大小不要小于原始分辨率，500 米产品输出 250 米不会增加信息，只会让文件大四倍。参数文件写好后，用-R传入文件清单、-O指定输出目录批量跑。

MRT_HOME=/opt/mrt/bin $MRT_HOME/resample -P albers.prm -R filelist.txt -O ./out_albers -V

-V打开详细日志，每景处理成功还是失败都会打印；批量几百景时把日志重定向到文件，失败景单独重跑。MRT 对文件名有要求，输入清单里写绝对路径更稳，输出文件名会自动加波段后缀，后期按规则批量改名即可。跑完先抽一景和原始正弦投影数据叠一下，看边界是不是吻合，再接下去做时间序列。

5. 避坑排查：下载、读取、投影五个翻车现场

5.1 日期差一天：UTC和儒略日叠在一起

现象：按“1 月 16 日”检索，下回来的数据文件名是A2003015，做日统计时这一天的值对不上。

原因：文件名里的日期是数据采集的 UTC 儒略日，检索页填的时间范围也是 UTC。北京时间 1 月 16 日 00:00 到 08:00 属于 UTC 1 月 15 日，这 8 小时的数据会落在A2003015里。如果只按文件名日期分组，就会把本地日期的一整天拆成两半。

解决：统一在脚本里做时区转换，把文件名解析成 UTC 时间后再加 8 小时得到本地日期，按本地日期分组；或者干脆全程用 UTC 做时间轴，出图时再标注本地时间。两种都行，关键是别混。

5.2 只下MOD02，打开没有地理定位

现象：MOD02 能打开，但找不到经纬度数据集，做裁剪和重投影时无从下手。

原因：MOD02 是 1B 辐亮度，几何定位信息在配套的 MOD03 文件里，两者是分开下载的。检索时只勾了 MOD02，MOD03 就漏了。

解决：检索清单里同时勾选同一天、同一轨的 MOD03，下载后按文件名里的时间戳配对。MRT 重投影时把 MOD03 作为定位输入，或者在 Python 里读 MOD03 的Latitude、Longitude做自定义校正。配对规则写进脚本，别靠人工比对。

5.3 用h5py读HDF4，报错还找不到原因

现象：h5py.File("MOD09GA...hdf")直接抛异常，提示不是 HDF5 文件。

原因：MODIS 标准产品是 HDF4，h5py只认 HDF5。两个格式同名不同代，API 完全不兼容。网上不少示例代码混着写，复制过来就报错。

解决：HDF4 用pyhdf，HDF5 用h5py，先看文件头再选库。GDAL 两类都能读，gdalinfo一跑就知道格式。批量处理前统一检查一遍文件头，避免一半文件用一种库、另一半用另一种库。

5.4 反射率值域不对：scale_factor没乘

现象：读出来的反射率在 0 到 10000 之间，做指数计算全是异常值。

原因：SDS 存的是整型，物理值要乘scale_factor再加add_offset。只读了get()没读属性，或者顺序反了先定标再掩膜，填充值被算成了有效值。

解决：固定写成“读原始数组 → 读属性 → 掩膜填充值 → 定标”四步，封装成函数复用。每次换产品先打印attributes()看一遍，MOD09、MOD11、MOD13 的 scale 和 fill 都不一样，不能一套代码硬套所有产品。

5.5 MRT输出接缝、错位与内存爆掉

现象：多景拼接后中间有一条错位带；或者跑到一半进程被系统杀掉。

原因：错位多半是先投影再拼接，或者输入清单里混了不同 Collection、不同分辨率的文件；内存爆掉通常是一次性把几十景全丢给 MRT，或者输出像元设得过细。

解决：确认 MRT 的处理顺序是先拼接后投影，输入清单按产品、版本、分辨率分组，不同组分开跑；批量时按瓦片或按年分批，每批 10 到 20 景，跑完检查接缝再继续。输出像元不小于原始分辨率，连续量用双线性、分类量用最近邻。跑之前先用两景做小样本测试，参数确认无误再放开批量。

6. 从单日瓦片到月度合成：QA筛云与结果验证的一个笨办法

单日反射率直接平均是最容易出问题的做法，云、云影、观测几何差异都会混进月度值里。我的习惯是先读state_1km做云筛，再按像元统计有效观测数，有效数低于阈值（比如月内少于 5 次）的像元直接置空，不参与后续分析。MOD09GA 的state_1km是位标志，低两位表示云状态，0 为晴空，1 为有云，2 为混合，3 为未定义，先掩膜再定标。

import numpy as np from pyhdf.SD import SD, SDC def read_clear_refl(path, band="sur_refl_b01"): hdf = SD(path, SDC.READ) raw = hdf.select(band).get().astype("float32") attr = hdf.select(band).attributes() qa = hdf.select("state_1km").get() hdf.end() refl = raw * attr.get("scale_factor", 1.0) + attr.get("add_offset", 0.0) refl[raw == attr.get("_FillValue", -28672)] = np.nan cloud = qa & 0b11 # 低两位：云状态 refl[(cloud != 0) | (qa == 1)] = np.nan # 非晴空和填充值都剔除 return refl

这段代码只做了一件事：把非晴空像元变成NaN，后续用np.nanmean做月度合成时自然跳过。跑完一个月后，别急着出图，先看两个数：有效像元比例和月度均值的直方图。有效比例低于 30% 的区域，说明这个月云太多，合成结果的可信度要打折扣；直方图出现双峰或者异常拖尾，多半是 QA 位定义用错了版本，或者把云影漏掉了。拿一景已知水体或城市核心区的像元值做对照，水体反射率在可见光波段应该低而平，城市建成区在短波红外明显偏高，对不上就回去查定标和掩膜。

从那以后我每次下完一批数据，都强制先跑一遍“文件头检查 + SDS 列表 + QA 有效比例”这三步，再谈分析和出图。这套流程不高级，但能把大部分返工挡在前期。文中这些参数、脚本和产品对照表，我整理成了同一份笔记，和那份 MODIS 资料放在一起，需要的按文末方式自取。希望帮到你。

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