做遥感的人应该都有同感：MODIS数据的下载和后处理，是你迟早要补的一门课。中分辨率成像光谱仪（MODIS）随Terra卫星在1999年发射、Aqua卫星在2002年发射，到现在已经积累了超过二十年的连续对地观测记录，36个波段覆盖从250米到1公里多种空间分辨率，产品线覆盖植被指数、地表温度、气溶胶、火烧迹地、积雪、初级生产力等，而且全部免费开放。也正因为免费加长时序，几乎所有的区域生态、农业、气候研究课题都会用到它。但免费并不等于省心，恰恰相反，数据下载和后处理才是真正的大头。

我见过不少同学卡在同一个地方：从平台搜到了数据，下载回来一堆HDF文件，然后对着命令行发愣——文件打不开、投影对不上、数值看着像噪声、拼接完有缝、遮云又不会处理。这篇文章就是把“MODIS数据下载及后处理”这条链路完整走一遍，从产品选型、批量下载、格式转换、重投影、拼接裁剪，到NDVI计算和质量控制，把每一步的原理和实操都讲清楚。适合刚开始接触MODIS的研究生、做生态和农业应用的项目人员，也包括想把时序分析流程跑通但一直被数据预处理卡住的开发者。

1. 内容整体设计与思路拆解

1.1 先想清楚你需要的到底是哪级产品

MODIS数据分为多个级别，最常接触的是L2级和L3级产品。L2级是经过辐射定标和大气校正后的逐轨数据，比如MOD021KM就是L2级的1公里辐射亮度数据；L3级则是把多轨数据按固定网格重新投影、合成后的产品，比如MOD13Q1、MOD11A2，这类产品空间位置整齐、时间尺度明确，绝大多数应用场景直接用L3级就够了。

你不需要自己去处理L2级的几何校正、大气校正那一整套流程，除非你要做非常高精度的定量反演。绝大多数情况下，选择L3级产品是正确的起点。以植被研究为例，直接下载MOD13Q1，拿NDVI/EVI做时序，比从MOD09GA逐日地表反射率自己合成16天NDVI要省事得多，而且官方产品在云检测、合成角度筛选上已经做了很多工作。

选产品之前，先回答三个问题：关注的地表参数是什么？研究区范围有多大？时间尺度是日、旬、月还是年？这三个答案直接决定你下哪些产品、多少景数据、用什么方式处理。

1.2 下载与后处理的完整链路

整个流程可以用一条链路说清楚：数据检索 → 批量下载 → 格式转换 → 重投影 → 拼接 → 裁剪 → 指数计算 → 质量过滤 → 时序分析。

每一步之间是有依赖关系的，顺序乱了容易白干活。比如先重投影再去拼接，和先拼接再重投影，结果基本一样，但计算效率差别很大——如果研究区跨多个瓦片，建议先拼接再裁剪，因为重投影本身有计算量，处理单张瓦片再拼，反而容易出现边缘变形。

还有一个经常被忽略的环节：数据质量标记（QA）处理。MODIS产品几乎都附带质量波段，里面用位掩码记录云、雪、传感器视角、气溶胶浓度等信息。很多人把NDVI算出来就直接用，忽略了QA，结果时序曲线上出现大量异常波动，其实是云污染没滤掉。所以在设计后处理流程时，一定要把QA波段处理放进主流程，而不是当作可选项。

2. 核心细节解析与实操要点

2.1 常用产品体系快览

把最常用的MODIS产品整理成一张表，后面选型时直接对照：

注意产品编号前缀，MOD表示Terra卫星数据，MYD表示Aqua卫星数据，MCD表示双星联合产品。比如MOD13Q1来自Terra，MYD13Q1来自Aqua。如果做时间序列分析，可以组合两星数据缩短重访周期，但必须注意两星过境时间不同，NDVI反射特性会有差异，合并前要先做一致性检查。

2.2 分辨率、尺度因子和填充值怎么读

