简介这是一份面向遥感初学者与地球科学方向研究者的MODIS数据学习文档内容覆盖数据背景、技术指标、下载方法和处理流程。文档首先介绍搭载于Terra与Aqua卫星上的中分辨率成像光谱仪说明其视场宽度两千三百三十公里、三十六个波段覆盖零点四零五至十四点三八五微米以及二百五十米、五百米、一千米等空间分辨率并结合波段设置解释其在土地覆被、自然灾害、气候变化等领域的应用价值。随后文档梳理NASA官网的数据获取步骤、HDF格式特点并给出使用ENVI、ArcGIS及开源软件开展预处理、校正、合成与分析的基本思路。资源包共包含一个doc文档大小五百零三KB内容紧凑便于系统阅读目前已有三百零四人学习下载。通过阅读读者能够快速建立对MODIS数据体系的整体认知掌握从原始数据到专题产品处理的常用流程适合作为课程作业、科研入门或项目预研的参考资料。1. 拿到 MODIS 数据的第一个小时先搞懂它是什么才敢往下做朋友接了个任务做某个流域近十年的植被指数变化分析他打开软件随手搜“MODIS数据下载”结果发现第一关不是下载而是产品多到眼花缭乱。MOD09GA、MOD13Q1、MOD11A1、MYD09……后缀不同、分辨率不同、时间频率不同连投影坐标都不是常见的经纬网格。MODIS 是当前地表遥感里覆盖面最广的中分辨率数据源两颗卫星交替过境、一天最多四次观测、产品按轨道和格网组织覆盖全球陆地。这篇笔记按“认识产品 → 拿到数据 → 处理成可用栅格”的顺序讲清楚适合刚接触遥感影像批量处理、准备把 MODIS 数据引入到自己统计或制图流程里的从业者。2. MODIS 数据产品怎么选从轨道到格网的三个层次2.1 两颗卫星的观测节奏Terra 与 Aqua 决定你拿到什么MODIS 同时装在两颗极轨卫星上这一点直接决定你下载时看到的一堆产品前缀。一颗在上午 10:30 左右由北向南穿过赤道另一颗在下午 13:30 左右穿过两颗都维持着稳定的对地观测轨道。上午那颗对应“MOD”前缀下午那颗对应“MYD”前缀当你看到“MCD”前缀时则表示该产品把双星数据做了合并处理。为什么要关心过境时间因为中分辨率光学影像最怕云。上午被云盖住的区域下午的轨道很可能还能补上一条清晰观测如果只认准单颗星的逐日产品时间序列上会出现大量空槽。像地表温度、火点这类需要日更新的应用我会优先选择双星产品而 NDVI 这类植被指数官方本身就有 16 天合成产品双星与否反而不用纠结。另一点容易被忽略的是扫描角。极轨卫星影像边缘的扫描角很大几何畸变和大气路径都变长同样是“逐日产品”轨道边缘的像元质量远不如星下点附近。所以下载时别只看日期也要关注目标区域在影像里的位置。2.2 产品分级的三级跳L1、L2、L3 分别适合什么MODIS 的 36 个波段覆盖 0.4 到 14.4 微米可见光、近红外、热红外一应俱全。原始观测经过辐射定标变成辐亮度数值对应 L1 级产品再经过大气校正、云掩膜等处理得到地表反射率、地表温度等 L2 级产品继续投影到固定全球格网并做多日合成就是 L3 级。绝大多数做工程分析的从业者不需要从原始辐亮度开始反演直接选用 L2、L3 能省掉一整套大气校正环节。我自己的选型习惯很直接做地表温度就取 L2 的 MOD11A1 逐日产品做植被长势就取 L3 的 MOD13Q1 或 MOD13A1 合成产品需要自己控制时间和空间范围、且对日期精度有要求时再退到 L2 的 MOD09GA 逐日反射率。这里有一条硬规则同一个分析里不同时期的文件产品集合版本必须一致比如都选版本 006 或版本 061混用会导致数值口径不一致后续任何变化检测都不可信。下表是常见产品的空间分辨率与适用场景方便你在下载前做第一轮筛滤空间分辨率代表波段常见产品适合场景250 米1、2 波段红光、近红外MOD09GQ、MOD13Q1高分辨率 NDVI、植被长势500 米1-7 波段MOD09GA、MOD14反射率订正、火点监测1000 米全部 36 波段MOD11A1、MOD21地表温度、云检测2.3 分幅与投影方式SIN 投影和 h/v 编码L3 级别的 MODIS 标准产品大多使用正弦投影Sinusoidal全球按 10 度 × 10 度划分成图幅每个图幅带一个 h水平和 v垂直编号。比如中国东部大部分区域落在 h25 到 h28、v04 到 v06 这个范围里下载时必须选对图幅编号否则拿到的是相邻区域的数据。图幅编号可以从数据分发平台的空间检索页面上直接读出也可以在本地用待分析区域的经纬度换算出来。为什么要用等面积性质的 SIN 投影因为 MODIS 的核心任务是地表参数反演NDVI、叶面积指数这些量都与面积相关。等面积投影保证了图幅内每个像元代表的真实地面面积大体一致区域统计才可比。但它的代价是经纬网和边界会“变形”所以处理时几乎总要再投影到目标坐标系比如 EPSG:4326 经纬度或 UTM 分带。对不熟悉投影概念的新手来说这一步经常成为第一个翻车点直接用原始 SIN 投影去做矢量裁剪结果边界错位、面积统计偏大。3. 把 MODIS 数据拿全注册、检索与批量下载脚本3.1 数据源与账户准备下载前的两件小事MODIS 数据通过公开的数据分发中心向全球开放下载前需要注册账户并申请下载凭据一串由字符组成的 token 或密钥就是后续脚本的“门禁”。常见的做法是把它写入环境变量而不是直接硬编码在脚本里这样既避免误提交泄露也方便在不同机器间迁移。另一个要注意的坑是产品版本号中心里同时存在新旧版本数据建议新任务直接用较新集合但历史任务务必保持版本一致否则不能做时间序列对比。本地目录结构我也建议一开始就固定下来。下面这种分层方式不算花哨但好用原始 HDF 文件按“产品名/日期/图幅”存中间产物单波段 TIFF单独放最终拼接裁剪结果再放一处。这个习惯能在批处理跑到第 200 个文件时帮你省掉大量排查时间。3.2 用检索接口按空间和时间找文件一条 Python 请求官方分发中心提供程序化检索接口返回结果是 JSON 列表。要点有三条一是用经纬度范围而不是图幅号去检索接口会返回命中的图幅列表二是在返回字段里再按 h/v 过滤避免把相邻图幅的重复数据全部下回来三是时间范围写成闭区间并明确起始日期。把这段逻辑封装成函数后面逐日、逐旬批量下载就只剩一个循环。import requests # 检索入口与凭据放配置避免散落在代码里 SEARCH_URL YOUR_DATA_PORTAL_SEARCH_API TOKEN YOUR_ARCHIVE_TOKEN def query_modis_files(product, h25, v5, start2020-01-01, end2020-03-31): params { short_name: product, bounding_box: 110,25,120,35, # 西,南,东,北 temporal: f{start}/{end}, } headers {Authorization: fBearer {TOKEN}} resp requests.get(SEARCH_URL, paramsparams, headersheaders, timeout30) resp.raise_for_status() items resp.json().get(items, []) urls [] for item in items: # 很多产品名携带图幅编号如 h25v05直接字符串过滤 if fh{h:02d}v{v:02d} in item.get(dataset_id, ): urls.append(item[download_url]) return urls files query_modis_files(MOD09GA) print(f命中 {len(files)} 个文件)这段代码里需要解释两个参数。bounding_box 的顺序是西经、南纬、东经、北纬写成 120,25,110,35 这种顺序会直接报错或返回空列表temporal 的起止格式必须补零到八位否则接口不会按你理解的“宽松匹配”处理。response 的字段名在实际环境中可能不同建议第一次调用时先 print 原始 JSON 再调整字段名不要照抄。函数里用字符串过滤 h/v 虽然粗暴但远比逐个判断经纬度交汇要可靠因为产品命名中的图幅编号是固定的两位编码。3.3 批量下载与断点续传把下载工具换成脚本单景 MODIS 文件从几十兆到几百兆不等网络中断是最常见的问题。直接用浏览器手工下载十几二十个文件还能忍受当任务扩展到上百个图幅、跨数月的时间范围时必须把下载交给脚本管理。我一般用 requests 的流式下载配合临时后缀实现断点续传先把文件写入 .part完整下载后才重命名为正式文件下次运行遇到已存在的正式文件就直接跳过。from pathlib import Path import time def download(url, dest_dir): dest Path(dest_dir) / url.split(/)[-1] if dest.exists() and dest.stat().st_size 0: print(f已存在跳过 {dest.name}) return dest tmp dest.with_suffix(dest.suffix .part) with requests.get(url, streamTrue, timeout60) as r: r.raise_for_status() with open(tmp, wb) as f: for chunk in r.iter_content(chunk_size1 20): f.write(chunk) tmp.rename(dest) return dest for url in files: try: download(url, ./data/MOD09GA) except Exception as exc: print(f下载失败 {url}: {exc}) time.sleep(10)这里有个容易被新手忽视的点断点续传的核心不是“接着下”而是“保证文件完整”。如果只按是否存在来跳过半截的 .part 文件会被误判为完整。所以上面逻辑里临时文件始终用独立后缀完整下载后才原子重命名下次运行看到正式文件就确认这份数据可用。另外短时间大量请求容易触发中心限流我在循环里加了 10 秒等待实际使用时如果文件数量极大建议把 sleep 提升到 20 至 30 秒或者采用官方推荐的多线程批量客户端而不要自己硬并发。4. 从 HDF 到可用栅格投影、裁剪、镶嵌与指数计算4.1 认识 HDF 里的数据层SDS 子数据集与比例因子MODIS 产品大多数是 HDF4-EOS 封装一个 .hdf 文件里并不是一张可以直接用软件打开的完整图片而是包含多个科学数据集SDS和全局属性的容器。比如 MOD09GA 一个文件里同时装着 250 米地表反射率波段 1 和 2、500 米反射率波段 3 到 7、质量评估波段 StateQA 等多个数据层。新手常见的错误是直接双击 .hdf 文件看到黑乎乎一块就以为数据坏了其实只是没选中正确的子数据集。规范化做法是先列出文件内部结构再决定提取哪一层。命令行工具可以立即看到每个 SDS 的完整路径和尺度信息gdalinfo MOD09GA.A2020081.h25v05.006.2020081.hdf输出里会列出所有 Subdatasets每一行都形如带前缀的完整路径。这一步很关键因为后续 gdal_translate 或 Python 打开时都必须引用这些完整路径而不是简单写文件名。同时在全局 Metadata 里你会看到比例因子scale_factor和有效值范围例如反射率存储为 int16 整数需要乘以 0.0001 才得到真实反射率。先看 gdalinfo 再写处理代码能省掉后面一半的调试时间。4.2 用 GDAL 做投影转换与重采样一条命令行把目标 SDS 转成 GeoTIFF 并完成投影转换用 gdalwarp 一条命令就能完成。但前提是你已经从 gdalinfo 里拿到了正确的 SDS 路径。下面示例把 MOD09GA 的第一个 250 米反射率波段转成经纬度坐标gdalwarp -t_srs EPSG:4326 \ -tr 0.0025 0.0025 \ -r bilinear \ HDF4_EOS:EOS_GRID:MOD09GA.A2020081.h25v05.006.2020081.hdf:250m Surface Reflectance Band 1 \ mod09_b1_wgs84.tif解释一下三个主要参数。-t_srs 指定输出坐标系EPSG:4326 是标准经纬度如果你的分析用 UTM 分带更合适把它换成对应的 EPSG 编号即可。-tr 控制输出像元分辨率0.0025 度约等于 250 米但如果只是做区域统计而非精确面积测量可以放大到 0.005 度以减少计算量。-r 指定重采样算法连续地表参数用 bilinear 能保持渐变分类或掩膜数据建议用 nearest 避免产生非整数类别。如果系统里 HDF4 驱动支持不完整可以先 gdal_translate 把 SDS 单独解包成 TIFF再对 TIFF 执行 gdalwarp效果完全等价。4.3 批量拼接与范围裁剪VRT 是最省的 IO 方案跨图幅的分析不可避免要把相邻的 h/v 图幅拼起来。先对每个文件单独做投影转换再拼接会在图幅接缝处留下明显的重采样痕迹和可能的错位更好的顺序是先用全部原始投影的 SDS 构建虚拟镶嵌VRT再由一次 gdalwarp 同时完成投影、拼接和裁剪。VRT 本身不拷贝数据只是记录各输入文件的空间范围和波段信息真正重采样发生在最终输出时IO 开销小且避免中间文件堆积。gdalbuildvrt mosaic.vrt band1_h25v04.tif band1_h25v05.tif \ band1_h26v04.tif band1_h26v05.tif gdalwarp -t_srs EPSG:4326 \ -cutline study_area.shp \ -crop_to_cutline \ -of GTiff \ -co COMPRESSDEFLATE \ mosaic.vrt scene.tif-crop_to_cutline 的作用是让输出范围严格贴合矢量边界不额外保留四周空白避免后续统计时把透明像元当成有效值。没有矢量边界时也可以直接用 -te 指定四至范围效果类似但少一步文件依赖。输出格式建议加上 COMPRESSDEFLATE遥感栅格的压缩率通常很可观而且读写速度未必变慢。做完这一步用 gdalinfo 或 Python 读一次统计值确认没有整块无效区域再进入指数计算。4.4 算 NDVI 的两种姿势直接读 SDS 还是先出 TIFF计算植被指数有两条路直接用 GDAL 读 HDF 里的反射率 SDS 做运算或者先把需要波段转成 TIFF 再算。前者省磁盘空间但对要反复测试、比较多套参数的场景每次运行都重新解析 HDF 会很慢我一般先转出少数字段波段后续所有计算都基于普通 TIFF。下面是一个把 MOD09GA 红光与近红外波段合成 NDVI 的示例from osgeo import gdal import numpy as np def read_band(path): ds gdal.Open(path) arr ds.ReadAsArray().astype(float32) geo ds.GetGeoTransform() proj ds.GetProjection() ds None return arr, geo, proj red, geo, proj read_band(mod09_b1_wgs84.tif) # 红光 nir, _, _ read_band(mod09_b2_wgs84.tif) # 近红外 red * 0.0001 # 乘比例因子还原真实反射率 nir * 0.0001 valid (red -0.1) (nir -0.1) (red 1.0) (nir 1.0) ndvi np.where(valid, (nir - red) / (nir red 1e-6), -9999) driver gdal.GetDriverByName(GTiff) out driver.Create(ndvi.tif, ndvi.shape[1], ndvi.shape[0], 1, gdal.GDT_Float32) out.SetGeoTransform(geo) out.SetProjection(proj) band out.GetRasterBand(1) band.WriteArray(ndvi) band.SetNoDataValue(-9999) out.FlushCache() out None这段代码的细节都在取值上。MODIS 反射率里存储的是 int16 整数真实反射率等于存储值乘 0.0001负的存储值表示无效或填充。valid 数组排除了极端值和无效值再在分母加上 1e-6 防止除零。写入时把 NoData 明确设为 -9999这样后端的统计、配色和裁剪工具都能识别无效像元。投影和地理变换信息从输入 TIFF 原样拷贝保证输出和输入严格套合。5. MODIS 处理常见问题排查5 个高频坑与解法5.1 现象QGIS 打开 HDF 只看到一片黑或整张单色很多人第一次用桌面软件打开 .hdf发现影像不是黑的就是一个颜色以为数据损坏。原因在于 HDF 是容器文件软件默认可能只渲染了第一个 SDS而这个 SDS 恰好不是你要的反射率层。解决方法是先在 gdalinfo 里确认子数据集列表再用 gdal_translate 转出目标 SDS 成 TIFF如果软件支持选择子数据集QGIS 里会提示也可以直接在下拉列表里选中目标 SDS但批量处理时还是建议统一走命令行避免交互操作不可复现。5.2 现象拼接后图幅接缝处出现明显错位或白线这是最让人头疼的一类问题因为它不报错。原因是相邻图幅之间在原始产品里本身就有一行像素的重叠如果逐图幅做完投影转换再叠加重采样会让重叠区数值略有差异拼接结果里就出现细微的双边或白缝。解决方法是把所有输入先放入 gdalbuildvrt让 gdalwarp 一次性完成投影、重采样和拼接重叠区交给插值算法统一决定输出值对 NDVI 这类连续指数即使有接缝数值差异也会被控制在很小的范围。另一个相关现象是裁剪后边缘出现黑色斜线那是 NoData 像元跨越了裁剪边界栅格统计前必须检查 NoData 设置。5.3 现象NDVI 时间序列里出现大量负峰或突变如果曲线在不应变绿的季节忽然掉到负数以下往往不是植被真的枯死而是云、云阴影或气溶胶污染没有被剔除。逐日反射率产品没有做时间合成任何一次云遮挡都会表现为异常低值。解决方法是使用质量评估波段 StateQA 做掩膜把云和云阴影对应的像元标记为无效再对时间序列做平滑或最大值合成。这个过程我在下一章给出具体位运算示例。记住一个原则逐日 MODIS 产品在进入时间序列分析前必须先做质量控制否则噪声会完全掩盖真实信号。5.4 现象下载后文件大小对但无法打开或转格式时报错文件下载完成、大小也和页面一致但 GDAL 一打开就报错通常是因为下载中途网络中断客户端保存了一个截断文件而恰好保留了原文件大小的部分信息。解决方法是采用我在下载章节里写的 .part 临时文件策略完整下载再改名同时看到 .part 残留时不要手动改名直接删除重下。另一种情况是产品号匹配错误拿着 MODIS 产品短名去请求了另一颗卫星的同名产品文件名里 A 与 M 前缀看错。下载后第一步就运行 gdalinfo能打开且 Metadata 里产品名与日期正确再进行后续流程。5.5 现象投影转换后经纬度偏差几十公里这个坑比较隐蔽但后果严重。常见原因有两种一是把图幅号选错又没做交叉验证比如实际想分析 h26v05检索返回了 h25v05 的文件也包含目标区域边缘转换后自然对不齐二是直接用软件内置的“自动投影”去解读 SIN 投影却忽略了 MODIS 产品的投影定义细节。解决方法是下载时按经纬度检索并固定 h/v 过滤转换后用矢量图层叠加目检一次我自己会在处理流程里留一个“抽样对比”步骤随机挑三五个点比对 gdalinfo 里的坐标与真实经纬度偏差超过一个像元就停住排查。6. 进阶把处理流程做成批处理的三个细节先处理一个高频需求能不能自动读取比例因子而不是每个脚本里硬编码 0.0001可以。用 GDAL 打开 SDS 后Metadata 里就带着 scale_factor 和 add_offset读取这两个键值替代常量代码在不同产品间迁移时不用改参数也避免多产品混用同一个比例因子的玄学错误。下面是从 SDS 元数据自动取 scale 的写法def read_scaled_band(sds_path): ds gdal.Open(sds_path) meta ds.GetMetadata() scale float(meta.get(scale_factor, 1.0)) offset float(meta.get(add_offset, 0.0)) arr ds.ReadAsArray().astype(float32) * scale offset return arr第二个细节是 QA 掩膜的位运算。StateQA 波段里每个像元值是二进制位标记位不同 bit 表示云、云阴影、气溶胶等状态。例如常见做法是取 bit 0-1 表示总体质量值 00 为优、01 为可用bit 2-3 表示云检测10 或 11 说明有云。想稳定地筛选晴空像元就按位与运算提取有效位再置无效qa ds.ReadAsArray().astype(uint16) # 取 bit0-1值为 0 或 1 时视为可用 quality_ok (qa 0b11) in (0, 1) # 取 bit10-11 等云掩膜位按产品文档定义 cloud_ok ((qa 10) 0b11) 0 valid_mask quality_ok cloud_ok第三个细节是批处理流程的断点设计。别把下载、转投影、算指数写在同一个巨大脚本里一次跑完而是每一步都输出到独立目录下一步从目录里扫描输入。这样某一步翻车时前面的中间文件还能复用我习惯用“日期-图幅-波段”的文件名模板比如 2020081_h25v05_b1.tif方便后续循环遍历和断点续跑。我第一次做多时相 MODIS 分析时全程没有加 QA 筛选结果 NDVI 时序曲线全是锯齿状噪声事后返工才发现是云污染在作祟从此流程里强制加一步质量评估宁肯丢失部分像元也不让坏值进统计。这三个细节看着不起眼但能把批处理从“能跑”推进到“能重复、可追溯”也希望帮到你。本文还有配套的精品资源点击获取