MODIS产品里存的通常不是直观的物理量，而是放大后的整型数值。以MOD13Q1为例，NDVI波段存储的是有符号16位整数，范围通常是-2000到10000，真实NDVI值要乘以0.0001的尺度因子（scale factor）。填充值（fill value）是-3000，代表该像元无有效数据。类似的，MOD11A2的地表温度存储值乘以0.02才是开尔文温度，不处理尺度因子直接把数值拿去用，结果会完全不对。

我习惯在拿到产品的第一件事，就用软件查看产品文档里的Scale Factor、Add Offset、Valid Range、Fill Value四项参数。每款产品略有差异，即使同一个产品不同版本（Collection 5与Collection 6）也可能不同。别凭经验记死参数，以产品附带元数据为准。

2.3 正弦曲线投影为什么非处理不可

MODIS的L3级产品采用的是一种全球正弦曲线投影（Sinusoidal投影），并把全球划分成36列18行的瓦片网格。文件命名里的h25v05，h对应水平方向列号，v对应垂直方向行号。这种投影的优势是全球面积变形相对均衡，但问题在于它不适合直接拿来和当地的地理坐标、矢量边界叠加，也不能直接对接大多数分析软件中常用的经纬度坐标。

所以后处理里几乎必做的一步就是把正弦投影转换为WGS84经纬度坐标，或者转换为研究区所在的高斯、UTM等投影。如果这一步不做，后续你想要裁剪到行政区边界，就会发现MODIS影像和矢量数据怎么都对不上，差着一整套坐标系的偏移。

2.4 QA波段为什么要专门处理

MODIS的每个科学数据集（SDS，Science Dataset）旁边，通常伴生一个质量评价数据集。MOD13Q1里有“250m 16 days pixel reliability”，用简单的数值分级表示像元质量：0表示理想数据、1表示边缘数据、2表示雪覆盖、3表示云覆盖。还有一种更复杂的“VI Quality”是16位位掩码，每一位都有具体含义。

实际使用时，对于大多数研究来说，pixel reliability已经够用。把等级大于等于2的像元标记为无效，NDVI时序曲线会干净很多。很多刚上手的人不懂这一点，觉得官方产品已经做过大气校正，就应该直接可用，但实际上云污染、传感器异常、冰雪覆盖仍然会混在数据里。

3. 实操过程与核心环节实现

3.1 官方数据平台的检索与下载

获取MODIS数据的官方渠道，主要是NASA的数据平台Earthdata，以及相对更直接的LAADS DAAC。第一步是在平台注册账号，这一步需要邮箱验证，部分情况下还要求设置双重验证。账号用于后续所有下载请求，不注册的话批量下载基本走不通。

以LAADS平台为例，选择Product（比如MOD13Q1），版本建议直接选最新版本（目前常见的是Collection 6.1），然后设定时间范围、空间范围或瓦片编号。平台支持按瓦片选择，比如我需要中国中部某区域，可以先确定覆盖该区域的瓦片编号，再按日期挑选数据。

下载方式有两种：直接在网页勾选文件下载，适合数据量小、偶尔下载的情况；另一种是用命令行批量下载，适合长时间序列、多瓦片的情况。网页方式就不多说了，重点说命令行。

Earthdata下载建议使用wget或curl配合.netrc文件实现自动化认证。先把账号信息写入用户目录下的.netrc文件，注意权限要设为600：

echo "machine urs.earthdata.nasa.gov login 你的用户名 password 你的密码" > ~/.netrc chmod 600 ~/.netrc

然后在产品对应的目录列表页找到URL。MODIS官网的产品目录通常按产品名、日期分目录组织，比如MOD13Q1的下载根路径类似于https://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.061/，下面是2018.01.01这样的日期目录。可以用wget递归抓取：

wget -e robots=off -m -np -nH --cut-dirs=3 \ "https://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.061/2018.01.01/" \ -P ./MOD13Q1_2018001

参数解释：-m开启镜像模式，-np不爬取上级目录，-nH不创建主机名目录，--cut-dirs=3截断前三级路径只保留文件和后续目录，-P指定保存目录。实测下来，用这个方式下载几十上百景数据很稳定，断线了重跑一次会跳过已下载文件，相当于有天然的断点续传。

提示：批量下载一定不要直接用浏览器一个个点。一是慢，二是容易漏。用命令行脚本，配合日期列表循环，可以一次把一整年的数据拉完。

3.2 AppEEARS和GEE：两种更省事的替代思路

如果不想下载整景影像，只想提取某个站点或小区域的时间序列曲线，有两个更省事的工具。

AppEEARS是NASA提供的在线提取工具，可以上传点坐标或矢量边界，选定产品和时间范围后，平台直接返回每个点/区域的像元值时间序列，以CSV或NetCDF格式输出。优点是不用处理投影、裁剪这些事，官方直接把数据提取好了。

Google Earth Engine（GEE）则是更灵活的云端处理平台。MODIS产品在GEE里以ImageCollection形式提供，比如MODIS/061/MOD13Q1。在GEE里筛选日期和区域，然后做reduce、chart，最后导出一个GeoTIFF或CSV。

var ndviCol = ee.ImageCollection("MODIS/061/MOD13Q1") .filterDate("2018-01-01", "2021-12-31") .filterBounds(roi) .map(function(img) { return img.select("NDVI").multiply(0.0001) .clip(roi).copyProperties(img, ["system:time_start"]); });

这两种替代思路适合不同需求：需要栅格原图做拼接裁剪研究的，老老实实下载；只需要点/区时间序列的，AppEEARS够用；如果需要按自己的算法批量处理，GEE更合适。三者不是互斥关系，很多项目里先下载部分数据做验证，再用GEE做全尺度分析，是很常见的组合。

3.3 HDF格式与子数据集定位

MODIS数据下载回来是HDF格式。别被这个格式吓到，HDF本质上就是一个文件里的多层数据结构。用GDAL打开时，你会看到这个文件里包含多个子数据集（Subdataset），每个子数据集对应一个科学数据集，比如NDVI、EVI、pixel_reliability等。单独打开文件并不会直接给你一个完整的栅格，你需要先选择要哪个子数据集。

用gdalinfo查看文件结构：

gdalinfo MOD13Q1.A2018001.h25v05.006.2018034103136.hdf

输出里会列出Subdatasets，形如：

SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD13Q1.A2018001.h25v05.006.2018034103136.hdf":MODIS_Grid_16DAY_250m_VI:"250m 16 days NDVI"

要读取某个子数据集，必须在文件名里带上这个完整的子数据集路径，直接用GDAL打开整个HDF文件通常会报错或者只能看到元数据。

注意：不同版本的MODIS产品，Grid名称和SDS名称可能有细微差别。不要凭网上老教程硬套，先gdalinfo看一眼实际名称，再写后续命令，这个习惯能避免大量低级错误。

3.4 用GDAL做格式转换和重投影

后处理最核心的两步是格式转换和重投影。用GDAL命令行最快。先把NDVI子数据集单独转成GeoTIFF：

gdal_translate -of GTiff \ "HDF4_EOS:EOS_GRID:\"MOD13Q1.A2018001.h25v05.006.2018034103136.hdf\":MODIS_Grid_16DAY_250m_VI:\"250m 16 days NDVI\"" \ NDVI_raw.tif

然后重投影到WGS84经纬度：

gdalwarp -t_srs EPSG:4326 -dstnodata -3000 NDVI_raw.tif NDVI_wgs84.tif

-dstnodata -3000是把填充值重写为-3000，这样后续计算时能统一识别无效像元。如果要输出UTM投影，把EPSG:4326换成对应的UTM带号代码就行。

如果研究区跨多个瓦片，可以先用gdalbuildvrt构建虚拟栅格，再做一次gdalwarp，这样只需要一次重投影，而不是每个瓦片单独转完再拼：

gdalbuildvrt NDVI.vrt NDVI_wgs84_*.tif gdalwarp -t_srs EPSG:4326 -cutline study_area.shp -crop_to_cutline NDVI.vrt NDVI_study.tif

-cutline配合-crop_to_cutline直接完成裁剪。就这么一条命令，把拼接和裁剪合并做了，效率很高。

3.5 NDVI计算和尺度因子还原

MODIS官方NDVI产品还需要乘以尺度因子才是真实NDVI。这一步我建议用Python的rasterio或GDAL的Python绑定来做，顺便把QA掩膜一起处理。

import numpy as np from osgeo import gdal ds_ndvi = gdal.Open("NDVI_wgs84.tif") ndvi_raw = ds_ndvi.ReadAsArray().astype(np.float32) # 还原真实NDVI ndvi = ndvi_raw * 0.0001 # 无效值处理 ndvi[ndvi_raw == -3000] = np.nan

如果你下载的是MOD09GA地表反射率产品，想自己算NDVI，那公式就是(NIR - Red) / (NIR + Red)，注意先各自乘以0.0001尺度因子，再相减相加，避免整数溢出。

3.6 QA波段的掩膜提取

以MOD13Q1的pixel reliability为例，用GDAL读取QA子数据集，然后生成一个掩膜：

ds_qa = gdal.Open("pixel_reliability_wgs84.tif") qa = ds_qa.ReadAsArray() # 生成质量掩膜：0和1保留，2（雪）、3（云）剔除 valid_mask = np.where((qa == 0) | (qa == 1), 1, 0) # 应用掩膜 ndvi_clean = np.where(valid_mask == 1, ndvi, np.nan)

这一步做完，NDVI时序曲线通常就会平滑很多。官方产品的应用标准里一般也会注明推荐使用哪些QA标志，建议在做任何分析之前先做这一步。

3.7 时间序列文件组织

做多年时序分析时，文件组织直接影响效率。我常用的目录结构是：

MOD13Q1/ 2018/ 001_2018-01-01_16day_NDVI.tif 017_2018-01-17_16day_NDVI.tif ... 2019/ ...

每个文件名里带上原始日期，后续写Python脚本处理时会非常方便，直接用glob排序即可。文件名里缺少日期信息会导致后期处理时序时还要去翻元数据，完全是自己给自己挖坑。

4. 常见问题与排查技巧实录

4.1 下载中断、认证失败的排查

下载是最先遇到的麻烦。常见问题包括：wget提示401认证失败、下载到一半连接被断开、下载的HTML文件伪装成HDF文件等。

认证失败通常是因为.netrc里的主机名和实际下载主机不一致。Earthdata有多个跳转域名，建议在.netrc里同时配置urs.earthdata.nasa.gov和实际下载域名，比如e4ftl01.cr.usgs.gov。

下载到一半断线，wget的-c参数可以断点续传，但注意和-m一起使用时行为略有不同，我更习惯用循环脚本逐个文件下载，配合-c和超时重试参数：

while read url; do wget -c -t 5 --timeout=60 "$url" done < download_list.txt

还有一种坑：wget抓下来的文件名是.hdf，但文件大小只有几KB，用file命令检查发现是HTML。这是平台返回了错误页面，通常是因为下载请求未携带有效认证，或者访问频率过高被限流。解决方法是放慢下载速度，增加请求间等待，并确认认证信息无误。

4.2 打开HDF文件提示无法识别

GDAL打开MODIS HDF报错，最常见原因是GDAL版本过老，不支持HDF4 EOS格式，或者编译时没有包含HDF4驱动。可以先执行gdalinfo --formats | grep HDF确认驱动是否存在。如果确实没有，建议更新到新版GDAL，或者通过Anaconda安装一个带全部驱动的版本，比如conda install gdal后默认就带HDF4驱动。

另外，很多人用ArcGIS直接打开MODIS HDF，只在目录里看到一个文件，拖进去后只有极少数子数据集能被识别。这不是数据坏了，而是GIS软件对HDF4 EOS容器的支持层次不同。用GDAL或MRT处理后再转成GeoTIFF，和GIS软件的兼容性最好。

4.3 重投影后影像范围偏移

重投影后影像偏离预期范围，常见原因是没有正确处理填充值。在新投影中，某些像元坐标无法被原数据覆盖，GDAL默认会填入0，有时这个0被误认为有效值参与计算，导致结果像元范围看起来“变了”。解决办法是在重投影时显式指定-dstnodata，确保新影像的无效区域和原数据填充值对应。

还有一种偏移是源数据本身的网格描述信息缺失，导致GDAL无法正确读取投影。这种情形较罕见，一般出现于非官方渠道下载的、处理过的数据。遇到这种情况，可以用MRT（MODIS Reprojection Tool）作为替代方案，它对MODIS的HDF元数据兼容性非常好，重投影结果稳定。

MRT的操作虽老但实用：通过图形界面打开HDF文件，选择需要的SDS，设置输出投影、分辨率、输出路径，然后生成批处理参数文件（.prm），后续用命令行批量处理：

modisTool -prm example.prm

MRT已经被官方停止维护，但至今仍能正常运行，作为GDAL之外的第二套后处理工具，在遇到批处理和重投影兼容性问题时可以互补使用。

4.4 数值范围异常

NDVI算出来大于1或小于-1，通常有三个原因：没乘尺度因子、没处理填充值、或者影像里含有云和雪。逐一排查顺序也是先乘尺度因子，再处理填充值和QA。MOD13Q1的NDVI理论上范围是-0.2到1.0左右，超过这个上界的应该怀疑是否含未滤除的噪声像元。

地表温度产品MOD11A2的常见错误则是把开尔文温度当摄氏度直接用。记住LST乘以0.02后得到的是开尔文值，减去273.15才是摄氏度。

4.5 拼接后出现接缝

跨瓦片拼接时，相邻瓦片间偶尔会出现明暗差异或错位。原因可能是两景数据采集时间不同，导致地表状态不一致，比如一景在云后、一景在云前。解决办法是检查输入数据的时间和QA，尽量使用同一天过境的数据；如果云造成差异，可以考虑用邻近期合成产品代替单日产品。

还有一种接缝并不明显，但在做植被指数后特别显眼：因为不同瓦片的太阳高度角、传感器视角不同，反射率本身就有差异。这类接缝可以通过后处理做匀色，但代价是改变原始物理量，定量研究中一般不建议。

4.6 内存不足和文件过大

MODIS数据单景不算大，但多年多瓦片拼接后，中间文件可能非常大。比如MOD09GA全中国一年每天500米反射率拼接，中间GeoTIFF动辄几个GB。排查思路是控制中间文件体积：用VRT虚拟栅格代替实际拼接GeoTIFF，只在最后分析时生成实际文件；用Float32而不是Float64存储中间结果；裁剪尽量提前，只保留研究区范围。

Python处理大文件时，建议用rasterio的窗口读取（window）分块处理，不要一次性ReadAsArray整个影像。当年我处理某区域的多年NDVI时，一次性读入一整个拼接文件直接把内存吃满，改成按块读取后速度反而更快了。

4.7 常见问题速查表

5. 自动化与批处理流程优化

5.1 一个相对完整的多景处理脚本

当手里的数据是几十景、上百景的时候，每次手动执行gdal_translate和gdalwarp效率太低。我习惯把整个流程写到Python脚本里，大致分为三部分。

第一步，遍历所有HDF文件，读取需要的SDS名称并转为GeoTIFF。第二步，统一重投影到目标坐标系并裁剪研究区。第三步，计算最终指数并叠加QA掩膜。

以MOD13Q1为例，简化版脚本结构如下：

import os import glob from osgeo import gdal hdf_list = sorted(glob.glob("MOD13Q1.A*.h25v05*.hdf")) for hdf_path in hdf_list: # 从文件名提取日期信息 basename = os.path.basename(hdf_path) doy = basename.split(".")[1][1:4] ndvi_sds = 'HDF4_EOS:EOS_GRID:"{}":MODIS_Grid_16DAY_250m_VI:"250m 16 days NDVI"'.format(hdf_path) qa_sds = 'HDF4_EOS:EOS_GRID:"{}":MODIS_Grid_16DAY_250m_VI:"250m 16 days pixel reliability"'.format(hdf_path) # 子数据集转GeoTIFF gdal.Translate("ndvi_{}.tif".format(doy), ndvi_sds) gdal.Translate("qa_{}.tif".format(doy), qa_sds) # 用gdal.Warp做重投影和裁剪 gdal.Warp("ndvi_wgs84_{}.tif".format(doy), "ndvi_{}.tif".format(doy), dstSRS="EPSG:4326", dstNodata=-3000, cutlineDSName="study_area.shp", cropToCutline=True) gdal.Warp("qa_wgs84_{}.tif".format(doy), "qa_{}.tif".format(doy), dstSRS="EPSG:4326", dstNodata=255, cutlineDSName="study_area.shp", cropToCutline=True)

这里的gdal.Translate对应命令行里的gdal_translate，gdal.Warp对应gdalwarp。脚本跑完后，每个日期都生成一对NDTI和QA tif文件。后续再做尺度还原和QA掩膜，输出一个整理好的时序影像立方体。

5.2 处理效率的取舍

处理MODIS数据时，效率问题往往不是单一环节造成的。文件数量、单个文件大小、处理步骤数量、中间文件是否落盘，每一项都会影响总耗时。实际项目里，我通常把中间GeoTIFF放在SSD上，避免机械硬盘写入成为瓶颈。

多景数据量大的时候，用concurrent.futures做简单的多进程并行也能显著提升效率。MODIS后处理每一步都是独立的，天然适合并行。但要小心同时处理几十个文件时磁盘IO吃满，反而拖慢整体速度。一般设置四到八个并行任务即可，不必贪多。

5.3 处理后数据的质量检查

任何自动化流程都建议在批处理跑完后做一次快速质量检查。我最常用的检查手段是把处理前后NDVI的直方图打印出来，看有没有异常的尖峰。比如大量像元集中在-0.2或1.0附近，大概率是无效值没处理好或者有异常像元没滤除。

另外，取一个研究区内的稳定点做时间序列曲线，是检验处理流程是否正确的有效方式。如果是草地或农田，NDVI时序应该有明显的季节波动；如果曲线是一条直线或毫无规律跳变，基本可以断定流程里某一步出了问题。

6. 后处理之后的进阶方向

数据下载和后处理本身不是终点。处理好的MODIS影像，通常会进入更深层的应用流程，比如物候提取、趋势分析、分类制图、模型模拟。这里简单说一下常见的衔接方向。

物候提取方面，处理好NDVI时序后，可以用分段拟合的方式提取生长季开始日期、结束日期和峰值时间。常用的方法包括动态阈值法、拐点法和Savitzky-Golay滤波。比如动态阈值法，给每个像元设定一个NDVI阈值比例，比如达到季节振幅的20%时视为生长季开始，这种方法在大量研究中被验证有效。

趋势分析方面，对多年时序做Mann-Kendall趋势检验和Theil-Sen斜率估计，是评估区域植被变化趋势的标准组合。研究区范围内的NDVI趋势图，可以直接服务于区域生态质量评价和植被恢复成效监测。

数据产品衔接方面，MODIS处理的成果也常作为输入送进其他模型。比如MOD15A2H的LAI产品直接是生态模型的输入，MOD17A2H的GPP产品用于碳循环估算。这些模型对数据的空间参考、时间连续性、质量标记都有要求，所以前期的下载和后处理质量直接决定模型输出可信度。

MODIS数据处理这套流程，我个人体会是，真正花时间的往往不是技术难点本身，而是对细节的控制。尺度因子、投影、QA这三个环节只要有一个疏忽，后续分析就会出现系统性偏差，而且这种偏差在统计结果里很难一眼看出来。建议每一个新项目开始，都用一幅已知区域、已知数值范围的影像做全流程测试，确认数值正确后再批量处理。我早期做某区域多年NDVI时，曾因为一个尺度因子疏忽，跑完一整批数据才发现所有NDVI都偏大，重新处理耽误了大半个月。把完整流程跑通一次并做好质量控制，比一开始急着处理所有数据更重要。另外，后处理脚本和中间文件也建议一并留存，MODIS产品版本会更新，半年后再想复现结论时，有这些记录能省很多事。